Filtration of low-qality cells for the 10x 8k PBMCs dataset

t-SNE with noise

plot_clusters <- setdiff(names(clusters_annotated), c(filtered_cbs, rescued_cbs))
gg <- PlotFiltrationResults(r_rescued, clusters_annotated[plot_clusters],
                            filtered.cbs=filtered_cbs, rescued.clusters=rescued_clusters, 
                            raster.width=4.28, raster.height=6)

gg <- Reduce(`+`, lapply(names(chull_list), function(n) 
  geom_polygon(data=chull_list[[n]], mapping=aes(x=V1, y=V2), alpha=0.3, 
               color=label_colors[n])), init=gg)

gg_repels <- lapply(1:3, function(i) 
  geom_label_repel(data=labels_df[i,], mapping=aes(x=V1, y=V2, label=Text), 
                   fill=alpha(labels_df$Color[i], 0.5), nudge_x=labels_df$OffsetX[i], 
                   nudge_y=labels_df$OffsetY[i], segment.size=1, size=3, point.padding=0,
                   arrow=arrow(length = unit(0.02, 'npc')), color='black', 
                   segment.color=labels_df$Color[i], fontface='bold'))

gg <- Reduce(`+`, gg_repels, init=gg)
gg_tsne_noise <- gg + theme_pdf(legend.pos=c(1, 1), show.ticks=F)
gg_tsne_noise

Error sources

bc_data <- PrepareLqCellsDataPipeline(holder, mit.chromosome.name='MT', scale=F)

real_cell_color <- '#0faf00'

gg <- function(y.max, filt.color) {
  ggplot() + scale_y_continuous(expand=c(0, 0), name='Fraction of cells', 
                                limits=c(0, y.max)) + 
    guides(fill=guide_legend(title=NULL)) + theme_pdf(legend.pos=c(1, 1)) +
    scale_fill_manual(values=as.vector(c(filt.color, real_cell_color)))
}

gg1 <- gg(0.9, label_colors['Mitochondrial']) + 
  geom_histogram(aes(x=bc_data[mitochondrial_inf$cbs, ]$MitochondrionFraction, 
                     y=..count.. / sum(..count..), fill = 'Mitochondrial'), 
                 color='black') +
  geom_histogram(aes(x=bc_data[intersect_cbs, ]$MitochondrionFraction, 
                     y=..count.. / sum(..count..), fill='Real'), alpha=0.5, 
                 color='black') +
  scale_x_continuous(limits=c(0, 1.01), expand=c(0, 0), 
                     name='Fraction of mitochondrial reads')

gg2 <- gg(0.51, label_colors['Intergenic']) + 
  geom_histogram(aes(x=bc_data[intergenic_inf$cbs, ]$IntergenicFrac, 
                     y=..count.. / sum(..count..), fill = 'Intergenic'), 
                 color='black') +
  geom_histogram(aes(x=bc_data[intersect_cbs, ]$IntergenicFrac, 
                     y=..count.. / sum(..count..), fill='Real'), 
                 alpha=0.5, color='black') +
  scale_x_continuous(limits=c(0, 1.01), expand=c(0, 0), 
                     name='Fraction of intergenic reads')
  
gg3 <- gg(0.6, label_colors['Low expression']) + 
  geom_histogram(aes(x=bc_data[low_rpu_inf$cbs, ]$UmiPerGene, 
                     y=..count.. / sum(..count..), fill = 'Low expression'), 
                 color='black') +
  geom_histogram(aes(x=bc_data[intersect_cbs, ]$UmiPerGene, 
                     y=..count.. / sum(..count..), fill='Real'), alpha=0.5, 
                 color='black') +
  scale_x_continuous(limits=c(1, 5.5), expand=c(0, 0), name='Mean #UMIs per gene')

gg_distributions <- plot_grid(gg1, gg2, gg3, nrow=3, align='v', labels=c('B', 'C', 'D'), 
                              label_x=0.025)

gg_distributions

Final figure

gg_noise_fig <- plot_grid(gg_tsne_noise, gg_distributions, ncol=2, rel_widths=c(1, 0.75), 
          labels=c('A', NULL))

gg_noise_fig

