# genomic-assembler

**Repository Path**: gkzhb/genomic-assembler

## Basic Information

- **Project Name**: genomic-assembler
- **Description**: No description available
- **Primary Language**: Unknown
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2020-11-21
- **Last Updated**: 2020-12-18

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# 宏基因组组装 - 算法 PJ

## 算法概述

利用错误率更低的短序列构建 De Brujin 图来推断原基因组。

具体来讲，先将每一条 read 都拆分为很多长度为 k 的片段 k-mers，然后对所有这些片段构造一个有向图，某个片段 A 有指向另一个片段 B 的边，当且仅当 A 的后 k-1 个碱基对序列与 B 的前 k-1 个碱基对序列完全相同。

在图中，我们取出尽可能长的路径将这些片段拼接起来就能得到推测的 DNA 序列。

### 算法

通过 `get_depth(id)` 对 id 节点递归寻找它的最大路径长度。
```
get_depth(id) {
  if (未访问过 id 节点) {
    标记访问过 id
    对所有 id 的子节点（出节点）递归调用 get_depth() 函数
    寻找所有子节点中最长路径长度 max_len
    id 的最长路径长度为 max_len + 1
  }
  return id 的最长路径长度
}
```

对所有 k-mers 节点调用以上 `get_depth()` 函数，找到其中最长的路径，得到 DNA 序列。

### 细节

构造完成 De Brujin 图后，我们通过递归搜索以某个节点为起点的最长路径。递归搜索过程中，优先递归搜索出现次数更多的节点。

在寻找到图中最长路径后，将它取出作为一条 DNA 序列。然后去除该序列中所有节点。

之后重复这一过程获得 DNA 序列。

## 代码实现

使用 Python 编写程序，分为两个文件，一个 `main.py` 处理数据的读取和输出，一个 `graph.py` 实现 De Brujin 图中的节点 `Vertex` 类和图 `Graph` 类。代码中原有的类 `Edge` 和 `Path` 都被简化掉了。

## 优化点

### 递归顺序

在读取 k-mers 节点时，对每个节点被访问的次数进行了计数（`weight` 属性）。 `weight` 属性的大小可以在某种程度上反映节点周边图结构的密度大小。所以对 `get_depth()` 函数中所有子节点的递归，按照 `weight` 从大到小的顺序递归，能够得到更好的结果。

另一方面，对于序列中出现的错误，这种递归顺序也能够较好地避免走到错误的序列上。

这个递归算法本质是一个随机算法。通过优化递归顺序，在理想情况下会得到更长更好的最长路径。

### 反向互补序列

对所有序列计算反向互补链并加入数据中，可以获得更完整、连通性更强的 De Brujin 图，更有利于构造最长路径。不过这会导致 Duplication ratio 提高，所以我最终没有添加反向互补序列（这样也可以节约内存和运行时间）。不过 `config` 中可以将 `rev_comp` 设为 `True` 启用这一项优化。

### 构造前导树

在构造 De Brujin 图的时候，可以构造一棵前导树，用来从 k-mers 序列起始开始匹配碱基对序列到对应 De Brujin 图中的节点。这样可以避免对所有 De Brujin 节点遍历匹配。

构造前导树的时候，我用 c++ 重新写了一个遍，代码见 `cpp` 文件夹中。每个 k-mers 节点都用 `Node` 类表示。前导树是一个四叉树，深度为 24 层（比 k-mers 长度小 1），所有节点都用 `TNode` 类表示。`TNode` 类的 `ch` 属性为指向子节点的指针数组，下标 0-3 分别对应 ATGC。树叶节点的子节点 `ch` 为指向 `Node` 类的指针。所以从根节点沿着前导树可以快速找到 k-mers 序列对应的节点。

同时，`Node` 类中有一个属性 `next` 为指向 `TNode` 叶子节点的指针，它表示以 k-mers 后 k-1 个序列为前缀指向的 `TNode` 叶节点，这样可以加快长序列切分成 k-mers 片段时寻找 k-mers 节点的速度。

但是从实际运行情况来看，构造这个前导树似乎是弊大于利的。k-mers 节点各不相同，字符串分布比较分散，所以在前导树中仍然有大量的节点，难以达到 O(log N) 的规模，而是更接近于 O(n)。虽然没有在 `TNode` 与 `Node` 中存储任何字符串（我本以为这样能够减少内存占用），但是由于节点众多，且每个节点存放多个指针，树状结构仍然无法有效地节约内存空间，甚至加剧了内存的使用。在实际运行时，即使用 data1 的数据构建 De Brujin 图也需要 6GB 以上的内存，导致程序直接被杀掉，甚至无法读完数据。所有前导树的优化方案最终被放弃了。

c++ 实现中，后面寻找 DNA 序列是通过两次广度优先搜索寻找最远路径实现，这个算法感觉可能存在一些问题，比如当路径分岔并最终合并的时候，广度优先搜索会选择更短的那一条分岔路径，而非更长的路径。

### Python 实现的优势

在 Python 中，字符串分片（slice）一般不会直接复制一份，而是类似于原字符串的指针。所以在 `Vertex` 类中存储大量的字符串子片段不会造成很大的内存占用。然后相比于构造前导树，利用 Hash 函数的 dict 字典(Map)数据结构更适合于此处的 string 到 id 的映射表示。因为 k-mers 字符串分布较为分散随机，不容易发生 Hash 碰撞，Hash 映射能够获得很好的性能。

## 运行

Python：

计算 `data1` 文件夹中的数据：
```bash
$ python main.py data1
```
结果在 `data1` 中的 `mycontig.fasta`

c++：

使用 cmake 编写项目，并且已附上 Makefile 文件，可通过 `make` 命令生成 `main.out` 可执行文件，通过 `./main.out data1` 同 Python 一样处理数据，不过没有验证能跑完，也没有输出结果（在读数据阶段就会因内存而被杀掉）。

## 参考

* [算法：SOAPdenovo的一个组装过程\_Krebs\_新浪博客](http://blog.sina.com.cn/s/blog_5d1edf6a0100w56l.html)
* [k-mer 如何影响宏基因组组装 ？](https://www.sohu.com/a/207040369_278730)
* [纯二代测序从头组装基因组(基础版） - 知乎](https://zhuanlan.zhihu.com/p/38317398)
* [DNA序列拼接中deBruijn图结构的研究 - 王东阳 - 百度文库](https://wenku.baidu.com/view/2f7e5c505727a5e9856a61a4.html)