程序安装部署
目前，由于github仓库被墙，单纯按照官网安装说明会存在各种问题。而使用conda安装，则只能安装 v0.4.0版本，而并非最新的v0.13.3版本。所以这个软件若能过成功安装，还是需要各种办法绕过墙。我已经在gitee上克隆了nanopolish及其依赖的各种包的镜像，程序编译依照如下脚本进行：
# 从我复制的gitee镜像拷贝代码
git clone https://gitee.com/wangshun1121/nanopolish.git
# 克隆nanopolish需要的各种库，并回滚到指定版本
cd ./nanopolish
# fast5 @ 18d6e34
rm -rf fast5
git clone https://gitee.com/huangyaolong/fast5.git
cd fast5
git reset --hard 18d6e346cadb9990ff8287e54447f90e76dd79bc # fast5 的 18d6e34 commit
cd ..
# htslib @ 3dc96c5
rm -rf htslib
git clone https://gitee.com/mirrors_samtools/htslib.git
cd htslib
git reset --hard 3dc96c5699a26c5a46e3ccca4092b3801e1874f4 # htslib 的 3dc96c5 commit
cd ..
# minimap2 @ d2de282
rm -rf minimap2
git clone https://gitee.com/W-zhennan/minimap2.git
cd minimap2
git reset --hard d2de282d21e9b621dc5f6cef293cb60b93cc553d # minimap2 的 d2de282 commit
cd ..
# slow5lib @ 3680e17
rm -rf slow5lib
git clone https://gitee.com/wangshun1121/slow5lib.git
cd slow5lib
git reset --hard 3680e17360366c8cbffb66eb1e07426756665454 # slow5lib 的 3680e17 commit
cd ..
# hdf5-1.8.14 这个文件下载可能比较慢
wget -c https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.14/src/hdf5-1.8.14.tar.gz
# eigen-3.3.7
wget -c https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2
# 编译，编译完成会在根目录出现可执行文件nanopolish
make
# ont-vbz-hdf-plugin，这个不部署，后面运行 nanopolish call-methylation时会报错
# 按照 https://github.com/nanoporetech/vbz_compression/issues/5#issuecomment-562476721说的办
# 拷贝这个plugin，以及环境变量永久化，可能需要高权限
wget -c https://github.com/nanoporetech/vbz_compression/releases/download/v1.0.1/ont-vbz-hdf-plugin-1.0.1-Linux-x86_64.tar.gz
tar xzvf ont-vbz-hdf-plugin-1.0.1-Linux-x86_64.tar.gz
cp -r ./ont-vbz-hdf-plugin-1.0.1-Linux/usr/local/hdf5/ /usr/local/.
echo "export HDF5_PLUGIN_PATH=/usr/local/hdf5/lib/plugin" >> /etc/profile # 永久化修改环境变量
这样，经过我的测试，后面一系列的分析能够参照官方甲基化分析说明正常走下去。
这里还有一点要说明：最新版本nanopolish是根据纳米孔R9版本建库试剂盒构建的模型，若测序使用的早期版本的试剂盒，则需要按官方说明下载早期版本的工具。截止到本文发布的时间，已经推出了R10试剂盒，但仅仅提供试用，尚未大规模商用。
甲基化分析
首先，待分析的数据存放于./BMK220118-AT448-01N0001中，数据目录如下：
./BMK220118-AT448-01N0001
├── cell_1
│ ├── fast5_0.fast5
│ ├── fast5_1.fast5
│ ├── fast5_2.fast5
│ ├── fast5_3.fast5
...
│ ├── fast5_19.fast5
│ ├── fastq_0.fastq
│ ├── fastq_1.fastq
│ ├── fastq_2.fastq
│ ├── fastq_3.fastq
...
│ ├── fastq_19.fastq
├── cell_2
│ ├── fast5_0.fast5
│ ├── fast5_1.fast5
│ ├── fast5_2.fast5
│ ├── fast5_3.fast5
...
│ ├── fast5_18.fast5
│ ├── fastq_0.fastq
│ ├── fastq_1.fastq
│ ├── fastq_2.fastq
│ ├── fastq_3.fastq
...
│ ├── fastq_18.fastq
├── cell_3
│ ├── fast5_0.fast5
│ ├── fast5_1.fast5
│ ├── fast5_2.fast5
│ ├── fast5_3.fast5
│ ├── fast5_4.fast5
│ ├── fast5_5.fast5
│ ├── fast5_6.fast5
│ ├── fastq_0.fastq
│ ├── fastq_1.fastq
│ ├── fastq_2.fastq
│ ├── fastq_3.fastq
│ ├── fastq_4.fastq
│ ├── fastq_5.fastq
│ ├── fastq_6.fastq
├── cell_4
│ ├── fast5_0.fast5
│ ├── fastq_0.fastq
│ └── sequencing_summary_0.txt
└── N0001.md5
├── cell1
│ └── data.md5
├── cell2
│ └── data.md5
├── cell3
│ └── data.md5
└── cell4
上述文件来自测序商。校验完文件，将所有fastq文件合并压缩：
cat ./BMK220118-AT448-01N0001/cell_*/*.fastq |gzip -c > map.N0001_clean.fq.gz
压缩后的fastq文件只有5G多，但配套的fast5文件则超过100G。
然后，使用 nanopolish index 将fastq文件与fast5文件进行关联：
./nanopolish/nanopolish index -v \
-d ./BMK220118-AT448-01N0001/cell_1 \
-d ./BMK220118-AT448-01N0001/cell_2 \
-d ./BMK220118-AT448-01N0001/cell_3 \
-d ./BMK220118-AT448-01N0001/cell_4 \
map.N0001_clean.fq.gz
nanopolish index 的-d参数可以给出多个。
于此同时，使用minimap2对序列进行比对，并使用sambamba按比对位置排序：
# 直接使用nanopolish目录中编译好的minimap2进行比对
./nanopolish/minimap2/minimap2 -t 48 -a -x map-ont ref.fa map.N0001_clean.fq.gz | \
sambamba view -t 4 -S -f bam /dev/stdin | \
sambamba sort -o output.sorted.bam -t 4 /dev/stdin
比对完成，产生的结果为output.sorted.bam，且对应的index也已经生成。
之后，使用nanopolish分析其中一小段的区间的甲基化：
# 这里再强调一遍，务必保证环境变量 HDF5_PLUGIN_PATH 对应的 hdf5 plugin设置正确！
# export HDF5_PLUGIN_PATH=/usr/local/hdf5/lib/plugin
./nanopolish/nanopolish call-methylation -t 48 \
-r map.N0001_clean.fq.gz \
-b output.sorted.bam \
-g ref.fa \
-w "1:500000-1000000" > methylation_calls.tsv
上面的-w参数，并未出现在软件自带的帮助信息中，只在官方说明的mannual里面才出现。
methylation_calls.tsv为每一条Reads上的甲基化位点，最后使用nanopolish附属的脚本calculate_methylation_frequency.py得到单个位点的甲基化水平：
./nanopolish/scripts/calculate_methylation_frequency.py methylation_calls.tsv > methylation_frequency.tsv
分析流程的封装与改进
试跑完，我发现流程有下面一些问题需要改进：
- 1、整个基因组
nanopolish call-methylation直接跑下来，速度非常慢，而且后续calculate_methylation_frequency.py必须将Reads上单碱基甲基化结果methylation_calls.tsv全部吞进去，才能够产生单个位点的甲基化水平。后者单线程运行，速度非常慢； - 2、单碱基甲基化水平，首先按Reads在基因组的位置排序，而并非甲基化位点在基因组的位置排序的，导致接下来数据压缩后无法直接读取特定位置的甲基化中间结果；
- 3、目前，
nanopolish call-methylation提供了4种甲基化模式，分别解释如下：
cpg
5'-----CG-----3' 3'-----GC-----5'gpc
5'-----GC-----3' 3'-----CG-----5'dam
5'-----GATC-----3' 3'-----CTAG-----5'dcm
5'-----CC[A/T]GG-----3' 3'-----GG[T/A]CC-----5'
若针对医学，则一个cpg模式已经足够。但分析植物，现有的几种模式仍不够；
- 4、
calculate_methylation_frequency.py分析单个位点甲基化，并没有区分正负链，所以得到的位点甲基化水平并非单个位点的结果，而是相邻反向互补的GC碱基甲基化水平按Reads覆盖程度的加权平均值，以cpg模式为例，若正链甲基化水平代表C碱基，则负链则是接下来G碱基匹配的那个C碱基的甲基化水平，这两者有必要区分开来。
根据工具自身的情况，我尝试进行如下改造。首先，尝试压缩methylation_frequency.tsv并构建索引：
awk -F "\t" '{print $1"\t"$3"\t"$4"\t"$5"\t"$11"\t"$2"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' methylation_frequency.tsv|head -n 1>compressed.methcall.txt
awk -F "\t" '{print $1"\t"$3"\t"$4"\t"$5"\t"$11"\t"$2"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' methylation_frequency.tsv| \
grep -v start | \
bedtools sort -i - >> compressed.methcall.txt
bgzip compressed.methcall.txt
tabix -p bed -S 1 compressed.methcall.txt.gz
上述操作，将正负链信息Strands放在第6列，这样文件结构会跟bed格式非常相似。同时保留表头，则即使行列顺序改变，仍然不影响calculate_methylation_frequency.py后续分析。最后，当调整后的文件用bgzip压缩，就能够成功构建index，方便提取指定位置的位点甲基化信息。
然后，若拆分正负链分别得到单个位点甲基化水平，则用下面的命令：
# 正链，使用调整列顺序压缩后的数据，正负链信息在第六列
gunzip -c compressed.methcall.txt.gz| \
awk -F "\t" '{if($6!="-"){print}}' | \
./nanopolish/scripts/calculate_methylation_frequency.py -s /dev/stdin >Plus.txt
# 负链
gunzip -c compressed.methcall.txt.gz| \
awk -F "\t" '{if($6!="+"){print}}' | \
./nanopolish/scripts/calculate_methylation_frequency.py -s /dev/stdin >Minus.txt
经过上述验证尝试，我编写了一份封装脚本Nanopolish_Methyl.pl，对nanopolish分析甲基化的过程进行了封装。只需将脚本Nanopolish_Methyl.pl放在nanopolish源文件的目录下，按照操作说明安装相应perl包，再按说明依次操作即可。
纳米孔数据分析碱基修饰
纳米孔数据可以完成多种碱基修饰，不仅仅是甲基化。而目前，使用nanopolish，并不能够将所有5mC甲基化的情况都能够分析。但是，甲基化碱基修饰分析的潜力仍有待挖掘，分析工具层出不穷。接下来简要介绍几篇文献：
Liu, Y., Rosikiewicz, W., Pan, Z. et al. DNA methylation-calling tools for Oxford Nanopore sequencing: a survey and human epigenome-wide evaluation. Genome Biol 22, 295 (2021). https://doi.org/10.1186/s13059-021-02510-z
7种针对纳米孔人类CpG甲基化分析的工具测评
Liu, Q., Fang, L., Yu, G. et al. Detection of DNA base modifications by deep recurrent neural network on Oxford Nanopore sequencing data. Nat Commun 10, 2449 (2019). https://doi.org/10.1038/s41467-019-10168-2
机器学习方法分析甲基化流程DeepSignal，王凯、肖传乐教授通讯作者
Ni, P., Huang, N., Nie, F. et al. Genome-wide detection of cytosine methylations in plant from Nanopore data using deep learning. Nat Commun 12, 5976 (2021). https://doi.org/10.1038/s41467-021-26278-9
针对纳米孔植物甲基化分析流程DeepSignal-plant，使用机器学习方法，可以分析CpG，CHG，CHH三种甲基化模式（对植物学很重要），肖传乐教授为通讯作者之一
Leger, A., Amaral, P.P., Pandolfini, L. et al. RNA modifications detection by comparative Nanopore direct RNA sequencing. Nat Commun 12, 7198 (2021). https://doi.org/10.1038/s41467-021-27393-3
使用纳米孔分析RNA 的 m6A修饰
