9. Gene Set Enrichment Analysis
基因集富集分析是生物信息學研究的重要環(huán)節(jié)甚疟。
在前面的步驟之后，已經(jīng)獲得了一些重要的DMPs或DMRs，因此我們可能想知道這些重要DMPs或dmr中涉及的基因是否因特定的生物學術語或通路而被富集揽祥。
要實現(xiàn)此分析，可以使用champ.GSEA()來進行GSEA分析府树。
champ.GSEA()將自動提取myDMP和myDMR中的基因(無論您使用什么方法生成DMR)。此外奄侠，如果您有自己的重要的CpG列表或基因列表計算從相同的甲基化陣列(450K或EPIC)其他方法不包括在ChAMP载矿。您還可以將它們格式化為list垄潮，并輸入要執(zhí)行GSEA的函數(shù)(請為您自己的CpG list或gene list中的每個元素指定一個名稱闷盔，否則可能會觸發(fā)錯誤)。champ.GSEA()會自動提取基因信息牡整，將CpG信息轉換為基因信息，然后對每個列表進行GSEA逃贝。在CpGs到基因的映射過程中，如果有多個CpGs映射到一個基因沐扳，則該基因只會被計數(shù)一次，以防計數(shù)過多歉闰。
做GSEA有三種方法。
在以前的版本中和敬，ChAMP使用了從MSigDB下載的通路信息戏阅。
然后利用Fisher精確試驗計算各通路的富集狀態(tài)昼弟。
經(jīng)過基因富集分析奕筐，champ.GSEA()函數(shù)會自動返回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基因數(shù)量不等的偏差芭逝，同時考慮到CpGs的顯著水平渊胸，可以將方法參數(shù)設置為ebayes旬盯，采用經(jīng)驗bayes方法翎猛。否則，你也可以用“gometh” method or “fisher” method 來做GSEA萨咳。
10. Empirical Bayes GSEA method
正如我們之前介紹的疫稿，ChAMP現(xiàn)在在ChAMP .GSEA()中加入了一種叫做ebayes的新方法培他，這種方法與大多數(shù)其他GSEA方法不同而克，因為它不需要DMP或DMR信息。
因此员萍，實際上用戶可以獨立于champ.GSEA()運行此方法。
我們在ChAMP中提供了ChAMP . ebaygsea()函數(shù)螃壤，用于那些希望直接從標準化的beta矩陣和表型進行GSEA的用戶。
像下面的
myebayGSEA <- champ.ebayGSEA(beta=myNorm,pheno=myLoad$pd$Sample_Group,arraytype="450K")
這是一種很有前途的GSEA方法奸晴，不受基因富集的CpGs數(shù)量或CpGs顯著水平的影響。
11.尋找作用通路網(wǎng)絡中的疾病關聯(lián)小網(wǎng)絡
這貨根本就沒有中文譯名逮光，上邊是我自己翻譯的墩划。簡而言之就是涕刚，人體內的作用網(wǎng)絡實在是太多了乙帮，幾百幾千吧，每一個都涉及了一系列基因察净。這就是過往的科學家研究出來的，比如某個通路會導致頭疼锈至，然后有幾百個蛋白（及其背后的基因）都是通過共同作用導致了這一場頭疼的。但是如果你頭疼裹赴，不見得是所有這些基因都出了問題诀浪，而是可能其中的某一部分延都，甚至于只有一兩個出了問題雷猪。所以你就可以基于已經(jīng)存在的那個網(wǎng)絡晰房，再集合數(shù)據(jù)，找出哪些網(wǎng)絡可能出問題了殊者？或者是這些大網(wǎng)絡中的哪些基因具體除了問題。
這個功能很重要摔刁，否則做完上邊幾步，用戶只會知道那些基因有問題共屈，至于他們之間有沒有關系，是不是會同時作用于某些網(wǎng)絡拗引，就沒法知道的了。但其實也不復雜矾削，只需要自己寫個程序，匹配一樣基因和網(wǎng)絡就完了欲间，只不過數(shù)據(jù)的準備啊挡逼，洗啊括改，匹配啊家坎，也是夠煩的，所以這個函數(shù)就提供了全套的分析虱疏。
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]])
