最新R语言因子实验设计和解释案例分析报告 附代码数据.docx
《最新R语言因子实验设计和解释案例分析报告 附代码数据.docx》由会员分享,可在线阅读,更多相关《最新R语言因子实验设计和解释案例分析报告 附代码数据.docx(22页珍藏版)》请在冰豆网上搜索。
![最新R语言因子实验设计和解释案例分析报告 附代码数据.docx](https://file1.bdocx.com/fileroot1/2022-11/24/ceb8cc4b-00a6-4a48-940e-94a309ed7d0c/ceb8cc4b-00a6-4a48-940e-94a309ed7d0c1.gif)
最新R语言因子实验设计和解释案例分析报告附代码数据
R语言因子实验设计和解释案例分析报告
·示例1:
两组比较
·示例2:
多个组
实例3:
两个条件,两个基因型,一个交互项
o野生型治疗效果(主效应)。
o突变体治疗的效果
o没有治疗的突变型和野生型之间有什么区别?
o通过治疗,突变型和野生型有什么区别?
o基因型的不同反应(相互作用项)
实例4:
两个条件,三个基因型
o基因型I的条件效应(主效应)
o基因型III的条件效应。
o基因型II的条件效应。
o在条件A下III与II的影响
o基因型III与基因型I的条件效应的相互作用项
○基因型III与基因型II的条件效应的相互作用项。
为了允许iDEP中的复杂模型(http:
//ge-lab.org/idep/),我尝试了解如何构建事实模型,并从DESeq2中提取期望的结果。
以下是基于DESeq2中resutls()函数的帮助文档,以及MikeLove对用户提问的回答。
我想要做的一个重点是,当研究设计涉及多个因素时(参见上面关于基因型+治疗实例的图),结果的解释是棘手的。
与R中的回归分析类似,分类因素的参考水平构成了我们的分歧的基础。
然而,默认情况下,它们是按字母顺序确定的。
选择每个因素的参考水平是至关重要的。
否则你的系数可能会有所不同,这取决于你如何进入DESeq2的实验设计。
这可以通过R中的relevel()函数完成。
参考级别是构成有意义比较基础的因素的基线级别。
在野生型与突变型实验中,“野生型”是参考水平。
在治疗与未治疗,参考水平显然是未经处理的。
例3中的更多细节。
例1:
两组比较
首先制作一些示例数据。
library(DESeq2)
dds<-makeExampleDESeqDataSet(n=10000,m=6)
assay(dds)[1:
10,]
##sample1sample2sample3sample4sample5sample6
##gene164111213
##gene291223131428
##gene35812117317811897
##gene4040383
##gene527369812
##gene648835382113
##gene7365061524422
##gene86816141819
##gene9214266419198157166
##gene1020121612162
这是一个非常简单的实验设计,有两个条件。
colData(dds)
##DataFramewith6rowsand1column
##condition
##
##sample1A
##sample2A
##sample3A
##sample4B
##sample5B
##sample6B
dds<-DESeq(dds)
resultsNames(dds)
##[1]"Intercept""condition_B_vs_A"
这显示了可用的结果。
请注意,默认情况下,R会根据字母顺序为因素选择一个参考级别。
这里A是参考水平。
折叠变化定义为B与A比较。
要更改参考级别,请尝试使用“同一个”()函数。
res<-results(dds,contrast=c("condition","B","A"))
res<-res[order(res$padj),]
library(knitr)
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene9056
360.168909
-2.045379
0.0000000
0.0001366
gene3087
43.897516
-2.203303
0.0000173
0.0858143
gene3763
72.409877
-1.834787
0.0000434
0.1434712
gene2054
322.494963
1.537408
0.0000681
0.1689463
gene4617
6.227415
6.125238
0.0002019
0.4008408
如果我们想用B作为控制,并用B作为基线定义倍数变化。
那我们可以这样做:
res<-results(dds,contrast=c("condition","A","B"))
ix=which.min(res$padj)
res<-res[order(res$padj),]
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene9056
360.168909
2.045379
0.0000000
0.0001366
gene3087
43.897516
2.203303
0.0000173
0.0858143
gene3763
72.409877
1.834787
0.0000434
0.1434712
gene2054
322.494963
-1.537408
0.0000681
0.1689463
gene4617
6.227415
-6.125238
0.0002019
0.4008408
正如你所看到的,折叠的方向是完全相反的。
这里我们展示最重要的基因。
barplot(assay(dds)[ix,],las=2,main=rownames(dds)[ix])
示例2:
多个组
假设我们有三个组A,B和C.
dds<-makeExampleDESeqDataSet(n=100,m=6)
dds$condition<-factor(c("A","A","B","B","C","C"))
dds<-DESeq(dds)
res=results(dds,contrast=c("condition","C","A"))
res<-res[order(res$padj),]
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene2
3.634986
-5.101773
0.0348679
0.5515088
gene20
4.678176
-4.490982
0.0445664
0.5515088
gene34
56.068672
-1.462155
0.0167820
0.5515088
gene35
537.847175
-1.177240
0.0087913
0.5515088
gene41
93.967810
1.064734
0.0412034
0.5515088
Example3:
twoconditions,twogenotypes,withaninteractionterm
Herewehavetwogenotypes,wild-type(WT),andmutant(MU).Twoconditions,control(Ctrl)andtreated(Trt).Weareinterestedintheresponsesofbothwild-typeandmutanttotreatment.Wearealsointerestedinthedifferencesinresponsebetweengenotypes,whichiscapturedbytheinteractionterminlinearmodels.
First,weconstructexampledata.Notethatwechangedsamplenamesfrom“sample1”to“Wt_Ctrl_1”,accordingtothetwofactors.
dds<-makeExampleDESeqDataSet(n=10000,m=12)
dds$condition<-factor(c(rep("Ctrl",6),rep("Trt",6)))
dds$genotype<-factor(rep(rep(c("WT","MU"),each=3),2))
colnames(dds)<-paste(as.character(dds$genotype),as.character(dds$condition),rownames(colData(dds)),sep="_")
colnames(dds)=gsub("sample","",colnames(dds))
kable(assay(dds)[1:
5,])
WT_Ctrl_1
WT_Ctrl_2
WT_Ctrl_3
MU_Ctrl_4
MU_Ctrl_5
MU_Ctrl_6
WT_Trt_7
WT_Trt_8
WT_Trt_9
MU_Trt_10
MU_Trt_11
MU_Trt_12
gene1
81
47
135
72
81
80
76
57
52
49
64
57
gene2
70
77
94
54
54
36
51
66
57
51
47
67
gene3
16
37
8
28
4
36
14
10
7
11
40
15
gene4
2
5
1
3
9
2
1
5
0
0
14
8
gene5
0
0
1
0
5
0
0
0
0
0
0
0
kable(colData(dds))
condition
genotype
WT_Ctrl_1
Ctrl
WT
WT_Ctrl_2
Ctrl
WT
WT_Ctrl_3
Ctrl
WT
MU_Ctrl_4
Ctrl
MU
MU_Ctrl_5
Ctrl
MU
MU_Ctrl_6
Ctrl
MU
WT_Trt_7
Trt
WT
WT_Trt_8
Trt
WT
WT_Trt_9
Trt
WT
MU_Trt_10
Trt
MU
MU_Trt_11
Trt
MU
MU_Trt_12
Trt
MU
Checkreferencelevels:
dds$condition
##[1]CtrlCtrlCtrlCtrlCtrlCtrlTrtTrtTrtTrtTrtTrt
##Levels:
CtrlTrt
Asyoucouldsee,“Ctrl”apearedfirstinthe2ndline,indicatingitisthereferencelevelforfactorcondition,aswecanexpectbasedonalphabeticalorder.Thisiswhatwewantandwedonotneedtodoanything.
dds$genotype
##[1]WTWTWTMUMUMUWTWTWTMUMUMU
##Levels:
MUWT
But“Mu”isthereferencelevelforgenotype,whichiswillgiveusresultsdifficulttointerpret.Weneedtochangeit.
dds$genotype=relevel(dds$genotype,"WT")
dds$genotype
##[1]WTWTWTMUMUMUWTWTWTMUMUMU
##Levels:
WTMU
Setupthemodel,andrunDESeq2:
design(dds)<-~genotype+condition+genotype:
condition
dds<-DESeq(dds)
resultsNames(dds)
##[1]"Intercept""genotype_MU_vs_WT"
##[3]"condition_Trt_vs_Ctrl""genotypeMU.conditionTrt"
Below,wearegoingtousethecombinationofthedifferentresults(“genotype_MU_vs_WT”,“condition_Trt_vs_Ctrl”,“genotypeMU.conditionTrt”)toderivebiologicallymeaningfulcomparisons.
Theeffectoftreatmentinwild-type(themaineffect).
res=results(dds,contrast=c("condition","Trt","Ctrl"))
ix=which.min(res$padj)#mostsignificant
res<-res[order(res$padj),]#sort
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene275
25.38836
2.607896
0.0000588
0.2684851
gene2744
101.02954
-1.670837
0.0001099
0.2684851
gene5441
58.67469
1.758921
0.0001261
0.2684851
gene7021
74.29359
1.610853
0.0001345
0.2684851
gene7795
326.43308
-1.624323
0.0000407
0.2684851
ThisisforWT,treatedcomparedwithuntreated.NotethatWTisnotmentioned,becauseitisthereferencelevel.Inotherwords,thisisthedifferencebetweensamplesNo.7-9,comparedwithsamplesNo.1-3.
Hereweshowthemostsignificantgene.
barplot(assay(dds)[ix,],las=2,main=rownames(dds)[ix])
Theeffectoftreatmentinmutant
Thisis,bydefinition,themaineffect plus theinteractionterm(theextraconditioneffectingenotypeMutantcomparedtogenotypeWT).
res<-results(dds,list(c("condition_Trt_vs_Ctrl","genotypeMU.conditionTrt")))
ix=which.min(res$padj)#mostsignificant
res<-res[order(res$padj),]#sort
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene5102
18.69017
4.757057
1.60e-06
0.0156910
gene7367
170.69259
1.481016
1.19e-05
0.0396834
gene8351
27.77227
3.033622
9.80e-06
0.0396834
gene8034
53.34342
-1.841023
1.62e-05
0.0403724
gene7272
92.82351
-1.414967
6.70e-05
0.1125808
Thismeasurestheeffectoftreatmentinmutant.Inotherwords,samplesNo.10-12comparedwithsamplesNo.4-6.
Hereweshowthemostsignificantgene,whichisdownregulatedexpressedinsamples10-12,thansamples4-6,asexpected.
barplot(assay(dds)[ix,],las=2,main=rownames(dds)[ix])
Whatisthedifferencebetweenmutantandwild-typewithouttreatment?
AsCtrlisthereferencelevel,wecanjustretrievethe“genotype_MU_vs_WT”.
res=results(dds,contrast=c("genotype","MU","WT"))
ix=which.min(res$padj)#mostsignificant
res<-res[order(res$padj),]#sort
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene5102
18.69017
-4.714519
0.0000020
0.0195594
gene2289
45.04711
1.763615
0.0000333
0.1218120
gene3388
122.85582
-1.529323
0.0000366
0.1218120
gene915
39.03250
2.104003
0.0001133
0.2374845
gene8351
27.77227
-2.653182
0.0001190
0.2374845
Inotherwords,thisisthesamplesNo.4-6comparedwithNo.1-3.Hereweshowthemostsignificantgene.
barplot(assay(dds)[ix,],las=2,main=rownames(dds)[ix])
Withtreatment,whatisthedifferencebetweenmutantandwild-type?
res=results(dds,list(c("genotype_MU_vs_WT","genotypeMU.conditionTrt")))
ix=which.min(res$padj)#mostsignificant
res<-res[order(res$padj),]#sort
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene4127
50.03677
1.925094
0.0000091
0.0910878
gene3879
34.67704
-2.133203
0.0000283
0.0940387
gene4799
24.11059
-2.459345
0.0000197
0.0940387
gene5441
58.67469
-1.744769
0.0001418
0.2357640
gene5926
161.16737
1.524163
0.0001339
0.2357640
ThisgivesusthedifferencebetweengenotypeMUandWT,underconditionTrt.Inotherwords,thisisthesamplessNo.10-12comparedwithsamples7-9.
Hereweshowthemostsignificantgene.
barplot(assay(dds)[ix,],las=2,main=rownames(dds)[ix])
Thedifferentresponseingenotypes(interactionterm)
Istheeffectoftreatment different acrossgenotypes?
Thisistheinteractionterm.
res=results(dds,name="genotypeMU.conditionTrt")
ix=which.min(res$padj)#mostsignificant
res<-res[order(res$padj),]#sort
kable(res[1:
5,-(3:
4)])
baseMean
log2FoldChange
pvalue
padj
gene5102
18.69017
7.395771
0.0000000
0.0000341
gene5441
58.67469
-3.297673
0.0000004
0.0019104
gene8034
53.34342
-2.564030
0.0000205
0.0682558
gene8858
14.68650
4.703936
0.0000297
0.0742066
gene4799
24.11059
-2.990101
0.0001234
0.2463639
Herewesho
wthemost
significantgene.
barplot(assay(dds)[ix,],las=2,main=rownames(dds)[ix])
Thedifferentresponseingenotypes(interactionterm)
Istheeffectoftreatment different acrossgenotypes?
Thisistheinteractionterm.
res=results(dds,name="genotypeMU.conditionTrt")
ix=which.min(res$padj)#mostsignifican