Copy number adjusted differential expression using edgeR and limma in a small region

[fpbarthel]

Hi all,

I would like to find differentially expressed genes within a region spanning approx 1 megabase while adjusting for the copy number status in that region. I'm looking for a way to do this using both limma (for microarray) and edgeR (for rnaseq). Reading the manuals I've come up with the basics but I'm stuck at defining the contrasts.

To start I've defined the following variables:
cases is a factor of length 50 with the different cases, all variables are of this length and in this order
genes is a list of genes in the 1 mb region
groups is a factor splitting cases in "Disease" (1) and "Normal" (0)
groups_cn is a factor of copy number states: "Normal" (0), "Amplified" (1) and "Deleted" (-1)

Limma

af_chr = af[rownames(af) %in% genes, , drop=F]
af_design = model.matrix(~ 0 + cn_groups + groups)
de_af = lmFit(af_chr, design=af_design)
fit = eBayes(de_af)
ans_af=topTable(fit, adjust="fdr", number=length(fit$coefficients), coef='groupsDisease')

When looking at ans_af, am I seeing the differentially expressed genes after adjusting for copy number? Or, am I seeing genes differentially expressed between disease and normal disregarding copy number?

edgeR

ge_chr = ge[rownames(ge) %in% genes, , drop=F]
ge_design = model.matrix(~ 0 + cn_groups + groups)
de_ge = DGEList(ge_chr, group = groups)
de_ge = estimateGLMCommonDisp(de_ge, design = ge_design)
de_ge = estimateGLMTagwiseDisp(de_ge, design = ge_design)
fit = glmFit(de_ge, design = ge_design) 
test = glmLRT(fit, coef = 'groupsDisease')
ans_ge = topTags(test)

Here the same question, when looking at ans_ge, am I seeing the differentially expressed genes after adjustment for copy number, or just the ones different between Disease and Normal?

Thanks

[Gordon Smyth]

It depends what you mean by "adjusting for copy number". Your analysis is adjusting for baseline differences between the copy number groups, but it is also assuming that the logFC between Disease and Normal is the same for all the copy number groups. Is that a reasonable assumption?

If you want to do the DE analysis between Disease and Normal separately for each copy number group, then that can be achieved by:

design <- model.matrix(~cn_groups + cn_groups:groups)