Last updated: 2025-06-03

Checks: 6 1

Knit directory: single-cell-jamboree/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(1)

The command set.seed(1) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: detected

The following chunks had caches available:
  • flashier-nmf

To ensure reproducibility of the results, delete the cache directory mcf7_cache and re-run the analysis. To have workflowr automatically delete the cache directory prior to building the file, set delete_cache = TRUE when running wflow_build() or wflow_publish().

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 1ac5a52

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 1ac5a52. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 


Untracked files:
    Untracked:  analysis/mcf7_cache/
    Untracked:  data/GSE132188_adata.h5ad.h5
    Untracked:  data/Immune_ALL_human.h5ad
    Untracked:  data/pancreas_endocrine.RData
    Untracked:  data/pancreas_endocrine_alldays.h5ad

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/mcf7.Rmd) and HTML (docs/mcf7.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 1ac5a52 Peter Carbonetto 2025-06-03 wflow_publish("mcf7.Rmd", verbose = TRUE, view = FALSE)
html 5dceacd Peter Carbonetto 2025-06-03 A few small adjustments to the topic modeling results in the mcf7 analysis.
Rmd 4499b61 Peter Carbonetto 2025-06-03 wflow_publish("mcf7.Rmd", verbose = TRUE, view = FALSE)
Rmd d3cb20f Peter Carbonetto 2025-06-03 Added topic model fit and structure plot to mcf7 analysis.
Rmd b191ab5 Peter Carbonetto 2025-06-03 Added PCA and GLM-PCA results to mcf7 analysis.
html a0b1bbe Peter Carbonetto 2025-06-03 First build of the mcf7 analysis.
Rmd a3d82ce Peter Carbonetto 2025-06-03 wflow_publish("mcf7.Rmd", verbose = TRUE, view = FALSE)

The MCF-7 data set from Sanford et al eLife paper is helpful for illustrating some of the basic ideas behind NMF analysis of single-cell data (even though it is actually a bulk RNA-seq data set!). This paper analyzed the transcriptional response of human MCF-7 cells to retinoic acid and TGF-\(\beta\), applied individually and in combination. These data were downloaded from GEO accession GSE152749.

The following abbreviations are used in analysis: EtOH = ethanol; RA = retinoic acid; TGFb = TGF-\(\beta\).

Load packages used to process the data, perform the analyses, and create the plots.

library(rsvd)
library(fastglmpca)
library(fastTopics)
library(flashier)
library(data.table)
library(GEOquery)
library(singlecelljamboreeR)
library(ggplot2)
library(cowplot)

Set the seed for reproducibility:

set.seed(1)

Prepare the data for analysis with fastTopics and flashier

Load the RNA-seq counts:

counts <- fread("../data/GSE152749_raw_counts_GRCh38.p13_NCBI.tsv.gz",
                sep = "\t",header = TRUE,stringsAsFactors = FALSE)
class(counts) <- "data.frame"
rownames(counts) <- counts$GeneID
counts <- counts[,-1]
counts <- as.matrix(counts)
storage.mode(counts) <- "double"
counts <- t(counts)
ids <- rownames(counts)

Load the gene information:

genes <- fread("../data/Human.GRCh38.p13.annot.tsv.gz",sep = "\t",
               header = TRUE,stringsAsFactors = FALSE)
class(genes) <- "data.frame"
genes <- genes[1:10]
genes <- transform(genes,
                   GeneType = factor(GeneType),
                   Status   = factor(Status))

Load the sample information:

geo <- getGEO(filename = "../data/GSE152749_family.soft.gz")
samples <- data.frame(id = names(GSMList(geo)),
                      treatment = sapply(GSMList(geo),
                                         function (x) Meta(x)$title))
samples <- samples[ids,]
rownames(samples) <- NULL
samples <- transform(samples,
                     EtOH = grepl("EtOH",treatment,fixed = TRUE),
                     RA   = grepl("RA",treatment,fixed = TRUE),
                     TGFb = grepl("TGFb",treatment,fixed = TRUE))
samples$label                 <- "EtOH"
samples[samples$RA,"label"]   <- "RA"
samples[samples$TGFb,"label"] <- "TGFb"
samples[with(samples,RA & TGFb),"label"] <- "RA+TGFb"
samples <- transform(samples,
                     label = factor(label,c("EtOH","RA","TGFb","RA+TGFb")))

Remove the non-protein-coding genes, and the genes that are expressed in fewer than 4 samples:

x <- colSums(counts > 0)
i <- which(x > 3 &
           genes$GeneType == "protein-coding" &
           genes$Status == "active")
genes  <- genes[i,]
counts <- counts[,i]

(Some background: it turns out that there is some structure in the non-coding RNA genes that is is picked up by the Poisson-based methods but not by the Gaussian-based methods, so to avoid this complication in the comparisons I have removed these genes from the data, which aren’t of interest anyhow.)

This is the dimension of the data set we will analyze:

dim(counts)
# [1]    41 16773

For the Gaussian-based analyses, we will need the shifted log counts:

a <- 1
s <- rowSums(counts)
s <- s/mean(s)
shifted_log_counts <- log1p(counts/(a*s))

PCA and GLM-PCA

Let’s see what happens when we apply PCA to the shifted log counts and GLM-PCA to the counts. As we will see, PCA and GLM-PCA both show a clear clustering of the data that corresponds to the different treatments. However, the PCs are not individually interpretable.

Run PCA on the shifted log counts, and plot the first 2 PCs:

pca <- rpca(shifted_log_counts,k = 2,center = TRUE,scale = FALSE)
colnames(pca$x) <- c("PC1","PC2")
pdat <- data.frame(samples,pca$x)
ggplot(pdat,aes(x = PC1,y = PC2,color = label)) +
  geom_point() +
  scale_color_manual(values = c("dodgerblue","tomato","darkblue",
                                "limegreen")) + 
  theme_cowplot(font_size = 10)

Version Author Date
5dceacd Peter Carbonetto 2025-06-03

GLM-PCA applied to the counts essentially produces the same result (aside from an arbitrary rotation):

fit_glmpca <- init_glmpca_pois(t(counts),K = 2)
fit_glmpca <- fit_glmpca_pois(t(counts),fit0 = fit_glmpca,verbose = FALSE,
                              control = list(maxiter = 50))
colnames(fit_glmpca$V) <- c("k1","k2")
pdat <- data.frame(samples,fit_glmpca$V)
ggplot(pdat,aes(x = k1,y = k2,color = label)) +
  geom_point() +
  scale_color_manual(values = c("dodgerblue","tomato","darkblue",
                                "limegreen")) + 
  theme_cowplot(font_size = 10)

Version Author Date
5dceacd Peter Carbonetto 2025-06-03

Topic model (fastTopics)

Fit a topic model with \(K = 3\) topics to the counts:

tm0 <- fit_poisson_nmf(counts,k = 3,init.method = "random",
                       numiter = 50,verbose = "none",
                       control = list(nc = 4,extrapolate = FALSE))
tm <- fit_poisson_nmf(counts,fit0 = tm0,numiter = 50,verbose = "none",
                      control = list(nc = 4,extrapolate = TRUE))

Structure plot comparing the topics to the clusters or treatment conditions:

topic_colors <- c("tomato","darkblue","dodgerblue")
n    <- nrow(tm$L)
L    <- poisson2multinom(tm)$L
rows <- order(pmax(L[,1],L[,2],L[,3]))
L    <- L[rows,]
structure_plot(L,grouping = samples$label[rows],topics = 1:3,
               loadings_order = 1:n,colors = topic_colors) +
  theme(axis.text.x = element_text(angle = 0,hjust = 0.5))

Version Author Date
5dceacd Peter Carbonetto 2025-06-03

EBNMF (flashier)

Next fit an NMF to the shifted log counts using flashier, with \(K = 4\):

n  <- nrow(counts)
x  <- rpois(1e7,1/n)
s1 <- sd(log(x + 1))
fl_nmf <- flash(shifted_log_counts,ebnm_fn = ebnm_point_exponential,
                var_type = 2, greedy_Kmax = 4,S = s1,backfit = TRUE,
                verbose = 0)

Warning: The above code chunk cached its results, but it won’t be re-run if previous chunks it depends on are updated. If you need to use caching, it is highly recommended to also set knitr::opts_chunk$set(autodep = TRUE) at the top of the file (in a chunk that is not cached). Alternatively, you can customize the option dependson for each individual chunk that is cached. Using either autodep or dependson will remove this warning. See the knitr cache options for more details.

EBNMF applied to the shifted log count decomposes the samples in a similar way as the topic model:

topic_colors <- c("olivedrab","dodgerblue","darkblue","tomato")
L    <- ldf(fl_nmf,type = "i")$L
rows <- order(pmax(L[,2],L[,3],L[,4]))
L    <- L[rows,]
structure_plot(L,grouping = samples[rows,"label"],topics = 4:1,
               loadings_order = 1:n,colors = topic_colors) +
  labs(y = "membership") +
  theme(axis.text.x = element_text(angle = 0,hjust = 0.5))


sessionInfo()
# R version 4.3.3 (2024-02-29)
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS 15.4.1
# 
# Matrix products: default
# BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
# LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# time zone: America/Chicago
# tzcode source: internal
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#  [1] cowplot_1.1.3             ggplot2_3.5.0            
#  [3] singlecelljamboreeR_0.1-3 GEOquery_2.70.0          
#  [5] Biobase_2.62.0            BiocGenerics_0.48.1      
#  [7] data.table_1.15.2         flashier_1.0.55          
#  [9] ebnm_1.1-34               fastTopics_0.7-25        
# [11] fastglmpca_0.1-108        rsvd_1.0.5               
# 
# loaded via a namespace (and not attached):
#   [1] distr_2.9.3          pbapply_1.7-2        rlang_1.1.5         
#   [4] magrittr_2.0.3       git2r_0.33.0         horseshoe_0.2.0     
#   [7] compiler_4.3.3       vctrs_0.6.5          reshape2_1.4.4      
#  [10] daarem_0.7           quadprog_1.5-8       stringr_1.5.1       
#  [13] pkgconfig_2.0.3      crayon_1.5.2         fastmap_1.1.1       
#  [16] labeling_0.4.3       utf8_1.2.4           promises_1.2.1      
#  [19] rmarkdown_2.26       tzdb_0.4.0           purrr_1.0.2         
#  [22] xfun_0.42            cachem_1.0.8         trust_0.1-8         
#  [25] jsonlite_1.8.8       progress_1.2.3       highr_0.10          
#  [28] later_1.3.2          irlba_2.3.5.1        parallel_4.3.3      
#  [31] prettyunits_1.2.0    R6_2.5.1             bslib_0.6.1         
#  [34] stringi_1.8.3        RColorBrewer_1.1-3   SQUAREM_2021.1      
#  [37] limma_3.58.1         jquerylib_0.1.4      Rcpp_1.0.12         
#  [40] knitr_1.45           R.utils_2.12.3       readr_2.1.5         
#  [43] httpuv_1.6.14        Matrix_1.6-5         splines_4.3.3       
#  [46] tidyselect_1.2.1     yaml_2.3.8           lattice_0.22-5      
#  [49] tibble_3.2.1         plyr_1.8.9           withr_3.0.2         
#  [52] evaluate_1.0.3       Rtsne_0.17           RcppParallel_5.1.10 
#  [55] startupmsg_0.9.6.1   xml2_1.3.6           pillar_1.9.0        
#  [58] whisker_0.4.1        plotly_4.10.4        softImpute_1.4-1    
#  [61] generics_0.1.3       rprojroot_2.0.4      invgamma_1.1        
#  [64] truncnorm_1.0-9      hms_1.1.3            munsell_0.5.0       
#  [67] scales_1.3.0         ashr_2.2-66          gtools_3.9.5        
#  [70] RhpcBLASctl_0.23-42  glue_1.8.0           scatterplot3d_0.3-44
#  [73] lazyeval_0.2.2       tools_4.3.3          fs_1.6.5            
#  [76] grid_4.3.3           tidyr_1.3.1          colorspace_2.1-0    
#  [79] sfsmisc_1.1-18       deconvolveR_1.2-1    cli_3.6.4           
#  [82] Polychrome_1.5.1     workflowr_1.7.1      fansi_1.0.6         
#  [85] mixsqp_0.3-54        viridisLite_0.4.2    dplyr_1.1.4         
#  [88] uwot_0.2.3           gtable_0.3.4         R.methodsS3_1.8.2   
#  [91] sass_0.4.9           digest_0.6.34        ggrepel_0.9.5       
#  [94] farver_2.1.1         htmlwidgets_1.6.4    R.oo_1.26.0         
#  [97] htmltools_0.5.8.1    lifecycle_1.0.4      httr_1.4.7          
# [100] statmod_1.5.0        MASS_7.3-60.0.1