可视化EWAS结果你可以这样做

今天给大家介绍一款用于甲基化关联分析（EWAS）的R包---coMET。coMET能够绘制CpG位点，DNA甲基化相关性图谱，同时可以添加ENSEMBL基因结构、ENCODE基因信息以及用户可以自定义的相关基因组注释信息。除了对甲基化数据进行可视化之外，该工具包还可扩展至GWAS和eQTL等数据。官网http://bioconductor.org/packages/release/bioc/html/coMET.html。

第一步：安装coMET并导入

if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("coMET")

此外，还需安装如下三个R包

install.packages("psych")
install.packages("corrplot")
install.packages("colortools")

library(coMET)

第二步：了解数据格式

1）info file文件格式

注意：此文件是必需的，而且必须是带有标题的表格格式。通过参数mydata.file指定，列数及其类型由config.file文件中的参数mydata.format定义。

第一列的名称必须以字母开头。该信息文件可以包含/不包含beta值（例如DNA甲基化水平）的CPG位点列表。

如果该信息文件是一个位点列表文件，那么它必须具有如下所示的4列，并且标题的顺序与Figure 1中的文件相同。beta值可以是第5列（可选），可以是数值（正值或负值）或仅是方向符号（+/-）。

Figure 1：info file位点列表文件格式

info file位点文件每一列所代表的信息如下所示：

(1) Name of omic feature

(2) Name of chromosome

(3) Position of omic feature

(4) P-value of omic feature

(5) Direction of association related to this omic feature. This can be the sign or an actual value of association effect size.【可选】

可通过如下代码查看info file示例位点文件：

extdata <- system.file("extdata", package="coMET",mustWork=TRUE)
infofile <- file.path(extdata, "cyp1b1_infofile.txt")
data_info <-read.csv(infofile, header = TRUE, sep = "\t", quote = "")
head(data_info)

如果该信息文件是一个基于区域的列表文件，则基于区域的信息文件必须具有5列（Figure 2），并按此顺序排列标题及相对应的各列信息。beta值或方向信息可以包括在第6列中（可选）。

Figure 2：info file区域列表文件格式

info file区域文件每一列所代表的信息如下所示：

(1) Name of omic feature

(2) Name of chromosome

(3) Start position of omic feature

(4) End position of omic feature

(5) P-value of omic feature

(6) Direction of association related to this omic feature. This can be the sign or an actual value of association effect size.【可选】

可通过如下代码查看info file示例区域文件：

extdata <- system.file("extdata", package="coMET",mustWork=TRUE)
infoexp <- file.path(extdata, "cyp1b1_infofile_exprGene_region.txt")
data_infoexp <-read.csv(infoexp, header = TRUE, sep = "\t", quote = "")
head(data_infoexp)

2）correlation matrix文件格式

注意：此文件是必需的，而且必须是带标题的表格格式。通过参数cormatrix.file指定，其数据格式通过config.file文件中的参数cormatrix.format指定（cormatrix.format= "cormatrix"/"raw"/"raw_rev"）。

该文件格式可以采用参数cormatrix.format中描述的3种格式：

(1)cormatrix格式：用户提供的预先计算好的相关系数矩阵；矩阵维数：【CpG_number】*【CpG_number】。文件中需要将CpG位点/区域按位置进行升序排列，并具有标题，其标题为CPG位点/区域的名称；

(2)raw格式：相关系数可以通过Spearman、Pearson、Kendall三种方法中的一种来计算，需要通过参数cormatrix.method指定。矩阵维数：【sample_size】*【CpG_number】。文件具有标题，其标题为CpG位点/区域的名称；该格式文件见Figure 3.

Figure 3：raw格式文件

可通过如下代码查看raw格式示例文件：

extdata <- system.file("extdata", package="coMET",mustWork=TRUE)
corfile <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt")
data_cor <-read.csv(corfile, header = TRUE, sep = "\t", quote = "")
data_cor[1:6,1:6]

(3)raw_rev格式：相关系数可以通过Spearman、Pearson、Kendall三种方法中的一种来计算，需要通过参数cormatrix.method指定。矩阵维数：【CpG_number】*【sample_size】。文件具有行名和标题，行名为CpG位点/区域的名称，标题为样本名称；

3）extra info file文件格式

注意：此文件是可选文件，如果提供，则必须是带标题的表格格式。第一列的名称必须以字母开头。通过参数mydata.large.file指定，其格式由参数mydata.large.format定义。可以使用多个extra info file文件，每个文件以逗号分隔。这可以是另一种类型info file文件（例如 expression or replication data），并且应遵循与info file文件相同的规则。

4）annotation file文件格式

注意：此文件是可选文件，通过biofeat.user.file指定，其格式为BED, GTF, and GFF3。

5）config.file文件

注意：文件中的每一行都是一个选项。选项的名称是小写字母，选项名称和其值通过“=”分隔，没有空格。如果有多个值，例如用于选项list.tracks或用于附加数据的选项，则需要用一个“逗号”将这些值分隔开，但不能有空格。如果想对绘图配置进行更改，可以下载配置文件，对其进行改动，然后将其上传到R中。该文件内容见Figure 4。该文件在绘图过程中占据非常重要的地位，可以通过其设置图表的标题，相关颜色，显著性P值，以及之前提到的cormatrix.format，data.format等。

