基础知识复习之比对

刘小泽写于18.11.30

一般我们得到原始数据、质控过滤后，要么进行比对，要么序列拼接。序列比对的话可以选择参考基因组、参考转录本，目的就是看看测序的reads分布在什么位置，然后根据这个去找变异或者看表达量多少；拼接的话可以拼接转录本或者构建基因组

序列比对

就是将测序reads重新定位到基因组/转录组上，又叫mapping

先说下特点

和平常我们常用的blast同源比对不同，这里的序列比对主要指高通量测序得到的短序列，这种序列主要有这几个特点：

• 均匀 覆盖全基因组
• 读长短：因此会出现一条reads比对到基因组许多位置，软件识别就会出现问题
• 有一定的错误率：reads中的错误会带到比对结果，产生噪音干扰，尤其是变异检测时，SNP和测序错误这两者需要区分
• 测序深度较高：目的就是解决上面读长短、测序错误率的问题，让一个位点多得到一些reads，帮助判断。一般测序量都是基因组的几十或者上百倍，目的就是一个：提高准确度！
• 双端数据pair-end（PE）：这是illumina的巧妙设计，不是测序读长短吗，那么它就一次测片段的两端（比如构建的500bp文库中有一个500bp的DNA片段，测序时两头各测100bp，中间空300）。这两条PE数据可是有相关关系的：它们是基因组同一区域片段上的两端，分别来自DNA的两条链，并且二者的物理距离（insert size）是500bp，且有方向性
• 相比于blast同源比对，测序数据比对的容错性会更低，因此体现在计算亲缘关系上，对于亲缘关系相对较远的序列比对来讲，blast计算的同源性为80%，测序比对可能就只有50%

比对的结果

测序是个麻烦活，得到的结果又各自千差万别，因此比对的结果也有好多种情况，先简单就一对PE reads的比对情况来了解：

• 最好的情况（Perfect match）：两条PE reads都没有错配地比对到了基因组唯一位置【1 vs 1 无错配】

• reads有错配地比对到了基因组唯一的位置，可能原因包括：测序错误；SNP和InDel 【1 vs 1 有错配】

• reads无错配地比对到基因组多个位置，可能原因：reads来自基因组上重复区域，由于序列长度短，软件无法准确判断具体来源的位置，只能都显示出来【1 vs 多 无错配】

• reads有错配比对到基因组多个位置，可能原因有很多：基因组重复区域的影响、测序错误或者突变

上面说的Pair end比对就是：两条reads同时比对到同一序列，当然，除了PE reads比对外，还有single end（SE）比对，包括了：

• 只有一条reads比对上
• 两条reads都比对上，但比对的是不同的序列
• 两条序列比对后的距离超过了insert size的长度

另外两条reads可能一条比对上的也没有，可能由于reads中错配太多或者两条reads同源性比较低

序列比对的应用

应用一：与自身拼接结果进行比对

比如自身的基因组、基因集

• 计算位点覆盖深度
• 计算参考序列覆盖比率

应用二：与参考进行比对

比如参考基因组、基因集、公共数据库等

• 变异检测
• 有参转录组

通过比对，我们可以得到一些具体的有用信息：

• reads利用率：比对上的reads/总reads。例如：总reads数是100w，其中PE比对上的有90w，那么PE比对率为90%，另外single比对上的有5w，则SE比对率为5%，因此总reads利用率为95%
  在应用一中，可以利用数据利用率评价序列拼接的可靠性；在应用二中，可以衡量样本与目标参考序列的同源性

• 覆盖深度（Coverage depth）/覆盖度/乘数：就是平常说的“测序测了10X或30X”，它表示每个碱基平均被测了多少次【测序量的首要衡量标准，如果公司测了100X基因组的数据，拿回来检测看到每个位点的覆盖度都在100左右，那么这个结果就是不错的；同时侧面反映了建库测序时随机打断的过程是不是均匀】
  覆盖比率（Coverage ratio）/覆盖率：被测序的碱基占全基因组大小的比率，它随覆盖度的升高而升高，同时受到测序偏差（bias）的影响【最直观的理解就是illumina测序会受到GC bias的影响】【全基因组测序理论上要覆盖所有的区域，即测序要饱和】

  它们虽然很像，但绝不一样！
  覆盖深度是可以想象是纵向的概念，而覆盖比率是横向的，例如：
  有一条1k长度的序列，通过序列比对，这1k个位点中有990个被测序到，那么覆盖率就是990/1000=99%；而覆盖深度则是将每个位点被覆盖的次数求和，然后除以基因组的长度
  

  覆盖深度、覆盖率

  有时分析比对结果发现，有的区域覆盖深度很高，是平均深度的几倍以上，称之为“高覆盖区”，这一般是基因组上的重复区域，因为来自不同区域的测序reads都可以mapping到这一块区域上；

  有高就有低，“低覆盖区” 可能属于GC存在偏差的区域，如高GC区域测序不均匀；或者基因组的复杂区域（杂合率较高或者简单重复区域）中拼接准确性比较低，导致mapping比率低

比对软件

常用算法

• 空位种子片段索引法，如Maq、ELAND，先将读段切分，并选其中一段或几段作为种子建立搜索索引，再查找索引并延展匹配来定位读段，通过轮换种子来定义允许的错配和各种可能的位置组合
• Burrows Wheeler转换【最常用】，如BWA、SOAP、Bowtie，利用BW转换将基因组序列按一定的规则进行压缩，并构建索引，再回头查找定位读段，通过碱基查找与替换定义允许的错配
• Smith-Waterman动态规划，如BFAST，利用迭代关系计算两个序列所有可能的比对分值，将结果存在一个矩阵中，再回头寻找最优比对结果

比对过程

• 目标序列fasta构建索引
  因为比对数据量很大，通过索引可以很快查找到比对的参考序列对应位置
• 短序列比对

最常用的BWA软件

顾名思义，就是采用bwt算法的aligner工具，输出sam/bam，目前比对最常用

• bwa index构建索引，其中注意-a是选择建立索引的方法（包括bwtsw、is、div三种，默认是is）其中bwtsw适用于比较大的参考基因组，如人，不能用于小于10M的基因组，如细菌； is不能用于大于2G的基因组
• bwa mem进行比对，如果下游需要用到gatk，就需要用-R指定类似这种"@RG\tID:$sample\tSM:$sample\tLB:WES\tPL:Illumina"的read group信息，用于区分不同的样本，其中ID每个group的唯一ID, SM表示样本名称， LB代表library,表示文库的名字，PL代表platform, 表示测序平台的名字，可选值有Illumina, Pacbio

再看看一些注意点

• 比对前都要构建索引，我们可以对基因组、基因集、数据库等构建索引（fasta格式），目标序列不要太短，不要有回车符（也就是不要直接将NCBI的一些碱基直接粘贴到记事本，再上传到linux服务器处理，因为从windows=》linux会自动加上回车符，如果发现可以用dos2unix命令去除）；另外选择正确的bwa构建方法，比如人要用bwtsw
• 比对的过程是资源消耗比较大的计算过程，对硬盘要求比较高。因此尽量用bam存储，或者利用管道直接跳过sam进行下一步分析
• 关于短序列与长序列比对：短序列一般考虑能不能比对上，而长序列考虑比对上多少；短序列一般设为5个gap，长序列相比能容许更多

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com