【5.2.3】将events匹配参考基因组
一、分析流程
需要提前安装好：
- nanopolish
- samtools
- minimap2
1.下载测试数据
wget http://s3.climb.ac.uk/nanopolish_tutorial/ecoli_2kb_region.tar.gz
tar -xvf ecoli_2kb_region.tar.gz
cd ecoli_2kb_region
curl -o ref.fa https://ftp.ncbi.nih.gov/genomes/archive/old_genbank/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid225/U00096.ffn # 参考基因组
样品信息：
- Sample : E. coli str. K-12 substr. MG1655
- Instrument : MinION sequencing R9.4 chemistry
- Basecaller : Albacore v2.0.1
- Region: “tig00000001:200000-202000”
- Note: Ligation-mediated PCR amplification performed
2.用minimap2 来比对
nanopolish也接受fastq的文件呢！
建立索引
nanopolish index -d fast5_files/ reads.fasta
比对：
minimap2 -ax map-ont -t 8 ref.fa reads.fasta | samtools sort -o reads-ref.sorted.bam -T reads.tmp;
samtools index reads-ref.sorted.bam
检查bam的完整性
samtools quickcheck reads-ref.sorted.bam
1.3 将nanopore 的events比对参考基因组
nanopolish eventalign --reads reads.fasta --bam reads-ref.sorted.bam --genome ref.fa --scale-events > reads-ref.eventalign.txt
1.4 生成结果
结果文件： reads-ref.eventalign.txt
contig position reference_kmer read_index strand event_index event_level_mean event_stdv event_length model_kmer model_mean model_stdv standardized_level
gi|545778205|gb|U00096.3|:c514859-514401 3 ATGGAG 0 t 16538 89.82 3.746 0.00100 CTCCAT 92.53 2.49 -0.88
gi|545778205|gb|U00096.3|:c514859-514401 3 ATGGAG 0 t 16537 88.89 2.185 0.00100 CTCCAT 92.53 2.49 -1.18
gi|545778205|gb|U00096.3|:c514859-514401 3 ATGGAG 0 t 16536 94.96 2.441 0.00125 CTCCAT 92.53 2.49 0.79
gi|545778205|gb|U00096.3|:c514859-514401 3 ATGGAG 0 t 16535 81.63 2.760 0.00150 NNNNNN 0.00 0.00 inf
gi|545778205|gb|U00096.3|:c514859-514401 7 AGTTAA 0 t 16534 78.96 2.278 0.00075 TTAACT 75.55 3.52 0.79
gi|545778205|gb|U00096.3|:c514859-514401 8 GTTAAT 0 t 16533 98.81 4.001 0.00100 ATTAAC 95.87 3.30 0.72
gi|545778205|gb|U00096.3|:c514859-514401 9 TTAATG 0 t 16532 96.92 1.506 0.00150 CATTAA 95.43 3.32 0.36
gi|545778205|gb|U00096.3|:c514859-514401 10 TAATGG 0 t 16531 70.86 0.402 0.00100 CCATTA 68.99 3.70 0.41
gi|545778205|gb|U00096.3|:c514859-514401 11 AATGGT 0 t 16530 91.24 4.256 0.00175 ACCATT 85.84 2.74 1.60
文献： https://www.nature.com/articles/nmeth.4184 图1，我们显示了选择 k-mers 的直方图，event_level_mean以展示甲基化如何改变观察到的电流。这些数字的数据是由 eventalign 生成的，我们随后使用 ggplot2 在 R 中绘制。
二、流程说明
纳米孔测序仪将单链 DNA 穿过嵌入膜中的孔。孔允许电流从膜的一侧流向另一侧。当 DNA 穿过这个孔时，它会部分阻断电流的流动，电流由仪器测量。在 Oxford Nanopore 的 MinION 系统中，测量的电流取决于进行测量时驻留在孔中的 5-mer。
MinION每秒对电流进行数千次采样；当 5-mers 滑过孔时，应该在多个样品中观察到它们。MinION 的事件检测软件处理这些样本并尝试检测当前电平变化的点。这些跳跃表明一个新的 5-mer 存在于孔中。为了帮助说明这一点，我从下面的预印本中复制了一个图：
这是一个理想化的纳米孔测序过程的模拟。黑点表示采样电流，红线表示构成检测事件的连续段。例如，前 0.5 秒的平均电流约为 60 皮安，加上一点噪音。然后电流在 0.1 秒内降至 40 pA，然后跳至 52 pA，依此类推。
事件检测软件将事件写入 HDF5 文件。原始 kHz 样本通常不会存储，因为输出文件会大得不切实际。这是此模拟的事件表：
为了帮助将事件转化为 DNA 序列，Oxford Nanopore 提供了一个孔模型（pore model），该模型描述了每个 5 聚体的预期电流信号。孔隙模型是一组 1024 个正态分布 - 示例可能如下所示：
这表明测得的电流预计来自ñ( 53.5 ,1.32)当 AAAAA 在孔隙中时等等。
2.2 Inferring Bases from Events
使用孔隙模型（pore model）和观察到的数据，我们可以解决许多推理问题。例如，我们可以推断通过孔的核苷酸序列。这是碱基调用（base calling）问题。我们还可以在给定一组重叠读数的情况下推断基因组的序列（ overlapping reads）。
这些推理问题因两个重要因素而变得复杂。
- 首先，5-mers 的正态分布重叠。有 1024 种不同的 5 聚体，但信号范围通常约为 40-70 pA。这使得很难推断出哪个 5-mer 产生了特定事件。数据结构部分缓解了这种情况；解决方案必须尊重 5-mer 之间的重叠，因此当我们查看后续事件时，难以解决的位置可能会变得清晰。
- 其次，事件检测是实时进行的，不可避免地会出错。如果某些事件太短或 DNA 链相邻 5 聚体的信号非常相似，则可能无法检测到它们。后一种情况的极端情况发生在对长均聚物进行测序时——在这里我们预计电流不会发生可检测的变化。相反的问题也会出现。由于系统中看起来像是电流变化的噪声，事件检测器可能会将应该是单个事件的事件拆分为多个事件。处理这些人工制品（artefacts）是准确推断产生事件的 DNA 序列的关键。
2.3 Aligning Events to a Reference
我们为共识问题设计（e consensus problem）的隐马尔可夫模型将提议的共识序列的 5 聚体作为 HMM 的主干，并具有额外的状态和转换来处理跳过/分裂伪影（artefacts）。
eventalign模块nanopolish将此功能公开为命令行工具。该程序接收一组在基础空间中与参考序列（或基因组组装草图）对齐的纳米孔读数，并在事件空间中重新对齐读数。
管道使用bwa mem路线作为指导。我们从正常的 bwa 工作流程开始：
bwa mem -x ont2d -t 8 ecoli_k12.fasta reads.fa | samtools view -Sb - | samtools sort - alignments.sorted
samtools index alignments.sorted.bam
然后我们可以使用 nanopolish 在事件空间中重新对齐：
nanopolish eventalign -r reads.fa -b alignments.sorted.bam -g ecoli_k12.fasta "gi|556503834|ref|NC_000913.3|:10000-20000" > eventalign.tsv
输出如下所示：
contig position reference_kmer read_index strand event_index event_level_mean event_length model_kmer model_mean model_stdv
gi|556503834|ref|NC_000913.3| 10000 ATTGC 1 c 27470 50.57 0.022 ATTGC 50.58 1.02
gi|556503834|ref|NC_000913.3| 10001 TTGCG 1 c 27471 52.31 0.023 TTGCG 51.68 0.73
gi|556503834|ref|NC_000913.3| 10001 TTGCG 1 c 27472 53.05 0.056 TTGCG 51.68 0.73
gi|556503834|ref|NC_000913.3| 10001 TTGCG 1 c 27473 54.56 0.011 TTGCG 51.68 0.73
gi|556503834|ref|NC_000913.3| 10002 TGCGC 1 c 27474 65.56 0.012 TGCGC 66.96 2.91
gi|556503834|ref|NC_000913.3| 10002 TGCGC 1 c 27475 69.97 0.071 TGCGC 66.96 2.91
gi|556503834|ref|NC_000913.3| 10003 GCGCT 1 c 27476 67.11 0.017 GCGCT 68.08 2.20
gi|556503834|ref|NC_000913.3| 10004 CGCTG 1 c 27477 69.47 0.052 CGCTG 69.84 1.89
这是从 Nick 的大肠杆菌数据中读取的二维纳米孔的互补链 (c)，与大肠杆菌 K12 对齐。列出的第一个事件（事件 27470）的测量电流水平为 50.57 pA。它与参考基因组第 10,000 位的参考 5-mer ATTGC 对齐。孔隙模型表明，针对 5-mer ATTGC 测量的事件应该来自ñ( 50.58 ,1.022)，与观测数据非常吻合。接下来的 3 个事件（27471、27472、27473）都与相同的参考 5-mer (TTGCG) 对齐，这表明事件检测器错误地调用了 3 个事件，而其中应该只发出一个事件。请注意，这 3 个事件的电流似乎都来自预期分布ñ( 51.68 ,0.732).
此输出对于每个事件都有一行。如果跳过参考 5-mer，则输出中将出现未观察到信号的间隙：
gi|556503834|ref|NC_000913.3| 10009 GCACC 1 c 27489 67.52 0.028 GCACC 66.83 2.46
gi|556503834|ref|NC_000913.3| 10011 ACCGC 1 c 27490 65.17 0.012 ACCGC 65.03 1.92
在这里，我们没有在 10010 位置观察到 5-mer 的事件。
该模块有望使处理信号级纳米孔数据变得更加容易，并有助于开发改进的模型。该eventalign模块可以在最新版本的nanopolish中找到
三、我的案例
参考资料
