转录组拼接踩过的那些坑-原理

前前后后做了一些转录组拼接的工作，而且目前还在持续进行中。期间遇到了各种各样莫名其妙的坑，也尝试了一些不同的方法和软件。做一个阶段性记录和小结。本篇是原理部分。

拼接基本原理

拼接可以分为基因组和转录组拼接，基因组拼接对数据量和测序深度要求更好，而转录组用平时的RNA-seq数据就可以。目前手里更多的是RNA-seq数据，所以做的也是和转录组拼接相关的内容。

无论拼基因组还是转录组，归根结底都是拼DNA序列。拼接最简单的逻辑就是一个由让由长变短再由短到长的过程，基因打断成一个个片段进行测序，生成测序数据，然后再用reads拼成contigs再拼成scaffolds。

从reads到contigs的过程中，需要进行多序列比对并将一致的 reads（consensus sequence）拼接起来来生成contigs，其中最大的问题是基因组上存在大量重复区域，会对拼接带来困扰。Scaffolds数据则是通过pair end reads 信息来判断contigs 的顺序、方向和相邻 contigs 之间的缺口（gap）大小来生成的。从contigs到scaffolds是一个排序和定向的过程。

关于拼接，Biostar大牛Istvan Albert在一个回答下面评论道：Assembly is somewhat of a "dark art"。

拼接常用算法

目前常用的拼接算法都是基于数学中的图论思想（Graph theory）产生,其中图论中的两个点表示两个read，而两点之间的连线表示两条read的重叠区域。拼接要做的事情就是在所有的路线中找最优解，类似于小时候玩过的一笔画问题。

如上图左所示，一个简易的基因组产生了ABCD……若干read，在理想的情况下我们可以根据所有reads彼此的重叠区域重构出圆圈所表示的基因组。如果简化成图来表示则应该是上图右的所示，所有黑连接的是正确的基因组，但实际情况是基因组有很多区域比较相似，以至于序列间会产生本没有的联系（如红线所示）。

在进行上述分析过程时，需要把所有的reads都进行比对，以便找到重叠区域，这个步骤非常耗费计算资源。Graph画出来之后的问题就是如何从中得到最优路径，即从有多种reads组合方式中找到从合适的一个形成contig。

下面是寻找最优路径的常用算法。

Greedy extension

贪婪算法（贪心图）是早期提出的拼接算法。首先选定初始read, 然后找和其重叠区域最高的read进行延伸，直到拼接后的read两端都不能再进行扩展为止。

每一次都是从最优匹配开始，然后次优匹配，到不能匹配时停止。这样一来，贪婪算法通常会得到局部最优解，而不是全局最优解。因此，这种算法在遇到重复序列时会出现比较大的问题。

A greedy assembler compares all pairs of read fragments continually and replaces any pair with sufficient overlap between the edges with a combined sequence. The assembly completes, when the edges of the remaining sequences do not have any significant overlap.

Overlap Layout Consensus

OLC图算法主要是用来针对长reads序列拼接，如一代测序数据（三代测序数据），简单理解就是把测序产生的长序列用彼此之间的overlap区域连接起来。对于数据量很大的数据或者全基因组数据来说，形成的olc图非常复杂，会消耗大量内存。

OLC算法共有三步：

• Overlap

  对所有reads计算任意两条之间的重叠区域，挑选出满足筛选条件的reads。这里区别与贪婪算法，会先把所有overlap都找到。这一步，通常会将一个reads分成若干个长度比较短的序列（kmer/seed/word），要求是每个片段序列之间至少有若干个碱基的重叠区域。

• layout garph

  简单化过程。对reads进行排序，确定reads之间的位置，建立overlap图，将重叠的reads组合成contig。

• Consensus

  在已经建好overlap图的基础上，将所有的read序列排列起来，找一条从起始节点到终止节点的最佳近似路径使得最终路径将会遍历一次重叠区域中的每个节点，相当于对初始的reads集合中全部序列进行重构得到目标基因序列。

de Bruijn graph

De Bruijn 图是目前最常用的二代测序拼接算法。比较流行的拼接软件如Velvet、Abyss和SOAP denovo都使用该算法。

与OLC不同之处在于，这个算法将已经非常短的reads再分割成更多个kmer短序列（k 小于reads 序列的长度），相邻的kmers序列通过（k-1）个碱基连接到一起（即每次只移动一个位置），进而降低算法计算重叠区域的复杂度，降低内存消耗。

这里的kmer首先不能太短，比如2个碱基肯定拼不出来基因组，它的长度既需要能够使其携带足够的基因组的信息，也要短到可以进行后续的错误矫正。除此之外，一个read中小的片段被分割之后还不会丢失原来reads 的前后位置信息。

