如何进行基因组序列比对？如何进行基因组序列比对？天空中最亮的星百家号阅读用时：全文共3小节，约6400字，约11分钟关键词：参考序列、比对软件、SAM文件拿到人基因组全外显子illumina下机数据fastq文件之后，如何进行后续的变异检测呢？首先要做的就是将测序得到的reads比对到人基因组参考序列上。1人基因组参考序列是如何获得的？随着人类基因组计划（Human Genome Project，HGP）的进行，International Human Genome Sequencing Consortium在2001年首次公布了人基因组序列的草图，2003年人类基因组计划宣布完成。HGP的样本来源是大量的欧洲匿名捐献者（收集血液或精子），然后从中选取少数样本进行提取DNA，为使捐献者的身份得到保护，研究人员也不知道谁的DNA进行了后续的测序(PMID：)。2004年发布了更加精确的基因组序列，版本号为NCBI Human Bulid 35，并于2004年在Natrue发表文章(PMID：)对该版本基因组的组装过程进行了详细描述。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:）, 2010年的BWA-SW算法（PMID：）和2013年的BWA-MEM算法（arXiv:[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.samSingle-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：）。SAM文件包括两部分内容：注释信息（header section）和比对结果部分（alignment section）。注释信息在此主要介绍三个信息：@SQ、@PG、@HD。@SQ参考序列的说明示例：@SQ SN:chr1 LN:SN: 人参考基因组中染色体的编号，如chr1,chr2LN: 参考序列中某序列的长度@PG使用的程序说明示例：@PG ID:bwa VN:0.7.13-r1126 CL: bwa mem hg19.fa read1.fq read2.fqID: 软件名称VN: 软件版本号CL: 命令行总结来说，注释部分主要是记录比对reads时用的参考序列的基本信息、对SAM文件使用过的程序、及reads的排序规则。接下来我们来详细讲解一下SAM文件中的比对结果部分。比对结果比对结果部分中的每一行表示一条read的比对信息 (如果是双端测序，则每一行记录的是单端read比对情况，在SAM/BAM文件中read称为Query template)，包括11个必填字段和其他可选字段，字段之间用tab键分割。11个必填字段，其顺序固定，详见如下：@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.sortedColFieldTypeBrief Description示例QNAMEStringfastq文件中read的sequence identifierE:H2GTJCCXY:4::72631FLAGIntCombination of bitwise FLAGs（详细介绍见下面）pPr1dRNAME该条read比对到参考序列中哪条染色体上chr14POS该条read比对到染色体上的位置，并且是1- based coordinate system，即参考序列的第一个碱基位置编号为1（0- based coordinate system，参考序列的第一个碱基位置编号为0）131985MAPQMAPping Quality is Phred-scaled：MapQ = -10 log10(P), 如MapQ=30,表示比对到该位置的概率为千分之一，相对MapQ=20来说，更加不是一个随机事件，该reads比对更加准确。486CIGAR简要比对信息表达式（详细介绍见下面）150M7RNEXT示例中该条read是Paired end测序中的read1,该信息是与其配对的read2比对上的染色体编号（若该条read是Paired end测序中的read2,则该信息是与其配对的read1比对上的染色体编号）：“=”表示与该条read一致；若不一致，则填写相应的染色体号。=8PNEXT在Paired end测序中与该read配对的另一条read比对到染色体上的位置。130199TLEN根据RNAME+POS与RNEXT+PNEXT计算该对reads对应的DNA文库片段的长度，如该示例中展示的：POS – RNEXT + CIGAR中的M+I个数 = +150 = 329。那为何是负329呢，因为该read1比对到read2的下游。-32910SEQ该条read的碱基序列CCTTAGCCTCTGT…11QUAL该条read碱基序列对应的质量值JJFJFA&JJJJJJ…第二列是FLAG，Combination of bitwise FLAGs为标识的加和，各种比对情况分别以不同的数字表示，每一个数字代表一种比对情况。将该query template的符合某些比对情况的数字分别相加，得到的数字即为FLAG。Bitwise FLAGs详见下面的表格：FLAG(十六进制)FLAG(十进制)Description0x0001pPaired end测序得到的单端read0x0002P0x01存在时，该Flag才有意义：Paired end测序中成对的read1和read2均比对到参考基因组上的合适位置上。见下面详解1：上面的图---成对的read1和read2比对到了同一条染色体上；而下面的图---read1和read2比对到了不同染色体上。详解1：0x0004u该条read自身没有比对到参考基因组上,见下面详解2的图示。0x0008U0x01存在时，该Flag才有意义：与其配对的read没有比对到参考基因组上，见下面详解2的图示。详解2：0x001016r该条read自身比对到参考基因组的负链上, 同时SAM文件中第10列的SEQ是fastq文件中碱基序列的反向互补序列，见详解3。0x002032R0x01存在时，该Flag才有意义：与其配对的read比对到参考基因组的负链上。详解3：该E:H2GTJCCXY:4::2381的read2比对到参考基因组的负链上啦，原fastq文件中序列为：CAAATCTGCATGCAGTCTGCCAAACTGA。。。。。。而在SAM文件中第十列SEQ为：。。。。。。TCAGTTTGGCAGACTGCATGCAGATTTG0x0040640x01存在时，该Flag才有意义：Paired-end测序中的read10x00801280x01存在时，该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/M I-I D- S-SH-H P-P N-Skipped region from the reference，示例如下：总结来说，对于Paired end测序，SAM文件的比对结果部分中每条信息不单详细的介绍了本条read的比对情况 (如FLAG、RNAME、POS、CIGAR)，还记录了与其匹配read的比对情况 (如FLAG、RNEXT、PNEXT)。大家是否对于比对过程中涉及到的人基因组参考序列的来源、比对软件和比对结果文件SAM/BAM有了更深入的了解？下期我们来一起聊聊如何进行变异检测。本文由百家号作者上传并发布，百家号仅提供信息发布平台。文章仅代表作者个人观点，不代表百度立场。未经作者许可，不得转载。天空中最亮的星百家号最近更新：简介:表现宇宙现象，普及科学知识。作者最新文章相关文章BWA比对详解 - 简书
BWA比对详解
bwa主要用于将低差异度的短序列（一般是同物种）与参考基因组进行比对。主要包含三种比对算法：backtrack、SW和MEM，第一种只支持短序列比对（&100bp），后两种支持长序列比对(70bp~1M)，并支持分割比对（split alignment）。MEM算法是最新的也是官方推荐的。
注：如果短序列是嵌合体，可能会输出两条结果。
注：参考基因组单挑染色体长度要小于2G
主要语法：
bwa index -a bwtsw ref.fasta
#对于大基因组建立FM-Index
bwa index -a is ref.fasta
#对小基因组建立index，速度快，内存消耗大
bwa mem ref.fa reads.fq & aln-se.sam
bwa mem ref.fa read1.fq read2.fq & aln-pe.sam
bwa aln ref.fa short_read.fq & aln_sa.sai
bwa samse ref.fa aln_sa.sai short_read.fq & aln-se.sam
bwa sampe ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq & aln-pe.sam
bwa bwasw ref.fa long_read.fq & aln.sam
内置程序包括：index，mem，aln samse，sampe，bwasw，
Index：对参考基因组建立索引
-p： prefix of output database
-a：Algorithm for constructing index
:速度快，内存消耗大，适合小基因组&2G
：速度慢，省内存，适合大参考基因组
MEM:maximal exact matches
先做seeding alignment，再用SW算法做延伸
主要参数：
Number of threads [1]
线程数，默认为1
Minimum seed length. Matches shorter than INT will be missed. The alignment speed is usually insensitive to this value unless it significantly deviates 20. [19]
最小种子匹配长度，影响比对速度，默认为19
Band width. Essentially, gaps longer than INT will not be found. Note that the maximum gap length is also affected by the scoring matrix and the hit length, not solely determined by this option. [100]
gap最大长度，超过后将不会被匹配，默认为100
Off-diagonal X-dropoff (Z-dropoff). Stop extension when the difference between the best and the current extension score is above |i-j|*A+INT, where i and j are the current positions of the query and reference, respectively, and A is the matching score. Z-dropoff is similar to BLAST’s X-dropoff except that it doesn’t penalize gaps in one of the sequences in the alignment. Z-dropoff not only avoids unnecessary extension, but also reduces poor alignments inside a long good alignment. [100]
影响停止延伸的一个参数。
Trigger re-seeding for a MEM longer than minSeedLen*FLOAT. This is a key heuristic parameter for tuning the performance. Larger value yields fewer seeds, which leads to faster alignment speed but lower accuracy. [1.5]
一个关键的启发式搜索参数，较大的数值意味着较少的种子，比对速度也更快，默认1.5，即1.5倍的最小种子长度
Discard a MEM if it has more than INT occurence in the genome. This is an insensitive parameter. [10000]
In the paired-end mode, perform SW to rescue missing hits only but do not try to find hits that fit a proper pair.
Matching score. [1]
Mismatch penalty. The sequence error rate is approximately: {.75 * exp[-log(4) * B/A]}. [4]
Gap open penalty. [6]
Gap extension penalty. A gap of length k costs O + k*E (i.e. -O is for opening a zero-length gap). [1]
Clipping penalty. When performing SW extension, BWA-MEM keeps track of the best score reaching the end of query. If this score is larger than the best SW score minus the clipping penalty, clipping will not be applied. Note that in this case, the SAM AS tag reports the best SW clipping penalty is not deducted. [5]
Penalty for an unpaired read pair. BWA-MEM scores an unpaired read pair as scoreRead1+scoreRead2-INT and scores a paired as scoreRead1+scoreRead2-insertPenalty. It compares these two scores to determine whether we should force pairing. [9]
Assume the first input query file is interleaved paired-end FASTA/Q. See the command description for details.
Complete read group header line. ’\t’ can be used in STR and will be converted to a TAB in the output SAM. The read group ID will be attached to every read in the output. An example is ’@RG\tID:foo\tSM:bar’. [null]
Don’t output alignment with score lower than INT. This option only affects output. [30]
Output all found alignments for single-end or unpaired paired-end reads. These alignments will be flagged as secondary alignments.
Append append FASTA/Q comment to SAM output. This option can be used to transfer read meta information (e.g. barcode) to the SAM output. Note that the FASTA/Q comment (the string after a space in the header line) must conform the SAM spec (e.g. BC:Z:CGTAC). Malformated comments lead to incorrect SAM output.
Use hard clipping ’H’ in the SAM output. This option may dramatically reduce the redundancy of output when mapping long contig or BAC sequences.
Mark shorter split hits as secondary (for Picard compatibility).
Control the verbose level of the output. This option has not been fully supported throughout BWA. Ideally, a value 0 for disabling all
1 for ou 2 for 3 for 4 or higher for debugging. When this option takes value 4, the output is not SAM. [3]
参考来源：http://starsyi.github.io//BWA-%E5%91%BD%E4%BB%A4%E8%AF%A6%E8%A7%A3/ 一、BWA 简介 BWA，即Burrows-Wheeler-Alignment Tool。BWA 的学习主要...
这篇文章很长，超过1万字，是本系列中最重要的一篇，因为我并非只是在简单地告诉大家几条硬邦邦的操作命令。对于新手而言不建议碎片时间阅读，对于有一定经验的老手来说，相信依然可以有所收获。在开始之前，我想先说一句：流程的具体形式其实是次要的，WGS本质上只是一个技术手段，重要的是...
首先是方法的选择。基于距离的方法有UPGMA、ME（Minimum Evolution，最小进化法）和NJ（Neighbor-Joining，邻接法）等。其他的几种方法包括MP（Maximum parsimony，最大简约法）、ML（Maximum likelihood，最...
欢迎来GitHub上fork，一起进步： https://github.com/xuzhougeng 比对软件很多，首先大家去收集一下，因为我们是带大家入门，请统一用hisat2，并且搞懂它的用法。直接去hisat2的主页下载index文件即可，然后把fastq格式的rea...
最抓不住的是时间，也许昨天还是艳阳高照的夏，今天就成了白雪纷飞的冬。 十二月，2016的最后一个月，就这样悄无声息的到来，我还没有缓过夏末的炎热，就到了初冬的寒冷。 有时候，短短几天，可能改变的就是一生。 就好像，突如其来的知道，他要去远方。 就好像，毫无征兆的听说，奶奶要...
丹麦是世界上最古老的君主国，丹麦王室则是欧洲最古老的王室。从丹麦的第一位君主高姆老国王即任之后的一千多年里，丹麦整个国家的历史都与丹麦王室的兴迭紧密地关联在一起。所以可以说，丹麦的历史就是丹麦国王们的历史。 高姆老国王的儿子哈拉尔国王可以说是丹麦国家的正式奠基人，被称为蓝牙...
《早安诗002o一个诗人的礼物》 文/刘汉皇 我住在诗歌的隔壁，隔壁的隔壁有一个花园清晨，花园里有鸟儿在唱歌那歌声比诗人的吟咏更动听，唱的整个三月都在枝头，田野舞动花园有一个池塘，有春光在游动鱼儿一样。我在墙头独坐闻着花园的幽香。像一个诗人在静静的夜借星光入酒，收月色入诗渲...
你知道《大宅门》中白景琦常唱的那几句：“看前面，黑洞洞，定是那贼巢穴；待俺赶上前去，杀他个干干净净！”是出自哪里的戏词呢？ 你知道“青山有幸埋忠骨，白铁无辜铸佞臣。”这句诗中“忠骨”指的是谁么？ 让我来告诉大家吧！ ——白景琦唱的戏词是出自戏曲名段《挑滑车》中的一段台词。剧...
只有正确可以原谅。SAM/BAM格式文件操作软件samtools――二代测序测序分析必备_邓飞龙博客
SAM/BAM格式文件操作软件samtools――二代测序测序分析必备
SAM和BAM是序列比对之后常用的输出格式，比如tophat输出BAM格式，bowtie和bwa等都采用了SAM格式。BAM格式其实就是SAM格式的二进制格式，占用存储空间更小。samtools由中国学者开发，专门用于sam/bam格式文件的各种操作。
1）查看BAM/SAM格式文件
& & $&samtools&view&[\bhuHS] [\t in.refList] [\o output] [\f reqFlag] [\F skipFlag] [\q minMapQ] [\l library] [\r read\Group] &in.bam&|&in.sam&
& & [region1 [...]]，提取/打印所有的或者部分序列文件。如果没有指定区域，所有的序列都将会被打印。否则，仅仅打印覆盖特定区域的序列才会被打印出来。如果一个序列覆盖多个区域，那么它将会被多次打印出来。一个区域可以以如下例格式呈现：`chr2' (the whole chr2), `chr2:; (region starting from 1,000,000bp) or `chr2:1,000,000\2,000,000' (region between 1,000,000 and 2,000,000bp including the end points).（即染色体名：起始打印位点..终止打印位点）位置是以1开始记录第一个碱基的方式定位（1-based coordinate）。
&&&&&&&& -b 以BAM格式输出，可以用于samtools的后续分析
&&&&&&&& -u 以未压缩的BAM格式输出，可以节约时间，一般在管道执行时使用
&&&&&&&& -h 在结果中包含头header
-H 只输出头
&&&&&&&& -S 输入文件为SAM格式，如果确实@SQ头，则需要-t选项
-t FILE 这个文件是以TAB分隔的，每行必须包含参考序列名字（如染色体名），一行只包含单个参考序列名字，其它区域被忽略。这个文件定义排序是的参考序列顺序。如果你运行samtools faidx &ref.fa&，得到的索引结果文件可以作为这里的FILE文件（见6）。
&&&&&&&& -o FILE 输出文件[stdout]，如果在该管道中可以使用command1|tee FILE|command2。
& & & & -f INT 只输出在FLAG中含有整个INT的序列，INT可以使用十六进制的/^0x[0\9A\F]+/ [默认0]
&&&&&&&&-F INT 跳过含有INT的序列
&&&&&&&&-q INT 跳过MAPQ（定位的质量）比INT小的序列[默认0]
&&&&&&&&-l STR 只输出STR库中的序列
&&&&&&&&-r STR 只输出在STR读段组中的序列
2）samtools&import&&in.ref_list& &in.sam& &out.bam&，从0.1.4起，它与samtools view \bt &in.ref_list& \o &out.bam& &in.sam&意义相同。
3）samtools&sort&[\n] [\m maxMem] &in.bam& &out.prefix&，根据左起位点对序列排序，将产生&out.prefix&.bam文件，当整个序列无法完全装载到内存（\m选项控制）时，可能还会产生&out.prefix&.%d.bam这个临时文件。
-n 设定排列方式为读段名称而非染色体位置
-m 设定内存使用量[](默认500M)
4）samtools&merge&[\h inh.sam] [\n] &out.bam& &in1.bam& &in2.bam& [...]，将多个排序后的序列文件合并。参考头文件列出所有的输入BAM文件和inh.sam中@SQ头，如果存在，必须全部指向同一系列的参考序列。参考头文件列表（除非被-h选项替代）和in1.bam中的&@&头将会被拷贝到out.bam中，其他的输入文件中的头将会被忽略。
-h FILE 使用FILE中的行作为@头拷贝到out.bam中，代替从in1.bam中拷贝的头。（FILE实际上是BAM格式的。但是其中的任何序列信息都会被忽略）。
-n 输入文件是以读段名称排序的，而非染色体定位（当sort中使用-n选项时，此处应该选-n）。
5）samtools&index&&aln.bam& 对排序后的序列索引，用于快速的随机处理，此处将产生&aln.bam&.bai文件。
6）samtools&faidx&&ref.fasta& [region1 [...]]，对FASTA格式的参考基因组序列建立索引或者提取已经索引的序列中的子序列（subsequence）。如果没有指定区域，它将会索引文件，并在磁盘上建立&ref.fasta&.fai，子序列将会被检索并且以FASTA格式打印在标准输出。输入文件可以压缩成RAZF格式。
7）samtools&pileup&[\f in.ref.fasta] [\t in.ref_list] [\l in.site_list] [\iscgS2] [\T theta] [\N nHap] [\r pairDiffRate] &in.bam&|&in.sam&，以pileup格式输出序列，在pileup格式中，每一行代表一个基因组位点，包括染色体名称、位点、参考碱基、读段质量和序列定位质量。匹配、错配、插入和缺失、链、mapping质量以及读段的开始和结束位点都在读段碱基一列。在这一列中，点代表与参考序列的正向链匹配，逗号表示与参考序列的反向链匹配，ACTGN代表正向链上的错配，而actgn代表反向链上的错配，+[0\9]+[ACGTNacgtn]+这一模式表示在参考序列之间存在插入序列，模式中，数字代表插入长度，其后为插入序列。\[0\9]+[ACGTNacgtn]+这一模式代表与参考序列相比存在缺失，后面以*代表被删除的碱基。在读段碱基列中，符号^标记一个读段片段的起始，它是读段中被N/S/H CIGAR操作符分隔开临近的子序列（subsequence）。^后面紧跟的ASCII码字符减去33给出mapping质量。$符号标记读段片段（read segment）的结尾。（注意-c选项不同的输出）
如果应用-c选项，一致碱基、Phred-scaled一致性质量、SNP质量（指与参考序列相同的Phred-scaled一致性的可能性）和包含该位点的读段mapping质量均方根（root mean square RMS）将会被插入到参考碱基和读段序列之间（即会多出几列）。插入缺失将占有额外一列。
每个插入缺失行包含染色体名称、位点、*、基因型、一致性质量、SNP质量、mapping质量均方根，#支持第一种等位序列的读段，#支持第二种等位序列的读段，#与前两种等位形式不同且包含插入缺失的读段。
&&&&&&&& -s 最后一行打印mapping质量。这一选项使得输出结果易于分析，但较耗费空间
-S 输入文件时SAM格式
-i 只输出含有插入缺失的pileup行
-f FILE FASTA格式的参考序列，如果缺少索引文件将产生FILE.fai
-M INT设定mapping质量上限（cap）为INT
-t FILE以import 命令描述的格式列出参考名称和序列长度列表。如果出现这一个选项，samtools将假设输入文件为SAM格式；否则假设为BAM格式。
-l(L) FILE指定pileup输出的位点，第一列为染色体，第二列为1-based位点。其他的列将会被忽略，推荐与-s选项同时使用，因为默认格式我们可能无法知道mapping质量
-c 使用MAQ一致性模型检测一致序列，在使用-g和-c选项时，仅-T、-N、-l和-r选项可用
-g 产生GLFv3格式的基因型相似性，这一选项抑制-c、-i和-s选项
-T FLOAT MAQ一致性检测模式的theta参数[0.85]（错误依赖系数）
-N INT样本中的单倍型数量[默认2，要求&=2]
-r FLOAT 每一对单倍体之间期望的片段差异值[0.001]
-I(i) INT 序列中插入缺失的Phred概率[40]
8）samtools&tview&&in.sorted.bam& [ref.fasta]，文本定位查看器，在查看器中，使用？取得帮助，按下g从一个区域开始位置检查定位（alignment），格式如chr10:10,000,000。
9）samtools&fixmate&&in.nameSrt.bam& &out.bam&，为以名称排序的定位alignment填入配对坐标，ISIZE（inferred insert size猜测的插入序列大小）和配对相应的标签（flag）
10）samtools&rmdup&&input.srt.bam& &out.bam&，移除潜在的PCR重复：如果多个读段含有相同的外部位点（external coordinates），只保留那些高mapping质量的碱基对。这一命令只适用于FR开始的并且要求正确设置ISIZE。
11）samtools&rmdupse&&input.srt.bam& &out.bam&，除去单端读段（single-ended reads）。这一命令将把所有的读段当作单端读段，即使它们实际上是配对的。
12）samtools&fillmd&[\e] &aln.bam& &ref.fasta&，产生MD标签，如果已经存在并且产生的MD标签不同，这个命令将给予警告。
-e 如果读段碱基与参考序列相同，将其转换成=，indel检测目前还不支持=作为碱基标志。
SAM是以TAB分隔的，除去以@开始的header行，每一定位行（alignment）包含以下项目：
|Col | Field | Description |
+\\\\+\\\\\\\+\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\+
| 1 | QNAME | Query (pair) NAME 读段对应模板名称|
| 2 | FLAG | bitwise FLAG 按位的标签（意义见后面）|
| 3 | RNAME | Reference sequence NAME 参考序列名称|
| 4 | POS | 1\based leftmost POSition/coordinate of clipped sequence 1-based左起位置/整齐的序列定位|
| 5 | MAPQ | MAPping Quality (Phred\scaled) 读段定位质量|
| 6 | CIAGR | extended CIGAR string 扩展的CIGAR字符串（详见SAM格式详解）|
| 7 | MRNM | Mate Reference sequence NaMe (`=' if same as RNAME) 配对的参考序列名称，=表示相同|
| 8 | MPOS | 1\based Mate POSistion 参考序列中1-based配对位置|
| 9 | ISIZE | Inferred insert SIZE 猜测的插入序列大小|
|10 | SEQ | query SEQuence on the same strand as the reference 与参考序列处于同一链的读段序列|
|11 | QUAL | query QUALity (ASCII\33 gives the Phred base quality) 度短序列的碱基质量ASCII码-33为碱基质量|
|12 | OPT | variable OPTional fields in the format TAG:VTYPE:VALUE 变量选项，格式TAG:VTYPE:VALUE |
+\\\\+\\\\\\\+\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\+
下表列出了FLAG列的含义：
+\\\\\\\+\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\+
| Flag | Description |
+\\\\\\\+\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\+
|0x0001 (1)| the read is paired in sequencing |读段序列是成对的
|0x0002 (2)| the read is mapped in a proper pair |读段定位在适当位置
|0x0004 (4)| the query sequence itself is unmapped |读段序列自身没有定位
|0x0008 (8)| the mate is unmapped |与其配对的读段为定位
|0x0010 (16)| strand of the query (1 for reverse) |读段对应链
|0x0020 (32)| strand of the mate |配对链
|0x0040 (64)| the read is the first read in a pair |读段是读段对的第一个
|0x)| the read is the second read in a pair |读段是读段对的第二个
|0x)| the alignment is not primary |定位不是最优选
|0x)| the read fails platform/vendor quality checks |读段质量未生成
|0x)| the read is either a PCR or an optical duplicate |读段是PCR或者光学重复
+\\\\\\\+\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\+
1、& 在bam_import.c, bam_endian.h, bam.c和bam_aux.c中的非定位单词
2、& CIGAR操作符P不能正确处理
3、& 合并时，输入文件需要有相同数目的参考基因序列。设备要求可以降低，另外，合并时不会自动重构头字典（header dictionary）。用户必须提供正确的头，Picard做合并表现更好。
4、& samtools的rmdup无法处理单端数据，不能移除染色体之间的重复。Picard表现更好。
Heng Li from the Sanger Institute wrote the C version of samtools. Bob
Handsaker from the Broad Institute implemented the BGZF library and Jue
Ruan from Beijing Genomics Institute wrote the RAZF library. Various
people in the 1000Genomes Project contributed to the SAM format specification.
Samtools website: &http://samtools.sourceforge.net&
主要转载自：http://www.plob.org//1700.html
[相关日志]
  ( 58 )
  ( 21 )
  ( 7 )
  ( 15 )
  ( 7 )