关于唯一比对和序列冗余

Unique mapping & Duplicates

Posted by Wenlong Shen on April 8, 2020

技术服务于思想

比对

高通量测序要么是为了拼接组装出一套完整的基因组,要么是为了得到感兴趣区域的序列信息,对于后一种情况,我们需要利用比对这一手段来确定序列在基因组的位置。二代测序常用的比对软件有BWA、Bowtie2、SOAP等。

基因组上存在大量序列重复的片段,它们在测序、比对时都会对数据产生干扰(你不知道测到的片段究竟在基因组的哪个位置),当然,目前有越来越多的实验和分析方法来解决这个问题,但如果不是特意要研究重复序列或相应的基因组区段,我们为了后续分析和生物学解释的方便,通常只选取唯一比对的序列。但相应地,何为唯一?序列比对算法本身是基于某种特定的打分规则来判断存在于基因组何处,因而,什么样的阈值能确定是唯一比对?

早些时候,Bowtie等比对软件会在sam文件中提供tag,如XS,表明还有另一个比对位置的分数也超过当前软件设定的阈值,我们可以利用该tag来区分是否唯一比对。但不同的比对软件未必使用同一种tag,设定的阈值也未必相同,使得这种方法在软件复用、批量处理中带来不便。目前,更多地是利用mapq值作为比对唯一性、可靠性的判断标准。而在WGS/WES的分析中,考虑到任何一条测序质量合格的reads都有可能带来有用的信息,很多算法并不对mapq进行硬过滤,如Strelka2会设置双重阈值,不同阈值的序列会被用于算法的不同步骤中。

冗余

在比对之后,我们往往还要进行一次去冗余的过程。这是因为,对于当前大多数二代测序来说,上机文库在构建过程中往往需要经过PCR扩增,由于PCR对不同的分子具有不同的扩增偏性,导致某些分子会“更多”地出现在测序文库中,这就是序列冗余的主要来源(要注意的是,虽然发生概率低,但测序仪也可能带来冗余,比如Illumina的HiSeq,相邻两个cluster可能混杂、光学解析可能错判等)。

这里我们引出一个问题:文库构建时通常会有多轮循环,比如10个循环就能将分子扩增1000倍,但为何实际测序的冗余并不会特别高?这是由于文库中unique的分子数要远高于测序仪flowcell上能够结合的cluster数,所以,冗余率其实和文库多样性有着密切联系,这也是目前Tn5在较少样本量下建库冗余率高的原因之一。

我们在后续分析中通常不会使用冗余序列,但考虑到冗余对考察样本整体的序列分布、质控依然有着一定意义,所以目前不同算法对待冗余序列的处理方法不尽相同,实际应用中需要有所注意。一般地,可利用Picard工具包的MarkDuplicates模块进行冗余去除。