Exploration of CD4 RNA-Seq Dataset

Using Dataset salmon_hg38.analysisSet_ensembl.85

Preliminary Setup

First we load the necessary libraries, along with a set of utility functions.

library(here)
library(stringr)
library(magrittr)
library(Matrix)
library(lme4)
library(SummarizedExperiment)
library(dplyr)
library(edgeR)
library(limma)
library(ggplot2)
library(scales)
library(GGally)
library(ggalt)
library(reshape2)
library(assertthat)
library(ggfortify)
library(broom)
library(codingMatrices)
library(variancePartition)

library(doParallel)
ncores <- getOption("mc.cores", default=parallel::detectCores(logical = FALSE))
registerDoParallel(cores=ncores)
library(BiocParallel)
BiocParallel:::.registry_init() # Workaround for https://github.com/Bioconductor/BiocParallel/issues/65
register(DoparParam())

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

We also set the output directory for all plot files.

plotdir <- here("plots/RNA-seq", params$dataset)
dir.create(plotdir, recursive = TRUE, showWarnings = FALSE)

We also set the random seed for reproducibility:

set.seed(1986)

Data Loading and Preprocessing

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

sexpfile <- here("saved_data", sprintf("SummarizedExperiment_rnaseq_%s.RDS", params$dataset))
sexp <- readRDS(sexpfile)
colnames(sexp) <- colData(sexp)$SampleName

We extract the sample metadata from the SummarizedExperiment. We also tell R to use a coding matrix for each factor that puts the intercept at the mean of all factor levels when incorporating it into a design matrix.

sample.table <- colData(sexp) %>%
    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=""))
for (i in names(sample.table)) {
    if (is.factor(sample.table[[i]]) && nlevels(sample.table[[i]]) > 1) {
        contrasts(sample.table[[i]]) <- code_control_named(levels(sample.table[[i]]))
    }
}

Next we extract the count matrix from the 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")
if (all(libtype.assayNames %in% assayNames(sexp))) {
    message("Selecting stranded counts for each sample")
    sample.table %<>% mutate(count_type=libtype.assayNames[libType])
    assay(sexp, "unstranded.counts") <- assay(sexp, "counts")
    assay(sexp, "counts") <- lapply(seq_len(nrow(sample.table)), function(i) {
        message("Using ", sample.table[i,]$count_type, " for ", colnames(sexp)[i])
        assay(sexp, 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.

if (all(libtype.assayNames %in% assayNames(sexp))) {
    total.counts <- sexp %>% 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))))
}

Exploratory Analysis

Normalization & Filtering

Now we create a DGEList from the counts.

## 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

Next we take care of the initial scaling normalization for sequencing depth and composition bias. We also discard any genes with all zero counts, since there is no meaningful analysis that can be done with these genes.

## Remove all genes with zero counts in all samples
nonzero <- rowSums(dge$counts) > 0
dge %<>% .[nonzero,]
dge %<>% calcNormFactors

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 Kallisto.

if ("length" %in% assayNames(sexp)) {
    normMat <- assay(sexp, "length") %>% divide_by(exp(rowMeans(log(.)))) %>%
        .[nonzero,]
    normCounts <- dge$counts/normMat
    lib.offsets <- log(calcNormFactors(normCounts)) + log(colSums(normCounts))
    dge$offset <- t(t(log(normMat)) + lib.offsets)
}

We plot the distribution of average log2 CPM values to verify that our chosen presence threshold is appropriate. The distribution is expected to be bimodal, with a low-abundance peak representing non-expressed genes and a high-abundance peak representing expressed genes. The chosen threshold should separate the two peaks of the bimodal distribution.

a <- aveLogCPMWithOffset(dge)
avelogcpm.presence.threshold <- -1

p <- list(
    Histogram=ggplot(data.frame(logCPM=a)) +
        aes(x=logCPM) +
        geom_histogram(aes(y=100*(..count..)/sum(..count..)), binwidth=0.25, boundary=0) +
        geom_vline(xintercept=avelogcpm.presence.threshold, color="red", linetype="dashed") +
        xlab("Average logCPM") + ylab("Percent of genes in bin") +
        coord_cartesian(xlim=quantile(a, c(0, 0.995)), ylim=c(0,10)) +
        labs(title="Average gene LogCPM distribution",
             subtitle="for genes with at least 1 read") +
        theme(plot.caption = element_text(hjust = 0)),
    ECDF=ggplot(fortify(ecdf(a))) +
        aes(x=x, y=y*100) +
        geom_step() +
        geom_vline(xintercept=avelogcpm.presence.threshold, color="red", linetype="dashed") +
        xlab("Average logCPM") + ylab("Percent of genes with smaller average logCPM") +
        coord_cartesian(xlim=quantile(a, c(0, 0.995))) +
        labs(title="Empirical Cumulative Distribution Function of gene LogCPM values",
             subtitle="for genes with at least 1 read") +
        theme(plot.caption = element_text(hjust = 0)))

