Last updated: 2025-06-18

Checks: 4 3

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: uncommitted changes

The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Environment: objects present

The global environment had objects present when the code in the R Markdown file was run. These objects 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. Use wflow_publish or wflow_build to ensure that the code is always run in an empty environment.

The following objects were defined in the global environment when these results were created:

Name Class Size
a numeric 56 bytes
celltype_topics character 504 bytes
clusters factor 13.2 Kb
counts dgCMatrix 148.6 Mb
de_le topic_model_de_analysis;list 25.1 Mb
de_vsnull topic_model_de_analysis;list 25.1 Mb
fl_nmf flash;list 164.8 Mb
genes data.frame 2.5 Mb
i integer 4.3 Kb
j integer 1.4 Mb
L matrix;array 221.5 Kb
n integer 56 bytes
other_topics character 432 bytes
out data.frame 283.4 Kb
outliers character 320 bytes
p1 gg;ggplot 227.3 Kb
p2 gg;ggplot 120.2 Kb
s numeric 294.1 Kb
s1 numeric 56 bytes
samples data.frame 370.5 Kb
shifted_log_counts dgCMatrix 148.4 Mb
tm poisson_nmf_fit;list 10.7 Mb
tm_merged multinom_topic_model_fit;list 9.7 Mb
x integer 38.1 Mb

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:
  • de-analysis
  • fit-topic-model
  • flashier-nmf
  • flashier-nmf-k-12

To ensure reproducibility of the results, delete the cache directory pancreas_cytokine_S1_factors_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: b7f4df5

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 b7f4df5. 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/lps_cache/
    Untracked:  analysis/mcf7_cache/
    Untracked:  analysis/pancreas_cytokine_S1_factors_cache/
    Untracked:  analysis/temp2.R
    Untracked:  data/GSE132188_adata.h5ad.h5
    Untracked:  data/GSE183010/
    Untracked:  data/Immune_ALL_human.h5ad
    Untracked:  data/pancreas_cytokine.RData
    Untracked:  data/pancreas_endocrine.RData
    Untracked:  data/pancreas_endocrine_alldays.h5ad

Unstaged changes:
    Modified:   analysis/pancreas_cytokine_S1_factors.Rmd

