Analysis of totals counts

Setup

library("cowplot")
library("dplyr")
library("ggplot2")
library("tidyr")
theme_set(theme_cowplot())
source("../code/functions.R")
library("Biobase")

Import data.

eset <- readRDS("../data/eset.rds")
anno <- pData(eset)

Remove samples with bad cell number or TRA-1-60.

eset_qc <- eset[, anno$cell_number == 1 & anno$tra1.60]
anno_qc <- pData(eset_qc)

Extract the columns relevant for this analysis.

totals <- anno_qc %>% select(experiment:well, raw:mol_hs)

What percentage of reads are mapped to the genome?

totals_long <- totals %>%
  gather(key = "type", value = "count", raw:mol_hs) %>%
  filter(type %in% c("raw", "umi", "mapped", "molecules")) %>%
  mutate(type = factor(type, levels = c("raw", "umi", "mapped", "molecules"),
                       labels = c("Raw reads", "Reads with UMI",
                                  "Reads mapped\nto genome",
                                  "Molecules mapped\nto genome")))
head(totals_long)

  experiment well      type   count
1   02192018  A01 Raw reads 4755105
2   02192018  A02 Raw reads 1955271
3   02192018  A03 Raw reads 4876011
4   02192018  A04 Raw reads 6770341
5   02192018  A05 Raw reads 8782029
6   02192018  A06 Raw reads 7043643

ggplot(totals_long, aes(x = type, y = count)) +
  geom_boxplot() +
  labs(x = "Classification", y = "Number of sequences",
       title = "Processing of raw reads")

How does the number of mapped reads vary by C1 chip?

mapped <- totals_long %>%
  filter(type == "Reads mapped\nto genome")
ggplot(mapped, aes(x = experiment, y = count)) +
  geom_boxplot() +
  labs(x = "C1 chip", y = "Number of mapped reads",
       title = "Mapped reads per C1 chip") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

mapped_sum <- mapped %>%
  group_by(experiment) %>%
  summarize(total = sum(count))
ggplot(mapped_sum, aes(x = experiment, y = total)) +
  geom_bar(stat= "identity") +
  labs(x = "C1 chip", y = "Total number of mapped reads",
       title = "Total mapped reads per C1 chip") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

Sorted by total number of mapped reads (in millions).

mapped_sum %>%
  arrange(total) %>%
  mutate(total = total / 10^6) %>%
  knitr::kable()

What is the conversion of reads to molecules?

convert_all <- ggplot(totals, aes(x = mapped, y = molecules)) +
  geom_point() +
  labs(x = "Number of reads mapped to any genome",
       y = "Number of molecules mapped to any genome",
       title = "Conversion of reads to molecules")
convert_hs <- ggplot(totals, aes(x = reads_hs, y = mol_hs)) +
  geom_point() +
  labs(x = "Number of reads mapped to human genome",
       y = "Number of molecules mapped to human genome",
       title = "Conversion of reads to molecules for human")
convert_dm <- ggplot(totals, aes(x = reads_dm, y = mol_dm)) +
  geom_point() +
  labs(x = "Number of reads mapped to fly genome",
       y = "Number of molecules mapped to fly genome",
       title = "Conversion of reads to molecules for fly")
convert_ce <- ggplot(totals, aes(x = reads_ce, y = mol_ce)) +
  geom_point() +
  labs(x = "Number of reads mapped to worm genome",
       y = "Number of molecules mapped to worm genome",
       title = "Conversion of reads to molecules for worm")

plot_grid(convert_all, convert_hs, convert_dm, convert_ce,
          labels = letters[1:4])

Drosophila spike-in

p_reads_dm <- ggplot(anno_qc, aes(x = as.factor(fly), y = reads_dm / mapped * 100)) +
  geom_boxplot() +
  labs(x = "Fly spike-in concentration",
       y = "Percentage of mapped reads that mapped to fly genome")
p_mol_dm <- ggplot(anno_qc, aes(x = as.factor(fly), y = mol_dm / molecules * 100)) +
  geom_boxplot() +
  labs(x = "Fly spike-in concentration",
       y = "Percentage of mapped molecules that mapped to fly genome")
plot_grid(p_reads_dm, p_mol_dm, labels = letters[1:2])

p_reads_dm_counts <- ggplot(anno_qc, aes(x = as.factor(fly), y = reads_dm)) +
  geom_boxplot() +
  labs(x = "Fly spike-in concentration",
       y = "Number of reads that mapped to fly genome")
p_mol_dm_counts <- ggplot(anno_qc, aes(x = as.factor(fly), y = mol_dm)) +
  geom_boxplot() +
  labs(x = "Fly spike-in concentration",
       y = "Number of molecules that mapped to fly genome")
plot_grid(p_reads_dm_counts, p_mol_dm_counts, labels = letters[1:2])

