Analysis of droplet data from Montoro et al. (using preprocessing)

Introduction and Code

I again fit 30 factors to the droplet dataset discussed in Montoro et al. This time, I follow pre-processing steps used by the authors (normalizing individual counts by counts per cell and only using “highly variable genes”). Let’s see whether it makes any difference.

Click “Code” to view the code used to pre-process the data.

library(Matrix)
library(mgcv)
library(ggplot2)
library(mixsqp)

# Functions from https://github.com/adamh-broad/single_cell_airway:

extract.field <- function(string, field = 1, delim = "_", fixed = TRUE) {
  return(strsplit(string,delim, fixed=fixed)[[1]][field])
}

info <- function(text, ...) {
  cat(sprintf(paste(Sys.time(),"INFO:", text,"\n")))
}

tpm <- function(counts, mult = 10000) {
  info("Running TPM normalisation")
  total.counts = Matrix::colSums(counts)
  scaled.counts = t(t(counts) / total.counts)
  scaled.counts * mult
}

get.variable.genes.umis <- function(umi.cts, residual.threshold = -0.25,
                                    use.spline = FALSE, batch = NULL,
                                    ret.plot = FALSE) {
  if(!is.null(batch)){
    v = as.vector(table(batch))
    total_transcripts = data.frame(as.matrix(t(Matrix.utils::aggregate.Matrix(t(umi.cts), groupings = batch, fun="sum"))))
    detection_frac = Matrix.utils::aggregate.Matrix(t(umi.cts > 0), groupings = batch, fun="sum")
    detection_frac = data.frame(as.matrix(t(detection_frac / v)))
    test_genes = rownames(detection_frac)[Matrix::rowSums(detection_frac > 0) == length(unique(batch))]
    detection_frac = detection_frac[test_genes, ]
    total_transcripts = total_transcripts[test_genes, ]
    detection_frac$gene = rownames(detection_frac)
    total_transcripts$gene = rownames(total_transcripts)
    detection_frac = melt(detection_frac, id.vars="gene")
    colnames(detection_frac) = c("gene", "batch", "alpha")
    total_transcripts = melt(total_transcripts, id.vars="gene")
    colnames(total_transcripts) = c("gene", "batch", "UMIs")
    z = cbind(total_transcripts, detection_frac)[, c("gene", "batch", "alpha", "UMIs")]
    info("Fitting logistic GLM (controlling for batch covariate)")
    model.logit = glm(data = z, formula = alpha ~ log10(UMIs) + batch, family = binomial)
    #model.logit = robust::glmRob(data = z, formula = alpha ~ log10(UMIs), family = binomial)
    info("Fitting spline quantile regression (controlling for batch covariate)")
    model.gam = quantreg::rq(data = z, formula = alpha ~ splines::ns(log10(UMIs), df=15) + batch, tau=0.8)
    #model.gam = mgcv::gam(data = z, formula = alpha ~ s(log10(UMIs)), method="REML")
  }else{
    info("Computing gene dection rates (alphas)..")
    z = data.frame(UMIs = Matrix::rowSums(umi.cts), alpha= Matrix::rowSums(umi.cts>0) / ncol(umi.cts))
    z = subset(z, UMIs > 0 | alpha > 0)
    info("Fitting GLMs..")
    model.logit = glm(data = z, formula = alpha ~ log10(UMIs), family = binomial)
    #model.logit = robust::glmRob(data = z, formula = alpha ~ log10(UMIs), family = binomial)
    model.gam = quantreg::rq(data = z, formula = alpha ~ splines::ns(log10(UMIs), df=15), tau=0.8)
    #model.gam = mgcv::gam(data = z, formula = alpha ~ s(log10(UMIs)), method="REML")
  }

  if(use.spline){
    info("use.spline is ON. Using GAM fit (blue), logit in red")
    z$predicted = predict(object = model.gam, z, type="response")
    z$predicted.alternate = predict(object = model.logit, z, type="response")
    z$residual = model.gam$residuals
  }else{
    info("use.spline is OFF. Using logit fit (blue), GAM in red")
    z$predicted = predict(model.logit, type="response") #predict(object = model.logit, z, type="response")
    z$predicted.alternate = predict(object = model.gam, z, type="response")
    z$residual = model.logit$residuals
  }
  if(is.null(batch)) {z$gene = rownames(z)}
  outliers = subset(z, residual < residual.threshold)
  g = ggplot(z, aes(x=log10(UMIs), y=alpha, label=gene)) + geom_point(color="grey50", size=0.5, stroke=0) +
    ylim(c(0,1)) + geom_line(aes(y=predicted), size=0.5, color="blue", linetype="dotted") +
    geom_line(aes(y=predicted.alternate), size=0.5, color="red", linetype="dotted") +
    geom_text(data=outliers, color="black", size=1.5, vjust=2)
  if(!is.null(batch)){g = g + facet_wrap(~batch)}
  if(!ret.plot){print(g)}
  rv = unique(unlist(lapply(rownames(outliers), extract.field, 1, delim="_")))
  if(ret.plot){return(list("var.genes"=rv, "plot"=g, "fit.data"=z, "logit"=model.logit))}
  rv
}