ggprint(p)

The red dashed line in each plot indicates the chosen presence threshold. We now subset the DGEList to only those genes above the threshold.

dge %<>% .[aveLogCPMWithOffset(.) >= avelogcpm.presence.threshold,]

Variance & Heteroskedasticity

Now we estimate the dispersions for each gene, to get an idea of what the variability of this data set is like. In order to evaluate the effect of empirical Bayes shrinkage on the dispersions, we estimate the gene dispersions in 3 different ways: once with no inter-gene information sharing, once with ordinary shrinkage, and once with robust shrinkage, which reduces the strength of shrinkage for outlier genes whose dispersion is farthest away from the trend.

design <- model.matrix(~0 + group + donor_id, sample.table)
colnames(design) %<>% str_replace("^group", "")

dge %<>% estimateDisp(design, robust=TRUE)

message("Common dispersion: ", dge$common.dispersion)

Common dispersion: 0.634358039982791

message("BCV: ", sqrt(dge$common.dispersion))

BCV: 0.796465969130377

dge.with.eBayes <- dge %>% estimateDisp(design, robust=FALSE)
dge.with.robust.eBayes <- dge %>% estimateDisp(design, robust=TRUE)
dge.without.eBayes <- dge %>% estimateDisp(design, prior.df=0)

We now plot all 3 dispersion estimates, along with the overall average and estimated trend. Each plot includes the points from the previous plots in lighter colors for comparison.

disptable <- data.frame(
    logCPM=dge.without.eBayes$AveLogCPM,
    CommonBCV=dge.with.eBayes$common.dispersion %>% sqrt,
    TrendBCV=dge.with.eBayes$trended.dispersion %>% sqrt,
    GeneWiseBCV=dge.without.eBayes$tagwise.dispersion %>% sqrt,
    eBayesBCV=dge.with.eBayes$tagwise.dispersion %>% sqrt,
    RobustBCV=dge.with.robust.eBayes$tagwise.dispersion %>% sqrt) %>%
    cbind(dge$genes)

## Reduce the number of points to plot for each line for performance
## reasons
npoints <- c(Common=2, Trend=500)
disp.line.table <-
    disptable %>%
    select(logCPM, TrendBCV, CommonBCV) %>%
    melt(id.vars="logCPM", variable.name="DispType", value.name = "BCV") %>%
    mutate(DispType=str_replace(DispType, "BCV$", "")) %>%
    group_by(DispType) %>%
    do({
        spline(x=.$logCPM, y=.$BCV, n=npoints[.$DispType[1]]) %>% data.frame(logCPM=.$x, BCV=.$y)
    })

baseplot <- ggplot(disptable) +
    aes(x=logCPM)
raw.disp.plot <- baseplot +
    geom_point(aes(y=GeneWiseBCV), size=0.1, color="black") +
    geom_line(data=disp.line.table, aes(x=logCPM, y=BCV, group=DispType), color="white", size=1.5, alpha=0.5) +
    geom_line(data=disp.line.table, aes(x=logCPM, y=BCV, linetype=DispType), color="darkred", size=0.5) +
    scale_linetype_manual(name="Dispersion Type", values=c(Trend="solid", Common="dashed")) +
    ylab("Biological coefficient of variation") +
    ggtitle("BCV plot (Raw dispersions)")

eBayes.disp.plot <- baseplot +
    geom_point(aes(y=GeneWiseBCV), size=0.4, color="gray") +
    geom_point(aes(y=eBayesBCV), size=0.1, color="darkblue") +
    geom_line(data=disp.line.table, aes(x=logCPM, y=BCV, group=DispType), color="white", size=1.5, alpha=0.5) +
    geom_line(data=disp.line.table, aes(x=logCPM, y=BCV, linetype=DispType), color="darkred", size=0.5) +
    scale_linetype_manual(name="Dispersion Type", values=c(Trend="solid", Common="dashed")) +
    ylab("Biological coefficient of variation") +
    ggtitle("BCV plot (Raw & squeezed dispersions)")

