Drop-seq human/mouse mixture analysis

Load data

Here bam file was filtered by realigning it with kallisto 0.43 separately on mouse and human genome. Only reads, which were aligned only on one of them were used in dropEst.

# holder <- readRDS('../../data/dropest/dropseq/thousand/est_2018_01_26_tophat/thousand.rds')
# saveRDS(holder$cm_raw, '../../data/dropest/dropseq/thousand/est_2018_01_26_tophat/cm.rds')
cm <- readRDS('../../data/dropest/dropseq/thousand/est_2018_01_26_tophat/cm.rds')
kPlotDir <- '../../output/figures/'

cell_number <- 1100

gene_species <- ifelse(substr(rownames(cm), 1, 2) == "HU", 'Human', 'Mouse') %>% 
  as.factor()
umi_by_species <- lapply(levels(gene_species), function(l) cm[gene_species == l,] %>% 
                           Matrix::colSums()) %>% as.data.frame() %>% 
  `colnames<-`(levels(gene_species)) %>% tibble::rownames_to_column('CB') %>% 
  as_tibble() %>% 
  mutate(Total = Human + Mouse, Organism=ifelse(Human > Mouse, "Human", "Mouse"), 
         IsReal=rank(Total) >= length(Total) - cell_number) %>%
  filter(Total > 20)

mit_genes <- rownames(cm)[grep("_MT:.+", rownames(cm))]
umi_by_species$MitochondrionFraction <- GetGenesetFraction(cm, mit_genes)[umi_by_species$CB]

umi_by_species$Type <- ifelse(umi_by_species$IsReal, umi_by_species$Organism, "Background")
umi_by_species$Type[umi_by_species$Mouse > 5e3 & umi_by_species$Human > 5e3] <- 'Dublets'

Common view

gg_template <- ggplot(umi_by_species, aes(x=Mouse, y=Human)) + 
  geom_abline(aes(slope=1, intercept=0), linetype='dashed', alpha=0.5) +
  scale_x_log10(limits=c(1, 2e5), name="#Mouse molecules") + 
  scale_y_log10(name="#Human molecules") + annotation_logticks() +
  theme_pdf(legend.pos=c(0.97, 0.05)) + theme(legend.margin=margin(l=3, r=3, unit="pt"))

gg_left <- gg_template + 
  geom_point(aes(color=IsReal), size=0.1, alpha=0.15) +
  guides(color=guide_legend(override.aes=list(size=1.5, alpha=1)))

gg_right <- gg_template + 
  geom_point(aes(color=MitochondrionFraction), size=0.1, alpha=0.15) +
  scale_color_gradientn(colours=c("#1200ba", "#347fff", "#cc4000", "#ff3333"), 
                        values=scales::rescale(c(0, 0.1, 0.3, 0.8))) +
  guides(color=guide_colorbar(direction="horizontal", title.position="top", 
                              title="Mitochondrial\nfraction", 
                              barwidth=unit(1.0, units="in")))

cowplot::plot_grid(gg_left, gg_right)

ggplot(umi_by_species) + 
  geom_point(aes(x=Total, y=pmin(Human, Mouse) / Total, color=Organism), size=0.1, 
             alpha=0.1) +
  scale_x_log10(name='#Real UMIs', limits=c(10, 2e5)) + annotation_logticks() + 
  ylab('Fraction of mixed UMIs') +
  guides(color=guide_legend(override.aes=list(size=1.5, alpha=1))) +
  theme_pdf(legend.pos=c(1, 1))

Check for constant background

Background cells have constant fraction of mouse and human reads. Though this fraction is not determined by the corresponding ratio in real cells:

mouse_frac <- umi_by_species %>% filter(IsReal) %>% 
  summarise(Mouse=sum(Mouse[Organism == 'Mouse']), Human=sum(Human[Organism == 'Human']), 
            MF=Mouse / (Mouse + Human)) %>% .$MF

ggplot(umi_by_species) + 
  geom_histogram(aes(x=Mouse / Total, y=..density.., fill=IsReal), binwidth=0.005, position="identity") + 
  geom_vline(xintercept=mouse_frac) +
  xlab("Fraction of mouse reads") +
  theme_pdf(legend.pos=c(1, 1))

Distribution of total number of molecules by background cells:

gg <- ggplot(umi_by_species %>% filter(!IsReal)) +
  geom_histogram(aes(x=Total), bins=100) +
  scale_x_continuous(limits=c(0, 1000), expand=c(0, 0), name="Total #UMIs") +
  scale_y_continuous(limits=c(0,4500), expand=c(0, 0), name="#Cells") +
  theme_pdf()

gg

Though, it looks better on logscale:

gg + 
  scale_x_log10(limits=c(19, 3000), expand=c(0, 0), name="Total #UMIs") +
  annotation_logticks(sides="b")

Figure

arrows_df <- umi_by_species %>% group_by(Type) %>% 
  summarise(MouseEnd=median(Mouse), HumanEnd=median(Human)) %>%
  mutate(Mouse=c(1e1, 6e4, 2e1, 6e4), 
         Human=c(1e3, 5e3, 7e4, 7e1))

fig_width <- 3.7
fig_height <- 4.4

gg_fig <- gg_template + 
  geom_point_rast(aes(color=MitochondrionFraction), size=0.2, alpha=0.15, 
                  width=fig_width, height=fig_height, dpi=200) +
  scale_color_gradientn(colours=c("#1200ba", "#347fff", "#cc4000", "#ff3333"), 
                        values=scales::rescale(c(0, 0.1, 0.3, 0.8)), 
                        breaks=seq(0, 1.0, 0.2)) +
  guides(color=guide_colorbar(direction="horizontal", title.position="top", 
                              title="Mitochondrial\nfraction", 
                              barwidth=unit(1.0, units="in"))) + 
  stat_ellipse(aes(group=Type), level=0.99) +
  geom_segment(aes(xend=MouseEnd, yend=HumanEnd, group=Type),  data=arrows_df, 
               arrow=arrow(length = unit(0.03, "npc"))) +
  geom_label(aes(label=Type),  data=arrows_df, fill=alpha('white', 1)) +
  theme(plot.margin=margin(1, 1, 1, 1))

try(invisible(dev.off()), silent=T)
ggsave(paste0(kPlotDir, 'fig_human_mouse.pdf'), gg_fig, 
       width=fig_width, height=fig_height)

gg_fig

Session information