Question: Why the result of pairwise comparison and Anova comparison of the same contrast is different ?
0
gravatar for star
5 months ago by
star10
star10 wrote:

I have a group line with 2 treatments and one control and I did DE using EdgeR. My design matrix is as below and I run contrast for those samples using pairwise and anova test. Although their LigFC are the same for DE genes, their FDR are different.

data
           Ensembl_ID Gene_name    H9_wt_B1  H9_wt_B2   H9_KO_B1     H9_KO_B2   H9_KI_B1       H9_KI_B2 
## 1 ENSG00000000003    TSPAN6     4274     4233          3755          4937       4681          4061  
## 2 ENSG00000000005      TNMD      116      117           106           154       105            86     
## 3 ENSG00000000419      DPM1     2443     2391          2389          2597       3026          2453


design
##                  samples
## H9_wt_B1           H9_wt
## H9_wt_B2           H9_wt
## H9_KO_B1           H9_KO
## H9_KO_B2           H9_KO
## H9_KI_B1           H9_KI
## H9_KI_B2           H9_KI

all(rownames(design) %in% colnames(data))
## [1] TRUE
all(rownames(design) == colnames(data))
## [1] TRUE

design.matrix_H9 <- model.matrix(~0+samples,design)
design.matrix_H9
##               samplesH9_KI samplesH9_KO samplesH9_wt
## H9_wt_B1                 0                 0            1
## H9_wt_B2                 0                 0            1
## H9_KO_B1                 0                 1            0
## H9_KO_B2                 0                 1            0
## H9_KI_B1                 1                 0            0
## H9_KI_B2                 1                 0            0
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$samples
## [1] "contr.treatment"


H9_list <- DGEList(counts = data[c(3:8)], group=design$samples, genes=data[1:2])
H9_filter <- filterByExpr(H9_list)
H9_keep<-H9_list[H9_filter, , keep.lib.sizes=FALSE]
H9_keep_norm <- calcNormFactors(H9_keep)
H9_dispersion<- estimateDisp(H9_keep_norm, design.matrix_H9, robust=TRUE)
H9_fited <- glmQLFit(H9_dispersion, design.matrix_H9, robust=TRUE)




contrasts_H9 <- makeContrasts(H9.KOvswt = samplesH9_KO - samplesH9_wt, 
                          H9.KIvswt = samplesH9_KI - samplesH9_wt, 
                          H9.KIvsKO = samplesH9_KI - samplesH9_KO, levels=design.matrix_H9)
contrasts_H9
##                    Contrasts
## Levels              H9.KOvswt H9.KIvswt H9.KIvsKO
##   samplesH9_KI         0         1         1
##   samplesH9_KO         1         0        -1
##   samplesH9_wt        -1        -1         0


qlf_H9.contrasts <- glmQLFTest(H9_fited, contrast=contrasts_H9)
topTags(qlf_H9.contrasts)
## Coefficient:  LR test on 2 degrees of freedom 
##            Ensembl_ID Gene_name logFC.H9.KOvswt logFC.H9.KIvswt
 ## 16448 ENSG00000187325     TAF9B     0.047290741      -8.5215156
## 2694  ENSG00000102710   SUPT20H    -0.009879572      -9.8090108
## 6131  ENSG00000129317     PUS7L     0.071584025     -13.0416646
## 7625  ENSG00000138161     CUZD1     1.072742182      -3.5320786
## 17959 ENSG00000198934    MAGEE1   -10.949547315       0.9180734

##       logFC.H9.KIvsKO   logCPM        F       PValue          FDR
## 16448       -8.568806 5.044505 955.5740 1.059827e-11 1.005574e-07
## 2694        -9.799131 5.032961 949.3324 1.093252e-11 1.005574e-07
## 6131       -13.113249 4.246053 834.5035 1.670160e-10 1.024142e-06
## 7625        -4.604821 3.376219 261.1062 4.800669e-09 2.207828e-05
## 17959       11.867621 2.659453 347.6069 6.771092e-09 2.491220e-05


 qlf_H9.KOvswt <- glmQLFTest(H9_fited, contrast=contrasts_H9[,"H9.KOvswt"])
 topTags(qlf_H9.KOvswt)
## Coefficient:  1*samplesH9_KO -1*samplesH9_wt 
##            Ensembl_ID  Gene_name      logFC    logCPM        F
## 50201 ENSG00000268658  LINC00664  -8.590411 2.3669917 404.3259
## 17959 ENSG00000198934     MAGEE1 -10.949547 2.6594526 516.2982
## 5095  ENSG00000120738       EGR1   4.377344 4.6205615 218.4309
## 15758 ENSG00000184515       BEX5  -4.938635 0.6839495 188.3181
## 27893 ENSG00000228065  LINC01515  -4.375104 0.9611039 187.7014

 ##             PValue          FDR
## 50201 4.027064e-09 5.809191e-05
## 17959 6.315711e-09 5.809191e-05
## 5095  8.535940e-08 5.107783e-04
## 15758 1.367709e-07 5.107783e-04
## 27893 1.388286e-07 5.107783e-04
  • I want to know why my result in qlf_H9.contrasts is different than qlf_H9.KOvswt? e.g the LogFC of MAGEE1 gene in the qlf_H9.contrasts and qlf_H9.KOvswt is the same but its FDR is different, and these kind of differences in FDR cause I lost lots of significant genes and had different list of genes when consider qlf_H9.contrasts or qlf_H9.KOvswt?

  • which one (Anova or pairwise) is better for my comparisons?

Is it because of my colnames of design.matrix_H9 or the Levels name of contrast_H9 are not the same as the colnames of data!

edger • 146 views
ADD COMMENTlink modified 5 months ago by Gordon Smyth39k • written 5 months ago by star10
1

This isn’t a DESeq2 question. Please don’t add unnecessary tags as it automatically emails package maintainers who are not relevant to your question.

ADD REPLYlink written 5 months ago by Michael Love25k
1

It's hard to answer your question because it isn't clear what you are asking. The anova test and the single contrast tests are different and test different hypotheses, as you already seem to understand. Obviously different tests will give different results. Why would you expect different tests to give the same FDRs?

What scientific comparison do you want to test and what do think you have tested?

ADD REPLYlink modified 5 months ago • written 5 months ago by Gordon Smyth39k

Thanks for your reply! In fact I would like to test which Genes differentially Express in H9.kI vs wt and H9.ko vs wt samples, then select common up-down regulated genes between those comparisons. I would like to know can I use Anova or I should consider single comparison, then find common Genes , as righ one?

ADD REPLYlink written 5 months ago by star10
Answer: Why the result of pairwise comparison and Anova comparison of the same contrast
0
gravatar for Gordon Smyth
5 months ago by
Gordon Smyth39k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth39k wrote:

Get lists of DE genes for KO vs wt and KI vs wt:

qlf_H9.KOvswt <- glmQLFTest(H9_fited, contrast=contrasts_H9[,"H9.KOvswt"])
KO <- topTags(qlf_H9.KOvswt, n=Inf, sort="none")$table
qlf_H9.KIvswt <- glmQLFTest(H9_fited, contrast=contrasts_H9[,"H9.KIvswt"])
KI <- topTags(qlf_H9.KIvswt, n=Inf, sort="none")$table

Find common genes:

i <- KO$FDR < 0.05 & KO$FDR < 0.05 & KO$logFC*KI$logFC > 0
sum(i)
CommonDE <- data.frame(KO[i,], KI[i,])
head(CommonDE)
ADD COMMENTlink modified 5 months ago • written 5 months ago by Gordon Smyth39k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 228 users visited in the last hour