测序数据比对

一、测序数据比对

高通量测序数据分析一共有测序数据分析主要有两条路径：一条是进行基因组拼接，得到基因组序列；另一条则是不经过拼接，直接与参考序列进行比对。由于拼接基因组需要消耗较多的计算资源，目前很多分析主要采用测序数据比对的方式。例如变异检测，RNAseq，甲基化检测，病原微生物鉴定等。

其实，即使是基因组拼接，其中也会大量用到测序数据比对的部分，例如二代测序拼接过程中，在得到 contig 之后，将 pairend reads 重新比对到 contig 上，寻找 contig 之间的关系，用户构建 scaffold。因此，测序数据比对是高通量测序分析中最核心的操作。

二、数据比对的意义

测序数据比对到参考序列上，得到一种“堆叠”的效果。这种效果是将测序数据比对到参考序列上。将测序数据信息，参考序列信息以及二者之间的比对信息都保存起来。

2.1 得到每个位点细节信息

首先，比对之后可以得到每个位点的细节信息，例如参考序列位点是否被覆盖到，有多少reads 覆盖到，每个位点的碱基序列是什么，这些信息将作为下一步序列分析的依据。

2.2 变异检测

测序数据比对一个重要的应用就是变异检测，例如当前的人基因组变异检测，检测 SNP，indel，SV，CNV 等都是通过该方法得到的。通过比对，可以比较参考序列位点与测序样品覆盖位点是否一致，如果不一致，则可能是一个潜在的突变位点，再结合覆盖深度，测序错误率，先验概率等信息综合评估。例如图中该位点参考序列为 A 碱基，测序数据中既包含 A 碱基，也包含 G 碱基，则该位点为一个潜在的 SNP。

利用比对结果进行变异检测

2.3 基因表达量计算

基因表达量属于表达分析，将 RNAseq 测序数据与参考序列进行比对，然后结合基因位置信息即可根据覆盖度情况判断每个基因的表达信息。不仅可以进行定性，例如判断基因是否有表达，还可以进行定量，比较不同基因表达量的差别。例如两个同样长度的基因 A 和 B，在相同测序深度情况下，A 基因测序深度为 100X，B 基因测序深度为 50X，则二者基因表达量差异为 2 倍。相比于 B 基因组，A 基因为高表达。

2.4 计算覆盖深度

根据计算每个位点比对的 reads 数目，即可计算覆盖深度。将全部比对数据除以基因组总长即可计算平均覆盖深度，例如基因组大小为 1M，全部比对上的碱基为 100M，则平均覆盖度为 100X。该值可以用来衡量测序数据覆盖情况。在组装过程中，覆盖深度用于识别重复序列；在变异检测过程中，该值用来识别缺失数据。在人类遗传疾病研究中，可以用于识别21 三体综合症；在植物数量性状研究中可以用于 QTL 定位。

2.5 计算覆盖比率

与计算覆盖深度类似，将参考序列被覆盖位点除以位点总数，即可计算覆盖比率，覆盖比率可以反映出测序数据与参考序列之间的亲缘关系。例如参考序列总长为 1M，其中 990K 被覆盖到，则覆盖比率为 99%，二者之间具有较高的亲缘关系。

2.6 计算 reads 利用率

与覆盖深度，覆盖比率类似，reads 利用率也是一个衡量测序数据比对的重要指标。只不过前两个指标是相对于参考序列的计算，而 reads 利用率相对于测序数据来计算。将比对上的reads 碱基数除以全部 reads 数，则得到 reads 利用率。reads 利用率可以反映出测序数据的比对效率，与参考序列之间的亲缘关系，测序数据的错误率，拼接结果准确性等。例如全部reads 有 100M 的碱基，比对上 90M 碱基，则 reads 利用率为 90%，表示有 90%的碱基可以比对到参考序列上。

2.7 可变剪切检验

将表达数据比对到参考序列上，如果同一条 reads 或者 pairend 的双末端 reads，同时比对到不同的外显子，则可以用于识别可变剪切。

通过序列比对识别可变剪切

2.8 基因融合检测

与识别可变剪切类似，如果测序数据比对到不同基因或者不同染色体上，则可以用于识别基因融合。

基因融合示意图

2.9 组装结果纠错

