當我們得到差異的探針或者差異的甲基化區(qū)域之后裸燎，通常都會分析這些差異區(qū)域?qū)?yīng)的基因是否在特定功能上有富集。在ChAMP中疙描，通過champ.GSEA函數(shù)來實現(xiàn)功能富集分析馆揉。
用法示例：
myNorm <- champ.norm()
myDMP <- champ.DMP()
myDMR <- champ.DMR()
myGSEA <- champ.GSEA()
在ChAMP中嗡贺，提供了兩種富集分析的方法:
fisher
gometh
champ.GSEA默認對差異CpG位點和差異甲基化區(qū)域?qū)?yīng)的基因做富集分析隐解，采用的方式是Fisher exact test， 分析的是Gene Set 來自MSigDB诫睬。
富集分析早已經(jīng)是研究基因功能的常用工具之一了煞茫，那么對于甲基化芯片的富集分析和傳統(tǒng)的富集分析有沒有不一樣的地方呢？
Gene Set摄凡，叫做基因集合续徽，本質(zhì)上是一系列具有相同功能的基因構(gòu)成的集合，比如某一條代謝通路亲澡，在該代謝通路中有很多的基因钦扭，位于同一條pathway下的基因就構(gòu)成了一個基因集合。
基因集合中最基本的元素是一個一個的基因床绪，而芯片中客情，我們直接得到的是差異的探針或者差異的區(qū)域，首先需要將探針或者區(qū)域映射到基因上会涎，在映射的過程中裹匙，我們必須考慮到一個因素瑞凑，基因和探針之間的關(guān)系末秃。大部分的基因具有多個CpG位點，會對應(yīng)多個探針I(yè)D籽御。比如基因A上有50個差異CpG位點练慕，基因B上具有2個CpG位點，很明顯二者是有很大差別的技掏，如果只考慮基因铃将，那么A和B就是相同的，都是差異探針對應(yīng)的基因哑梳。所以需要將基因?qū)?yīng)的CpG位點考慮進來劲阎，gometh算法將基因覆蓋的CpG位點個數(shù)作為基因的長度，用來矯正P值鸠真。
GSEA 結(jié)果如下：
str(myGSEA)
List of 2
$ DMP:’data.frame’: ? ?666 obs. of ?9 variables:
?..$ Gene_List: Factor w/ 8338 levels “3_5_CYCLIC_NUCLEOTIDE_PHOSPHODIESTERASE_ACTIVITY”,..: 355 822 1359 7228 3732 14 364 350 2564 3656 …
?..$ nList ? ?: num [1:666] 1118 2485 1972 1426 102 …
?..$ nRep ? ? : num [1:666] 1026 2267 1768 1172 98 …
?..$ fRep ? ? : num [1:666] 0.918 0.912 0.897 0.822 0.961 …
?..$ nOVLAP ? : int [1:666] 190 329 268 194 39 251 157 161 28 165 …
?..$ OR ? ? ? : num [1:666] 2.37 1.81 1.88 2.05 6.59 …
?..$ P.value ?: num [1:666] 2.55e-21 4.41e-18 4.35e-17 2.76e-16 4.91e-16 …
?..$ adjPval ?: num [1:666] 2.13e-17 1.84e-14 1.21e-13 5.76e-13 8.18e-13 …
?..$ Genes ? ?: Factor w/ 6617 levels “A2M”,”A2M ACVRL1 ACVR1 IL12RB2 TNFRSF1B IL2RA”,..: 1604 5352 2139 2554 4011 5150 2711 519 2708 6157 …
$ DMR:’data.frame’: ? ?115 obs. of ?9 variables:
?..$ Gene_List: Factor w/ 8338 levels “3_5_CYCLIC_NUCLEOTIDE_PHOSPHODIESTERASE_ACTIVITY”,..: 350 7089 246 2061 355 1206 2564 361 3653 7228 …
?..$ nList ? ?: num [1:115] 1062 148 212 201 1118 …
?..$ nRep ? ? : num [1:115] 964 137 194 183 1026 …
?..$ fRep ? ? : num [1:115] 0.908 0.926 0.915 0.91 0.918 …
?..$ nOVLAP ? : int [1:115] 44 21 23 22 43 22 15 29 19 37 …
?..$ OR ? ? ? : num [1:115] 7.69 25.32 19.04 19.22 6.94 …
?..$ P.value ?: num [1:115] 1.77e-21 4.04e-21 1.99e-20 1.04e-19 1.43e-19 …
?..$ adjPval ?: num [1:115] 1.47e-17 1.69e-17 5.53e-17 2.16e-16 2.39e-16 …
?..$ Genes ? ?: Factor w/ 1432 levels “ADCYAP1”,”ADCYAP1 BOLL VAX1 NKX2-3 PRDM14 DBX1 ALX1 PRRT1 RFX4 HOXD4 TNXB LBX1”,..: 446 985 348 984 239 117 387 1321 39 1294 ..
默認對DMP和DMR對應(yīng)的基因都是富集分析悯仙，所以結(jié)果是一個長度為2的列表，第一個列表是DMP富集分析的結(jié)果吠卷，第二個列表是DMR富集分析的結(jié)果锡垄，每個富集結(jié)果是一個data.frame對象。
每列的含義如下：
*Gene_listMSigDB數(shù)據(jù)庫中定義的基因集合
nList
每個基因集合包括的基因個數(shù)nRep
基因集合的基因與所有輸入的gene list 中overlap的基因個數(shù)fRep
overla的基因的比例nOVLAP
位于該基因集合下的基因與輸入的gene list 中overlap的個數(shù)OR
費舍爾精確檢驗的odds ratio-
Pvalue
單尾fisher exact test檢驗的p值祭隔，具體代碼如下listPV.v_2 <- t(apply(fisher.lm_2,1,function(x) unlist(fisher.test(matrix(x,2,2),alternative=”greater”)[c(1,3)])))
adjPval
多重假設(shè)檢驗校正之后的P值货岭，默認采用”BH”方法Genes
gene symbol， 個數(shù)和nOVLAP相同
需要注意一點，對于fRep < 0.6 的基因集合千贯，會過濾掉屯仗。官方是這樣解釋的 remove lists with less than 60% representation on array。
ChAMP中提供的富集分析并不是我們常說的GO/KEGG 富集分析丈牢，有很多的R包祭钉，比如clusterProfiler， 都可以用來做GO/KEGG富集分析己沛，但是都不會考慮基因?qū)?yīng)的CpG位點個數(shù)慌核，這也算是一個小的遺憾。如果要做GO/KEGG 富集分析申尼，同時又要考慮CpG位點個數(shù)垮卓，就必須考慮`missMethyl`包。