totals_percent <- anno_qc %>%
  mutate(sample = paste(experiment, well, sep = "-"),
         sample = factor(sample, levels = sample[order(reads_dm / raw)])) %>%
  arrange(reads_dm / raw) %>%
  group_by(experiment, sample, fly, worm, ERCC) %>%
  summarize(p_reads_ce = reads_ce / raw * 100,
            p_reads_dm = reads_dm / raw * 100,
            p_reads_ercc = reads_ercc / raw * 100,
            p_reads_hs = reads_hs / raw * 100,
            p_unmapped = unmapped / raw * 100,
            p_no_umi = (raw - unmapped - mapped) / raw * 100) %>%
  gather(key = "type", value = "percent", p_reads_ce:p_no_umi) %>%
  ungroup() %>%
  mutate(type = factor(type, levels = c("p_no_umi", "p_unmapped", "p_reads_hs",
                                        "p_reads_ercc", "p_reads_ce", "p_reads_dm"),
                       labels = c("No UMI", "Unmapped", "H. sapiens", "ERCC",
                                  "C. elegans", "D. melanogaster")),
         fly = factor(fly, levels = c(5000, 50000), labels = c("5 pg fly spike-on",
                                                               "50 pg fly spike-on")))

ggplot(totals_percent, aes(x = sample, y = percent, fill = type)) +
  geom_col() +
  facet_wrap(~fly) +
  labs(x = "Sample", y = "Percent", title = "Read classification per sample") +
  theme(legend.title = element_blank())

C. elegans spike-in

p_reads_ce <- ggplot(anno_qc, aes(x = as.factor(worm), y = reads_ce / mapped * 100)) +
  geom_boxplot() +
  labs(x = "Worm spike-in concentration",
       y = "Percentage of mapped reads that mapped to worm genome")
p_mol_ce <- ggplot(anno_qc, aes(x = as.factor(worm), y = mol_ce / molecules * 100)) +
  geom_boxplot() +
  labs(x = "Worm spike-in concentration",
       y = "Percentage of mapped molecules that mapped to worm genome")
plot_grid(p_reads_ce, p_mol_ce, labels = letters[1:2])

p_reads_ce_counts <- ggplot(anno_qc, aes(x = as.factor(worm), y = reads_ce)) +
  geom_boxplot() +
  labs(x = "Worm spike-in concentration",
       y = "Number of mapped reads that mapped to worm genome")
p_mol_ce_counts <- ggplot(anno_qc, aes(x = as.factor(worm), y = mol_ce)) +
  geom_boxplot() +
  labs(x = "Worm spike-in concentration",
       y = "Number of mapped molecules that mapped to worm genome")
plot_grid(p_reads_ce_counts, p_mol_ce_counts, labels = letters[1:2])

ERCC spike-in

p_reads_ercc <- ggplot(anno_qc, aes(x = as.factor(ERCC), y = reads_ercc / mapped * 100)) +
  geom_boxplot() +
  labs(x = "ERCC spike-in concentration",
       y = "Percentage of mapped reads that mapped to ERCCs")
p_mol_ercc <- ggplot(anno_qc, aes(x = as.factor(ERCC), y = mol_ercc / molecules * 100)) +
  geom_boxplot() +
  labs(x = "ERCC spike-in concentration",
       y = "Percentage of mapped molecules that mapped to ERCCs")
plot_grid(p_reads_ercc, p_mol_ercc, labels = letters[1:2])

p_reads_ercc_counts <- ggplot(anno_qc, aes(x = as.factor(ERCC), y = reads_ercc)) +
  geom_boxplot() +
  labs(x = "ERCC spike-in concentration",
       y = "Number of reads that mapped to ERCCs")
p_mol_ercc_counts <- ggplot(anno_qc, aes(x = as.factor(ERCC), y = mol_ercc)) +
  geom_boxplot() +
  labs(x = "ERCC spike-in concentration",
       y = "Number of molecules that mapped to ERCCs")
plot_grid(p_reads_ercc_counts, p_mol_ercc_counts, labels = letters[1:2])

Session information

sessionInfo()

R version 3.4.1 (2017-06-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS: /project2/gilad/jdblischak/miniconda3/envs/scqtl/lib/R/lib/libRblas.so
LAPACK: /project2/gilad/jdblischak/miniconda3/envs/scqtl/lib/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  methods   stats     graphics  grDevices utils     datasets 
[8] base     

other attached packages:
[1] bindrcpp_0.2        Biobase_2.38.0      BiocGenerics_0.24.0
[4] tidyr_0.7.1         dplyr_0.7.4         cowplot_0.9.1      
[7] ggplot2_2.2.1      

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.13     compiler_3.4.1   git2r_0.19.0     plyr_1.8.4      
 [5] highr_0.6        bindr_0.1        tools_3.4.1      digest_0.6.12   
 [9] evaluate_0.10.1  tibble_1.3.3     gtable_0.2.0     pkgconfig_2.0.1 
[13] rlang_0.1.2      yaml_2.1.14      stringr_1.2.0    knitr_1.20      
[17] rprojroot_1.2    grid_3.4.1       tidyselect_0.2.3 glue_1.1.1      
[21] R6_2.2.0         rmarkdown_1.8    purrr_0.2.2      magrittr_1.5    
[25] backports_1.0.5  scales_0.5.0     htmltools_0.3.6  assertthat_0.1  
[29] colorspace_1.3-2 labeling_0.3     stringi_1.1.2    lazyeval_0.2.0  
[33] munsell_0.4.3