robust.eBayes.disp.plot <- baseplot +
    geom_point(aes(y=GeneWiseBCV), size=0.4, color="gray") +
    geom_point(aes(y=eBayesBCV), size=0.4, color="deepskyblue") +
    geom_point(aes(y=RobustBCV), size=0.1, color="darkgreen") +
    geom_line(data=disp.line.table, aes(x=logCPM, y=BCV, group=DispType), color="white", size=1.5, alpha=0.5) +
    geom_line(data=disp.line.table, aes(x=logCPM, y=BCV, linetype=DispType), color="darkred", size=0.5) +
    scale_linetype_manual(name="Dispersion Type", values=c(Trend="solid", Common="dashed")) +
    ylab("Biological coefficient of variation") +
    ggtitle("BCV plot (Raw & squeezed & robust dispersions)")

p <- list(raw.disp.plot, eBayes.disp.plot, robust.eBayes.disp.plot)
ggprint(p)

Next, we use limma’s sample weight calculating methods to investigate possible quality issues. To confirm our results, we also split the samples into treatment groups and use edgeR to estimate the dispersion within each group. This is a crude method, and these group-specific dispersion estimate would be too unstable for use in a differential expression analysis, but comparing the overall mean dispersion for each sample to the sample quality weights determined by limma provides a useful sanity check.

elist.w <- voomWithQualityWeightsAndOffset(dge, design)
dbg <- estimateDispByGroup(dge, sample.table$group, sample.table$batch)

To see whether the weights are correlated with specific experimental factors, we create a boxplot of the weights and group dispersions against each relevant covariate.

covars <- sample.table %>% dplyr::select(group, time_point, donor_id, batch, cell_type)
qcmetrics <- data.frame(Weight=elist.w$sample.weights,
                        GroupBCV=sapply(dbg, `[[`, "common.dispersion")[as.character(covars$group)])

p <- ggduo.dataXY(covars, qcmetrics %>% transmute(Log2Weight=log2(Weight), GroupBCV)) +
  ggtitle("Weights and dispersions by group") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=0.5))
ggprint(p)

We also make plots of the individual weights against each covariate, and compute the ANOVA p-value for the relationship between the weights and each covariate separately.

(Note for later: It would be good to try all covariates in the same model as random effects.)

awdf <- data.frame(covars, qcmetrics)
anovas <- lapply(colnames(covars), function(x) {
    formula <- as.formula(str_c("log2(Weight) ~ ", x))
    lm(formula, awdf) %>% aov %>% tidy %>% filter(term == x)
}) %>%
    do.call(what=rbind) %>%
    mutate(padj=p.adjust(p.value, method="BH"))
pvals <- anovas %$% setNames(p.value, term)
aw.plot.base <- ggplot(awdf) +
    aes(y=Weight) +
    scale_y_continuous(trans=log_trans(2)) +
    geom_dotplot(binaxis = "y", stackdir="center", binwidth=0.1)
p <- lapply(colnames(covars), function(x) {
    pretty.covar.name <- x %>% str_replace_all("_", " ") %>% str_to_title
    aw.plot.base + aes_string(x=x) +
        labs(title=str_c("Array weights by ", pretty.covar.name),
             subtitle=sprintf("ANOVA p-value = %0.3g", pvals[x]))
})
ggprint(p)

Batch correction

Next we perform several methods of batch correction on the voom-transformed data. First, we try direct batch subtraction using limma’s removeBatchEffect function, which fits a linear model with the batches as sum-to-zero coefficients and then subtracts the batch coefficients from the data.

elist.bc <- elist.w
# Batch is confounded with time point, so we leave it out of the design for this step
design.NoTime <- model.matrix(~ cell_type + donor_id, sample.table)
elist.bc$E %<>%
    removeBatchEffect(batch=sample.table %$% batch:donor_id, design=design.NoTime,
                      weights=elist.bc$weights)

Coefficients not estimable: batch5 batch6 batch7 

Warning: Partial NA coefficients for 17029 probe(s)

Second, we use ComBat, which performs empirical Bayes shrinkage of the batch correction parameters. First, we run ComBat in parametric prior mode in order to produce a diagnostic plot of the priors. Note that we use ComBat here to subtract both the batch and donor effects, so that the MDS plots will be more reflective of the desired biological effects.

# Right now we're just running ComBat to produce the plot, so we discard the output.
invisible(capture.output(ComBat(elist.w$E, batch=sample.table$batch, mod=design.NoTime, par.prior=TRUE, prior.plots = TRUE)))

Found2batches

