Analysis of scRNA data from Montoro et al.

Introduction and Code

I fit 30 factors to the drop-seq dataset discussed in Montoro et al. I removed all genes with zero total counts and did a log1p transform of the counts. I used nonnegative priors for gene loadings and normal-mixture priors for cell loadings.

# Add 30 factors with rough backfits after every 5 factors.
fl <- flashier(data, greedy.Kmax = 30, var.type = 1, 
               prior.type = c("nonnegative", "normal.mix"), 
               ash.param = list(optmethod = "mixSQP"), 
               backfit.every = 5, final.backfit = TRUE, 
               backfit.order = "montaigne", warmstart.backfits = FALSE)
# Refine by backfitting.
fl <- flashier(data, flash.init = fl, backfit = "only",
               backfit.order = "dropout", backfit.maxiter = 200)

The following code is used to produce the boxplots and tables below.

suppressMessages({
  library(ggplot2)
  library(topGO)
  library(org.Mm.eg.db)
})
#> Warning: package 'ggplot2' was built under R version 3.4.4

fl <- readRDS("~/Downloads/DropSeq30_backfit.rds")
# Remove large data object to free up memory.
fl$fit$Y <- NULL

# Data frame containing cell type and loadings for each factor.
PSdf <- data.frame(fl$loadings$normalized.loadings[[2]])
cell.names <- rownames(fl$loadings$normalized.loadings[[2]])
cell.types <- as.factor(sapply(strsplit(cell.names, "_"), `[`, 3))
levels(cell.types)
#> [1] "Basal"          "Ciliated"       "Club"           "Goblet"        
#> [5] "Ionocyte"       "Neuroendocrine" "Tuft"
levels(cell.types) <- c("Bas", "Cil", "Clb", "Gob", "Ion", "Nec", "Tft")
PSdf$cell.type <- cell.types

# Need to select signficant genes for topGO.
#   Scale gene loadings by scale constant * maximum cell loading.
s <- fl$loadings$scale.constant * 
  apply(abs(fl$loadings$normalized.loadings[[2]]), 2, max)
gene.loadings <- fl$loadings$normalized.loadings[[1]]
gene.loadings <- gene.loadings * rep(s, each = nrow(gene.loadings))
#   Get a pseudo-t statistic by dividing by the residual SE.
gene.t <- gene.loadings * sqrt(fl$fit$tau)
#   Convert to a p-value.
gene.p <- 2 * (1 - pnorm(abs(gene.t)))
#   Select significant genes using Benjamini-Hochberg.
BH <- function(k, alpha = 0.01) {
  pvals <- gene.p[, k]
  selected <- rep(0, length(pvals))
  names(selected) <- names(pvals)
  
  n <- length(pvals)
  sorted.pvals <- sort(pvals)
  BH <- sorted.pvals < alpha / (n - 0:(n - 1))
  cutoff <- min(which(!BH))
  selected[pvals < sorted.pvals[cutoff]] <- 1
  
  return(selected)
}
gene.sig <- sapply(1:ncol(gene.p), BH)

# How many significant genes per factor are there?
colSums(gene.sig)
#>  [1] 286 512  98  74 163 183 131 198 168 541 364 365 104 558 535 499  31
#> [18] 504 406 436   3 183   0  96 448 326 361 404  45 156

# Omit factors with PVE below a certain (empirically determined) threshold.
kset <- (1:ncol(gene.p))[fl$pve > 0.0001 & colSums(gene.sig) > 9]

# I also want to see which factors contain genes mentioned in the paper.
paper.genes <- data.frame(
  gene = c("Nfia", "Ascl1", "Ascl2", "Ascl3", "Foxq1", "Cdhr3", "Rgs13",
           "Muc5b", "Notch2", "Il13ra1", "Krt4", "Krt13", "Krt8", "Ecm1",
           "S100a11", "Cldn3", "Lgals3", "Anxa1", "Il25", "Tslp", "Alox5ap",
           "Ptprc", "Pou2f3", "Gfi1b", "Spib", "Sox9", "Gp2", "Tff1",
           "Tff2", "Lman1l", "P2rx4", "Muc5ac", "Dcpp1", "Dcpp2", "Dcpp3",
           "Atp6v1c2", "Atp6v0d2", "Cftr", "Aqp3", "Krt5", "Dapl1", "Hspa1a",
           "Trp63", "Scgb1a1", "Krt15", "Cyp2f2", "Lypd2", "Cbr2", "Foxj1",
           "Ccdc153", "Ccdc113", "Mlf1", "Lztfl1", "Chga", "Dclk1"),
  type = c("Club", "NEC", "Tuft", "Ion", "Goblet", "Cil", "Tuft",
           "Club", "Club", "Club", "Hill", "Hill", "Basal", "Hill",
           "Hill", "Hill", "Hill", "Hill", "Tuft", "Tuft", "Tuft2",
           "Tuft2", "Tuft1", "Tuft2", "Tuft2", "Tuft2", "Goblet", "Gob1",
           "Gob1", "Gob1", "Gob1", "Gob1", "Gob2", "Gob2", "Gob2",
           "Ion", "Ion", "Ion", "BM", "BM", "BM", "BM",
           "BM", "ClbM", "ClbM", "ClbM", "ClbM", "ClbM", "CilM",
           "CilM", "CilM", "CilM", "CilM", "NECM", "TftM"))
