SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Question about model design in edgeR user guide case study example 4.5 MBWatson Bioinformatics 2 10-08-2014 11:25 AM
EdgeR or GLM model design question younko Bioinformatics 10 09-01-2014 12:22 AM
Help with multi factorial model fitting in edgeR/DEseq HTnoob Bioinformatics 0 01-28-2014 10:00 PM
EdgeR cpm values after batch effects danielfortin86 Bioinformatics 4 08-25-2013 10:12 PM
edgeR choose the best model of GLM eszter.ari RNA Sequencing 1 02-23-2013 07:49 PM

Reply
 
Thread Tools
Old 01-13-2015, 04:34 PM   #1
rhcr56
Junior Member
 
Location: Columbia, MO

Join Date: Aug 2011
Posts: 7
Default edgeR GLM model - factorial design with batch effects

Hello,

I am analyzing some RNA-Seq data with edgeR, and I have a question about how to set up the GLM model.matrix. My experiment is a 2X2 factorial design where both WT and KO animals are subjected to either of two treatments. The 4 groups are WT-Oil, WT-E2, KO-Oil, and KO-E2. I am interested in determining which genes are DE in the WT-E2 group (our phenotype after treating with E2 is rescued in similarly treated KO animals).

Here is the code that I used

> data.frame(Sample=colnames(a),genotype2,trt2)
Sample genotype2 trt2
1 X1_WT_Oil_140710 WT OIL
2 X5_WT_Oil_140710 WT OIL
3 X1_WT_Oil_141109 WT OIL
4 X2_WT_Oil_141109 WT OIL
5 X2_WT_E2_140710 WT E2
6 X6_WT_E2_140710 WT E2
7 X3_WT_E2_141109 WT E2
8 X4_WT_E2_141109 WT E2
9 X3_KO_Oil_140710 KO OIL
10 X8_KO_Oil_140710 KO OIL
11 X4_KO_E2_140710 KO E2
12 X5_KO_E2_141109 KO E2
13 X6_KO_E2_141109 KO E2
14 X7_KO_E2_141109 KO E2

> design2 <- model.matrix(~0+genotype2:trt2)
> design2
genotype2KO:trt2E2 genotype2WT:trt2E2 genotype2KO:trt2OIL genotype2WT:trt2OIL
1 0 0 0 1
2 0 0 0 1
3 0 0 0 1
4 0 0 0 1
5 0 1 0 0
6 0 1 0 0
7 0 1 0 0
8 0 1 0 0
9 0 0 1 0
10 0 0 1 0
11 1 0 0 0
12 1 0 0 0
13 1 0 0 0
14 1 0 0 0
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$genotype2
[1] "contr.treatment"

attr(,"contrasts")$trt2
[1] "contr.treatment"

> b <- estimateGLMCommonDisp(b,design2,verbose=TRUE)
Disp = 0.20378 , BCV = 0.4514
> fit <- glmFit(b,design2)
> lrt <- glmLRT(fit,contrast=c(-1,3,-1,-1))

This worked just fine, but now I would like to reanalyze this data taking batch effects into account. The last part of the sample IDs (140710 and 141109) indicate the two separate dates in which the samples were collected/sequenced. The new data frame would look like this:

> data.frame(Sample=colnames(a),genotype2,trt2,chip2)
Sample genotype2 trt2 chip2
1 X1_WT_Oil_140710 WT OIL 1
2 X5_WT_Oil_140710 WT OIL 1
3 X1_WT_Oil_141109 WT OIL 2
4 X2_WT_Oil_141109 WT OIL 2
5 X2_WT_E2_140710 WT E2 1
6 X6_WT_E2_140710 WT E2 1
7 X3_WT_E2_141109 WT E2 2
8 X4_WT_E2_141109 WT E2 2
9 X3_KO_Oil_140710 KO OIL 1
10 X8_KO_Oil_140710 KO OIL 1
11 X4_KO_E2_140710 KO E2 1
12 X5_KO_E2_141109 KO E2 2
13 X6_KO_E2_141109 KO E2 2
14 X7_KO_E2_141109 KO E2 2

The manual doesn't describe how to set up the model.matrix when there are two factors plus batch effects, and I'm struggling setting this up.

Could somebody guide me as to how to properly define the model.matrix as well as properly define the correct contrasts in the glmLRT function? Thanks!

Last edited by rhcr56; 01-13-2015 at 04:38 PM. Reason: highlight important changes
rhcr56 is offline   Reply With Quote
Old 01-14-2015, 01:02 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

I never got the fixation with treating a simple factorial design as if it's not. You can skip the contrasts in your first design and simply use:
Code:
design2 <- model.matrix(~genotype2*trt2)
You then don't need any contrasts:
Code:
fit <- glmFit(b,design2)
lrt <- glmLRT(fit) #You'll get the last coefficient, which will be the interaction, by default.
BTW, just make sure that KO and OIL are the base levels of their respective factors.

This also lends itself nicely to adding a batch (not that it's significantly different using the "collapse the groups and make contrasts" example, where you will still contrast the same groups):
Code:
design = model.matrix(~batch+genotype2*trt2)
dpryan is offline   Reply With Quote
Old 01-14-2015, 07:45 AM   #3
rhcr56
Junior Member
 
Location: Columbia, MO

Join Date: Aug 2011
Posts: 7
Default

Thanks for the reply. This was my first inclination, however I'm confused about the last coefficient. The code you provided was used below

> design3 <- model.matrix(~genotype2*trt2)
> fit <- glmFit(b,design3)
> lrt2 <- glmLRT(fit)
> topTags(lrt2)
Coefficient: genotype2WT:trt2OIL
logFC logCPM LR PValue FDR
1600002K03Rik -8.515051 2.068219 37.67803 8.343907e-10 8.445590e-06
Mettl21c 9.151695 2.339162 36.03935 1.933725e-09 8.445590e-06
Il10ra 9.719101 2.009320 35.98865 1.984707e-09 8.445590e-06
Mpo 5.802613 3.754349 34.41918 4.443238e-09 1.265384e-05
Cd177 7.962947 2.954040 34.20659 4.956072e-09 1.265384e-05
Hspa1a -4.576308 4.852461 33.63713 6.641285e-09 1.413044e-05
Tpte 7.864546 2.116211 32.57592 1.146264e-08 2.090457e-05
Bace2 -10.589050 3.180579 31.81367 1.696943e-08 2.707897e-05
Hspa1b -4.314080 4.794562 30.50166 3.335810e-08 4.731661e-05
Fuom -8.057760 1.793066 30.06869 4.170094e-08 5.323542e-05

If I am interpreting this correctly, this will identify genes that are DE in the WT-Oil group compared to every other group. I am interested in those DE genes in the WT-E2 group. Do I need to use a different coefficient? When I used the contrasts, I got the appropriate test.

> topTags(lrt,n=343)
Coefficient: -1*genotype2KO:trt2E2 3*genotype2WT:trt2E2 -1*genotype2KO:trt2OIL -1*genotype2WT:trt2OIL

Thanks
rhcr56 is offline   Reply With Quote
Old 01-14-2015, 07:52 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

I'll reiterate my previous "BTW":
Quote:
BTW, just make sure that KO and OIL are the base levels of their respective factors.
dpryan is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:43 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO