如何进行基因组序列比对？
阅读用时：全文共3小节，约6400字，约11分钟
关键词：参考序列、比对软件、SAM文件
拿到人基因组全外显子illumina下机数据fastq文件之后，如何进行后续的变异检测呢？首先要做的就是将测序得到的reads比对到人基因组参考序列上。
1
人基因组参考序列是如何获得的？
随着人类基因组计划（Human Genome Project，HGP）的进行，International Human Genome Sequencing Consortium在2001年首次公布了人基因组序列的草图，2003年人类基因组计划宣布完成。HGP的样本来源是大量的欧洲匿名捐献者（收集血液或精子），然后从中选取少数样本进行提取DNA，为使捐献者的身份得到保护，研究人员也不知道谁的DNA进行了后续的测序(PMID：11230172)。
2004年发布了更加精确的基因组序列，版本号为NCBI Human Bulid 35，并于2004年在Natrue发表文章(PMID：15496913)对该版本基因组的组装过程进行了详细描述。Bulid 35版本的基因组序列包括2.85 billion核酸，仅存在341个gaps（gaps主要是由于重复片段导致），它覆盖了99%的基因组序列。
随后，GRC（Genome Reference Consortium）在NCBI Human Bulid 35版本的基础上于2009年发布了GRCh37（hg19）版本，总长3,137,144,693；又于2013年发布了GRCh38版本，总长3,209,286,105。
目前很多关于人基因组的数据库，如dbSNP、Clivar等数据库均更新到了GRCh38版本。至于具体分析过程中使用哪个版本的参考基因组序列，可根据自己需求进行选择。
各版本人基因组信息详见：
https://www.ncbi.nlm.nih.gov/grc/human/data
各版本下载地址详见：
http://hgdownload.cse.ucsc.edu/downloads.html#human）
2
如何利用BWA进行序列比对呢？
人基因组参考序列的来源、详细信息及下载地址已经知道啦，那么我们来看看利用什么软件或算法来将测序数据比对到31亿碱基序列上呢？
目前应用最为广泛的软件非BWA莫属了。BWA软件除可用于最为常见的illumina平台测序数据，还可用于SOLiD, 454, Sanger reads, PacBio reads and assemblycontigs等。
目前Heng Li发布的BWA（Burrows-Wheeler Aligner）软件已包含三个算法，依次为2009年的BWA-backtrack算法（BWA-ALN；PMID:19451168）, 2010年的BWA-SW算法（PMID：20080505）和2013年的BWA-MEM算法（arXiv:1303.3997v2[q-bio.GN]）：
1）BWA-ALN算法可用于reads长度小于等于100 bp的测序数据；
2）BWA-SW算法可用于reads长度在70 bp~1 Mb的测序数据；
3）BWA-MEM算法可用于reads长度在70 bp~1 Mb的测序数据；而相对BWA-SW，该算法更快更准确；相对于BWA-ALN算法，在70 bp~100 bp的reads比对中，该算法也有更好的性能。
使用三种算法之前，均需要先对参考序列构建FM-index (Full-text index Minute space)。 FM-index是基于Burrows-Wheeler transform进行全文压缩和构建索引的算法。
bwa index hg19.fa
构建FM-index后，就可选择三种算法之一进行比对啦：
Single-end alignment --- aln算法：
bwa aln hg19.fa reads.fastq> aln_sa.sai
bwa samse hg19.fa aln_sa.sai reads.fastq > aln.sam
Single-end alignment --- bwasw算法：
bwa bwasw hg19.fa long_read.fq > aln.sam
Single-end alignment --- mem算法：
bwa mem hg19.fa reads.fastq> aln.sam
Paired-end alignment --- aln算法：
bwa aln hg19.fa reads1.fq> aln1.sai
bwa aln hg19.fa reads2.fq> aln2.sai
bwa sampe hg19.fa aln1.sai aln2.sai reads1.fq reads2.fq > aln.sam
Paired-end alignment --- bwasw算法：
bwa bwasw hg19.fa long_read1.fq long_read2.fa> aln.sam
Paired-end alignment --- mem算法：
bwa mem hg19.fa read1.fq read2.fq > aln.sam
比对之后得到的信息，就存储在SAM文件里。别看SAM文件只有11列必填字段，其包含的信息量可是非常丰富的，接下来我们就来看看该文件都存储了哪些信息吧。
3
SAM文件格式简介
SAM文件格式是Heng Li于2009年提出的用于存储比对结果信息的文本文件（以tab键分割），并同时发布了一款用于处理SAM文件的软件SAMtools（PMID：19505943）。
SAM文件包括两部分内容：注释信息（header section）和比对结果部分（alignment section）。
注释信息
@SQ | 参考序列的说明 |
示例：@SQ SN:chr1 LN:249250621 | |
SN: 人参考基因组中染色体的编号，如chr1,chr2 LN: 参考序列中某序列的长度 |
@PG | 使用的程序说明 |
示例：@PG ID:bwa VN:0.7.13-r1126 CL: bwa mem hg19.fa read1.fq read2.fq | |
ID: 软件名称 VN: 软件版本号 CL: 命令行 |
@HD | 用SO来记录比对后的reads其排列顺序 |
示例： SO:unsorted SO:queryname SO:coordinate | |
unsorted或者sam文件中没有@HD这行信息：比对之后默认的排序顺序，与比对时输入的fastq文件中的Sequence identifier一致。 那如何对比对后的文件进行排序呢？用samtools软件中sort模块，该模块需要输入bam格式： samtools view -Sb aln.sam > aln.bam
下面以图示的形式给大家展示下不同排序方式的sam文件： 先来看看fastq文件中Sequence identifier的排序序列 fastq文件：less -SN read1.fq 比对之后默认的sam文件，与比对时输入的fastq文件中的Sequence identifier一致，如下图： samtoolsview -X aln.bam | less -SN 备注：大家看到下面的文件第二列是字符串，与下面的几个截图不一致，是因为在samtools view时添加了-X参数，可以将第二列的FLAG转为字符形式，FLAG具体会在比对结果部分解说，在这主要是看FLAG这列最后的数字，“1”表示Paired end测序中的read1,“2”表示Paired end测序中的read2。
以query name排序的sam文件，以Sequence identifier（即比对结果部分的第一列）从小到大排序：samtools sort -n aln.bam aln.sorted 以coordinate排序的sam文件，以比对到参考序列上的位置（即比对结果部分的第三列和第四列）从小到大排序：samtools sort aln.bam aln.sorted |
总结来说，注释部分主要是记录比对reads时用的参考序列的基本信息、对SAM文件使用过的程序、及reads的排序规则。
接下来我们来详细讲解一下SAM文件中的比对结果部分。
比对结果
比对结果部分中的每一行表示一条read的比对信息 (如果是双端测序，则每一行记录的是单端read比对情况，在SAM/BAM文件中read称为Query template)，包括11个必填字段和其他可选字段，字段之间用tab键分割。11个必填字段，其顺序固定，详见如下：
Col | Field | Type | Brief Description | 示例 |
1 | QNAME | String | fastq文件中read的sequence identifier | E00513:140:H2GTJCCXY:4:1223:25723:72631 |
2 | FLAG | Int | Combination of bitwise FLAGs（详细介绍见下面） | pPr1d |
3 | RNAME | String | 该条read比对到参考序列中哪条染色体上 | chr1 |
4 | POS | Int | 该条read比对到染色体上的位置，并且是1- based coordinate system，即参考序列的第一个碱基位置编号为1（0- based coordinate system，参考序列的第一个碱基位置编号为0） | 13198 |
5 | MAPQ | Int | MAPping Quality is Phred-scaled：MapQ = -10 log10(P), 如MapQ=30,表示比对到该位置的概率为千分之一，相对MapQ=20来说，更加不是一个随机事件，该reads比对更加准确。 | 48 |
6 | CIGAR | String | 简要比对信息表达式（详细介绍见下面） | 150M |
7 | RNEXT | String | 示例中该条read是Paired end测序中的read1,该信息是与其配对的read2比对上的染色体编号（若该条read是Paired end测序中的read2,则该信息是与其配对的read1比对上的染色体编号）：“=”表示与该条read一致；若不一致，则填写相应的染色体号。 | = |
8 | PNEXT | Int | 在Paired end测序中与该read配对的另一条read比对到染色体上的位置。 | 13019 |
9 | TLEN | Int | 根据RNAME+POS与RNEXT+PNEXT计算该对reads对应的DNA文库片段的长度，如该示例中展示的：POS – RNEXT + CIGAR中的M+I个数 = 13198-13019+150 = 329。那为何是负329呢，因为该read1比对到read2的下游。 | -329 |
10 | SEQ | String | 该条read的碱基序列 | CCTTAGCCTCTGT… |
11 | QUAL | String | 该条read碱基序列对应的质量值 | JJFJFA<JJJJJJ… |
第二列是FLAG，Combination of bitwise FLAGs为标识的加和，各种比对情况分别以不同的数字表示，每一个数字代表一种比对情况。将该query template的符合某些比对情况的数字分别相加，得到的数字即为FLAG。Bitwise FLAGs详见下面的表格：
FLAG(十六进制) | FLAG(十进制) | String | Description |
0x0001 | 1 | p | Paired end测序得到的单端read |
0x0002 | 2 | P | 0x01存在时，该Flag才有意义：Paired end测序中成对的read1和read2均比对到参考基因组上的合适位置上。见下面详解1：上面的图---成对的read1和read2比对到了同一条染色体上；而下面的图---read1和read2比对到了不同染色体上。 |
详解1： | |||
0x0004 | 4 | u | 该条read自身没有比对到参考基因组上,见下面详解2的图示。 |
0x0008 | 8 | U | 0x01存在时，该Flag才有意义：与其配对的read没有比对到参考基因组上，见下面详解2的图示。 |
详解2： | |||
0x0010 | 16 | r | 该条read自身比对到参考基因组的负链上, 同时SAM文件中第10列的SEQ是fastq文件中碱基序列的反向互补序列，见详解3。 |
0x0020 | 32 | R | 0x01存在时，该Flag才有意义：与其配对的read比对到参考基因组的负链上。 |
详解3： 该E00513:140:H2GTJCCXY:4:1101:11475:2381的read2比对到参考基因组的负链上啦，原fastq文件中序列为： CAAATCTGCATGCAGTCTGCCAAACTGA。。。。。。 而在SAM文件中第十列SEQ为： 。。。。。。TCAGTTTGGCAGACTGCATGCAGATTTG | |||
0x0040 | 64 | 1 | 0x01存在时，该Flag才有意义：Paired-end测序中的read1 |
0x0080 | 128 | 2 | 0x01存在时，该Flag才有意义: Paired-end测序中的read2 |
0x0100 | 256 | s | 这条比对信息不是该read的最优比对位置，见详解4的图示。 |
详解4： | |||
0x0200 | 512 | f | the read fails platform/vendor quality checks |
0x0400 | 1024 | d | 该read是PCR或光学 (optical) 导致的duplication read，即至少存在另外一条read碱基序列与该条read一致。PCR过程就是DNA模板的复制，如果测序数据量足够多时，就有可能测到两条碱基序列一致的reads；光学重复简单的来理解就是在测序拍照获取信号过程中本身是来自于一个cluster的信号，却识别成两个cluster，导致出现了两条碱基序列一致的reads。 |
如此多的FLAG，如何快速提取出某指定FLAG的reads呢？只需要在samtools view后添加-f参数即可，如提取出没有比对到参考基因组上的reads：
-f参数指定该FLAG对应的十六进制值： samtools view -X -f 0x0004 aln.bam | less -SN
-f参数指定该FLAG对应的十进制值：samtools view -X -f 4 aln.bam | less -SN
第六列CIGAR，是简要比对信息表达式 (Compact Idiosyncratic Gapped Alignment Report)，其以参考序列为基础，使用数字加字母表示比对结果，其中M-Match/Mismatch; I-Insertion; D-deletion; S-Soft clip;H-Hard clip; P-Padding; N-Skipped region from the reference，示例如下：
总结来说，对于Paired end测序，SAM文件的比对结果部分中每条信息不单详细的介绍了本条read的比对情况 (如FLAG、RNAME、POS、CIGAR)，还记录了与其匹配read的比对情况 (如FLAG、RNEXT、PNEXT)。
大家是否对于比对过程中涉及到的人基因组参考序列的来源、比对软件和比对结果文件SAM/BAM有了更深入的了解？
下期我们来一起聊聊如何进行变异检测。
【完】
作者原创作品，未经授权禁止转载！
扫码关注，获取更多精彩内容
我
是
彩
蛋
关注公众号后：
回复文字：好好学习，收听喜马拉雅FM电台栏目《一分钟听懂NGS基础概念》，让生信分析不再遥不可及。
回复文字：果然科学，给你看一篇好玩的科普文章。