Adjusting for4covariate(s) or covariate level(s)

Fitting L/S model and finding priors

Finding parametric adjustments

Adjusting the Data

cairo_pdf(file.path(plotdir, "rnaseq-ComBat-qc.pdf"))
# Run the same thing again to output the plot to the PDF file
invisible(ComBat(elist.w$E, batch=sample.table$batch, mod=design.NoTime, par.prior=TRUE, prior.plots = TRUE))

Standardizing Data across genes

invisible(dev.off())

Since the parametric fit for the variance is not a good match for the empirical distribution, we now perform the actual batch correction using the non-parametric mode of ComBat.

# Now perform the actual batch correction using non-parametric prior
design.cb <- model.matrix(~cell_type, sample.table)
elist.cb <- elist.w
elist.cb$E %<>% ComBat(batch=sample.table$batch, mod=design.NoTime, par.prior=FALSE)

Found2batches

Adjusting for4covariate(s) or covariate level(s)

Standardizing Data across genes

Fitting L/S model and finding priors

Finding nonparametric adjustments

Adjusting the Data

For each of the batch correction methods, we compute the sample distance matrix using multidimensional scaling and plot the first 3 principal coordinates. We reflect all the principal coordinates so as to make the mean of the Naive D0 samples negative in all dimensions, so that MDS plots have a greater chance of being oriented consistently with each other.

naive.d0.samples <- sample.table$group == "NaiveD0"
dmat <- suppressPlot(plotMDS(elist.w)$distance.matrix) %>% as.dist
mds <- cmdscale(dmat, k=attr(dmat, "Size") - 1, eig=TRUE)
mds$points %<>% scale(center=FALSE, scale=-sign(colMeans(.[naive.d0.samples,]))) %>%
    add.numbered.colnames("Dim") %>% data.frame(sample.table, .)
dmat.bc <- suppressPlot(plotMDS(elist.bc)$distance.matrix) %>% as.dist
mds.bc <- cmdscale(dmat.bc, k=attr(dmat, "Size") - 1, eig=TRUE)
mds.bc$points %<>% scale(center=FALSE, scale=sign(colMeans(.[naive.d0.samples,]))) %>%
    add.numbered.colnames("Dim") %>% data.frame(sample.table, .)
dmat.cb <- suppressPlot(plotMDS(elist.cb)$distance.matrix) %>% as.dist
mds.cb <- cmdscale(dmat.cb, k=attr(dmat, "Size") - 1, eig=TRUE)
mds.cb$points %<>% scale(center=FALSE, scale=sign(colMeans(.[naive.d0.samples,]))) %>%
     add.numbered.colnames("Dim") %>% data.frame(sample.table, .)

ggmdsbatch <- function(dat, dims=1:2) {
    if (length(dims) == 1) {
        dims <- dims + c(0,1)
    }
    assert_that(length(dims) == 2)
    ggplot(dat) +
        aes_string(x=str_c("Dim", dims[1]), y=str_c("Dim", dims[2])) +
        aes(color=batch, label=SampleName) +
        geom_text() +
        scale_x_continuous(expand=c(0.15, 0)) +
        coord_equal()
}

p <- list(
    ggmdsbatch(mds$points) +
        labs(title="limma voom Principal Coordinates 1 & 2",
             subtitle="No batch correction"),
    ggmdsbatch(mds.bc$points) +
        labs(title="limma voom Principal Coordinates 1 & 2",
             subtitle="After naive batch subtraction"),
    ggmdsbatch(mds.cb$points) +
        ggtitle("limma voom Principal Coordinates 1 & 2",
             subtitle="After ComBat batch correction"),
    ggmdsbatch(mds$points, dims=2:3) +
        ggtitle("limma voom Principal Coordinates 2 & 3",
             subtitle="No batch correction"),
    ggmdsbatch(mds.bc$points, dims=2:3) +
        ggtitle("limma voom Principal Coordinates 2 & 3",
             subtitle="After naive batch subtraction"),
    ggmdsbatch(mds.cb$points, dims=2:3) +
        ggtitle("limma voom Principal Coordinates 2 & 3",
             subtitle="After ComBat batch correction"))
ggprint(p)

The naive batch subtraction seems to leave one donor as an outlier for unclear reasons, while ComBat avoids this artifact, instead yielding abiologically more plausible MDS plot.

MDS Plots

Choosing ComBat as the best-looking batch correction, we make more MDS plots for this data, this time plotting the first 5 PCs and adding the experimental information to the plot as colors and shapes.