# Load and preprocess data as in Montoro et al.
counts <- read.table("~/Downloads/GSE103354_Trachea_droplet_UMIcounts.txt")
counts <- Matrix::Matrix(as.matrix(counts))
var.genes <- get.variable.genes.umis(counts)
counts <- tpm(counts)
counts <- log1p(counts)
counts <- counts[var.genes, ]

# Add 30 factors with rough backfits after every 5 factors.
fl <- flashier(counts, 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,
               verbose.lvl = 3)

saveRDS(fl, "~/Downloads/ppdrop30.rds")

# Refine by backfitting.
fl <- flashier(counts, flash.init = fl, backfit = "only",
               backfit.order = "dropout", backfit.maxiter = 200,
               verbose.lvl = 3)

saveRDS(fl, "~/Downloads/ppdrop30_backfit.rds")

The following code, used to produce the boxplots and tables below, is in all important respects identical to the code used in my previous analysis of this dataset.

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

fl <- readRDS("~/Downloads/ppdrop30_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] 206   0  17  45   7  40  88  34   7  43  19  22  18   3   7  42  22
#> [18]   0   1  22   0   2  37   5   2  26   7   5   3   2

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

# 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")
paper.genes <- paper.genes[paper.genes$gene %in% rownames(gene.sig), ]

# 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 .....
#>  ( 2323 GO terms found. )
#> 
#> Build GO DAG topology ..........
#>  ( 5695 GO terms and 13135 relations. )
#> 
#> Annotating nodes ...............
#>  ( 593 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.1336)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 2 (PVE: 0.0039)

Gene Ontology terms:

#> Warning in .local(object, test.stat, ...): No enrichment can pe performed -
#> there are no feasible GO terms!

Factor 3 (PVE: 0.2234)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 4 (PVE: 0.0036)

Gene Ontology terms:

Factor 5 (PVE: 0.0752)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 6 (PVE: 0.0107)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 7 (PVE: 0.0187)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 8 (PVE: 0.0022)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 9 (PVE: 0.1554)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 10 (PVE: 0.0041)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 11 (PVE: 0.0011)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 14 (PVE: 0.0016)

Gene Ontology terms:

Factor 15 (PVE: 0.0121)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 16 (PVE: 0.0359)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 17 (PVE: 0.0004)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 18 (PVE: 0.0191)

Gene Ontology terms:

#> Warning in .local(object, test.stat, ...): No enrichment can pe performed -
#> there are no feasible GO terms!

Factor 19 (PVE: 0.0007)

Gene Ontology terms:

Factor 21 (PVE: 0.0087)

Gene Ontology terms:

#> Warning in .local(object, test.stat, ...): No enrichment can pe performed -
#> there are no feasible GO terms!

Factor 22 (PVE: 0.0049)

Gene Ontology terms:

#> Warning in .local(object, test.stat, ...): No enrichment can pe performed -
#> there are no feasible GO terms!

Factor 25 (PVE: 0.0004)

Gene Ontology terms:

Factor 26 (PVE: 0.0076)

Significant genes that are also mentioned in the paper:

Gene Ontology terms:

Factor 30 (PVE: 0.0026)

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.52.2 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.20.22      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.10     R6_2.3.0          
#> [49] git2r_0.21.0       compiler_3.4.3