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.tabletopTagslrt.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}}