9. Gene Set Enrichment Analysis
基因集富集分析是生物信息学研究的重要环节。
在前面的步骤之后,已经获得了一些重要的DMPs或DMRs,因此我们可能想知道这些重要DMPs或dmr中涉及的基因是否因特定的生物学术语或通路而被富集。
要实现此分析,可以使用champ.GSEA()来进行GSEA分析。
champ.GSEA()将自动提取myDMP和myDMR中的基因(无论您使用什么方法生成DMR)。此外,如果您有自己的重要的CpG列表或基因列表计算从相同的甲基化阵列(450K或EPIC)其他方法不包括在ChAMP。您还可以将它们格式化为list,并输入要执行GSEA的函数(请为您自己的CpG list或gene list中的每个元素指定一个名称,否则可能会触发错误)。champ.GSEA()会自动提取基因信息,将CpG信息转换为基因信息,然后对每个列表进行GSEA。在CpGs到基因的映射过程中,如果有多个CpGs映射到一个基因
,则该基因只会被计数一次
,以防计数过多
。
做GSEA有三种方法。
在以前的版本中,ChAMP使用了从MSigDB下载的通路信息。
然后利用Fisher精确试验计算各通路的富集状态。
经过基因富集分析,champ.GSEA()函数会自动返回p值小于adjPval的通路。
myGSEA <- champ.GSEA(beta=myNorm,DMP=myDMP[[1]], DMR=myDMR, arraytype="450K",adjPval=0.05, method="fisher")
# myDMP and myDMR could (not must) be used directly.
myGSEA2 <- champ.GSEA(beta=myNorm,DMP=myDMP[[1]], DMR=myDMR, arraytype="450K",adjPval=0.05, method="ebayes")
# myDMP and myDMR could (not must) be used directly.
注:如果要修正CpG基因数量不等的偏差,同时考虑到CpGs的显著水平,可以将方法参数设置为ebayes,采用经验bayes方法。否则,你也可以用“gometh” method or “fisher” method 来做GSEA。
10. Empirical Bayes GSEA method
正如我们之前介绍的,ChAMP现在在ChAMP .GSEA()中加入了一种叫做ebayes的新方法,这种方法与大多数其他GSEA方法不同,因为它不需要DMP或DMR信息。
因此,实际上用户可以独立于champ.GSEA()运行此方法。
我们在ChAMP中提供了ChAMP . ebaygsea()函数,用于那些希望直接从标准化的beta矩阵和表型进行GSEA的用户。
像下面的
myebayGSEA <- champ.ebayGSEA(beta=myNorm,pheno=myLoad$pd$Sample_Group,arraytype="450K")
这是一种很有前途的GSEA方法,不受基因富集的CpGs数量或CpGs显著水平的影响。
11.寻找作用通路网络中的疾病关联小网络
这货根本就没有中文译名,上边是我自己翻译的。简而言之就是,人体内的作用网络实在是太多了,几百几千吧,每一个都涉及了一系列基因。这就是过往的科学家研究出来的,比如某个通路会导致头疼,然后有几百个蛋白(及其背后的基因)都是通过共同作用导致了这一场头疼的。但是如果你头疼,不见得是所有这些基因都出了问题,而是可能其中的某一部分,甚至于只有一两个出了问题。所以你就可以基于已经存在的那个网络,再集合数据,找出哪些网络可能出问题了?或者是这些大网络中的哪些基因具体除了问题。
这个功能很重要,否则做完上边几步,用户只会知道那些基因有问题,至于他们之间有没有关系,是不是会同时作用于某些网络,就没法知道的了。但其实也不复杂,只需要自己写个程序,匹配一样基因和网络就完了,只不过数据的准备啊,洗啊,匹配啊,也是够烦的,所以这个函数就提供了全套的分析。
myEpiMod <- champ.EpiMod(beta=myNorm,pheno=myLoad$pd$Sample_Group)
library("ChAMP")
testDir=system.file("extdata",package="ChAMPdata")
myLoad <- champ.load(testDir,arraytype="450K",method ="ChAMP")
str(myLoad)
myLoad2=champ.load(testDir,arraytype="450K",method ="minfi")
str(myLoad2)
mydata=champ.import(testDir)
beta=mydata$beta
myfilter <- champ.filter(mydata$beta,pd=NULL,au)
# Or you may separate about code as champ.import(testDir) + champ.filter()
myLoad$pd
myLoad2$pd
myfilter$pd
CpG.GUI(CpG=rownames(myLoad$beta),arraytype="450K")
champ.QC()
CpG.GUI()
champ.QC(Feature.sel = "SVD") # Alternatively: QC.GUI()
QC.GUI()
myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
QC.GUI(beta=myNorm)
kk=myLoad$pd[,c(1,3)]
str(myLoad$pd)
champ.SVD(beta=myNorm,pd=kk)
myLoad$pd=myLoad$pd[,c(1,3)]
champ.SVD(beta=myNorm,pd=myLoad$pd)
champ.SVD()
myDMP <- champ.DMP(beta = myNorm,pheno=myLoad$pd$Sample_Group)
head(myDMP[[1]])
myNorm[rownames(myNorm)=="cg06822689",]
DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=myLoad$pd$Sample_Group)
# myDMP is a list now, each data frame is stored as myDMP[[1]], myDMP[[2]], myDMP[[3]]...
myDMR <- champ.DMR(beta=myNorm,pheno=myLoad$pd$Sample_Group,method="Bumphunter")
DMR.GUI()
myBlock <- champ.Block()
Block.GUI()
myGSEA <- champ.GSEA(beta=myNorm,DMP=myDMP[[1]], DMR=myDMR, arraytype="450K",adjPval=0.05, method="fisher")
myGSEA2 <- champ.GSEA(beta=myNorm,DMP=myDMP[[1]], DMR=myDMR, arraytype="450K",adjPval=0.05, method="ebayes")
# myDMP and myDMR could (not must) be used directly.
myGSEA_DMP<-champ.GSEA(beta=myNorm,DMP=myDMP[[1]], arraytype="450K",adjPval=0.05, method="fisher")
myebayGSEA <- champ.ebayGSEA(beta=myNorm,pheno=myLoad$pd$Sample_Group,arraytype="450K")
kk=myGSEA_DMP$DMP
myEpiMod <- champ.EpiMod(beta=myNorm,pheno=myLoad$pd$Sample_Group)
myebayGSEA <- champ.ebayGSEA(beta=myNorm,pheno=myLoad$pd$Sample_Group,arraytype="450K")
myEpiMod <- champ.EpiMod()
TCGA
Medat=data.table::fread("G:/小张聊科研项目/甲基化/HumanMethylation450 (1)/HumanMethylation450")
head(Medat)[,1:10]
Medat=data.frame(Medat)
row.names(Medat)=Medat[,1]
Medat2=Medat[,-1]
myLoad=Medat2
head(Medat2)
names(pd)
grou=data.frame("Sample_Name"=names(Medat2),"Sample_Group"=c(rep("C",100),rep("T",55)))
myfilter <- champ.filter(as.matrix(Medat2) ,pd=grou)
data=as.data.frame(myfilter$beta)
myfilter$beta=na.omit(myfilter$beta)
# Or you may separate about code as champ.import(testDir) + champ.filter()
CpG.GUI(CpG=rownames(myfilter$beta),arraytype="450K")
myNorm <- champ.norm(beta=myfilter$beta,arraytype="450K",cores=5)
QC.GUI(beta=myNorm)
myDMP <- champ.DMP(beta = myNorm,pheno=myfilter$pd$Sample_Group)
head(myDMP[[1]])