Staged changes:
    New:        data/pancreas_cytokine_S1_tm_k=13_clusters.csv

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/pancreas_cytokine_S1_factors.Rmd) and HTML (docs/pancreas_cytokine_S1_factors.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
html b7f4df5 Peter Carbonetto 2025-06-18 Ran wflow_publish("pancreas_cytokine_S1_factors.Rmd").
Rmd 928cbd3 Peter Carbonetto 2025-06-18 Split the Structure plot for flashier into two plots in the pancreas_cytokine_S1_factors analysis.
Rmd 710d6b6 Peter Carbonetto 2025-06-17 Added de_analysis calls to the pancreas_cytokine_S1_factors analysis.
Rmd d584b40 Peter Carbonetto 2025-06-17 Small fix to the structure plot.
Rmd bd05725 Peter Carbonetto 2025-06-17 Added flashier NMF analysis to pancreas_cytokine_S1_factors, with k=9.
Rmd 28b597c Peter Carbonetto 2025-06-17 Reworking the topic modeling results with k=13 in the pancreas_cytokine_S1_factors analysis.
Rmd 9d95144 Peter Carbonetto 2025-06-16 Working on changes to the pancreas_cytokine_S1_factors analysis (still a work-in-progress).
Rmd 8988553 Peter Carbonetto 2025-06-13 Working on a bunch of changes to the pancreas_cytokine_S1_factors analysis.
Rmd 7442af1 Peter Carbonetto 2025-06-12 A few fixes to the code for the k=13 fits in the pancreas_cytokine_S1_factors analysis.
Rmd 8789250 Peter Carbonetto 2025-06-12 Added k=13 fits to the pancreas_cytokine_S1_factors analysis.
Rmd 18a86f3 Peter Carbonetto 2025-06-11 A couple small changes to pancreas_cytokine_S1_factors.Rmd.
Rmd 89d3f1e Peter Carbonetto 2025-06-11 Added a link to the pancreas_cytokine_S1_factors analysis.
html 1507be2 Peter Carbonetto 2025-06-11 Fixed up structure plots and added annotation heatmaps to the
Rmd 2153b30 Peter Carbonetto 2025-06-11 wflow_publish("pancreas_cytokine_S1_factors.Rmd", verbose = TRUE)
Rmd 980e670 Peter Carbonetto 2025-06-11 Fixed the clustering for the pancreas_cytokine data slightly.
Rmd d1fdbe9 Peter Carbonetto 2025-06-11 Made a few improvements to the pancreas_cytokine_S1_factors analysis.
Rmd ce314bb Peter Carbonetto 2025-06-09 First try at running fastTopics and flashier on the pancreas_cytokine data, for mouse = S1 only; from this analysis I learned that I need to remove the mt and rp genes.
Rmd 422c8ed Peter Carbonetto 2025-06-09 Added steps to the pancreas_cytokine_S1_factors analysis to prepare the data for fastTopics and flashier.
Rmd 46ba21a Peter Carbonetto 2025-06-06 Started new analysis in pancreas_cytokine_S1_factors.Rmd.

Here we perform a NMF analyses of the “pancreas cytokine” data set, focussing on the scRNA-seq data from untreated mouse only.

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

library(Matrix)
library(fastTopics)
library(flashier)
library(singlecelljamboreeR)
library(ggplot2)
library(cowplot)

Set the seed for reproducibility:

set.seed(1)

Load the prepared data set:

load("../data/pancreas_cytokine.RData")

Here we will analyze the cells from the untreated mouse only:

i       <- which(samples$mouse == "S1")
samples <- samples[i,]
counts  <- counts[i,]

Remove three cells that appear to be outliers (one of them appears to be an acinar cell based on Eric’s analysis):

outliers <- c("TTTGTTGTCGTTAGTG-1","TTTGTTGGTAGAGCTG-1","CCCAACTCACTCATAG-1")
i        <- which(!is.element(samples$barcode,outliers))
samples  <- samples[i,]
counts   <- counts[i,]

Remove genes that are expressed in fewer than 5 cells:

j      <- which(colSums(counts > 0) > 4)
genes  <- genes[j,]
counts <- counts[,j]

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

dim(counts)
# [1]  3136 16359

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

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

Topic model (fastTopics)

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

set.seed(1)
tm <- fit_poisson_nmf(counts,k = 13,init.method = "random",method = "em",
                      numiter = 40,verbose = "none",
                      control = list(numiter = 4,nc = 8,extrapolate = FALSE))
tm <- fit_poisson_nmf(counts,fit0 = tm,method = "scd",numiter = 40,
                      control = list(numiter = 4,nc = 8,extrapolate = TRUE),
                      verbose = "none")

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.

Structure plots comparing the topics to the clusters (some of which are inferred from the topics):

set.seed(1)
celltype_topics <- paste0("k",c(2,3,5,7:9,13))
other_topics <- paste0("k",c(4,1,6,10:12))
L <- poisson2multinom(tm)$L
clusters <- as.character(samples$cluster)
clusters[clusters == "islet"]                  <- "beta"
clusters[clusters == "beta" & L[,"k3"] > 0.25] <- "alpha"
clusters[clusters == "beta" & L[,"k8"] > 0.25] <- "delta+epsilon"
clusters[clusters == "beta" & L[,"k9"] > 0.25] <- "gamma"
clusters <- factor(clusters,c("beta","alpha","delta+epsilon","gamma","duct",
                              "endothelial-mesenchymal","macrophage"))
i <- c(sample(which(clusters == "beta"),400),
       which(clusters != "beta"))
p1 <- structure_plot(L[i,],grouping = clusters[i],topics = celltype_topics,
                     gap = 10,n = Inf) +
  labs(fill = "")
p2 <- structure_plot(L[i,],grouping = clusters[i],topics = other_topics,
                     gap = 10,n = Inf) +
  labs(fill = "")
plot_grid(p1,p2,nrow = 2,ncol = 1)

Version Author Date
b7f4df5 Peter Carbonetto 2025-06-18
1507be2 Peter Carbonetto 2025-06-11

Save this topic-model-based clustering because it could be useful elsewhere:

out <- data.frame(barcode = samples$barcode,
                  cluster = clusters,
                  stringsAsFactors = FALSE)
write.csv(out,"pancreas_cytokine_S1_tm_k=13_clusters.csv",quote = FALSE,
          row.names = FALSE)

Select genes in two ways: (a) genes that are different from the mean expression level, (b) genes that are “distinctive”.

set.seed(1)
tm_merged <- merge_topics(poisson2multinom(tm),paste0("k",c(1,2,4,11)))
de_vsnull <- de_analysis(tm_merged,counts,pseudocount = 0.1,lfc.stat = "vsnull",
                         verbose = FALSE,control = list(ns = 1e4,nc = 8))
de_le <- de_analysis(tm_merged,counts,pseudocount = 0.1,lfc.stat = "le",
                     verbose = FALSE,control = list(ns = 1e4,nc = 8))

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.

Based on the estimated \(\mathbf{F}\), we have the following potential interpretation of the topics:

scale_rows <- function (A)
  A / apply(A,1,max)
marker_genes <- c("Ins1","Ins2","Mafa","Gcg","Mafb","Sst","Ghrl",
                  "Ppy","Chga","Iapp","Krt19","Ccr5","Pecam1","Esam",
                  "Col1a1","Ghrl")
j <- match(marker_genes,genes$symbol)
F <- poisson2multinom(tm)$F
F <- F[j,]
F <- scale_rows(F)
rownames(F) <- marker_genes
topics <- paste0("k",c(5,4,6,7,2,3))
annotation_heatmap(F[,topics],select_features = "all",verbose = FALSE)

EBNMF (flashier)

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

set.seed(1)
n  <- nrow(samples)
x  <- rpois(1e7,1/n)
s1 <- sd(log(x + 1))
fl_nmf <- flash(shifted_log_counts,S = s1,ebnm_fn = ebnm_point_exponential,
                var_type = 2,greedy_Kmax = 9,backfit = FALSE,
                nullcheck = FALSE,verbose = 0)
fl_nmf <- flash_backfit(fl_nmf,extrapolate = FALSE,maxiter = 40,verbose = 0)
fl_nmf <- flash_backfit(fl_nmf,extrapolate = TRUE,maxiter = 80,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.

(Note that I tried setting greedy_Kmax > 9, but this didn’t result in more interesting fits—the additional factors loaded on only a single cell.)

Structure plot comparing the factors to the clusters:

set.seed(1)
L <- ldf(fl_nmf,type = "i")$L
colnames(L) <- paste0("k",1:9)
i <- c(sample(which(clusters == "beta"),400),
       which(clusters != "beta"))
p1 <- structure_plot(L[i,],grouping = clusters[i],topics = c(2:6,8,9),
                    gap = 10,n = Inf) +
  labs(y = "membership",fill = "")
p2 <- structure_plot(L[i,],grouping = clusters[i],topics = c(1,7),
                    gap = 10,n = Inf) +
  labs(y = "membership",fill = "")
plot_grid(p1,p2,nrow = 2,ncol = 1)

Version Author Date
b7f4df5 Peter Carbonetto 2025-06-18
1507be2 Peter Carbonetto 2025-06-11

Possible interpretation of the factors:

scale_cols <- function (A) {
  b <- apply(A,2,max)
  return(t(t(A) * b))
}
marker_genes <- c("Ins1","Ins2","Mafa","Gcg","Mafb","Sst","Ghrl",
                  "Ppy","Chga","Iapp","Krt19",
                  "Ccr5","Pecam1","Esam","Col1a1")
j <- match(marker_genes,genes$symbol)
F <- ldf(fl_nmf,type = "i")$F
F <- scale_cols(F)
F <- F[j,]
rownames(F) <- marker_genes
colnames(F) <- paste0("k",1:7)
factors <- paste0("k",c(4,7,3,5,2,6))
annotation_heatmap(F[,factors],select_eatures = "all",verbose = FALSE)

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 flashier_1.0.55          
# [5] ebnm_1.1-34               fastTopics_0.7-25        
# [7] Matrix_1.6-5              workflowr_1.7.1          
# 
# loaded via a namespace (and not attached):
#  [1] tidyselect_1.2.1     viridisLite_0.4.2    farver_2.1.1        
#  [4] dplyr_1.1.4          fastmap_1.1.1        lazyeval_0.2.2      
#  [7] promises_1.2.1       digest_0.6.34        lifecycle_1.0.4     
# [10] processx_3.8.3       invgamma_1.1         magrittr_2.0.3      
# [13] compiler_4.3.3       rlang_1.1.5          sass_0.4.9          
# [16] progress_1.2.3       tools_4.3.3          utf8_1.2.4          
# [19] yaml_2.3.8           data.table_1.17.4    knitr_1.45          
# [22] labeling_0.4.3       prettyunits_1.2.0    htmlwidgets_1.6.4   
# [25] scatterplot3d_0.3-44 RColorBrewer_1.1-3   plyr_1.8.9          
# [28] Rtsne_0.17           withr_3.0.2          purrr_1.0.2         
# [31] grid_4.3.3           fansi_1.0.6          git2r_0.33.0        
# [34] colorspace_2.1-0     scales_1.3.0         gtools_3.9.5        
# [37] cli_3.6.4            rmarkdown_2.26       crayon_1.5.2        
# [40] generics_0.1.3       RcppParallel_5.1.10  rstudioapi_0.15.0   
# [43] httr_1.4.7           reshape2_1.4.4       pbapply_1.7-2       
# [46] cachem_1.0.8         stringr_1.5.1        splines_4.3.3       
# [49] parallel_4.3.3       softImpute_1.4-1     vctrs_0.6.5         
# [52] jsonlite_1.8.8       callr_3.7.5          hms_1.1.3           
# [55] mixsqp_0.3-54        ggrepel_0.9.5        irlba_2.3.5.1       
# [58] horseshoe_0.2.0      trust_0.1-8          plotly_4.10.4       
# [61] jquerylib_0.1.4      tidyr_1.3.1          glue_1.8.0          
# [64] codetools_0.2-19     ps_1.7.6             uwot_0.2.3          
# [67] stringi_1.8.3        Polychrome_1.5.1     gtable_0.3.4        
# [70] later_1.3.2          quadprog_1.5-8       munsell_0.5.0       
# [73] tibble_3.2.1         pillar_1.9.0         htmltools_0.5.8.1   
# [76] truncnorm_1.0-9      R6_2.5.1             rprojroot_2.0.4     
# [79] evaluate_1.0.3       lattice_0.22-5       highr_0.10          
# [82] RhpcBLASctl_0.23-42  SQUAREM_2021.1       ashr_2.2-66         
# [85] httpuv_1.6.14        bslib_0.6.1          Rcpp_1.0.12         
# [88] deconvolveR_1.2-1    whisker_0.4.1        xfun_0.42           
# [91] fs_1.6.5             getPass_0.2-4        pkgconfig_2.0.3