总体而言，该算法将reads打断成长度为K的核酸片段，再用Kmer间的overlap关系构建DBG，最后通过DBG得到基因组序列。

拼接步骤通常包括：

• 构建DBG图，将read分割为一系列连续kmer

如下图，1和2两个序列会拼接成一条长的序列，在引入第三条序列后会出现一个圆环。

下图是一个最简单的DBG拓扑结构，两球一线（一进一出）代表相邻的两个kmer，圆圈则代表有多种连接方式。

下图是三种常见的DBG结构，分别代表了拼接过程中的不同情况。

• 合并DBG图

合并路径中出度入度唯一（one incoming and one outcoming ）的节点，去除段末端，低覆盖度节点和泡状结构。

• 构建contig

寻找最优路径（经过每个节点且仅经过一次），最优路径对应的碱基序列构成一个contig

• 构建scaffold

通过PE reads 位置信息确定contig之间的相对位置和方向，组装contig，填充contig之间的gap，得到scaffold序列。

几个注意事项

1. 当把双链信息考虑进来之后，可能的连接情况就会增加，因为任何一个节点都可能和某个节点的反向互补序列相连。

2. 有重复区域时的情况如下图所示，可以发现，一个确定的genome只可能有一个DBG，但反过来一个DBG不一定找到的是唯一的潜在基因组。

kmer和内存

在拼接相关的文章中，kmer是出现频率非常高的一个词。而kmer在整个生物信息分析过程中的用处也是非常之多。这里有一份资料进行了比较详细的介绍

k值越大可辨别更多的小重复序列，越容易把DBG转换为唯一的序列，但得到的拼接过程含有更多的gaps；小的k值对应的DBG能够得到较好的连通性，但是算法的复杂度会提高，repeats序列处理会更复杂，增加了错拼的可能性。

在拼接数据预处理软件khmer的文献中有这样一段关于kmer和内存大小与处理结果关系的描述：

The interaction between these three parameters and the filtering process is complex and depends on the data set being processed, but higher coverage levels and longer k-mer sizes result in less data being removed. Lower memory allocation increases the rate at which reads are removed due to erroneous estimates of their abundance, but this process is very robust in practice

kmer越大需要的内存就越多，所以计算机的内存大小也会限制kmer的取值。这里需要说明的是，输入数据的多少不会影响memory用量，但是输入数据的错误越少，占用的内存也就越少，假设所有测序数据都没有任何错误，那么DBG的大小并不会因为测序深度的增加，因为不需要将因为几个碱基不一致的kmer存入到DBG中（下一部分会具体提到）。至于需要多大的RAM则取决于DBG的大小和组装基因组的大小。

另外，在拼接的过程中尽量避免使用偶数kmer，否则容易是kmer产生回文序列，特别是在链特意性的数据中。

在平时分析中，一般会设置一个kmer的梯度（21，23，25，27,2931），来解决DBG算法loss of read coherence的问题。然后从中选择最好的结果。另外，还有一种说法是在进行拼接过程时，kmer应该选择read长度的1/2到2/3大小，否则可能拼接出过多的Contig。这一点，也可能是我们平常使用trinity拼接时拼出Contig 过多的原因，trinity的默认拼接大小是25。上限是32？（有待确定）。如果kmer有上线，是否也可以考虑在预处理的时候，处理的力度大一点，把序列截短一些？

拼接的干扰因素

在实际情况中，拼接往往是在覆盖度不均匀且含噪声的数据中进行，这为拼接带来了三个方面的困难：

1. 增加了大量假kmer从而提高了对缓存的要求
2. 错误的reads通过增加tips，bubbles和corss-links等改变了DBG的结构
3. 不统一的read覆盖度使得拼接参数对拼接结果有非常大的影响

假kmer

在一次测序得到的数据中，kmer matches 的数量和测序深度以及read长度相关。假设在完全没有测序错误的数据中，read长度是100，测序深度是50X，选取kmer值为21，那么一个只匹配基因组一次的kmer 出现的次数应该是$(100-21+1)*50/100=40$（基因组的每个位置被测了50次，100bp的read有80个21bpkmer）,如果匹配基因组两次的kmer应该出现80次。因此，峰值应该在40，80的位置有一个小峰。

但是，当存在测序错误的时候，会在匹配1次的位置出现大量的kmer,就是由于测序的误差导致的。为什么说是测序错误呢，因为在50x的测序中只出现了一次，如果一个read中有一个碱基错了，那么这个read就会产生21个错误的21kmer。更大的问题是，随着测序深度的增加，这样错的kmer数量也会增加，

改变结构

错误的reads 会为DBG引入三种类型的错误

tips

