Filtration of low-qality cells for 10x Frozen BMMCs dataset

Number of rescued cells per cluster

rescued_table <- TableOfRescuedCells(clusters_annotated_resc[c(intersect_cbs, rescued_cbs)], 
                                     rescued_cbs)
write.csv(rescued_table, paste0(kOutputFolder, "tables/rescued_cbc_bmmc1.csv"), row.names=F)
rescued_table

Seurat

Parameters are the same as in the demonstration.

presented_cbs <- intersect(colnames(r_cm_rescued), names(clusters_annotated))
seurat_cm <- r_cm_rescued[, presented_cbs]
seurat_cm <- seurat_cm[Matrix::rowSums(seurat_cm) > 200, ]

srt <- CreateSeuratObject(raw.data = seurat_cm, project = "bmmc1", display.progress=F)
srt <- NormalizeData(object = srt, normalization.method = "LogNormalize",
                     scale.factor = 10000, display.progress=F)
srt <- FindVariableGenes(object = srt, mean.function = ExpMean,
                         dispersion.function = LogVMR, x.low.cutoff = 0.0125,
                         x.high.cutoff = 3, y.cutoff = 1, do.plot=F, display.progress=F)
srt <- ScaleData(object = srt, vars.to.regress = "nUMI", display.progress=F)

Find genes to compare cell types:

srt@ident <- as.factor(clusters_annotated[colnames(srt@raw.data)])
names(srt@ident) <- colnames(srt@raw.data)
compared_clusters <- unique(srt@ident) %>% as.character()
cluster_markers <- mclapply(compared_clusters, function(i)
  mclapply(setdiff(compared_clusters, i), FindClusterMarkers, i, srt, mc.cores=4),
  mc.cores=11)

de_genes <- GetOverexpressedGenes(srt, compared_clusters, cluster_markers)

321 differentially expressed genes found.

Plots

10x CellRanger used wrong threshold:

scores <- ScorePipelineCells(holder, mit.chromosome.name='MT', 
                             predict.all=T)[names(umis_per_cell)]

smoothScatter(scores[1:6000], bandwidth=c(60, 0.015), xlab='Cell rank',
              ylab='Quality score')
abline(v=ncol(cm_10x), col='#bc2727', lty=2, lw=2.5)
abline(v=est_cell_num$expected, col='#0a6607', lty=2, lw=2.5)
arrows(x0=c(1000, 4000), y0=c(0.6, 0.7), 
       x1=c(ncol(cm_10x) - 100, est_cell_num$expected + 100), y1=c(0.42, 0.52),
       lw=2)
text(x=c(900, 4100), y=c(0.67, 0.75),
     labels=c("10x threshold", "dropEst threshold"), cex=1.3)

tested_clusts <- clusters_annotated[presented_cbs]

separation <- c(setNames(rep('rescued', length(rescued_cbs)), rescued_cbs),
                setNames(rep('real', length(intersect_cbs)), intersect_cbs))

umis_per_cb_subset <- log10(Matrix::colSums(r_cm_rescued[, names(tested_clusts)]))
tested_clusts <- tested_clusts[order(tested_clusts, -umis_per_cb_subset)]

Prepare heatmaps:

plot_df <- ExpressionMatrixToDataFrame(r_rescued$counts[names(tested_clusts), de_genes],
                                       umis_per_cb_subset, tested_clusts,
                                       filtration.type=separation)

plot_df$Cluster <- as.character(plot_df$Cluster)
plot_df$Cluster[plot_df$Cluster == "Non-dividing Pro B cells"] <- "Non-dividing\nPro B cells"

plot_df <- plot_df %>% filter(UmisPerCb < 3.4)
plot_dfs <- split(plot_df, plot_df$FiltrationType)

ggs <- lapply(plot_dfs, HeatmapAnnotGG, umi.per.cell.limits=range(plot_df$UmisPerCb),
              raster.width=3, raster.height=3, raster.dpi=100)

legend_guides <- list(HeatmapLegendGuide('Expression'),
                      HeatmapLegendGuide('Cell type', guide=guide_legend, ncol=3),
                      HeatmapLegendGuide('log10(#molecules)'))
gg_legends <- mapply(`+`, ggs$real, legend_guides, SIMPLIFY=F) %>%
  lapply(`+`, theme(legend.margin=margin(l=4, r=4, unit='pt'))) %>% lapply(get_legend)

ggs$real$heatmap <- ggs$real$heatmap + rremove('xlab') + ylab('Cells')
ggs$rescued$heatmap <- ggs$rescued$heatmap + labs(x = 'Genes', y = 'Cells')
ggs_annot <- lapply(ggs, function(gg) cowplot::plot_grid(
  plotlist=lapply(gg, `+`, theme(legend.position="none", plot.margin=margin())),
  nrow=1, rel_widths=c(1.5, 0.1, 0.1), align='h'))

gg_legends_plot <- cowplot::plot_grid(plotlist=gg_legends, nrow=3, align='v')

Compile plot parts:

gg_left <- cowplot::plot_grid(ggs_annot$real, ggs_annot$rescued, nrow=2,
                              labels=c('B', 'C'))
gg_right <- gg_tsne + theme(plot.margin=margin(l=0.1, unit='in'),
                            axis.text=element_blank(), axis.ticks=element_blank())
gg_bottom <- cowplot::plot_grid(plotlist=gg_legends[c(1, 3, 2)], ncol=3,
                                rel_widths=c(1, 1, 2.6))

gg_filtration <- cowplot::plot_grid(gg_left, gg_right, labels=c('', 'D'), ncol=2) %>%
  cowplot::plot_grid(gg_bottom, nrow=2, rel_heights=c(1, 0.27), align='v')

coords <- list(optimal=list(x_l=est_cell_num$expected,
                            x_t=4000, y_l=3.2e5, y_t=5e5),
               `10x`=list(x_l=ncol(cm_10x),
                          x_t=1300, y_l=2.5e5, y_t=4e5))

gg_repel <- function(coords, label) {
  coords$label <- label
  gg <- ggrepel::geom_label_repel(
    data=as.data.frame(coords),
    mapping=aes(x=x_l, y=y_l, label=label), nudge_x=coords$x_t - coords$x_l,
    nudge_y=coords$y_t - coords$y_l, size=5.5, segment.size=0.7, force=0,
    arrow=ggplot2::arrow(length = unit(0.03, 'npc')), fill=alpha("white", 0.7)
  )
  return(gg)
}

gg_cell_number <- PlotCellsNumberLine(Matrix::colSums(holder$cm_raw)) +
  geom_vline(aes(xintercept=coords$`10x`$x_l), linetype='dashed',
             color='#bc2727', size=1) +
  geom_vline(aes(xintercept=coords$optimal$x_l), linetype='dashed',
             color='#0a6607', size=1) +
  gg_repel(coords$`10x`, label="10x threshold") +
  gg_repel(coords$optimal, label="dropEst threshold") +
  scale_x_continuous(limits=c(0, 5750), expand=c(0, 0)) +
  theme_pdf() +
  theme(axis.ticks=element_blank(), axis.text.y=element_blank(),
        panel.grid.major.y=element_blank(), panel.grid.minor.y=element_blank())

Final figure

gg_fig <- cowplot::plot_grid(gg_cell_number + theme(plot.margin=margin(b=0.1, unit="in")),
                             gg_filtration, nrow=2, rel_heights=c(1.2, 3),
                             labels=c('A', '')) +
  theme(plot.margin=margin(1, 1, 1, 1))

ggsave(paste0(kOutputFolder, 'figures/fig_bmmc_filtration.pdf'), gg_fig, width=7.5, height=7)

gg_fig