Comparison of Quantification methods for CD4 RNA-Seq Dataset

Preliminary Setup

First we load the necessary libraries.

library(here)
library(stringr)
library(magrittr)
library(openxlsx)
library(SummarizedExperiment)
library(dplyr)
library(edgeR)
library(limma)
library(sva)
library(ggplot2)
library(scales)
library(GGally)
library(ggalt)
library(reshape2)
library(assertthat)
library(ggfortify)
library(broom)

library(doParallel)
ncores <- getOption("mc.cores", default=parallel::detectCores(logical = FALSE))
registerDoParallel(cores=ncores)
library(BiocParallel)
register(DoparParam())

source(here("scripts/utilities.R"))

Data Loading and Preprocessing

Now we’ll load the RNA-seq data set from an RDS file containing a SummarizedExperiment object, and modify it to use the sample names as column names.

datasets <- c("hisat2_grch38_snp_tran_ensembl.85",
              "hisat2_grch38_snp_tran_knownGene",
              "kallisto_hg38.analysisSet_ensembl.85",
              "kallisto_hg38.analysisSet_knownGene",
              "salmon_hg38.analysisSet_ensembl.85",
              "salmon_hg38.analysisSet_knownGene",
              "shoal_hg38.analysisSet_ensembl.85",
              "shoal_hg38.analysisSet_knownGene",
              "star_hg38.analysisSet_ensembl.85",
              "star_hg38.analysisSet_knownGene")
datasetmeta <- datasets %>% str_match("^([^_]+)_(.+)_([^_]+)$") %>%
    as_data_frame %>% setNames(c("dataset", "program", "genome", "transcriptome"))
assert_that(all(table(datasetmeta$program) == 2))

[1] TRUE

assert_that(all(table(datasetmeta$transcriptome) == 5))

[1] TRUE

sexp.files <- setNames(here("saved_data",
                            sprintf("SummarizedExperiment_rnaseq_%s.RDS", datasets)),
                       datasets)
sexps <- bplapply(sexp.files, readRDS)
for (dset in names(sexps)) {
    sexps[[dset]] %<>% set_colnames(.$SampleName)
}

We extract the sample metadata from the SummarizedExperiment (we only do this once since it’s the same for all datasets).

sample.table <- colData(sexps[[1]]) %>%
    as.data.frame %>% autoFactorize %>%
    rename(batch=technical_batch) %>%
    mutate(time_point=factor(days_after_activation) %>% `levels<-`(sprintf("D%s", levels(.))),
           group=interaction(cell_type, time_point, sep=""))

Next we extract the count matrix from each SummarizedExperiment. This is made more complicated than usual by the fact that half of the samples were sequenced with a different protocol than the other half, and the two protocols produce reads with opposite strand orientations. Hence, we need the sense counts for half of the samples and the antisense counts for the other half. The appropriate strand for each sample is documented in the libType column of the sample metadata, using the library type abbreviations established by Salmon.

For SummarizedExperiments generated using tximport, this step is skipped, since the quantification tool has already been told which strand to use and only provides counts for that strand.

libtype.assayNames <- c(SF="sense.counts", SR="antisense.counts")
for (dset in seq_along(sexps)) {
    if (all(libtype.assayNames %in% assayNames(sexps[[dset]]))) {
        sample.table %<>% mutate(count_type=libtype.assayNames[libType])
        assay(sexps[[dset]], "unstranded.counts") <- assay(sexps[[dset]], "counts")
        assay(sexps[[dset]], "counts") <- lapply(seq_len(nrow(sample.table)), function(i) {
            assay(sexps[[dset]], sample.table[i,]$count_type %>% as.character)[,i]
        }) %>% do.call(what=cbind)
    }
}

As a sanity check, we make sure that we selected the strand sense with the higher total count for each sample.

for (dset in seq_along(sexps)) {
    if (all(libtype.assayNames %in% assayNames(sexps[[dset]]))) {
        total.counts <- sexps[[dset]] %>% assays %>% sapply(colSums) %>% data.frame %>%
            mutate(SampleName=row.names(.)) %>%
            inner_join(sample.table, by="SampleName")
        total.counts %$% invisible(assert_that(all(counts == pmax(sense.counts, antisense.counts))))
    }
}

