Hello,

I am analysing RNA-seq data to find differentially expressed genes between a wild-type (Ctl) and mutant (Mut). The mutant has a large DNA duplication encompassing several hundred genes, so they always come up as ~2x upregulated but this is not biologically interesting. Hence the H0 (null) for these particular genes is that Mut = 2 x Ctl, while H0 for the rest of the genes is Mut = Ctl.

My workaround for this was to do two runs of contrast testing on the limma.voom output, one run where H0 is Mut = Ctl, the second where H0 is Mut = 2 x Ctl. Then I would get 2 topTable outputs, one from which I'd select rows for the non-duplicated genes and the other from which I'd select the duplicated genes. However, I am having difficulty specifying (or perhaps even understanding) the contrasts for these comparisons and the results are not making sense to me. In particular, I would have expected the log-FC values from both comparisons to be well-correlated, but they are actually very different, especially for genes with high expression. So is this even a valid approach?

Here is an example using the zebrafish dataset:

library(zebrafishRNASeq)
library(limma)
library(edgeR)
library(ggplot2)
data(zfGenes)
filter <- apply(zfGenes, 1, function(x) length(x[x>10]) > 1)
cmat <- zfGenes[filter,]
design <- matrix(c(rep(1,3), rep(0,6), rep(1,3)), ncol=2, dimnames = list(c(),c('ctl','mut')))
nf <- calcNormFactors(cmat)
y <- voom(cmat, design, lib.size=colSums(cmat)*nf)
fit <- lmFit(y,design)

## make a function to run eBayes etc. with different null hypotheses ##
run_cf_eb <- function(H0 = 1, fit){
  contrasts = paste0('mut-', H0, '*ctl')
  cont.matrix <- makeContrasts(contrasts = contrasts, levels=c('ctl','mut'))
  print(cont.matrix)
  fit_b <- contrasts.fit(fit, cont.matrix)
  fit_b <- eBayes(fit_b)
  return( topTable(fit_b, n=Inf, sort.by='none') )
}

tt_1 <- run_cf_eb(H0 = 1, fit) ## H0 is 'ctl = 1 x mut' ##
tt_2 <- run_cf_eb(H0 = 2, fit) ## H0 is 'ctl = 2 x mut' ##
tt_all <- merge(tt_1, tt_2, by = 'row.names', suffixes = c('.H0_1','.H0_2'))

## Plot logFC from H0('ctl=mut') vs logFC from H0('ctl=2mut'): ##
ggplot(tt_all) + geom_point(aes(x=logFC.H0_1, y=logFC.H0_2, colour = AveExpr.H0_1))

Output:

      Contrasts
Levels mut-1*ctl
   ctl        -1
   mut         1

      Contrasts
Levels mut-2*ctl
   ctl        -2
   mut         1

I realize that I will probably have to tweak the calcNormFactors bit to account for the duplication and also re-run the p-value adjustment after the analysis, however my question right now is just about how to approach the contrasts.

Thanks,

Stuart

You are forgetting that limma is working with expression values on the log2-scale, so a 2-fold increase in Mutant vs Control is represented by Mut - Ctl = log2(2) = 1. If you conduct a standard DE analysis by:

design <- cbind(Intercept=1, MutvsCtl=c(0,0,0,1,1,1))
y <- voom(cmat, design, lib.size=colSums(cmat)*nf)
fit <- lmFit(y, design)
fit <- eBayes(fit)
topTable(fit, coef=2)

then you would see a lot of logFC values about 1 in the top-table for genes in the duplicated region. Is this what you actually see?

If 'InDup' is TRUE for genes in the duplicated region, then look at:

plotMD(fit, col=2, status=InDup)

to see whether the logFC are doing what you expect.

To formally test whether the true logFC is equal to 1, you can do this:

t.stat <- (fit$coefficients[,2] - 1) / fit$stdev.unscaled[,2] / sqrt(fit$s2.post)
p.value <- 2*pt(-abs(t.stat), df=fit$df.total)
FDR <- p.adjust(p.value, method="BH")

To formally test whether genes in the duplicated region tend to be up-regulated overall you can conduct a gene set test:

roast(y, index=InDup, design=design)

or even better:

roast(y, index=InDup, design=design, gene.weights=pmax(fit$Amean,0))

to weight genes according to their expression levels.

Note: this answer was edited after Stuart's comment.