xlims <- range(unlist(mds.cb$points[c("Dim1", "Dim2")]))
ylims <- range(unlist(mds.cb$points[c("Dim2", "Dim3")]))
pbase <- ggplot(mds.cb$points) +
    aes(x=Dim1, y=Dim2, label=SampleName, color=batch, fill=time_point, shape=cell_type, linetype=donor_id, group=cell_type:donor_id) +
    geom_encircle(aes(group=time_point:cell_type, color=NULL), s_shape=0.75, expand=0.05, color=NA, alpha=0.2) +
    geom_path(color=hcl(c=0, l=45), aes(color=NULL)) +
    geom_point(size=4) +
    scale_shape_manual(values=c(Naive=21, Memory=24)) +
    scale_color_manual(values=col2hcl(c(B1="green", B2="magenta"), l=80)) +
    scale_fill_hue(l=55) +
    scale_linetype_manual(values=c("solid", "dashed", "dotdash", "twodash")) +
    guides(colour = guide_legend(override.aes = list(shape = 21)),
           fill = guide_legend(override.aes = list(shape = 21)),
           shape = guide_legend(override.aes = list(color=hcl(c=0, l=80), fill=hcl(c=0, l=55)))) +
    labs(title="limma voom Principal Coordinates 1 & 2",
         subtitle=" (after ComBat batch correction)") +
    coord_equal(xlim=xlims, ylim=ylims)
p <- list(PC12=pbase,
          PC23=pbase +
              aes(x=Dim2, y=Dim3) +
              labs(title="limma voom Principal Coordinates 2 & 3"),
          PC34=pbase +
              aes(x=Dim3, y=Dim4) +
              labs(title="limma voom Principal Coordinates 3 & 4"),
          PC45=pbase +
              aes(x=Dim4, y=Dim5) +
              labs(title="limma voom Principal Coordinates 4 & 5"))
ggprint(p)

Variance Partitioning analysis

To further investigate the sources of variance within the data, we can use the variancePartition package. We fit the model to the uncorrected data and the two corrected data sets (simple batch subtraction and ComBat) so that we can see how the percent of variance explained is affected by batch correction. We use random effects for the factors in the model since batch and time point are confounded.

# Fit all factors as random effects, since group and batch are confounded
vp.formula <- ~ (1|group) + (1|donor_id) + (1|batch)
elists <- list(NoBC=elist.w, BatSub=elist.bc, ComBat=elist.cb)
# Function is already parallelized, so don't call it in parallel
varParts <- mapply(function(...) try(fitExtractVarPartModel(...)), 
                   exprObj = elists, 
                   MoreArgs = list(
                       formula = vp.formula, 
                       data=sample.table),
                   SIMPLIFY = FALSE)

Projected run time: ~ 191 min 
Projected run time: ~ 4 min 
Projected run time: ~ 4 min 

# Collapse SVs to a single column
varTables <- list()
for (i in names(varParts)) {
    if (is(varParts[[i]], "try-error")) {
        message("Could not run variancePartition for ", i, ", probably due to collinearity of covariates.")
    } else {
        assert_that(is(varParts[[i]], "varPartResults"))
        x <- as(varParts[[i]], "data.frame")
        x.sv <- x %>% select(dplyr::matches("^SV\\d+$"))
        if (ncol(x.sv) > 0) {
            x.nosv <- x[setdiff(colnames(x), colnames(x.sv))]
            x <- data.frame(x.nosv, SV=rowSums(x.sv))
        }
        varTables[[i]] <- cbind(select(x, -Residuals), select(x, Residuals))
    }
}

The variancePartition function failed to run for the model that included both known covariates and surrogate variables. This is expected if the surrogate variables are highly correlated with any of the known covariates, which we have already observed is the case for donor ID.

p <- list()
for (i in names(varTables)) {
    incl <- str_replace_all(i, "_and_", " + ")
    p[[i]] <- plotVarPart(varTables[[i]]) +
        labs(title=str_c("Variance Partitions, ", incl))
}
ggprint(p)

The results seen here are consistent with what can be seen qualitatively from the MDS plots. With no batch correction, batch effects are a significant contributor to overall variance, which is problematic since batch and time point are confounded. With simple batch subtraction. recall that the MDS plots showed an outlier donor. This is reflected here in the very large contribution of donor ID to the variance of many genes. Lastly, the variance partition plot for ComBat looks quite reasonable, with most of the variance for a majority of genes explained by group (i.e. the combination of time point and cell type) and some smaller part explained by donor ID.