Now we create a DGEList from each count matrix, and perform the initial scaling normalization. In addition, if there is a length assay, we also use that to derive an offset matrix that corrects for sample-specific biases detected by Salmon or shoal.

dges <- bplapply(sexps, function(sexp) {
    ## Extract gene metadata and colapse lists
    all.gene.meta <- mcols(sexp) %>% as.data.frame
    # Convert list columns to character vectors
    all.gene.meta[] %<>% lapply(function(x) if (is.list(x)) sapply(x, str_c, collapse=",") else x)
    dge <- DGEList(counts=assay(sexp, "counts"))
    dge$genes <- all.gene.meta
    rownames(dge$genes) <- rownames(dge)
    dge %<>% calcNormFactors
    if ("length" %in% assayNames(sexp)) {
        normMat <- assay(sexp, "length") %>% divide_by(exp(rowMeans(log(.))))
        normCounts <- dge$counts/normMat
        lib.offsets <- log(calcNormFactors(normCounts)) + log(colSums(normCounts))
        dge$offset <- t(t(log(normMat)) + lib.offsets)
    }
    dge
})

Finally, we compute the normalized logCPM matrix for each dataset, and convert it into a tidy format. We also compute the normalized logCPM without taking into account the offsets derived from gene length normalization.

logcpm.matrices <- bplapply(dges, cpmWithOffset, log=TRUE)
logcpm.matrices.nooffset <- bplapply(dges, cpm, log=TRUE)
logcpm <- bplapply(dges, function(dge) {
    cpm.mat <- cpmWithOffset(dge, log=TRUE)
    cpm.df <- melt(cpm.mat, varnames=c("GeneID", "Sample"), value.name = "logCPM") %>%
        mutate(GeneID=as.character(GeneID), Sample=as.character(Sample))
    genes <- dge$genes %>% mutate(GeneID=rownames(.))
    inner_join(genes[c("GeneID", "ENTREZID", "ENSEMBL")], cpm.df, by="GeneID") %>%
        select(-GeneID)
})
logcpm.nooffset <- bplapply(dges, function(dge) {
    cpm.mat <- cpm(dge, log=TRUE)
    cpm.df <- melt(cpm.mat, varnames=c("GeneID", "Sample"), value.name = "logCPM") %>%
        mutate(GeneID=as.character(GeneID), Sample=as.character(Sample))
    genes <- dge$genes %>% mutate(GeneID=rownames(.))
    inner_join(genes[c("GeneID", "ENTREZID", "ENSEMBL")], cpm.df, by="GeneID") %>%
        select(-GeneID)
})

Comparative Analysis

Now we can begin comparing the different methods of quantification. First, we will compare Ensembl vs UCSC annotations for each quantification method.

Ensembl vs UCSC

First, we format the logCPM tables into tidy format, and merge the tables of both datasets (Ensembl and UCSC knownGene) for each alignment/quantification program. The merging is performed only for genes with a one-to-one relationship between Ensembl and Entrez gene IDs.

logcpm.tables <- datasetmeta %>%
    arrange(transcriptome) %>%
    group_by(program) %>% do({
        df <- .
        assert_that(nrow(df) == 2)
        dsetA <- df$dataset[1]
        dsetB <- df$dataset[2]
        txA <- df$transcriptome[1]
        txB <- df$transcriptome[2]
        logcpm_A <- logcpm[[dsetA]] %>% 
            rename(logCPM_A=logCPM) %>%
            filter(lengths(ENTREZID) == 1, lengths(ENSEMBL) == 1) %>%
            mutate(ENTREZID=unlist(ENTREZID),
                   ENSEMBL=unlist(ENSEMBL))
        logcpm_B <- logcpm[[dsetB]] %>%
            rename(logCPM_B=logCPM) %>%
            filter(lengths(ENTREZID) == 1, lengths(ENSEMBL) == 1) %>%
            mutate(ENTREZID=unlist(ENTREZID),
                   ENSEMBL=unlist(ENSEMBL))
        logcpm_merged <- inner_join(logcpm_A, logcpm_B)
        data_frame(program=df$program[1], 
                   genome=df$genome[1],
                   datasetA=dsetA, datasetB=dsetB,
                   transcriptomeA=txA, transcriptomeB=txB,
                   logcpm=list(logcpm_merged))
    })
