TCGA甲基化数据质量控制和差异分析（使用ChAMP包）

参考文章：
1.甲基化芯片入门学习-ChAMP包（二）
2.TCGA数据库的癌症甲基化芯片数据重分析
3.TCGA甲基化芯片数据质控和过滤
4.甲基化芯片数据的差异分析
5.甲基化芯片注释中的CpG shores, open sea 是什么

上一篇文章，写了如何下载以及整理我们需要的临床样品。得到了一个ChAMP对象，里面包含了甲基化信号beta值，以及样品的信息（肿瘤/正常）。上一篇里的过滤只是过滤了甲基化信号里的NA值，现在要用ChAMP包来进行进一步的过滤以及质量控制。这个包的官方使用说明PDF版：here

> library(ChAMP)
> library(dplyr)
> library(tibble)
> load('OSCC_29paired_methydata_ChAMPfiltered.Rdata')

（一）标准化前的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.]

根据运行的提示，你会在当前文件夹里得到一个新文件夹，里面有3张QC的图。mdsPlot图 (Multidimensional Scaling Plot)主要看看不同分组样本是否分开；densityPlot图，每个样本的beta值的分布图，主要看看有无异常的样本；dendrogram，样本的聚类图。下面是我得到的3张图：

肿瘤和正常组织的分群

beta值看起来没有异常值

样品聚类图

（二）数据标准化（归一化）

#这里我设置的是4个核，跑了有一小会儿，不过电脑有些热。这一步比较费内存。
> myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=4)
> dim(myNorm)
[1] 412481     58
#归一化后会产生很多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个样品里都没有NA，有3个样品里产生了很多NA,接下来把这个样品和相对应的配对样品删掉
> 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

（三）标准化后的QC

你可以使用ChAMP包的champ.QC功能再次验证一下是否需要去掉异常值（代码见上面）；你也可以根据参考文章2和3，进行常规的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）相关关系矩阵热图

组内的样本的相似性应该是要高于组间的：

>pheatmap::pheatmap(cor(myNorm[cg,]),
                   annotation_col = ac,
                   show_rownames = F,
                   show_colnames = F)

根据上面的QC结果看，倒是没有什么离群值，所以不用去除样品了。

（四）去除异常值

上面的关系热图看出有两个样品不太好，要把它们去掉，以及它们对应的配对样品：

#去除异常值（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）甲基化位点差异分析

比较常见的差异甲基化分析是DMP（Differential Methylation Probe），用于找出差异的甲基化位点，ChAMP包里包含分析过程；然后根据你的分组信息，获得差异甲基化位点；最后用DMP.GUI查看下结果。并且分析得到的结果数据里自带了注释，可以拿去做富集分析。

参考文章：ChAMP 包分析甲基化数据

champ.DMP() 实现了 limma包中利用linear model计算差异甲基化位点的p-value。

> library(ChAMP)
> library(tibble)
> x <- pd$sample_type
#先把样品分类的名称改一下，中间有空格下面运行会报错
> x[which(x=="Primary Tumor")] <- "Tumor"
> x[which(x=="Solid Tissue Normal")] <- "Normal"
> pd$sample_type = x
> group_list <- pd$sample_type

#利用ChAMP找出差异甲基化位点
> myDMP <- champ.DMP(beta = myNorm,pheno=group_list)
> head(myDMP$Tumor_to_Normal) #差异甲基化位点结果，你会发现这个结果里基因名都给你标好了
  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")

参考文章:甲基化芯片数据的差异分析
此处设置的logFC和p值的阈值是与原文一致的，由于甲基化的beta值不同于表达量，实际上用logFC也不是很对。推荐用deltabeta值替代logFC，就是甲基化信号值的差值。

火山图

> 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)

除了上面的常规的火山图和热图，你还可以用ChAMP包里GUI来查看结果，里面也包含了可视化结果：

#你可以利用下面一行代码查看分析结果，也可以用head()来看
> DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=group_list)

比如热图：

或者探针位置的比列：

可以看到大部分的探针落在了CpG岛上。Shores 指的是位于CpG island上下游2kb 以内的区域；Shelves指的是位于shores 上下游2kb以内的区域；open sea指的是除了CpG islands, CpG shores, CpG shelves之外的其他区域。

看探针在基因组上的分布：

整个基因组分为Promoter, Body, 3UTR, Intergenic 4种区域，其中Promoter区又分为TSS200, TSS1500, 5UTR, ‘1stExon’ 4个区域。

（2）甲基化片段差异分析

DMRs(Differential Methylation Regions)主要指一连串的CpG都出现明显差异的区域。

#还可以利用ChAMP找出差异甲基化片段（DMR,Differential Methylation Regions),但目前只支持两组样品进行分析
> 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查看结果：

> DMR.GUI(DMR=myDMR,beta=myNorm,pheno=group_list)

当然你选择不同的pvalue-cutoff会得到不一样的表格和图

（六）GSEA分析

有了差异甲基化位点/区域，就可以进行GSEA分析了。这里ChAMP包里也有自己的GSEA分析的功能。
参考文章：生信修炼手册：ChAMP分析甲基化芯片数据-GSEA篇

champ.GSEA默认对差异CpG位点和差异甲基化区域对应的基因做富集分析。
芯片中，我们直接得到的是差异的探针或者差异的区域，首先需要将探针或者区域映射到基因上，在映射的过程中，我们必须考虑到一个因素，基因和探针之间的关系。大部分的基因具有多个CpG位点，会对应多个探针ID。比如基因A上有50个差异CpG位点，基因B上具有2个CpG位点，很明显二者是有很大差别的，如果只考虑基因，那么A和B就是相同的，都是差异探针对应的基因。所以需要将基因对应的CpG位点考虑进来，gometh算法将基因覆盖的CpG位点个数作为基因的长度，用来矫正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)

#查看结果
> str(myGSEA)
List of 2
 $ DMP:'data.frame':    604 obs. of  9 variables:
  ..$ Gene_List（MSigDB数据库中定义的基因集合）: Factor w/ 8328 levels "3_5_CYCLIC_NUCLEOTIDE_PHOSPHODIESTERASE_ACTIVITY",..: 355 364 3728 822 361 3645 3547 8328 350 3733 ...
  ..$ nList （每个基因集包括的基因个数）: num [1:604] 1118 1038 590 2485 652 ...
  ..$ nRep （基因集合的基因与所有输入的gene list 中overlap的基因个数） : 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的个数）: int [1:604] 963 844 547 1973 561 942 719 1233 865 320 ...
  ..$ OR （费舍尔精确检验的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 （个数和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：对于fRep < 0.6 的基因集合，会过滤掉。官方是这样解释的 remove lists with less than 60% representation on array。

最后保存：

> save(myDMP,myDMR,myGSEA,file = "ChAMP_DMP_DMR_GSEA.Rdata")