将同一样品的测序数据与该样品的拼接结果进行比对，然后根据比对情况进行纠错。如果是测序数据与参考序列进行比对则是找突变，与自身数据比对，则是进行纠错。

2.10 微生物鉴定

得到测序数据之后，可以不进行拼接，直接与物种分类数据库进行比对，用于鉴定微生物。该方法不仅可以用于定性，同样也可以进行定量研究，可以得到一个样品中物种的组成以及丰度。

2.11 基因组成环鉴定

如果测序数据可以同时比对上基因组的首位两端，则认为基因组成环。成环是细菌基因组是否为完成图的重要标志。可以使用 pairend reads，例如两条 reads 分别比对到基因组的首位两侧，或者使用长读长数据，直接跨过基因组首位两端来进行判断。

2.12 overlap 关系

三代测序 reads 与 reads 之间直接进行比对，可以得到 reads 之间的 overlap 关系，用于基因组拼接。

2.13 计算测序错误率

将已知基因组测序数据与基因组进行比对，可以用于计算测序错误率。例如 lambda DNA序列已知，将其测序数据与 lambda DNA 基因组进行比对，可以用于计算测序错误率。

2.14 生成一致性序列

与组装结果纠错类似，如果是测序数据与近源参考序列进行比对，则可以根据每个位点的比对情况生成一条一致性序列，该方法用于病毒基因组的拼接中。

三、短序列比对

最早的高通量测序数据读长都比较短，所以测序数据的比对，直接就称为短序列比对。随着三代长读长测序的兴起，目前有越来越多的长读长测序数据。三代测序数据长度更长，有更多的测序错误，包括一些小的插入缺失，之前的二代测序算法已经不再使用三代测序的比对，这就需要新的比对软件及算法。

3.1 短读长 reads 比对

随着以 illumina 为主的高通量测序技术的流行，短序列比对成为二代测序分析中的核心分析，短序列比对也就是将测序得到的短片段再回帖到基因组上，这个过程也被称为 mapping。mapping 的结果包含了所有的信息，后续所有的分析都是基于这个 mapping 的结果。

二代高通量测序具有以下特点：

1.测序覆盖全基因组

2.测序数据读长短

3.测序数据具有一定的错误率

4.测序数据深度高

5.测序数据具有 Pair-end 的关系

因此短序列比对具有其自身的特点。

1、需要测序 reads（单端或双端均可），fastq 格式与参考序列 fasta 格式；

2、参考序列可以是基因组也可以是基因集，只能是核酸序列；

3、需要对参考序列就建立索引；

4、比对是整条比上或者比不上，不能像 blast 比对，分开比对；

5、比对仅能容许一定数目的错配和空位；

6、序列太短，会出现一条序列比对到多个位置的情况；

7、数据量较大，比对比较耗时。

3.2 比对算法

短序列比对有很多比对软件，例如 bwa，soap，bowtie2，hisat2，subread 等，在众多的短序列比对软件中，BWA 几乎已经成为默认的行业标准。

随着测序技术的提高，BWA 软件也在不断地进化。使用的算法分别是 BWA-backtrack，BWA-SW 和 BWA-MEM 以及最新的 BWA-MEM2。最早的短序列比对主要使用的是 bwt 算法，(Burrows–Wheeler transform），bwa 软件的名字也来源于利用 bwt 算法的 aligner。

首先，BWA-backtrack，软件中的 aln，samse，sampe 比对都使用 bwa-backtack 算法中，主要适合比较短的 reads，尤其是 100bp 以下的，不过目前很少了；BWA-SW 中的 SW 表示affine-gap Smith-Waterman，BWA-MEM 中的 MEM 表示 maximal exact matches。BWA-SW 和 BWA-MEM 都支持长 read 适合 split 比对，reads 长度可以从 70bp 到 1M，BWA-MEM 是更新的算法，更快也更准确，是最推荐使用的算，也是使用范围最广泛的，Illumina, 454, Ion Torrent，Sanger，拼接的 contigs 或者 BAC 都可以使用这种算法，对于包含很多 gap 的情况，bwa-sw 会更好些。所以很明显，目前的版本我们使用 bwa-mem 算法就行了。

目前 bwa 软件又额外发布了 BWA-MEM2 比对软件，该软件比对速度更快，精确度更高。

bwa-mem2 官网：https://github.com/bwa-mem2/bwa-mem2

3.3 比对结果