# Print the correlations
logcpm.tables %>% group_by_(names(.)[1:5]) %>% do({
    data.frame(Correlation=.$logcpm[[1]] %$% cor(logCPM_A, logCPM_B))
})

Now we plot the logCPM for Ensembl vs UCSC, for each quantification method (program + transcriptome).

# Use the same x and y limits for all plots
limits <- logcpm.tables$logcpm %>% lapply(. %$% c(range(logCPM_A), range(logCPM_B))) %>% unlist %>% range
logcpm.plots <- logcpm.tables %>% group_by(program) %>% do({
    df <- .
    assert_that(nrow(df) == 1)
    xlabel <- str_c("Gene logCPM in ", df$datasetA)
    ylabel <- str_c("Gene logCPM in ", df$datasetB)
    ptitle <- sprintf("logCPM, %s vs %s for %s on %s genome", df$transcriptomeA, df$transcriptomeB, df$program, df$genome)
    tab <- df$logcpm[[1]] %>% filter(logCPM_A != min(logCPM_A) | logCPM_B != min(logCPM_B))
    p <- ggplot(tab) + 
        aes(x=logCPM_A, y=logCPM_B) +
        geom_point(size=0.2, alpha=0.2) +
        geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
        coord_fixed(xlim=limits, ylim=limits) +
        labs(title=ptitle, 
             subtitle="Genes matched by unique Entrez & Ensembl IDs; identity line (red) plotted red for reference",
             x=xlabel, y=ylabel)
    cbind(.[c("program", "genome", "datasetA", "datasetB", "transcriptomeA", "transcriptomeB")],
          data_frame(plot=list(p)))
})
ggprint(logcpm.plots$plot)

Overall concordance is very high for all methods. The alignment-based methods (STAR and HISAT2) have slightly higher correlation, which is matched by the apparently tighter spread of points about the identity line in the plots. This is expected since these methods discard ambiguous reads, and are therefore biased toward unambiguously-defined genes, which are the same genes more likely to be annotated consistently between Entrez and Ensembl.

Interestingly, the correlation for Shoal is lower than for Salmon, despite the fact that the Shoal plot looks like it has a tighter spread and the fact that Shoal is meant to regularize Salmon’s quantification results, which would be expected to increase the correlation. Perhaps Shoal decreases the overall correlation by spreading many points away from the bottom left corner, where previously all the genes with zero counts were exactly overlaid on a single point. (This kind of effect is why pearson correlations on gene abundances is rarely informative.)

STAR vs HISAT2

The two count-based methods differ only in the aligner that was used, so any discrepancy between them represents a difference in which reads they could align unambiguously.

logcpm_A <- logcpm$hisat2_grch38_snp_tran_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_A=logCPM)
logcpm_B <- logcpm$star_hg38.analysisSet_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_B=logCPM)
logcpm.merged <- inner_join(logcpm_A, logcpm_B) %>%
    filter(logCPM_A != min(logCPM_A) | logCPM_B != min(logCPM_B))
limits <- logcpm.merged %>% select(logCPM_A, logCPM_B) %>% lapply(range) %>% unlist %>% range
p <- ggplot(logcpm.merged) +
    aes(x=logCPM_A, y=logCPM_B) +
    geom_point(size=0.2, alpha=0.2) +
    geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
    coord_fixed(xlim=limits, ylim=limits) +
    labs(title="logCPM, STAR vs HISAT2", 
         subtitle="Identity line (red) plotted red for reference",
         x="logCPM by HISAT2", y="logCPM by STAR")
ggprint(p)

Effect of gene length normalization

We now compare Salmon, Shoal, and Kallisto with and without normalizing for gene length. This will show how much or how little the gene length normalization affects the abundance estimation.

