coMETR包可视化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 <- 38290160end <- 38303219gen <- "hg19"strand <- "*"BROWSER.SESSION="UCSC"mySession <- browserSession(BROWSER.SESSION)genome(mySession) <- gengenetrack <-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。通过点击式的方式就可以绘制出高端大气上档次的图片,为不会编程的同学们提供大展身手的机会。使用的参数不同或注释文件不同,所绘制出的图片就会达到不同的效果,感兴趣的小伙伴赶快动手试试吧!!!

生物信息学

评估肿瘤纯度的方法(三): 基于拷贝数变异 ABSOLUTE和DoAbsolute

2020-8-28 4:05:50

生物信息学

visNetworkR包:交互式网络可视化

2020-8-28 4:10:48

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索