paper.genes$gene <- as.character(paper.genes$gene)
levels(paper.genes$type) <- c("Basal", "Basal (Marker)", "Ciliated", 
                              "Ciliated (Marker)", "Club (Marker)", "Club",
                              "Goblet-1", "Goblet-2", "Goblet", "Hillock",
                              "Ionocyte", "NEC", "NEC (Marker)",
                              "Tuft (Marker)", "Tuft", "Tuft-1", "Tuft-2")

# Set up topGOdata object.
GO.list <- as.factor(gene.sig[, 1])
GOdata <- new("topGOdata", ontology = "BP", allGenes = GO.list,
              annot = annFUN.org, mapping = "org.Mm.eg", ID = "symbol")
#> 
#> Building most specific GOs .....
#>  ( 11286 GO terms found. )
#> 
#> Build GO DAG topology ..........
#>  ( 15215 GO terms and 35696 relations. )
#> 
#> Annotating nodes ...............
#>  ( 15591 genes annotated to the GO terms. )

# Loop over kset.
for (k in kset) {
  cat(paste0("## Factor ", k, " (PVE: ",
             formatC(fl$pve[k], format = "f", digits = 4), ")\n"))
  fctr <- paste0("X", k)
  plot(ggplot(PSdf, aes_string(x = "cell.type", y = fctr)) +
         geom_boxplot(outlier.shape = NA, color = "red") +
         geom_jitter(position = position_jitter(0.2), cex = 0.1) +
         labs(x = "Cell Type", y = paste("Factor", k, "value")))
  cat("\n")
  
  mentions <- which(gene.sig[paper.genes$gene, k] == 1)
  if (length(mentions) > 0) {
    mentions.df <- paper.genes[mentions, ]
    mentions.df$t.val <- gene.t[names(mentions), k]
    cat("Significant genes that are also mentioned in the paper:\n")
    print(knitr::kable(mentions.df))
    cat("\n")
  }
  
  cat("Gene Ontology terms:\n")
  GOdata@allScores <- as.factor(gene.sig[, k])
  result <- suppressMessages(
    runTest(GOdata, algorithm = "classic", statistic = "fisher")
  )
  allRes <- GenTable(GOdata, classic = result, topNodes = 10)
  print(knitr::kable(allRes))
  cat("\n")
}

Factor 1 (PVE: 0.5246)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 2 (PVE: 0.0189)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 3 (PVE: 0.0387)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 4 (PVE: 0.0581)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 5 (PVE: 0.0016)

Gene Ontology terms:

Factor 6 (PVE: 0.0076)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 7 (PVE: 0.0185)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 8 (PVE: 0.0028)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 9 (PVE: 0.0012)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 13 (PVE: 0.0027)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 17 (PVE: 0.0003)

Gene Ontology terms:

Factor 24 (PVE: 0.0006)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 29 (PVE: 0.0001)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Session information

sessionInfo()
#> R version 3.4.3 (2017-11-30)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS High Sierra 10.13.6
#> 
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] org.Mm.eg.db_3.5.0   topGO_2.30.1         SparseM_1.77        
#>  [4] GO.db_3.5.0          AnnotationDbi_1.40.0 IRanges_2.12.0      
#>  [7] S4Vectors_0.16.0     Biobase_2.38.0       graph_1.56.0        
#> [10] BiocGenerics_0.24.0  ggplot2_3.1.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] xfun_0.4           lattice_0.20-35    colorspace_1.3-2  
#>  [4] htmltools_0.3.6    yaml_2.2.0         blob_1.1.0        
#>  [7] rlang_0.3.0.1      R.oo_1.21.0        pillar_1.2.1      
#> [10] glue_1.3.0         withr_2.1.2.9000   DBI_0.7           
#> [13] R.utils_2.6.0      bit64_0.9-7        bindrcpp_0.2      
#> [16] matrixStats_0.54.0 bindr_0.1          plyr_1.8.4        
#> [19] stringr_1.3.1      munsell_0.5.0      gtable_0.2.0      
#> [22] workflowr_1.0.1    flashier_0.1.0     R.methodsS3_1.7.1 
#> [25] evaluate_0.12      memoise_1.1.0      labeling_0.3      
#> [28] knitr_1.21.6       highr_0.7          Rcpp_1.0.0        
#> [31] scales_1.0.0       backports_1.1.2    bit_1.1-12        
#> [34] digest_0.6.18      stringi_1.2.4      dplyr_0.7.4       
#> [37] grid_3.4.3         rprojroot_1.3-2    tools_3.4.3       
#> [40] magrittr_1.5       lazyeval_0.2.1     tibble_1.4.2      
#> [43] RSQLite_2.0        whisker_0.3-2      pkgconfig_2.0.1   
#> [46] assertthat_0.2.0   rmarkdown_1.11     R6_2.3.0          
#> [49] git2r_0.21.0       compiler_3.4.3