阅读用时：全文共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.saibwa samse hg19.fa aln_sa.sai reads.fastq > aln.samSingle-end alignment --- bwasw算法：bwa bwasw hg19.fa long_read.fq > aln.sam Single-end alignment --- mem算法：bwa mem hg19.fa reads.fastq> aln.samPaired-end alignment --- aln算法：bwa aln hg19.fa reads1.fq> aln1.saibwa aln hg19.fa reads2.fq> aln2.saibwa sampe hg19.fa aln1.sai aln2.sai reads1.fq reads2.fq > aln.samPaired-end alignment --- bwasw算法：bwa bwasw hg19.fa long_read1.fq long_read2.fa> aln.samPaired-end alignment --- mem算法：bwa mem hg19.fa read1.fq read2.fq > aln.sam比对之后得到的信息，就存储在SAM文件里
别看SAM文件只有11列必填字段，其包含的信息量可是非常丰富的，接下来我们就来看看该文件都存储了哪些信息吧
3SAM文件格式简介 SAM文件格式是Heng Li于2009年提出的用于存储比对结果信息的文本文件（以tab键分割），并同时发布了一款用于处理SAM文件的软件SAMtools（PMID：19505943）
SAM文件包括两部分内容：注释信息（header section）和比对结果部分（alignment section）
注释信息在此主要介绍三个信息：@SQ、@PG、@HD
@SQ参考序列的说明示例：@SQ SN:chr1 LN:249250621SN: 人参考基因组中染色体的编号，如chr1,chr2LN: 参考序列中某序列的长度@PG使用的程序说明示例：@PG ID:bwa VN:0.7.13-r1126 CL: bwa mem hg19.fa read1.fq read2.fqID: 软件名称VN: 软件版本号CL: 命令行@HD用SO来记录比对后的reads其排列顺序示例：SO:unsortedSO:querynameSO:coordinateunsorted或者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个必填字段，其顺序固定，详见如下：ColFieldTypeBrief Description示例1QNAMEStringfastq文件中read的sequence identifierE00513:140:H2GTJCCXY:4:1223:25723:726312FLAGIntCombination of bitwise FLAGs（详细介绍见下面）pPr1d3RNAMEString该条read比对到参考序列中哪条染色体上chr14POSInt该条read比对到染色体上的位置，并且是1- based coordinate system，即参考序列的第一个碱基位置编号为1（0- based coordinate system，参考序列的第一个碱基位置编号为0）131985MAPQIntMAPping Quality is Phred-scaled：MapQ = -10 log10(P), 如MapQ=30,表示比对到该位置的概率为千分之一，相对MapQ=20来说，更加不是一个随机事件，该reads比对更加准确
486CIGARString简要比对信息表达式（详细介绍见下面）150M7RNEXTString示例中该条read是Paired end测序中的read1,该信息是与其配对的read2比对上的染色体编号（若该条read是Paired end测序中的read2,则该信息是与其配对的read1比对上的染色体编号）：“=”表示与该条read一致；若不一致，则填写相应的染色体号
=8PNEXTInt在Paired end测序中与该read配对的另一条read比对到染色体上的位置
130199TLENInt根据RNAME+POS与RNEXT+PNEXT计算该对reads对应的DNA文库片段的长度，如该示例中展示的：POS – RNEXT + CIGAR中的M+I个数 = 13198-13019+150 = 329
那为何是负329呢，因为该read1比对到read2的下游
-32910SEQString该条read的碱基序列CCTTAGCCTCTGT…11QUALString该条read碱基序列对应的质量值JJFJFA<JJJJJJ…第二列是FLAG，Combination of bitwise FLAGs为标识的加和，各种比对情况分别以不同的数字表示，每一个数字代表一种比对情况
将该query template的符合某些比对情况的数字分别相加，得到的数字即为FLAG
Bitwise FLAGs详见下面的表格：FLAG(十六进制)FLAG(十进制)StringDescription0x00011pPaired end测序得到的单端read0x00022P0x01存在时，该Flag才有意义：Paired end测序中成对的read1和read2均比对到参考基因组上的合适位置上
见下面详解1：上面的图---成对的read1和read2比对到了同一条染色体上；而下面的图---read1和read2比对到了不同染色体上
详解1：0x00044u该条read自身没有比对到参考基因组上,见下面详解2的图示
0x00088U0x01存在时，该Flag才有意义：与其配对的read没有比对到参考基因组上，见下面详解2的图示
详解2：0x001016r该条read自身比对到参考基因组的负链上, 同时SAM文件中第10列的SEQ是fastq文件中碱基序列的反向互补序列，见详解3
0x002032R0x01存在时，该Flag才有意义：与其配对的read比对到参考基因组的负链上
详解3：该E00513:140:H2GTJCCXY:4:1101:11475:2381的read2比对到参考基因组的负链上啦，原fastq文件中序列为：CAAATCTGCATGCAGTCTGCCAAACTGA





而在SAM文件中第十列SEQ为：





TCAGTTTGGCAGACTGCATGCAGATTTG0x00406410x01存在时，该Flag才有意义：Paired-end测序中的read10x008012820x01存在时，该Flag才有意义: Paired-end测序中的read20x0100256s这条比对信息不是该read的最优比对位置，见详解4的图示
详解4：0x0200512fthe read fails platform/vendor quality checks0x04001024d该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有了更深入的了解？下期我们来一起聊聊如何进行变异检测