logcpm.plots <- datasetmeta %>%
    filter(transcriptome == "ensembl.85", program %in% c("salmon", "shoal", "kallisto")) %>% 
    group_by(dataset, program) %>% do({
        assert_that(nrow(.) == 1)
        dset <- .$dataset
        program <- .$program
        logcpm_A <- logcpm.nooffset[[dset]] %>% 
            select(ENSEMBL, Sample, logCPM) %>%
            rename(logCPM_A=logCPM)
        logcpm_B <- logcpm[[dset]] %>% 
            select(ENSEMBL, Sample, logCPM) %>%
            rename(logCPM_B=logCPM)
        logcpm.merged <- inner_join(logcpm_A, logcpm_B) %>%
            filter(logCPM_A != min(logCPM_A) | logCPM_B != min(logCPM_B))
        limits <- logcpm.merged %>% select(logCPM_A, logCPM_B) %>% lapply(range) %>% unlist %>% range
        p <- ggplot(logcpm.merged) +
            aes(x=logCPM_A, y=logCPM_B) +
            geom_point(size=0.2, alpha=0.2) +
            geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
            coord_fixed(xlim=limits, ylim=limits) +
            labs(title=sprintf("%s logCPM, length-normalized vs non-length-normalized", program), 
                 subtitle="Identity line (red) plotted red for reference",
                 x="logCPM without gene length normalization", y="logCPM with gene length normalization")
        data_frame(plot=list(p))
    })
ggprint(logcpm.plots$plot)

STAR vs Salmon

We now compare STAR, a count-based method, to Salmon, an alignment-free method. We use the Salmon results with no gene length normalization for an apples-to-apples comparison.

logcpm_A <- logcpm$star_hg38.analysisSet_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_A=logCPM)
logcpm_B <- logcpm.nooffset$salmon_hg38.analysisSet_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_B=logCPM)
logcpm.merged <- inner_join(logcpm_A, logcpm_B) %>%
    filter(logCPM_A != min(logCPM_A) | logCPM_B != min(logCPM_B))
limits <- logcpm.merged %>% select(logCPM_A, logCPM_B) %>% lapply(range) %>% unlist %>% range
p <- ggplot(logcpm.merged) +
        aes(x=logCPM_A, y=logCPM_B) +
        geom_point(size=0.2, alpha=0.2) +
        geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
        coord_fixed(xlim=limits, ylim=limits) +
        labs(title="logCPM, Salmon vs STAR", 
             subtitle="Identity line (red) plotted red for reference",
             x="logCPM by STAR", y="logCPM by Salmon")
ggprint(p)

Salmon vs Kallisto

The next comparison is between Salmon and Kallisto, two conceptually very similar pseudoalignment-based methods.

logcpm_A <- logcpm$kallisto_hg38.analysisSet_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_A=logCPM)
logcpm_B <- logcpm$salmon_hg38.analysisSet_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_B=logCPM)
logcpm.merged <- inner_join(logcpm_A, logcpm_B) %>%
    filter(logCPM_A != min(logCPM_A) | logCPM_B != min(logCPM_B))
limits <- logcpm.merged %>% select(logCPM_A, logCPM_B) %>% lapply(range) %>% unlist %>% range
p <- ggplot(logcpm.merged) +
        aes(x=logCPM_A, y=logCPM_B) +
        geom_point(size=0.2, alpha=0.2) +
        geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
        coord_fixed(xlim=limits, ylim=limits) +
        labs(title="logCPM, Salmon vs Kallisto", 
             subtitle="Identity line (red) plotted red for reference",
             x="logCPM by Kallisto", y="logCPM by Salmon")
ggprint(p)

Salmon vs Shoal

Last, we compare the basic Salmon results to the same results regularized by Shoal.

logcpm_A <- logcpm$salmon_hg38.analysisSet_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_A=logCPM)
logcpm_B <- logcpm$shoal_hg38.analysisSet_ensembl.85 %>% 
    select(ENSEMBL, Sample, logCPM) %>%
    rename(logCPM_B=logCPM)
logcpm.merged <- inner_join(logcpm_A, logcpm_B) %>%
    filter(logCPM_A != min(logCPM_A) | logCPM_B != min(logCPM_B))
limits <- logcpm.merged %>% select(logCPM_A, logCPM_B) %>% lapply(range) %>% unlist %>% range
p <- ggplot(logcpm.merged) +
        aes(x=logCPM_A, y=logCPM_B) +
        geom_point(size=0.2, alpha=0.2) +
        geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
        coord_fixed(xlim=limits, ylim=limits) +
        labs(title="logCPM, Shoal vs Salmon", 
             subtitle="Identity line (red) plotted red for reference",
             x="logCPM by Salmon", y="logCPM by Shoal")
ggprint(p)