NMF analysis of the “pancreas cytokine” data set (untreated mouse only)

Last updated: 2025-07-05

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:

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: 4e1bb00

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 4e1bb00. 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:  .DS_Store
    Untracked:  analysis/#pancreas_cytokine_S1_factors.Rmd#
    Untracked:  analysis/.#pancreas_cytokine_S1_factors.Rmd
    Untracked:  analysis/.DS_Store
    Untracked:  analysis/lps_cache/
    Untracked:  analysis/mcf7_cache/
    Untracked:  analysis/pancreas_cytokine_S1_factors_cache/
    Untracked:  data/GSE132188_adata.h5ad.h5
    Untracked:  data/GSE183010/
    Untracked:  data/Immune_ALL_human.h5ad
    Untracked:  data/pancreas_cytokine.RData
    Untracked:  data/pancreas_cytokine_lsa.RData
    Untracked:  data/pancreas_endocrine.RData
    Untracked:  data/pancreas_endocrine_alldays.h5ad
    Untracked:  output/panc_cyto_lsa_res/

Unstaged changes:
    Modified:   analysis/temp2.R

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
Rmd	4e1bb00	Peter Carbonetto	2025-07-05	wflow_publish("pancreas_cytokine_S1_factors.Rmd", verbose = T,
Rmd	8e87671	Peter Carbonetto	2025-07-05	A few more slight improvements to the pancreas_cytokine_S1_factors analysis (still very much a work-in-progress).
Rmd	69e6613	Peter Carbonetto	2025-07-04	Working on improvements/changes to the pancreas_cytokine_S1_factors analysis; these improvements are still a bit rough and need to refine them.
Rmd	65bedc4	Peter Carbonetto	2025-07-04	Updated code for running flashier on the pancreas_cytokine data.
Rmd	9147de4	Peter Carbonetto	2025-07-04	Fixed a bug in lps.Rmd.
Rmd	14db985	Peter Carbonetto	2025-07-04	Added note to pancreas_cytokine_S1_factors.Rmd.
Rmd	a8de13f	Peter Carbonetto	2025-07-04	Updated interpretation of topics in pancreas_cytokine_S1_factors analysis.
html	3d8a2aa	Peter Carbonetto	2025-07-04	Added lps_gsea_fl_nmf.csv output.
Rmd	7476abf	Peter Carbonetto	2025-06-20	Small edit to the pancreas_cytokine_S1_factors analysis.
Rmd	4ac7821	Peter Carbonetto	2025-06-18	Created pancreas_cytokine_S1_tm_k=13_clusters.csv containing the topic-model based clustering for the pancreas cytokine data (untreated mouse only).
html	4ac7821	Peter Carbonetto	2025-06-18	Created pancreas_cytokine_S1_tm_k=13_clusters.csv containing the topic-model based clustering for the pancreas cytokine data (untreated mouse only).
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(NNLM)
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

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

topic 2: $\beta$ cells (Ins1, Ins2, Mafa)
topic 3: $\alpha$ cells (Gcg, Mafb)
topic 5: duct cells (Krt19)
topic 6: replicate effect
topic 7: macrophages (Ccr5)
topic 8: $\delta$ cells (Sst)
topic 9: $\gamma$ cells (Ppy)
topic 12: cell cycle
topic 13: endothelial (Pecam1) and mesenchymal cells (Col1a1).

Topic 10 is clearly capturing a technical difference in the two replicates:

pdat <- cbind(samples,L)
ggplot(pdat,aes(x = replicate,y = k10)) +
  geom_boxplot() +
  theme_cowplot(font_size = 10)

EBNMF (flashier)

Next we will fit an NMF to the shifted log counts using flashier, with $K = 13$. Since the greedy initialization does not seem to work well in this example, I’ll use a different initialization strategy: obtain a “good” initialization using the NNLM package, then use this initialization to fit a NMF using flashier. This approach is implemented in the following function:

flashier_nmf <- function (X, k = 3, n.threads = 1) {
  n <- nrow(X)
  m <- ncol(X)
  Y <- as.matrix(X)
  nmf0 <- nnmf(Y,k = 1,loss = "mse",method = "scd",max.iter = 10,
               verbose = 2,n.threads = n.threads)
  W0 <- nmf0$W
  H0 <- nmf0$H
  W0 <- cbind(W0,matrix(runif(n*(k-1)),n,k-1))
  H0 <- rbind(H0,matrix(runif(m*(k-1)),k-1,m))
  nmf <- nnmf(Y,k,init = list(W = W0,H = H0),loss = "mse",method = "scd",
              max.iter = 10,verbose = 2,n.threads = n.threads)
  x  <- rpois(1e7,1/n)
  s1 <- sd(log(x + 1))
  out <- flash_init(X,var_type = 2,S = s1)
  out <- flash_factors_init(out,list(nmf$W,t(nmf$H)),ebnm_point_exponential)
  out <- flash_backfit(out,extrapolate = FALSE,maxiter = 100,verbose = 2)
  return(flash_backfit(out,extrapolate = TRUE,maxiter = 100,verbose = 2))
}

Now fit an NMF to the shifted log counts, with $K = 13$:

set.seed(1)
fl_nmf <- flashier_nmf(shifted_log_counts,k = 13,n.threads = 8)
# 
#  Iteration |        MSE |        MKL |     Target |  Rel. Err.
# --------------------------------------------------------------
#          1 |     0.1441 |     0.1960 |     0.0721 |          2
#          3 |     0.1426 |     0.1967 |     0.0713 |       0.01
#          5 |     0.1426 |     0.1967 |     0.0713 |      1e-08
# --------------------------------------------------------------
#  Iteration |        MSE |        MKL |     Target |  Rel. Err.
# 
# 
#  Iteration |        MSE |        MKL |     Target |  Rel. Err.
# --------------------------------------------------------------
#          1 |     0.1356 |     0.1913 |     0.0678 |          2
#          3 |     0.1248 |     0.1804 |     0.0624 |       0.08
#          5 |     0.1226 |     0.1779 |     0.0613 |       0.02
#          7 |     0.1218 |     0.1768 |     0.0609 |      0.007
#          9 |     0.1214 |     0.1763 |     0.0607 |      0.003
#         11 |     0.1213 |     0.1762 |     0.0606 |      0.001
# --------------------------------------------------------------
#  Iteration |        MSE |        MKL |     Target |  Rel. Err.
# 
# Backfitting 13 factors (tolerance: 7.64e-01)...
#   Difference between iterations is within 1.0e+03...
#   Difference between iterations is within 1.0e+02...
#   --Maximum number of iterations reached!
# Wrapping up...
# Done.
# Backfitting 13 factors (tolerance: 7.64e-01)...
#   Difference between iterations is within 1.0e+02...
#   --Maximum number of iterations reached!
#   Backfit complete. Objective: 4345377.189 
# Wrapping up...
# Done.

Structure plot comparing the factors to the clusters:

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

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

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

factor 1: $\beta$ cells (Ins1, Ins2, Mafa)
factors 2, 9, 11: endothelial cells (Pecam1, Esam)
factor 3: duct cells (Krt19)
factor 4: macrophages (Ccr5)
factor 5: $\alpha$ (Gcg, Mafb) and $\gamma$ cells (Ppy)
factor 6: replicate effect
factors 7, 12: cell cycle
factor 10: mesenchymal cells (Col1a1)
factor 13: $\delta$ cells (Sst)

Factor 6 almost perfectly captures the technical difference in the two replicates:

pdat <- cbind(samples,L)
r <- cor(as.numeric(samples$replicate),L[,"k6"])
ggplot(pdat,aes(x = replicate,y = k6)) +
  geom_boxplot() +
  ggtitle(paste("cor =",round(r,digits = 3))) +
  theme_cowplot(font_size = 12)

Factors 7 and 12 seem to capture cell cycle:

cell_cycle <- 
c("Aaas", "Abl1", "Abraxas1", "Acd", "Actr1a", "Adrm1", "Ahctf1",
"Ajuba", "Akap9", "Akt1", "Akt2", "Akt3", "Alms1", "Anapc1",
"Anapc10", "Anapc11", "Anapc15", "Anapc16", "Anapc2", "Anapc4",
"Anapc5", "Anapc7", "Ankle2", "Ankrd28", "Arpp19", "Atm", "Atrip",
"Atrx", "Aurka", "Aurkb", "B9d2", "Babam1", "Babam2", "Banf1",
"Bard1", "Blm", "Blzf1", "Bora", "Brca1", "Brcc3", "Brip1", "Bub1b",
"Bub3", "Cables1", "Cc2d1b", "Ccna2", "Ccnb1", "Ccnb2", "Ccnd1",
"Ccnd2", "Ccnd3", "Ccne1", "Ccne2", "Ccnh", "Ccp110", "Cdc14a",
"Cdc16", "Cdc20", "Cdc23", "Cdc25a", "Cdc25b", "Cdc26", "Cdc27",
"Cdc45", "Cdc6", "Cdc7", "Cdca5", "Cdca8", "Cdk1", "Cdk11b",
"Cdk2", "Cdk4", "Cdk5rap2", "Cdk6", "Cdk7", "Cdkn1a", "Cdkn1b",
"Cdkn1c", "Cdkn2b", "Cdt1", "Cenpa", "Cenpc1", "Cenpe", "Cenpf",
"Cenph", "Cenpi", "Cenpj", "Cenpk", "Cenpl", "Cenpm", "Cenpn",
"Cenpo", "Cenpp", "Cenpq", "Cenps", "Cenpt", "Cenpu", "Cenpw",
"Cenpx", "Cep131", "Cep135", "Cep152", "Cep164", "Cep192", "Cep250",
"Cep290", "Cep41", "Cep57", "Cep63", "Cep70", "Cep72", "Cep76",
"Cep78", "Cetn2", "Chek1", "Chek2", "Chmp2a", "Chmp2b", "Chmp3",
"Chmp4b", "Chmp4c", "Chmp6", "Chmp7", "Chtf18", "Chtf8", "Ckap5",
"Cks1b", "Clasp1", "Clasp2", "Clip1", "Clspn", "Cnep1r1", "Cop1",
"Csnk1d", "Csnk1e", "Csnk2a1", "Csnk2a2", "Csnk2b", "Ctc1", "Ctdnep1",
"Cul1", "Daxx", "Dbf4", "Dctn1", "Dctn2", "Dctn3", "Dkc1", "Dmc1",
"Dna2", "Dscc1", "Dsn1", "Dync1h1", "Dync1i1", "Dync1i2", "Dync1li1",
"Dync1li2", "Dynll1", "Dynll2", "Dyrk1a", "E2f1", "E2f2", "E2f3",
"E2f4", "E2f5", "Emd", "Eml4", "Ercc6l", "Esco1", "Esco2", "Espl1",
"Exo1", "Fbxl18", "Fbxl7", "Fbxo5", "Fbxw11", "Fen1", "Fignl1",
"Fkbpl", "Foxm1", "Fzr1", "Gar1", "Gins1", "Gins2", "Gins3",
"Gins4", "Gmnn", "Golga2", "Gorasp1", "Gorasp2", "Gtse1", "H3f3a",
"H3f3b", "Haus1", "Haus2", "Haus3", "Haus4", "Haus5", "Haus6",
"Haus7", "Haus8", "Hdac1", "Hdac8", "Herc2", "Hjurp", "Hmmr",
"Hsp90aa1", "Hsp90ab1", "Hus1", "Incenp", "Ist1", "Itgb3bp",
"Jak2", "Kat5", "Kif18a", "Kif20a", "Kif23", "Kif2a", "Kif2c",
"Kntc1", "Kpnb1", "Lbr", "Lcmt1", "Lig1", "Lin37", "Lin52", "Lin54",
"Lin9", "Lmna", "Lmnb1", "Lpin2", "Lpin3", "Lyn", "Mad1l1", "Mad2l1",
"Mapk1", "Mapk3", "Mapre1", "Mastl", "Mau2", "Mcm10", "Mcm2",
"Mcm3", "Mcm4", "Mcm5", "Mcm6", "Mcm7", "Mcm8", "Mcph1", "Mdc1",
"Mdm2", "Mdm4", "Mis12", "Mis18a", "Mis18bp1", "Mnat1", "Mre11a",
"Mzt1", "Mzt2", "Nbn", "Ncapd2", "Ncapd3", "Ncapg", "Ncapg2",
"Ncaph", "Ncaph2", "Ndc1", "Ndc80", "Nde1", "Ndel1", "Nedd1",
"Nek2", "Nek9", "Nhp2", "Ninl", "Nipbl", "Nme7", "Nop10", "Npm1",
"Nsd2", "Nsl1", "Nudc", "Nuf2", "Numa1", "Nup107", "Nup133",
"Nup153", "Nup155", "Nup160", "Nup188", "Nup205", "Nup210", "Nup214",
"Nup35", "Nup37", "Nup43", "Nup50", "Nup54", "Nup62", "Nup85",
"Nup88", "Nup93", "Nup98", "Odf2", "Ofd1", "Oip5", "Optn", "Orc2",
"Orc3", "Orc4", "Orc5", "Orc6", "Pafah1b1", "Pcm1", "Pcna", "Pds5a",
"Pds5b", "Phf20", "Phlda1", "Pias4", "Pkmyt1", "Plk1", "Plk4",
"Pmf1", "Pola1", "Pola2", "Pold1", "Pold2", "Pold3", "Pold4",
"Pole", "Pole2", "Pole3", "Pole4", "Pom121", "Pot1a", "Ppme1",
"Ppp1cb", "Ppp1cc", "Ppp1r12a", "Ppp1r12b", "Ppp2ca", "Ppp2cb",
"Ppp2r1a", "Ppp2r1b", "Ppp2r2a", "Ppp2r3d", "Ppp2r5a", "Ppp2r5b",
"Ppp2r5c", "Ppp2r5d", "Ppp2r5e", "Ppp6c", "Ppp6r3", "Prim1",
"Prim2", "Prkaca", "Prkcb", "Psma1", "Psma2", "Psma3", "Psma4",
"Psma5", "Psma6", "Psma7", "Psmb1", "Psmb2", "Psmb3", "Psmb4",
"Psmb5", "Psmb6", "Psmb7", "Psmc1", "Psmc2", "Psmc3", "Psmc4",
"Psmc5", "Psmc6", "Psmd1", "Psmd11", "Psmd12", "Psmd13", "Psmd14",
"Psmd2", "Psmd3", "Psmd6", "Psmd7", "Psmd8", "Pttg1", "Rab1a",
"Rab1b", "Rab2a", "Rab8a", "Rad1", "Rad17", "Rad21", "Rad50",
"Rad51", "Rad9a", "Rad9b", "Rae1", "Ran", "Ranbp2", "Rangap1",
"Rb1", "Rbbp4", "Rbbp7", "Rbbp8", "Rbl1", "Rbl2", "Rbm39", "Rbx1",
"Rcc1", "Rcc2", "Rfc1", "Rfc2", "Rfc3", "Rfc4", "Rfc5", "Rhno1",
"Rmi1", "Rmi2", "Rnf168", "Rnf8", "Rpa1", "Rpa2", "Rpa3", "Rsf1",
"Rtel1", "Ruvbl1", "Sdccag8", "Sec13", "Seh1l", "Set", "Sfi1",
"Sfn", "Sgo1", "Sgo2a", "Shq1", "Sirt2", "Ska2", "Skp2", "Smarca5",
"Smc1a", "Smc2", "Smc3", "Smc4", "Spast", "Spc24", "Spc25", "Spdl1",
"Src", "Ssna1", "Stag1", "Stag2", "Stn1", "Sumo1", "Taok1", "Ten1",
"Terf1", "Terf2", "Terf2ip", "Tert", "Tfdp1", "Top3a", "Topbp1",
"Tpr", "Tpx2", "Trp53", "Trp53bp1", "Tuba1a", "Tuba1b", "Tuba1c",
"Tuba4a", "Tuba8", "Tubb2a", "Tubb2b", "Tubb3", "Tubb4a", "Tubb4b",
"Tubb5", "Tubb6", "Tubg1", "Tubg2", "Tubgcp2", "Tubgcp3", "Tubgcp4",
"Tubgcp5", "Tubgcp6", "Uba52", "Ubb", "Ubc", "Ube2c", "Ube2d1",
"Ube2e1", "Ube2i", "Ube2n", "Ube2s", "Ube2v2", "Uimc1", "Vps4a",
"Vrk1", "Vrk2", "Wapl", "Wee1", "Wrap53", "Wrn", "Xpo1", "Ywhab",
"Ywhae", "Ywhag", "Ywhah", "Ywhaq", "Ywhaz", "Zfp385a", "Zw10",
"Zwilch", "Zwint")
F <- ldf(fl_nmf,type = "i")$F
rownames(F) <- genes$symbol
colnames(F) <- paste0("k",1:13)
sort(colSums(F[cell_cycle,]))
#        k6       k13        k8        k5        k1       k11        k3        k4 
#  4.572620  5.697999  7.211619  8.244211  9.451072 11.196086 11.865956 13.022190 
#        k2        k9       k10        k7       k12 
# 17.587821 21.151715 22.159905 73.938631 78.021005

Closer examination of the topics and factors

TO DO:

Run de_analysis with lfc.stat = “vsnull” to “shrink” the topics for comparison with the factors in scatterplots.
Compare topics and corresponding factors in scatterplots. See temp2.R.
Show that the topics and factors are “cluster-like” in that they are strongly correlated with each other.
Identify “distinctive” genes using de_analysis for topic model and with a custom function from the singlecelljamboreeR package for the factors.
The cell cycle factors (k7, k12) are other examples where the factors are largely independent from the others.
Create annotation plots for the topic model and the EBNMF model.

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

NMF analysis of the “pancreas cytokine” data set (untreated mouse only)

Peter Carbonetto

Topic model (fastTopics)

EBNMF (flashier)

Closer examination of the topics and factors