estimate the biological coefficient of variation across biological replicates

dispersion.R

#!/usr/bin/Rscript

library(edgeR)
library(limma)

## The BCV is the relative variability of expression between biological replicates
## The square root of the negative binomial dispersion for a gene 
## is the biological coefficient of variation (BCV) across replicates (stdev/mean)

## input: ge (raw count table filtered for lowly expressed genes, eg CPM > 1)
## condition=vector of treatment (eg. KO, WT)
exp.all <- data.frame(samples=colnames(ge), condition=condition)
dge<-DGEList(counts=ge, genes=rownames(ge), group=exp.all$condition)

cond <- factor(exp.all$condition)

model <- model.matrix(~0+cond)
colnames(model) <- c(levels(cond))
rownames(model) <- exp.all$samples

## TMM Normalization
dge <- calcNormFactors(dge, method="TMM")

## global(common) dispersion estimate averaged over all genes, a
## trended dispersion where the dispersion of a gene is predicted from its abundance
dge <- estimateDisp(dge, model, robust=TRUE)

estimateTrendedDisp(dge, verbose=TRUE)
estimateCommonDisp(dge, verbose=TRUE)
estimateTagwiseDisp(dge, verbose=TRUE)

## BCV plot: tagwise empirical Bayes moderated dispersion for each individual gene
edgeR::plotBCV(dge)

## bcv v. log2(CPM + 0.5)
v <- voom(dge, model, plot=TRUE)
vfit <- lmFit(v, model)

contr.matrix <- makeContrasts(
  KOvsWT = KO - WT,
  levels = colnames(model))

vfit.c <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit.c)

plotSA(efit) #variance no longer dependent on mean expression level as before precision weights were applied

## logFC v. logCPM, uses common dispersion
limma::plotMA(efit)