所谓tips指一个小分支，下图所示，我们有5条10bp的reads，其中第5条有一个碱基测序错误。使用7kmer会产生11个节点。如果用前4个read来拼是正常的，一旦引入第五条就会因为一个错误的碱基出现一条错误的有三个节点的分支。

在大量拼接的过程中，大量的read会匹配到正确的位置，一小部分会比配到错误的位置，因此可以对错误的tip进行清除。

bubbles

当read和kmer相比足够长时，错误可能出现read中间。此时会出现bubble的情况

需要注意的是，除了测序错误以外，可变剪切和snp以及插入缺失等也会导致tips的出现，因此增加了拼接的难度。因此，在进行拼接之间有必要对原始数据进行预处理。

覆盖度不均一

在实际拼接过程中，会去除一些低频率的kmer，这一操作在删除了大量错误kmer的同时，我们也不可避免的删除了很多是因为覆盖度低所以出现次数低的kmer。

综上，如果调低cutoff，高覆盖区域出错的可能就高，但是低覆盖区域的质量提升。kmer长度增加会使低覆盖区域更加分散，因为kmer的覆盖度会因为低于设置的cutoff被删除。

随着kmer的增加，分支会逐渐减少，DBG会越来越趋向于线性。对于一个很大的kmer,可能就完全线性，同时按照染色体分开。拼接质量通常会随着k值的增加先变好然后再变坏，因为这个过程中存在两种竞争性的过程。一方面，kmer的增加可以更好的处理重复区域，但是另一方面，由于覆盖度的原因，kmer的增加在一些区域使得他们出现的次数越来越低直到低于筛选的阈值。所以，kmer 的选择非常重要。

基因组拼接和转录组拼接

Problem	Structure	Coverage
Heterozygous genome	Combines two dissimilar genomic graphs	Equal coverage of the two
Metagenome	Combines hundreds of similar or dissiminar genomic graphs	Widely varying coverage
Transcriptome	Combines thousands of graphs from genes	Widely varying coverage

不同的拼接内容需要不同的拼接策略，其原因如上图所示，即不同的数据产生的DBG结构和覆盖度不同。

对于转录组拼接来说，如果假设一个转录组只有两个基因而且这两个基因没有重复区域，那只需要构建两个没有关联的DBG就可以了。

对于实际数据来说，一个转录组包含了成千上万个基因。且大多数基因没有overlap，我们构建的DBG实际上是众多不相关grphs的集合，每一个图都代表一个基因或者一个基因家族。下面的示意图显示构建了五个基因的图并且也展示了其覆盖度的差别，其中基因2和3可能是两个可变剪切或者一个家族非常类似的两个基因。

通过上图也可以发现，转录组拼接最大的问题在于kmer数目的高低很大程度上取决于某个基因的表达量高低。如果把kmer cutoff 设置的高一些，那么低表达的基因就可能拼不起来。

在基因组拼接过程中，kmer coverage 通常是单峰或者是双峰，但是转录组拼接不同。转录组拼接中，每一个基因的kmer分布可能是一个单峰，peak的位置取决于基因表达量的多少。而整体的分布则是所有这些单基因分布的总和。

转录组拼接时应该注意这样量个问题。首先，转录组拼接已经不满足DBG算法中覆盖度均一的假设，表达量非常高的基因可能主宰拼接的结果，因此需要能调整这种覆盖度的差异。另外，转录组拼接不用太担心低复杂度区域，不会有很多重复区域出现。

转录组拼接常用软件：trinity；SOAPdenovo-trans(华大)；Trans-ABySS；bridger

拼接质量的评估

基因组层面的拼接质量，一般会比较看重长度相关的指标。比如最大长度、平均长度、拼接后的总长度和 contig N50长度。

Contig N50 指 reads 拼接后获得一些不同长度的 contigs，将所有的 contigs 长度相加得到总长度。然后将所有的 contigs 按照长度从长到短进行排序，再将 contigs 按照这个顺序相加，当长度等于 contigs 总长度的一半时，最后一个加上的 contig 的长度称为 contig N50。对于总长度不同的两个拼接数据，直接对比N50 的数值没有什么意义。

对于转录组拼接而言，并不是越长越好，我们更在意的是拼接的质量，方向和回帖率等等信息。如果我们在转录组拼接过程中使用了kmer=25这个参数，在拼接好后应该用拼接的fastq文件mapping回拼接好的转录组，测试mapping效率，这里推荐使用salmon软件，需要注意salmon中的kmer应该和拼接时采用的kmer保持一致。

另外，transrate是一个专门用来评价拼接质量的软件。在下文的实际应用部分会有介绍。

---

加入靠谱熊基地，和大家一起交流