[BioC] Contrast Problem

[Gordon K Smyth]

Dear Aditi,

Thanks for going to quite a bit of effort to describe your problem, but 
I'm afraid that I still don't follow entirely what you're trying to do. 
I wonder how you have arrived at the contrasts you have defined.

When you say different expression (DE) between genomes for all the 
accessions, do you mean DE between genomes for each of the accessions 
separately?  That is DE genes for A.F1 vs D.F1, for A.Tom vs D.Tom, for 
A.Tx vs D.Tx and for A.Mx vs D.Mx?

When you say DE between accessions, is that for each genome separately? 
In other words, do you expect the differences between the accessions to be 
relatively the same for each genome, or will the differences between 
accessions be genome-specific?

Your comments about zero expression and not detecting significantly DE 
genes don't make sense to me.  I won't try to respond to those comments 
because I think sorting out the above questions will probably solve other 
perceived problems as well.

Best wishes
Gordon

> Date: Thu, 28 Jun 2012 12:01:50 -0700
> From: Aditi Rambani <aditi_rambani at yahoo.com>
> To: "bioconductor at stat.math.ethz.ch" <bioconductor at stat.math.ethz.ch>
> Subject: [BioC] Contrast Problem
>
> Hello, 
>
> I am a graduate student at Brigham Young University working on polyploid 
> cotton RNA seq data. Our study design has two explanatory variables, one 
> is 'accession' with four levels (F1,Tom,Tx,Mx) and other one is 'genome' 
> with two levels (A genome or D genome). We want to detect differential 
> expression of genes between 'genomes' from all the accessions and also 
> find genes that are differentially expressed between accessions. We 
> built a (Accession*Genome) model and did a contrast for two levels of 
> 'genomes'. In contrast results we see that many genes with zero 
> expression (0 RPKM) have significant FDRs and some significantly 
> differentially expressed genes are not detected. We are not sure why 
> this is happening, any help will be greatly appreciated.
>
> Thanks
>
> Aditi 
>
> We are using following script to do our analysis but 
>
>
> library("edgeR")
>
> counts <- read.table(INFILE, header=T, row.names=1)
> groups <- factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8))
> accessions <- factor(c(1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4))
> genomes <- factor(c(1,1,1,2,2,2,1,1,1,2,2,2,1,1,1,2,2,2,1,1,1,2,2,2))
>
>
> design2 <- model.matrix(~accessions*genomes)
>
> dge <- DGEList(counts=counts, group=groups)
> dge <- calcNormFactors(dge)
> dge2 <- estimateGLMCommonDisp(dge, design2)
> dge2 <- estimateGLMTrendedDisp(dge2, design2)
> dge2 <- estimateGLMTagwiseDisp(dge2, design2)
>
> fit2 <- glmFit(dge2, design2)
>
> lrt.acc1 <- glmLRT(dge2, fit2, contrast=c(0,0,1,1,-1,-1,0,0))
> lrt.acc2 <- glmLRT(dge2, fit2, contrast=c(0,0,1,1,1,1,-2,-2))
> lrt.acc3 <- glmLRT(dge2, fit2, contrast=c(-3,-3,1,1,1,1,1,1))
> lrt.F1 <- glmLRT(dge2, fit2, contrast=c(1,-1,0,0,0,0,0,0))
> lrt.Mx <- glmLRT(dge2, fit2, contrast=c(0,0,1,-1,0,0,0,0))
> lrt.Tx <- glmLRT(dge2, fit2, contrast=c(0,0,0,0,1,-1,0,0))
> lrt.Tom <- glmLRT(dge2, fit2, contrast=c(0,0,0,0,0,0,1,-1))
>
>
> write.table(topTags(lrt.acc1, n=10000), file="acc1.results", sep="\t", quote=F)
> write.table(topTags(lrt.acc2, n=10000), file="acc2.results", sep="\t", quote=F)
> write.table(topTags(lrt.acc3, n=10000), file="acc3.results", sep="\t", quote=F)
> write.table(topTags(lrt.F1, n=10000), file="F1.results", sep="\t", quote=F)
> write.table(topTags(lrt.Mx, n=10000), file="Mx.results", sep="\t", quote=F)
> write.table(topTags(lrt.Tx, n=10000), file="Tx.results", sep="\t", quote=F)
> write.table(topTags(lrt.Tom, n=10000), file="Tom.results", sep="\t", quote=F)


______________________________________________________________________
The information in this email is confidential and intend...{{dropped:5}}