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)
	[[alternative HTML version deleted]]