Figure 4：config.file文件信息

可通过如下代码查看config.file文件信息

extdata <- system.file("extdata", package="coMET",mustWork=TRUE)
configfile <- file.path(extdata, "config_cyp1b1_zoom_4webserver_Grch38.txt")
data_config <-read.csv(configfile, quote = "", sep="\t", header=FALSE)
data_config

第三步：通过comet.web绘图

comet.web是预先定制的功能，也就是说很多参数及相关注释信息已经被设置好，允许用户快速可视化相关甲基化信息，注意此模块只能可视化人类的结果。绘图代码如下（Figure 5）：

extdata <- system.file("extdata", package="coMET",mustWork=TRUE)
myinfofile <- file.path(extdata, "cyp1b1_infofile_Grch38.txt")
myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region_Grch38.txt")
mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt")
configfile <- file.path(extdata, "config_cyp1b1_zoom_4webserver_Grch38.txt")
comet.web(config.file=configfile, mydata.file=myinfofile,cormatrix.file=mycorrelation ,mydata.large.file=myexpressfile,print.image=FALSE,verbose=FALSE)

Figure 5：comet.web绘制甲基化信息相关图示

Figure 5主要包含三部分的信息：

(1)曼哈顿图（上面部分）：显示相关EWAS关联信号的强度和范围。涵盖相关P值及Beta值信息。

(2)annotation tracks（中间部分）：该部分主要通过相关数据库信息，展示已知的gene结构或者SNP等相关信息，主要通过list.tracks参数设置（例：该图中的相关设置为list.tracks=geneENSEMBL,ChromHMM,DNAse,RegENSEMBL）。目前所支持的相关数据库信息如下：

geneENSEMBL, CGI, ChromHMM, DNAse, RegENSEMBL, SNP, transcriptENSEMBL, SNPstoma, SNPstru, SNPstrustoma, GAD, ClinVar, GeneReviews, GWAS, ClinVarCNV, GCcontent, genesUCSC, xenogenesUCSC, metQTL, eQTL, BindingMotifs-Biomart, chromHMM_RoadMap, miRNATargetRegionsBiomart, Other-RegulatoryRegionsBiomart, RegulatoryEvidenceBiomart and segmentalDupsUCSC。

(3)热图（下面部分）：该部分主要利用热图（heatmap）显示基因组区域中选定的CpG位点之间的相关性。

第三步：通过comet绘图

comet与comet.web类似，但是其支持用户自定义相关注释信息。绘图实例代码如下：

library("Gviz")
library("rtracklayer")
extdata <- system.file("extdata", package="coMET",mustWork=TRUE)
configfile <- file.path(extdata, "config_cyp1b1_zoom_4comet.txt")
myinfofile <- file.path(extdata, "cyp1b1_infofile.txt")
myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region.txt")
mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt")
chrom <- "chr2"
start <- 38290160
end <- 38303219
gen <- "hg19"
strand <- "*"
BROWSER.SESSION="UCSC"
mySession <- browserSession(BROWSER.SESSION)
genome(mySession) <- gen
genetrack <-genes_ENSEMBL(gen,chrom,start,end,showId=TRUE)
snptrack <- snpBiomart_ENSEMBL(gen,chrom, start, end,dataset="hsapiens_snp_som",showId=FALSE)
cpgIstrack <- cpgIslands_UCSC(gen,chrom,start,end)
prombedFilePath <- file.path(extdata, "/RoadMap/regions_prom_E063.bed")
promRMtrackE063<- DNaseI_RoadMap(gen,chrom,start, end, prombedFilePath,
featureDisplay='promotor', type_stacking="squish")
bedFilePath <- file.path(extdata, "RoadMap/E063_15_coreMarks_mnemonics.bed")
chromHMM_RoadMapAllE063 <- chromHMM_RoadMap(gen,chrom,start, end,bedFilePath, featureDisplay = "all",colorcase='roadmap15' )
listgviz <- list(genetrack,snptrack,cpgIstrack,promRMtrackE063,chromHMM_RoadMapAllE063)
comet(config.file=configfile, mydata.file=myinfofile, mydata.type="file",cormatrix.file=mycorrelation, cormatrix.type="listfile",mydata.large.file=myexpressfile, mydata.large.type="listfile",tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE)

Figure 6：comet绘制甲基化信息相关图示

Figure 7：comet绘制甲基化信息默认参数

Figure 8：网页版绘图界面

可以看出图中Figure 6展示的信息与Figure 5类似，但是添加了更多自定义的相关内容。与comet.web函数相比，在绘制图片时comet主要通过修改Figure 7中的相关参数达到个性化绘制图片的目的。此外，该R包还非常人性化的提供了网页版的绘制工具（http://epigen.kcl.ac.uk/comet/upload.html），绘图页面见Figure 8。通过点击式的方式就可以绘制出高端大气上档次的图片，为不会编程的同学们提供大展身手的机会。使用的参数不同或注释文件不同，所绘制出的图片就会达到不同的效果，感兴趣的小伙伴赶快动手试试吧！！！