ggsave(paste0(kPlotsFolder, 'fig_pbmc8k_lq_cells.pdf'), gg_noise_fig, width=7.5, height=6)

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_pbmc8k.csv"), row.names=F)
rescued_table

Seurat

seurat_clusters <- tested_clusts
seurat_cm <- r_cm_rescued[, names(seurat_clusters)]
seurat_cm <- seurat_cm[Matrix::rowSums(seurat_cm) > 200, ]

srt <- CreateSeuratObject(raw.data = seurat_cm, project = "pbmc8k", 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)
save(srt, file="../../output/pbmc8k_srt.Rda")
# load("../../output/pbmc8k_srt.Rda")

srt@ident <- as.factor(seurat_clusters[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=5), 
  mc.cores=6)

overexpressed_genes <- GetOverexpressedGenes(srt, compared_clusters, cluster_markers, 
                                             expression.threshold=0.4)
length(overexpressed_genes)

[1] 177

Heatmaps

tested_clusts <- seurat_clusters

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)]

de_genes <- overexpressed_genes

raster_width <-  3
raster_height <- 3
raster_dpi <- 100

plot_df <- ExpressionMatrixToDataFrame(r_rescued$counts[names(tested_clusts), de_genes],
                                       umis_per_cb_subset, tested_clusts,
                                       filtration.type=separation)
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=raster_width, raster.height=raster_height, raster.dpi=raster_dpi)

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)
  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 <- plot_grid(plotlist=gg_legends, nrow=3, align='v')

gg_left <- 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 <- plot_grid(plotlist=gg_legends[c(1, 3, 2)], ncol=3, rel_widths=c(1, 1, 2.6))

gg_10x <- plot_grid(gg_left, gg_right, labels=c('', 'D'), ncol=2) %>%
  plot_grid(gg_bottom, nrow=2, rel_heights=c(1, 0.21), align='v')

Supplementary figure

gg_10x

ggsave(paste0(kPlotsFolder, 'supp_comparison_with_10x.pdf'), gg_10x, width=7.5, height=6)

Initial labels and scores

kEmbeddingType <- 'tSNE'
# kEmbeddingType <- 'largeVis'

cbs_full <- names(scores)[1:20000]
r_cm_full <- holder$cm_raw[, cbs_full]
r_full <- GetPagoda(r_cm_full, n.cores=1, embeding.type=kEmbeddingType)

20000 cells, 14590 genes; normalizing ... using plain model winsorizing ... log scale ... done.
calculating variance fit ... using gam 2460 overdispersed genes ... 2460 persisting ... done.
running PCA using 1000 OD genes .... done
calculating distance ... pearson ...running tSNE using 1 cores:
 - point 10000 of 20000
 - point 20000 of 20000

rast_width <- 7.5 / 2
rast_height <- 5
rast_dpi <- 150
initial_labels <- dropestr:::EstimateCellsQuality(Matrix::colSums(r_cm_full))
gg1 <- PlotPagodaEmbeding(r_full, clusters=initial_labels, mark.clusters=F, show.legend=T, alpha=0.6, embeding.type=kEmbeddingType,
                          size=0.5, raster=T, raster.width=rast_width, raster.height=rast_height, raster.dpi=rast_dpi) +
  scale_color_manual(values=scales::hue_pal(l=55)(3)[c(2, 1, 3)], name='Initial label') +
  guides(color=guide_legend(override.aes=list(size=1.2))) +
  theme_pdf(legend.pos=c(1, 1), show.ticks=F)

gg2 <- PlotPagodaEmbeding(r_full, colors=scores, mark.clusters=F, show.legend=T, alpha=0.6, embeding.type=kEmbeddingType, 
                          size=0.5, raster=T, raster.width=rast_width, raster.height=rast_height, raster.dpi=rast_dpi) +
  scale_color_distiller(palette='RdYlBu', direction=1, limits=c(0, 1), name='Score') +
  theme_pdf(legend.pos=c(1, 1), show.ticks=F)

gg_fig <- plot_grid(gg1, gg2 + rremove("ylab"), labels='AUTO', label_x=c(0, -0.05))

gg_fig

ggsave(paste0(kPlotsFolder, 'supp_pbmc8k_initial_labeling.pdf'), gg_fig, width=7.5, height=5)