參考文章：
1.甲基化芯片入門學(xué)習(xí)-ChAMP包（二）
2.TCGA數(shù)據(jù)庫(kù)的癌癥甲基化芯片數(shù)據(jù)重分析
3.TCGA甲基化芯片數(shù)據(jù)質(zhì)控和過(guò)濾
4.甲基化芯片數(shù)據(jù)的差異分析
5.甲基化芯片注釋中的CpG shores, open sea 是什么
上一篇文章，寫(xiě)了如何下載以及整理我們需要的臨床樣品。得到了一個(gè)ChAMP對(duì)象送丰，里面包含了甲基化信號(hào)beta值，以及樣品的信息（腫瘤/正常）炉媒。上一篇里的過(guò)濾只是過(guò)濾了甲基化信號(hào)里的NA值，現(xiàn)在要用ChAMP包來(lái)進(jìn)行進(jìn)一步的過(guò)濾以及質(zhì)量控制昆烁。這個(gè)包的官方使用說(shuō)明PDF版：here
> library(ChAMP)
> library(dplyr)
> library(tibble)
> load('OSCC_29paired_methydata_ChAMPfiltered.Rdata')
（一）標(biāo)準(zhǔn)化前的QC
> QC = champ.QC(beta = myLoad$beta, pheno = myLoad$pd$sample_type)
[===========================]
[<<<<< ChAMP.QC START >>>>>>]
-----------------------------
champ.QC Results will be saved in ./CHAMP_QCimages/
[QC plots will be proceed with 412481 probes and 58 samples.]
<< Prepare Data Over. >>
<< plot mdsPlot Done. >>
<< Plot densityPlot Done. >>
< Dendrogram Plot Feature Selection Method >: No Selection, directly use all CpGs to calculate distance matrix.
<< Plot dendrogram Done. >>
[<<<<<< ChAMP.QC END >>>>>>>]
[===========================]
[You may want to process champ.norm() next.]
根據(jù)運(yùn)行的提示吊骤，你會(huì)在當(dāng)前文件夾里得到一個(gè)新文件夾，里面有3張QC的圖静尼。mdsPlot圖 (Multidimensional Scaling Plot)主要看看不同分組樣本是否分開(kāi)白粉；densityPlot圖，每個(gè)樣本的beta值的分布圖鼠渺，主要看看有無(wú)異常的樣本鸭巴；dendrogram，樣本的聚類圖拦盹。下面是我得到的3張圖：
腫瘤和正常組織的分群
beta值看起來(lái)沒(méi)有異常值
樣品聚類圖
（二）數(shù)據(jù)標(biāo)準(zhǔn)化（歸一化）
#這里我設(shè)置的是4個(gè)核鹃祖，跑了有一小會(huì)兒，不過(guò)電腦有些熱普舆。這一步比較費(fèi)內(nèi)存恬口。
> myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=4)
> dim(myNorm)
[1] 412481 58
#歸一化后會(huì)產(chǎn)生很多NA,看一下哪些樣品里含有NA
> num.na <- apply(myNorm,2,function(x)(sum(is.na(x))))
> table(num.na)
num.na
0 257463 259559 263148
55 1 1 1
#這里顯示55個(gè)樣品里都沒(méi)有NA校读，有3個(gè)樣品里產(chǎn)生了很多NA,接下來(lái)把這個(gè)樣品和相對(duì)應(yīng)的配對(duì)樣品刪掉
> library(stringr)
> names(num.na) = colnames(myNorm)
> dt = names(num.na[num.na>0])
> dn = str_replace(dt,"-01","-11")
> keep = setdiff(colnames(myNorm),c(dt,dn)) #只取colnames(myNorm)里的元素
> myNorm = myNorm[,keep]
> dim(myNorm)
[1] 412481 52
> pd = myLoad$pd
> pd <- pd[pd$sample_submitter_id %in% colnames(myNorm),]
> dim(pd)
[1] 52 4
（三）標(biāo)準(zhǔn)化后的QC
你可以使用ChAMP包的champ.QC功能再次驗(yàn)證一下是否需要去掉異常值（代碼見(jiàn)上面）；你也可以根據(jù)參考文章2和3祖能，進(jìn)行常規(guī)的PCA和熱圖的檢查地熄，如下；
（1）主成分分析
> library(FactoMineR)
> library(factoextra)
> dat <- t(myNorm)
> group_list=pd$sample_type
> table(group_list)
group_list
Primary Tumor Solid Tissue Normal
26 26
> dat.pca <- PCA(dat, graph = FALSE)
> fviz_pca_ind(dat.pca,
geom.ind = "point",
col.ind = group_list,
addEllipses = TRUE,
legend.title = "Groups")
（2）熱圖
> cg=names(tail(sort(apply(myNorm,1,sd)),1000))
> library(pheatmap)
> ac=data.frame(group=group_list)
> rownames(ac)=colnames(myNorm)
> pheatmap(myNorm[cg,],show_colnames =F,show_rownames = F,
annotation_col=ac)
（3）相關(guān)關(guān)系矩陣熱圖
組內(nèi)的樣本的相似性應(yīng)該是要高于組間的：
>pheatmap::pheatmap(cor(myNorm[cg,]),
annotation_col = ac,
show_rownames = F,
show_colnames = F)
根據(jù)上面的QC結(jié)果看芯杀，倒是沒(méi)有什么離群值，所以不用去除樣品了雅潭。
（四）去除異常值
上面的關(guān)系熱圖看出有兩個(gè)樣品不太好揭厚，要把它們?nèi)サ簦约八鼈儗?duì)應(yīng)的配對(duì)樣品：
#去除異常值（5971扶供，6953）
> pn = c("TCGA-CV-5971-01","TCGA-CV-6953-11")
> drop = str_sub(colnames(myNorm),1,12) %in% str_sub(pn,1,12)
> table(drop)
#drop
#FALSE TRUE
# 48 4
> myNorm = myNorm[,!drop]
> dim(myNorm)
#[1] 412481 48
> pd = pd[!(pd$case_submitter_id %in% str_sub(pn,1,12)),]
> save(pd,myNorm,file = "OSCC_26paired_after_ChAMP_norm.Rdata")
（五）差異分析
（1）甲基化位點(diǎn)差異分析
比較常見(jiàn)的差異甲基化分析是DMP（Differential Methylation Probe）筛圆，用于找出差異的甲基化位點(diǎn)，ChAMP包里包含分析過(guò)程椿浓；然后根據(jù)你的分組信息太援，獲得差異甲基化位點(diǎn)；最后用DMP.GUI查看下結(jié)果扳碍。并且分析得到的結(jié)果數(shù)據(jù)里自帶了注釋提岔，可以拿去做富集分析。
champ.DMP() 實(shí)現(xiàn)了 limma包中利用linear model計(jì)算差異甲基化位點(diǎn)的p-value笋敞。
> library(ChAMP)
> library(tibble)
> x <- pd$sample_type
#先把樣品分類的名稱改一下碱蒙，中間有空格下面運(yùn)行會(huì)報(bào)錯(cuò)
> x[which(x=="Primary Tumor")] <- "Tumor"
> x[which(x=="Solid Tissue Normal")] <- "Normal"
> pd$sample_type = x
> group_list <- pd$sample_type
#利用ChAMP找出差異甲基化位點(diǎn)
> myDMP <- champ.DMP(beta = myNorm,pheno=group_list)
> head(myDMP$Tumor_to_Normal) #差異甲基化位點(diǎn)結(jié)果，你會(huì)發(fā)現(xiàn)這個(gè)結(jié)果里基因名都給你標(biāo)好了
logFC AveExpr t P.Value adj.P.Val B Tumor_AVG
cg12825070 0.6511792 0.3994829 33.31000 2.972469e-34 1.226087e-28 67.77294 0.7250725
cg09385093 0.6304708 0.3598621 30.87567 8.916709e-33 1.838987e-27 64.47550 0.6750975
cg14416371 0.6703583 0.3644433 29.56171 6.196156e-32 8.519323e-27 62.58603 0.6996225
cg15811515 0.6126625 0.4631829 28.27707 4.451309e-31 4.370805e-26 60.65743 0.7695142
cg18233405 0.4102125 0.2374704 28.06051 6.255518e-31 4.370805e-26 60.32398 0.4425767
cg07792478 0.6299208 0.4680037 28.05022 6.357828e-31 4.370805e-26 60.30808 0.7829642
Normal_AVG deltaBeta CHR MAPINFO Strand Type gene feature cgi
cg12825070 0.07389333 -0.6511792 5 148033708 F I HTR4 1stExon island
cg09385093 0.04462667 -0.6304708 2 119607885 F I IGR island
cg14416371 0.02926417 -0.6703583 11 43602847 R I MIR129-2 TSS200 island
cg15811515 0.15685167 -0.6126625 16 31580989 F I CSDAP1 TSS200 island
cg18233405 0.03236417 -0.4102125 8 98290148 F I TSPYL5 1stExon island
cg07792478 0.15304333 -0.6299208 8 65290484 F I MIR124-2 TSS1500 island
feat.cgi UCSC_CpG_Islands_Name DHS Enhancer
cg12825070 1stExon-island chr5:148033472-148034080 NA NA
cg09385093 IGR-island chr2:119607378-119607910 NA NA
cg14416371 TSS200-island chr11:43602545-43603215 NA TRUE
cg15811515 TSS200-island chr16:31580559-31581023 NA NA
cg18233405 1stExon-island chr8:98289604-98290404 TRUE NA
cg07792478 TSS1500-island chr8:65290108-65290946 NA NA
Phantom Probe_SNPs Probe_SNPs_10
cg12825070
cg09385093 rs71422594
cg14416371
cg15811515
cg18233405 high-CpG:98359303-98359424
cg07792478
> df_DMP <- myDMP$Tumor_to_Normal
#取基因名不為空白的行
> df_DMP=df_DMP[df_DMP$gene!="",]
> logFC_t <- 0.45
> P.Value_t <- 10^-15
> df_DMP$change <- ifelse(df_DMP$adj.P.Val < P.Value_t & abs(df_DMP$logFC) > logFC_t,
ifelse(df_DMP$logFC > logFC_t ,'UP','DOWN'),'NOT')
> table(df_DMP$change)
DOWN NOT UP
258 104980 619
> save(df_DMP,file = "OSCC_DF_methy.Rdata")
參考文章:甲基化芯片數(shù)據(jù)的差異分析
此處設(shè)置的logFC和p值的閾值是與原文一致的夯巷，由于甲基化的beta值不同于表達(dá)量赛惩，實(shí)際上用logFC也不是很對(duì)。推薦用deltabeta值替代logFC趁餐，就是甲基化信號(hào)值的差值喷兼。
火山圖
> library(dplyr)
> library(ggplot2)
> dat = rownames_to_column(df_DMP)
> for_label <- dat%>% head(3)
> p <- ggplot(data = dat,
aes(x = logFC,
y = -log10(adj.P.Val))) +
geom_point(alpha=0.4, size=3.5,
aes(color=change)) +
ylab("-log10(Pvalue)")+
scale_color_manual(values=c("green", "grey","red"))+
geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
theme_bw()
> p
熱圖
> cg <- rownames(df_DMP[df_DMP$change != "NOT",])
> plot_matrix <- myNorm[cg,]
> nnotation_col <- data.frame(Sample=pd$sample_type)
> rownames(annotation_col) <- colnames(plot_matrix)
> ann_colors = list(Sample = c(Normal="#4DAF4A", Tumor="#E41A1C"))
> library(pheatmap)
> pheatmap(plot_matrix,show_colnames = T,
annotation_col = annotation_col,
border_color=NA,
color = colorRampPalette(colors = c("white","navy"))(50),
annotation_colors = ann_colors,show_rownames = F)
除了上面的常規(guī)的火山圖和熱圖，你還可以用ChAMP包里GUI來(lái)查看結(jié)果后雷，里面也包含了可視化結(jié)果：
#你可以利用下面一行代碼查看分析結(jié)果季惯，也可以用head()來(lái)看
> DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=group_list)
比如熱圖：
或者探針位置的比列：
可以看到大部分的探針落在了CpG島上。Shores 指的是位于CpG island上下游2kb 以內(nèi)的區(qū)域喷面；Shelves指的是位于shores 上下游2kb以內(nèi)的區(qū)域星瘾；open sea指的是除了CpG islands, CpG shores, CpG shelves之外的其他區(qū)域。
看探針在基因組上的分布：
整個(gè)基因組分為Promoter, Body, 3UTR, Intergenic 4種區(qū)域惧辈，其中Promoter區(qū)又分為TSS200, TSS1500, 5UTR, ‘1stExon’ 4個(gè)區(qū)域琳状。
（2）甲基化片段差異分析
DMRs(Differential Methylation Regions)主要指一連串的CpG都出現(xiàn)明顯差異的區(qū)域。
#還可以利用ChAMP找出差異甲基化片段（DMR,Differential Methylation Regions),但目前只支持兩組樣品進(jìn)行分析
> myDMR <- champ.DMR(beta=myNorm,pheno=group_list,method="Bumphunter")
> head(myDMR$BumphunterDMR)
seqnames start end width strand value area cluster indexStart
DMR_1 chr6 28602543 28603437 894 * 3.663917 128.2371 168092 78333
DMR_2 chr6 29520965 29521803 838 * 3.373399 121.4424 168249 78769
DMR_3 chr6 31695903 31698226 2323 * 1.841876 112.3544 169009 81605
DMR_4 chr13 78492568 78494171 1603 * 2.666064 111.9747 57435 28411
DMR_5 chr7 94284258 94285887 1629 * 2.138407 102.6435 184882 91745
DMR_6 chr11 14993378 14995770 2392 * 2.385313 102.5685 33948 16158
indexEnd L clusterL p.value fwer p.valueArea fwerArea
DMR_1 78367 35 35 0 0 0 0
DMR_2 78804 36 39 0 0 0 0
DMR_3 81665 61 61 0 0 0 0
DMR_4 28452 42 49 0 0 0 0
DMR_5 91792 48 90 0 0 0 0
DMR_6 16200 43 43 0 0 0 0
使用GUI查看結(jié)果：
> DMR.GUI(DMR=myDMR,beta=myNorm,pheno=group_list)
當(dāng)然你選擇不同的pvalue-cutoff會(huì)得到不一樣的表格和圖
（六）GSEA分析
有了差異甲基化位點(diǎn)/區(qū)域盒齿，就可以進(jìn)行GSEA分析了念逞。這里ChAMP包里也有自己的GSEA分析的功能困食。
參考文章：生信修煉手冊(cè)：ChAMP分析甲基化芯片數(shù)據(jù)-GSEA篇
champ.GSEA默認(rèn)對(duì)差異CpG位點(diǎn)和差異甲基化區(qū)域?qū)?yīng)的基因做富集分析。
芯片中翎承，我們直接得到的是差異的探針或者差異的區(qū)域硕盹，首先需要將探針或者區(qū)域映射到基因上，在映射的過(guò)程中叨咖，我們必須考慮到一個(gè)因素瘩例，基因和探針之間的關(guān)系。大部分的基因具有多個(gè)CpG位點(diǎn)甸各，會(huì)對(duì)應(yīng)多個(gè)探針I(yè)D垛贤。比如基因A上有50個(gè)差異CpG位點(diǎn)，基因B上具有2個(gè)CpG位點(diǎn)趣倾，很明顯二者是有很大差別的聘惦，如果只考慮基因，那么A和B就是相同的儒恋，都是差異探針對(duì)應(yīng)的基因善绎。所以需要將基因?qū)?yīng)的CpG位點(diǎn)考慮進(jìn)來(lái)，gometh算法將基因覆蓋的CpG位點(diǎn)個(gè)數(shù)作為基因的長(zhǎng)度诫尽，用來(lái)矯正P值禀酱。
> myGSEA <- champ.GSEA(beta=myNorm,
DMP=myDMP[[1]],
DMR=myDMR,
CpGlist=NULL,
Genelist=NULL,
pheno=group_list,
method="fisher",
arraytype="450K",
Rplot=TRUE,
adjPval=0.01)
#查看結(jié)果
> str(myGSEA)
List of 2
$ DMP:'data.frame': 604 obs. of 9 variables:
..$ Gene_List（MSigDB數(shù)據(jù)庫(kù)中定義的基因集合）: Factor w/ 8328 levels "3_5_CYCLIC_NUCLEOTIDE_PHOSPHODIESTERASE_ACTIVITY",..: 355 364 3728 822 361 3645 3547 8328 350 3733 ...
..$ nList （每個(gè)基因集包括的基因個(gè)數(shù)）: num [1:604] 1118 1038 590 2485 652 ...
..$ nRep （基因集合的基因與所有輸入的gene list 中overlap的基因個(gè)數(shù)） : num [1:604] 1020 900 564 2266 586 ...
..$ fRep （overla的基因的比例） : num [1:604] 0.912 0.867 0.956 0.912 0.899 ...
..$ nOVLAP（位于該基因集合下的基因與輸入的gene list 中overlap的個(gè)數(shù)）: int [1:604] 963 844 547 1973 561 942 719 1233 865 320 ...
..$ OR （費(fèi)舍爾精確檢驗(yàn)的odds ratio）: num [1:604] 5.32 4.7 9.92 2.16 6.92 ...
..$ P.value : num [1:604] 4.10e-54 9.41e-44 2.87e-42 1.36e-37 3.09e-37 ...
..$ adjPval : num [1:604] 3.41e-50 3.92e-40 7.97e-39 2.83e-34 5.15e-34 ...
..$ Genes （個(gè)數(shù)和nOVLAP相同）: Factor w/ 8309 levels "A1BG YWHAZ PPP1R13B INPP5D SOS1 CYTH3 YWHAE AKT1 RPS6KA2 AKT3 YWHAH SHC1 RPS6KB1 FOXO1 GSK3B AKT2 CD55 IGFBP1 G"| __truncated__,..: 4531 48 6038 4767 4543 5201 2142 2988 7084 3334 ...
$ DMR:'data.frame': 387 obs. of 9 variables:
..$ Gene_List: Factor w/ 8328 levels "3_5_CYCLIC_NUCLEOTIDE_PHOSPHODIESTERASE_ACTIVITY",..: 364 355 350 361 3728 3733 3642 2417 3721 3652 ...
..$ nList : num [1:387] 1038 1118 1062 652 590 ...
..$ nRep : num [1:387] 900 1020 956 586 564 327 258 364 98 335 ...
..$ fRep : num [1:387] 0.867 0.912 0.9 0.899 0.956 ...
..$ nOVLAP : int [1:387] 267 282 269 196 155 101 87 96 48 88 ...
..$ OR : num [1:387] 7.16 6.54 6.64 8.1 5.88 ...
..$ P.value : num [1:387] 6.45e-106 7.58e-104 1.08e-100 8.40e-87 3.60e-55 ...
..$ adjPval : num [1:387] 5.37e-102 3.16e-100 2.98e-97 1.75e-83 6.00e-52 ...
..$ Genes : Factor w/ 5211 levels "ABCB1","ABCB1 ABCB4",..: 18 1977 4419 4871 241 2451 3292 4981 4109 2908 ...
NOTE：對(duì)于fRep < 0.6 的基因集合，會(huì)過(guò)濾掉牧嫉。官方是這樣解釋的 remove lists with less than 60% representation on array比勉。
最后保存：
> save(myDMP,myDMR,myGSEA,file = "ChAMP_DMP_DMR_GSEA.Rdata")
