甲基化数据可视化软件seqmonk不能预测CpG island

近期一直在做表观遗传学分析的开发，正好借机会分享一些方法和心得。

本次研究的物种是山羊（Capra hircus，vARS1），采用的是全基因组甲基化测序（WGBS）。

在分析过程中遇到2个难题：

1、目前山羊物种没有已知注释的CpG island信息，所以无法分析N_shore、N_shelf、island、S_shore及S_shelf的甲基化水平。

2、由于是WGBS，无法做差异甲基化CpG位点（DMC，Differentially methylated CpG sites ）分析。

现有方法是针对甲基化芯片数据分析的，采用illumina human methylation 450的甲基化芯片，数据来源是TCGA。而且芯片数据是提供N_shore、N_shelf、island、S_shore及S_shelf等具体结果的，如果利用WGBS分析数据，需要根据注释文件得到CpG island区域，自己写程序处理，得到shelf、shore区域的结果。

图 1.

The HumanMethylation450 BeadChip offers broad coverage across gene

regions, as well as CpG islands/CPG island regions, shelves, and shores for

the most comprehensive view of methylation state.

针对甲基化芯片数据分析DMC和DMR的软件是COHCAP，而COHCAP软件只能处理芯片和靶向测序数据。

所以我们迫切想找到一款软件可以分析山羊的CpG island。

通过查阅文献，我们找到seqmonk软件，其能够对甲基化数据可视化分析，和IGV差不多，而IGV主要是可视化DNA数据，通过bam看位点、断裂点等附近reads覆盖情况。

我们的假设是可以通过seqmonk找到CpG island。下面简单介绍下seqmonk的使用说明：

seqmonk

下载：http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/

安装：

linux：java环境变量配置1.8，环境变量配置：

chmod 755 seqmonk

export JAVA_HOME=/share/public/software/jdk1.8.0_11

export JRE_HOME=$/jre

export CLASSPATH=.:$/lib:$/lib

export PATH=$/bin:$PATH

## open X11

./seqmonk

windows：java 1.8（电脑是64bit）、R

注意事项：

1、电脑需要有java 1.8版本。如果是64位电脑需要下载x64的jdk，官网没有，去万能的百度搜吧。

2、需要安装R软件，我安装的最新R 3.5版本。但是安装R依赖包的时候出现报错。没关系，退出seqmonk再重新打开就可以了。

图 2. seqmonk软件在windows安装R依赖包出现报错

运行：

如果在linux上运行seqmonk会占很多资源，windows操作相对方便，但需要提前下载好参考基因组和gff注释文件。

新建参考基因组，导入genome和gff。但导入时出现问题：

1、chr书写不规范。山羊参考基因组没有具体chr1这样的染色体名称而是NC_022293.1，同样gff的染色体名seqmonk也是报错的。

2、scaffold太多。vARS1版本的山羊基因组共7w+个scaffold，所以染色体数目太多导致seqmonk报错，报错信息是java out of index。

图 3. seqmonk加载山羊基因组时java报错信息

发现问题后给seqmonk软件的负责人Simon Andrews写邮件。团队反馈很及时，他们立刻建立了seqmonk山羊的参考基因组，调整基因组和注释文件的染色体名称为1、2、3……X、Y，并去掉scaffold染色体，保留1-29号常染色体的基因组信息。

图 4. seqmonk在自己的服务器上及时更新了山羊基因组信息，便于用户导入标准基因组

找CG island（CGI）：

导入genome和gff后并没有找到CGI信息，导入bismark甲基化结果也不行。于是用小鼠参考基因组做测试，因为小鼠物种能够在UCSC上找到CGI信息。导入小鼠基因组、gff和cpgIslandExt.txt后能够看到CGI的区域信息了。

图 5. 以小鼠为例，只有导入CpG island注释文件才能看到CGI信息

结论：

通过以上分析，我们的结论是seqmonk需要导入CpG island的annotation信息，仅可视化，不可以通过seqmonk预测CpG island。

并且给seqmonk团队写邮件也证实了我的结论。simon andrews回复如下：

For CpG islands, there isn’t a way to infer these from within seqmonk. You would either need to find an existing annotation source for that, or run one of the detection programs against the genome sequence to find them for you.

不过他提到可以用EMBOSS这个软件试试。simon人真的超nice，立刻圈粉有木有。不过他善意的提醒一句，检测准确性因物种而异。

The one which most sites seem to use is newcpgreport from the EMBOSS suite (http://emboss.sourceforge.net/). You should try to validate though that standard CpG island detection thresholds work well in your species - the detection programs were designed to work well in humans, and don’t always detect correctly in species where the CpG islands aren’t as strong as they are there.

所以今天介绍的seqmonk并没有解决找CpG岛的问题，但后续我的想法是：

预测CpG island方法：EMBOSS

DMC方法：应用WGBS的methykit

如何利用WGBS数据找CpG island、DMC？大家有什么好方法欢迎后台留言交流哟~期待与你碰撞出知识的小火花

参考文献：

1、

Identification of DNAmethylation associated gene signatures in endometrial cancer via integratedanalysis of DNA methylation and gene expression systematically. J Gynecol Oncol. 2017 Nov; 28(6): e83.

2、

DNA methylation changesinendometrium and correlationwith gene expression during thetransition frompre-receptive to receptive phase. Scientific Reportsvolume 7, Article number: 3916 (2017)

3、

http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/

4、

Genome-Wide Epigenetic Characterization of Tissues from Three Germ Layers Isolated from Sheep Fetuses. Frontiers in Genetics [04 Sep 2017, 8:115]

我是yyt，爱做白日梦但又付诸行动的双子座

有科研梦想，做项目、看paper

信息分析工程师一枚

专注人类遗传病、肿瘤及表观遗传学分析等

喜欢碎碎念，分享职场小白心得