Thank you Yunshun!
But there are some confusion remains:
I did a similar comparison: finding DEGs in disease2 vs. control: while
classic EdgeR gives some DEGs, GLM approach gives 0 DEG, which is strange.
I believe that when using GLM approach, we also consider the patients
information and thus, the DEGs list should be larger than disregarding such
information (edgeR classic).


Also, when I want to find the response to hormone in any disease group, it
yields exceptions:

Also, the
#Genes that response to hormone in any disease group
lrt <- glmLRT(fit, coef=10:12)
topTags(lrt)
summary(de <- decideTestsDGE(lrt))
pdf("hormoneEffectInAnyThing.pdf")
detags <- rownames(y)[as.logical(de)]
plotSmear(lrt, de.tags=detags)
abline(h=c(-1, 1), col="blue")
dev.off()


 It yields an exception:

> plotSmear(lrt, de.tags=detags)
Error in plotSmear(lrt, de.tags = detags) :
  table$logFC slot in DGELRT object is NULL. We cannot produce an MA
(smear) plot if more than one coefficient from the GLM is being tested in
the likelihood ratio test as this results in more one logFC value per
gene---one for each coefficient.
> abline(h=c(-1, 1), col="blue")

Any suspect? The design is absolutely the same as the one in the example.
Son.


On Tue, Mar 11, 2014 at 6:54 PM, Yunshun Chen <yuchen@wehi.edu.au> wrote:

> Hi Son,
>
>
>
> First of all, under the classic edgeR you estimated the dispersions only
> from a subset of the data.
>
> The dispersion estimates were not as reliable as those estimated under the
> GLM approach (since the information of the other 12 samples were totally
> ignored).
>
>
>
> Secondly, your classic edgeR analysis did not account for the patient
> effect at all.
>
> These 6 samples are actually paired (by patient) and cannot be treated as
> independent samples.
>
> This is the main reason why you didn't get any DE genes under the classic
> edgeR.
>
>
>
> As stated in the users' guide, the classic edgeR is only applicable on
> datasets with a single factor design.
>
> Hence, you have to use the GLM approach for your data.
>
>
>
> Cheers,
>
> Yunshun
>
>
>
>
>
>
>
>
> --------------------------------------------------------------------------------------------------------------------------------------------
>
>
>
>
>
> Hi All,
>
> I have an experiment that is exactly the same as the one described in
>
> section 3.5 (Comparisons Both Between and Within Subjects)  in the
>
> manual<
> http://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
> >
>
> :
>
> The experiment has 18 RNA samples collected from 9 subjects. The samples
>
> correspond to cells from 3 healthy patients, either treated or not with a
>
> hormone; cells from 3 patients with disease 1, either treated or not with
>
> the hormone; and cells from 3 patients with disease 2, either treated or
>
> not with the hormone.
>
>
>
> I now find  genes responding to the hormone in disease1 patients using two
>
> approachs glm and edgeR classic. The results are dramatically different !!!
>
> glm gives hundreds of differentially expressed genes, while edgeR classic
>
> gives 0.
>
> Could anyone tell me that is what one can expect or I'm totally wrong !
>
>
>
>
>
> data.frame(Disease,Patient,Treatment)
>
> design <- model.matrix(~Disease+Disease:Patient+Disease:Treatment)
>
> colnames(design)
>
> y <- estimateGLMCommonDisp(y,design)
>
> y <- estimateGLMTrendedDisp(y,design)
>
> y <- estimateGLMTagwiseDisp(y,design)
>
> fit <- glmFit(y, design)
>
>
>
> lrt <- glmLRT(fit, coef="DiseaseDisease1:TreatmentHormone")
>
> topTags(lrt)
>
> detags <- rownames(topTags(lrt)$table)
>
> cpm(y)[detags, order(y$samples$group)]
>
> summary(de <- decideTestsDGE(lrt))
>
> This gives  612 (-1) and 1129 (+1) genes.
>
>
>
> I now try to just use 3 patients with Disease1, treated and non treated
>
> with hormone as 2 groups and apply the classic edgeR to find the effect of
>
> hormon to disease1:
>
>
>
> targets <-readTargets()
>
> x <-read.delim("new.counts", row.names="id")
>
> x_c_n <- x[c("d1_1_hormone", "d1_2_hormone", "d1_3_hormone", "d_1", "d_2",
>
> "d_3")]
>
> g <- factor(c(1,1,1,2,2,2))
>
> y <- DGEList(counts=x_c_n, group=g)
>
> y <- calcNormFactors(y)
>
> y <- estimateCommonDisp(y)
>
> y <- estimateTagwiseDisp(y)
>
> et <- exactTest(y)
>
> topTags(et)
>
> summary(de <- decideTestsDGE(et))
>
> This gives 0 (-1) and also 0 (+1) gene !!!!
>
>
>
> How can there be such a significant difference between these two
>
> approaches? Or I'm wrong?
>
>
>
> Thank you,
>
>
>
> Son.
>
>
>
> ______________________________________________________________________
> The information in this email is confidential and inte...{{dropped:10}}

