nanopolish-methyl

介绍

封装nanopolish call-methylation

安装

需要安装Parallel::ForkManager，以实现并行计算。命令如下：

cpan install Parallel::ForkManager

或者通过cpanm安装(cpanm使用方法看这里)：

cpanm --mirror http://mirrors.163.com/cpan Parallel::ForkManager

为保证程序正确运行，需要安装其他依赖的工具诸如samtools，tabix，bgzip（只要成功安装htslib套装即可），bedtools。不必拘泥于特定版本，使用conda即可方便安装。

安装最新版的nanopolish工具，按照我的简书说明进行编译。

接下来克隆代码：

git clone https://gitee.com/wangshun1121/nanopolish-methyl.git

若有必要，找到脚本Nanopolish_Methyl.pl的第19行，提供hdf5 plugin的正确路径：

$ENV{'HDF5_PLUGIN_PATH'} = '/usr/local/hdf5/lib/plugin'; # 这是一个坑，我在简书中会介绍

最后，将脚本拷贝到nanopolish源代码的根目录下。

使用说明

使用工具之前，请先用nanopolish index将reads和fast5文件进行关联。具体操作按照简书说明进行。

工具帮助信息如下：

perl ./nanopolish/Nanopolish_Methyl.pl
  封装 nanopolish call-methylation ，实现并行化
    -r|ref              参考基因组
    -i|reads            输入ONT的长reads，务必与其对应的fast5用 nanopolish index 做好关联
    -b|bam              输入的bam
    -o|out              输出的甲基化结果文件前缀（默认results)
    -t|target           待检测的基因组靶向区间（默认整个基因组）

    -cpu                使用的总线程数（默认节点上全部48个线程）
    -j|jobs             同时并行多少个任务（默认12)
    -w|window           每一个任务，基因组区间大小（默认200000）
    -tmp                临时文件路径（默认输出文件前缀.tmp）
    -s|split-group      -s表示启用calculate_methylation_frequency.py 的 --split-groups 参数。
    -k|keep             -k表示程序运行结束，不删除临时文件
    -options            nanopolish call-methylation的其他运行参数

    -h|help             Show this message
  关于本脚本配套的一系列操作，可参考下面：
    简书：https://www.jianshu.com/p/60d620a166fa
    gitee：https://gitee.com/wangshun1121/nanopolish-methyl

例如下面代码：

perl ./nanopolish/Nanopore_Methyl.pl -r ref.fa -i reads.fq.gz -b output.bam -o results -keep -split-group -jobs 6 -options "--methylation dam"

输出结果

输出结果包括两部分内容：

results.meth_call.tsv.gz(.tbi)

这是Reads甲基化位点信息，格式如下：

chromosome      start   end     read_name       sequence        strand  log_lik_ratio   log_lik_methylated      log_lik_unmethylated    num_calling_strands     num_motifs
1       356     356     3c3f7ae8-a8c7-426d-89e4-22384531a4ee    TATACCGATGT     +       -1.63   -94.18  -92.55  1       1
1       356     356     3fb519af-eaed-45b5-9b0a-e2ce197c6882    TATACCGATGT     +       0.13    -214.03 -214.16 1       1
1       356     356     b630ba83-8c85-4dae-b85a-6364dc8171e1    TATACCGATGT     +       2.39    -137.88 -140.27 1       1
1       356     356     29e0aaaa-4817-492b-89d6-7d75dd06a33a    TATACCGATGT     +       2.58    -140.40 -142.98 1       1
1       356     356     d4745cc3-c3b9-4e1c-a2be-a296040d273c    TATACCGATGT     +       -2.65   -180.61 -177.96 1       1
1       356     356     816db5a1-7c13-4064-a2a5-e4cae8844dae    TATACCGATGT     -       4.34    -149.16 -153.50 1       1
1       356     356     e1586034-dfdb-4d97-8674-775349b81d1e    TATACCGATGT     +       -0.11   -143.00 -142.89 1       1
1       356     356     98596332-ba2e-4244-b90c-f7fc25eabcb6    TATACCGATGT     -       -4.94   -143.14 -138.20 1       1
1       356     356     f076f7c4-7e38-41d7-a850-0e0a9a4af51d    TATACCGATGT     -       -2.11   -125.62 -123.51 1       1

列的顺序与官网默认的顺序不同，我为了方便使用tabix提取特定区间的甲基化信息，将列顺序调整为类似bed的格式，即将strand信息摆在第6列。若提取特定区间的甲基化信息，可使用下面命令：

gunzip -c results.meth_call.tsv.gz|head -n 1 > meth_call.chr20:10000-20000.tsv # 表头必须要！！！
tabix -p bed results.meth_call.tsv.gz chr20:10000-20000 >> meth_call.chr20:10000-20000.tsv

得到的文件meth_call.chr20:10000-20000.tsv仍然可以用nanopolish自带的工具calculate_methylation_frequency.py进行后续分析。

results.meth_freq.tsv.gz.(tbi)

基因组上特定位点的甲基化信息，格式如下：

chromosome      start   end     num_motifs_in_group     called_sites    called_sites_methylated methylated_frequency    group_sequence  called_sites.plus       called_sites_methylated.plus    methylated_frequency.minus       called_sites.minus
1       356     356     1       10      7       0.700   TATACCGATGT     6       5       0.833   4       2       0.500
1       378     378     1       7       1       0.143   split-group     6       0       0.000   1       1       1.000
1       380     380     1       7       1       0.143   split-group     6       0       0.000   1       1       1.000
1       471     471     1       26      8       0.308   AGTGGCGTTCC     13      5       0.385   13      3       0.231
1       513     513     1       29      6       0.207   split-group     13      3       0.231   16      3       0.188
1       515     515     1       29      6       0.207   split-group     13      3       0.231   16      3       0.188

我的格式比官方多了6列，分别记录基因组正链和负链上，C碱基的甲基化水平。前面原有的几列则与官方格式一致。

顺 / nanopolish-methyl

nanopolish-methyl

介绍

安装

使用说明

输出结果

简介

发行版

贡献者

近期动态

顺 / nanopolish-methyl .gitee-modal { width: 500px !important; }

nanopolish-methyl

介绍

安装

使用说明

输出结果

简介

发行版

贡献者

近期动态

搜索帮助

顺 / nanopolish-methyl