算法 - Project: DNA 测序错误检测 - Task 2: 有噪音分段比对
Algorithms (H) @ Fudan University, spring 2021.
题目简介
解题报告
1 解题思路
在 Task 1 的基础上,本题的主要难点分为两部分:
- 如何将
long.fasta
中的 $\textrm{read}$ 片段快速匹配到参考字符串 $\textrm{ref}$ 上。 - 如何在 $\textrm{read}$ 片段平均 $15\%$ 的噪音干扰下,准确找到 SV 片段。
这里我们参考了 Minimap2 的论文思路,根据实际情况进行了简化和调整。算法核心分为以下三个步骤:
- 对参考字符串 $\textrm{ref}$ 利用最小哈希建立索引。
- 利用 $\textrm{ref}$ 的索引,将 $\textrm{read}$ 片段匹配到 $\textrm{ref}$ 上的对应区域。
- 比较 $\textrm{read}$ 片段和 $\textrm{ref}$ 上匹配到的区域,查找 SV 片段。
1.1 建立索引
由于噪音和 SV 片段的存在,我们无法简单地在 $\textrm{ref}$ 字符串中直接查找一个 $\textrm{read}$ 子字符串片段,因此我们需要建立索引。
如何对参考字符串 $\textrm{ref}$ 建立索引?简单来说,就是将 $\textrm{ref}$ 中每个长度为 $k$(即 HASH_SIZE
,默认为 $15$,所有参数均可在配置文件 src/utils/config.cpp
中调整)的子字符串,利用一个哈希函数转换成一个整数。具体来说,我们利用一个 2 位二进制整数来表示一个 DNA 碱基:$\mathtt{00}$ 表示 $\textrm{A}$、$\mathtt{01}$ 表示 $\textrm{T}$、$\mathtt{10}$ 表示 $\textrm{C}$、$\mathtt{11}$ 表示 $\textrm{G}$。对于一条 DNA 链,其哈希值就是将所有碱基的二进制数表示连接起来,例如 $\textrm{GCTA}$ 的哈希值就是 $\mathtt{11100100}$。通过这种方式,我们可以将每个连续的 $k$ 位子字符串转换成一个 $2k$ 位二进制数作为其哈希值 $\textrm{hash}$。我们称这样一个从 $\textrm{hash}$ 到这个子字符串在 $\textrm{ref}$ 中位置 $[i,i+k)$ 的映射为一个 $\textrm{k-mer}$。
对于一条长度为 $N$ 的 DNA 链,就包含了 $N-k+1$ 个这样的 $\textrm{k-mer}$。为了减少 $\textrm{k-mer}$ 的数量,以减少使用的内存空间,我们维护一个长度为 $w$(即 WINDOW_SIZE
,默认为 $10$)的滑动窗口,只有滑动窗口中哈希值最小的 $\textrm{k-mer}$ 才会被保存作为索引。
具体代码可参见 src/common/dna.cpp
中函数 Dna::CreateIndex
的实现。为了方便重复使用,我们提供了函数 Dna::PrintIndex
用于将索引导出成文件,以及函数 Dna::ImportIndex
用于从文件中读取索引。
1.2 字符串片段的模糊匹配
1.2.1 生成 minimizer
建立完索引后,我们就可以在每个 $\textrm{read}$ 片段中遍历所有长度为 $k$ 的子字符串,根据其哈希值查找是否有相同哈希值的 $\textrm{k-mer}$。同时,我们对于 $\textrm{read}$ 片段的反向互补序列 $\textrm{read’}$ 也进行同样的操作。对于每个找到的 $\textrm{k-mer}$,我们保存一个这样的结构:$\{\textrm{range}_\textrm{ref},\ \textrm{key}_\textrm{read},\ \textrm{range}_\textrm{read}\}$,我们称其为一个 $\textrm{minimizer}$。其中,$\textrm{range}_\textrm{ref}$ 表示 $\textrm{k-mer}$ 映射到 $\textrm{ref}$ 上的位置 $[i,i+k)$,$\textrm{key}_\textrm{read}$ 表示 $\textrm{read}$ 的编号(例如 $\textrm{S1}_1$),$\textrm{range}_\textrm{read}$ 表示这个子字符串在 $\textrm{read}$(或 $\textrm{read’}$)上的位置 $[j,j+k)$。同时,在每个 $\textrm{range}$ 中还额外保存了一个原字符串($\textrm{ref}$ 或 $\textrm{read}$)的指针,用于之后读取及合并这个子字符串的值。在 $\textrm{range}_\textrm{read}$ 中还额外保存了 mode
字段和 unknown
字段,分别用于指示当前 $\textrm{read}$ 的模式(是否是反向互补序列),以及是否包含一定数量的未知字符 $\textrm{N}$(在合并时用于提高效率,不关键)。
随后,我们根据 $\textrm{read}$ 和 $\textrm{read’}$ 片段上 $\textrm{minimizer}$ 的数量决定是否对 $\textrm{read}$ 进行反向互补操作。即如果 $\textrm{read’}$ 上的 $\textrm{minimizer}$ 较多,则进行反向互补操作,反之则不进行。
具体代码可参见 src/common/dna.cpp
中函数 Dna::FindOverlaps
的实现。同时,我们提供了函数 Dna::PrintOverlaps
用于将 $\textrm{minimizer}$ 导出成文件,以及函数 Dna::ImportOverlaps
用于从文件中读取 $\textrm{minimizer}$。
1.2.2 合并 minimizer
生成 $\textrm{minimizer}$ 后,我们需要对它们进行过滤及合并。其中,过滤指的是将错误匹配的 $\textrm{minimizer}$ 移除,合并指的是将两个 $\textrm{minimizer}$ 根据其 $\textrm{range}_\textrm{ref}$ 的范围 $[i_1,i_1+k)$ 和 $[i_2,i_2+k)$ 进行合并。
具体来说,在与一个聚类合并时,对于每一个 $\textrm{minimizer}$,我们比较此次合并后 $\textrm{range}_\textrm{ref}$ 和 $\textrm{range}_\textrm{read}$ 表示范围的增量 $\Delta_\textrm{ref}$ 和 $\Delta_\textrm{read}$。如果它们的差距不大,则将这个 $\textrm{minimizer}$ 归并到当前聚类,同时此聚类的计数器加 $1$。反之则尝试合并到下一个聚类,如果没有可合并的聚类,则将其单独分到一个新的聚类。
合并后,新的 $\textrm{minimizer}$ 的 $\textrm{range}_\textrm{ref}$ 为 $[\min\{i_1,i_2\},\ \max\{i_1,i_2\}+k)$,$\textrm{key}_\textrm{read}$ 为空字符串,$\textrm{range}_\textrm{read}$ 为 $[0,l)$,其中 $l$ 为合并后新生成的 $\textrm{read}$ 字符串的长度,同时 $\textrm{range}_\textrm{read}$ 中保存的指针指向这个新字符串。
于是,我们就得到了若干 $\textrm{minimizer}$ 聚类。我们将其中计数器值较小或者范围较小的 $\textrm{minimizer}$ 过滤。
具体代码可参见 src/common/dna_overlap.cpp
中函数 DnaOverlap::Merge
的实现。
1.2.3 匹配 ref 链和 sv 链
由于 long.fasta
中同时包含了多条 $\textrm{sv}$ 链的采样,我们需要从中找到与 $\textrm{ref}$ 链匹配的 $\textrm{sv}$ 链。这里我们根据 $\textrm{minimizer}$ 在 $\textrm{ref}$ 上的覆盖率,选择覆盖率最高的 $\textrm{sv}$ 链与 $\textrm{ref}$ 链相匹配。
在 src/utils/config.cpp
中调整 LOG_LEVEL
为 DEBUG
,即可在日志 logs/output.log
中看到 $\textrm{minimizer}$ 的覆盖率(搜索 cover rate
)。对于本题的正式数据,我们的覆盖率分别达到了:
NC_010513.1
: $99.01\%$ ($\textrm{S1}$)NC_014752.1
: $98.07\%$ ($\textrm{S2}$)NC_017999.1
: $99.05\%$ ($\textrm{S3}$)
具体代码可参见 src/common/dna_overlap.cpp
中函数 DnaOverlap::SelectChain
和 DnaOverlap::CheckCoverage
的实现。
1.3 查找 SV 片段
最终,我们将问题化归到了类似于 Task 1 的情形。根据每个 $\textrm{minimizer}$ 中保存的 $\textrm{range}_\textrm{ref}$ 和 $\textrm{range}_\textrm{read}$,我们可以得到两个需要比较的字符串,接下来只需复用函数 Dna::FindDeltasChunk
的逻辑即可。
但是,由于 Task 2 的数据含有一定量的噪声,原先对 Task 1 的数据处理方式不再适用于 Task 2,我们需要重新研究如何处理通过 Dna::FindDeltasChunk
函数得到的 SV。
具体来说,由于噪声的存在,SV 变得更加零散,同时我们难以区分一个 SV 是真正的 SV 还是噪声。我尝试过利用 SV 的间隔来判断一个 SV 是否是噪声,也尝试过魔改 Myers 差分算法来消除部分噪声,但效果都不理想。最后经过助教的提示,我调整了方案,使用一定范围内 SV 的密度来估计 SV 可能存在的范围。这是因为如果是纯噪声,SV 的密度大约会在 $15\%$ 左右,而对于真实的 SV,其密度往往在 $50\%$ 以上。通过观察 SV 密度的变化,就有可能判断 SV 的位置。参见 src/common/dna_delta.cpp
中函数 DnaDelta::GetDensity
的实现。
具体代码可参见 src/common/dna.cpp
中函数 Dna::FindDeltasFromSegments
的实现。
接下来就是调参的工作了,在配置文件 src/utils/config.cpp
中有大量可以调整的参数,其中比较重要的参数有 SIGNAL_RATE
, DENSITY_WINDOW_SIZE
, DELTA_MIN_LEN
, SNAKE_MIN_LEN
, GAP_MIN_DIFF
等。由于时间关系,没有很多时间用来调参了,因此最终的实验结果尚不理想。
2 运行代码
本项目使用 C++17 编写,环境要求:
- GCC 9.0 或以上
- GNU Make 4.0 或以上
可用的 Make 指令如下:
make
:构建项目1make help
:显示参数及其用法make run
:完整启动算法程序make index
:只创建 $\textrm{ref}$ 的索引make minimizer
:只生成 $\textrm{minimizer}$make start
:只查找 SV 片段make clean
:清理构建文件
使用 Task 2 测试数据生成的结果位于 tests/test_2/sv.bed
。
3 测试环境
- Ubuntu 18.04.5 LTS (WSL2 5.4.72-microsoft-standard)
- GCC 10.3.0
- GNU Make 4.1