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
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
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
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
Code:
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')
edgeR
Code:
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)
Thanks
Comment