目录

1. samtools和picard的排序问题
2. SAM文件中FLAG值的理解
3. SAM文件中那些未比对的reads
4. 为什么一条read会有多条比对记录？
5. 同一份BAM文件samtools depth和samtools mpileup的结果不同？
6. 为什么samtools flagstat与hisat2的mapping rate不同？

1. samtools和picard的排序问题

samtools和picard都有对SAM/BAM文件进行排序的功能，一般都是基于坐标排序（还提供了-n选项来设定用reads名进行排序），先是对chromosome/contig进行排序，再在chromosome/contig内部基于start site从小到大排序，对start site排序很好理解，可是对chromosome/contig排序的时候是基于什么标准呢？

基于你提供的ref.fa文件中的chromosome/contig的顺序。当你使用比对工具将fastq文件中的reads比对上参考基因组后会生成SAM文件，SAM文件包含头信息，其中有以@SQ开头的头信息记录，reference中有多少条chromosome/contig就会有多少条这样的记录，而且它们的顺序与ref.fa是一致的。

SAM/BAM文件的头信息：

@HD     VN:1.3  SO:coordinate
@SQ     SN:chr1 LN:195471971
@SQ     SN:chr2 LN:182113224
@SQ     SN:chr3 LN:160039680
@SQ     SN:chr4 LN:156508116
@SQ     SN:chr5 LN:151834684
@SQ     SN:chr6 LN:149736546
@SQ     SN:chr7 LN:145441459
@SQ     SN:chr8 LN:129401213
@SQ     SN:chr9 LN:124595110
@SQ     SN:chr10        LN:130694993
@SQ     SN:chr11        LN:122082543
@SQ     SN:chr12        LN:120129022
@SQ     SN:chr13        LN:120421639
@SQ     SN:chr14        LN:124902244
@SQ     SN:chr15        LN:104043685
@SQ     SN:chr16        LN:98207768
@SQ     SN:chr17        LN:94987271
@SQ     SN:chr18        LN:90702639
@SQ     SN:chr19        LN:61431566
@SQ     SN:chrX LN:171031299
@SQ     SN:chrY LN:91744698
@SQ     SN:chrM LN:16299
@RG     ID:ERR144849    LB:ERR144849    SM:A_J  PL:ILLUMINA

ref.fa中chromosome/contig的排列顺序：

>chr1
>chr2
>chr3
>chr4
>chr5
>chr6
>chr7
>chr8
>chr9
>chr10
>chr11
>chr12
>chr13
>chr14
>chr15
>chr16
>chr17
>chr18
>chr19
>chrX
>chrY
>chrM

它们的顺序一致

当使用samtools或picard对SAM/BAM文件进行排序时，这些工具就会读取头信息，按照头信息指定的顺序来排chromosome/contig。所以进行排序时需要提供包含头信息的SAM/BAM文件。

那么普通情况下我们的chromosome/contig排序情况是什么样的?

一般情况下我们获取参考基因组序列文件的来源有三个：

• NCBI
• ENSEMBEL
• UCSC Genome Browser

这里以UCSC FTP下载源为例：

这是一个压缩文件，使用tar zxvf chromFa.tar.gz解压后，会得到多个fasta文件，每条chromosome/contig一个fasta文件：chr1.fa, chr2.fa ...

之后我们会将它们用cat *.fa >ref.fa合并成一个包含多条chromosome/contig的物种参考基因组序列文件

用grep ">" ref.fa可以查看合并后发ref.fa文件中染色体的排列顺序为：

>chr10
>chr11
>chr12
>chr13
>chr14
>chr15
>chr16
>chr17
>chr18
>chr19
>chr1
>chr1_GL456210_random
>chr1_GL456211_random
>chr1_GL456212_random
>chr1_GL456213_random
>chr1_GL456221_random
>chr2
>chr3
>chr4

这和我们平时想象的染色体的排列顺序是不是有一些出入？难道不应该是从chr1开始到chr22，最后是chrX和chrY这样的顺序吗？

想象归想象，实际上它是按照字符顺序进行的，chr11就应该排在chr2前面

一般情况下在进行SAM文件的排序时，染色体的排序到底是按照哪种规则进行排序的，不是一个很重要的问题，也不会对后续的分析产生影响，但是在执行GATK流程时，GATK对染色体的排序是有要求的，必须按照从chr1开始到chr22，最后是chrX和chrY这样的顺序，否则会报错

面对这样变态的要求，我们怎么解决？

在构造ref.fa文件时，让它按照从chr1开始到chr22，最后是chrX和chrY这样的顺序进行组织就可以了：

for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg19.fasta;
done

2. SAM文件中FLAG值的理解

FLAG列在SAM文件的第二列，这是一个很重要的列，包含了很多mapping过程中的有用信息，但很多初学者在学习SAM文件格式的介绍时，遇到FLAG列的说明，常常会一头雾水

what?还二进制，这也太反人类的设计了吧！

不过如果你站在开发者的角度去思考这个问题，就会豁然开朗

在mapping过程中，我们想记录一条read的mapping的信息包括：

• 这条read是read1 (forward-read) 还是read2 (reverse-read)？
• 这条read比对上了吗？与它对应的另一头read比对上了吗？
• ...

这些信息总结起来总共包括以下12项：

而每一项又只有两种情况，是或否，那么我可以用一个12位的二进制数来记录所有的信息，每一位表示某一项的情况，这就是原始FLAG信息的由来，但是二进制数适合给计算机看，不适合人看，需要转换成对应的十进制数，也就有了我们在SAM文件中看到的FLAG值

但是FLAG值所包含信息的解读还是要转换为12位的二进制数

3. SAM文件中那些未比对的reads

SAM格式文件的第3和第7列，可以用来判断某条reads是否比对成功到了基因组的染色体，左右两条reads是否比对到同一条染色体

有两个方法可以提取未比对成功的测序数据：

• SAM文件的第3列是*的(如果是PE数据，需要考虑第3,7列)

$ samtools view sample.bam | perl -alne '{print if $F[2] eq "*" or $F[6] eq "*" }' sample.unmapped.sam

• 或者SAM文件的flag标签包含0x4的

# 小写的f是提取，大写的F是过滤
$ samtools view -f4 sample.bam sample.unmapped.sam

虽然上面两个方法得到的结果是一模一样的，但是这个perl脚本运行速度远远比不上上面的samtools自带的参数

对于PE数据，在未比对成功的测序数据可以分成3类：

• 仅reads1没有比对成功

该提取条件包括：

• 该read是read1，对应于二进制FLAG的第7位，该位取1，其十进制值为64；
• 该read未成功比对到参考基因组，对应于二进制FLAG的第3位，该位取1，其十进制值为4；
• 另一配对read成功比对到参考基因组，对应于二进制FLAG的第4位，该位取0，其十进制值为8；

# 对于取1的位点采用提取的策略，用-f参数，值设为64+4=68
# 对于取0的位点采取过滤的策略，用-F参数，值设为8
$ samtools view -u -f 68 -F 8 alignments.bam >read1_unmap.bam

• 仅reads2没有比对成功

该提取条件包括：

• 该read是read2，对应于二进制FLAG的第8位，该位取1，其十进制值为128；
• 该read未成功比对到参考基因组，对应于二进制FLAG的第3位，该位取1，其十进制值为4；
• 另一配对read成功比对到参考基因组，对应于二进制FLAG的第4位，该位取0，其十进制值为8；

# 对于取1的位点采用提取的策略，用-f参数，值设为128+4=132
# 对于取0的位点采取过滤的策略，用-F参数，值设为8
$ samtools view -u -f 132 -F 8 alignments.bam >read2_unmap.bam

• 两端reads都没有比对成功

该提取条件包括：

• 该read未成功比对到参考基因组，对应于二进制FLAG的第3位，该位取1，其十进制值为4；
• 另一配对read未成功比对到参考基因组，对应于二进制FLAG的第4位，该位取1，其十进制值为8；

# 对于取1的位点采用提取的策略，用-f参数，值设为4+8=12
$ samtools view -u -f 12 alignments.bam >pairs_unmap.bam

看完这一部分，是不是有一个感觉：FLAG玩得溜，SAM文件可以处理得出神入化

4. 为什么一条read会有多条比对记录？

首先，思考一个问题：对于PE数据，一条测序片段（fragment）有read1和read2两条测序片段，它们俩的名字相同，那么对于这一条测序片段，对它进行mapping之后得到的SAM文件中会出现几条记录呢？

先声明以下只对BWA比对得到的SAM文件进行讨论，对其他比对工具输出的SAM文件可能不适用

首先，基于经验积累告诉我，它会得到有且只有两条记录，原因在于，BWA在对每条read执行比对时只会给出一个hit，若这条read是multiple mapping的情况，它会从中选择MAPQ值最高的那个hit作为输出，若存在多个hits的MAPQ值相等且最高，那么BWA会从中随机选择一个作为输出

对于我的这个假设可以用以下的方法进行验证：

# 将SAM文件的第一列提出来，排序去重，同时统计每个QNAME出现的次数
$ samtools view alignment.bam | cut -f1 | sort | uniq -c >record.count

得到的统计结果如下：
      2 ERR144849.1
      2 ERR144849.10
      2 ERR144849.100
      2 ERR144849.1000
      2 ERR144849.10000
      2 ERR144849.100000
      2 ERR144849.1000000
      2 ERR144849.10000000
      2 ERR144849.10000001
      2 ERR144849.10000002
      2 ERR144849.10000003
      2 ERR144849.10000004
      2 ERR144849.10000005
      ...

上面的测试结果与我们的假设吻合

但是在一次处理三代测序数据（三代测序数据是Single-End）中发现了不同：

在输出中出现了一些不太和谐的结果：有极少部分的QNAME对应2条以上的记录，这意味着存在一条read会有多条比对记录的情况，why？

对这个与预期不完全相符的结果，尝试去寻找里面的原因，其间进行了各种各样的推理、假设、验证，最终在 李恒的github 中找到了答案

2. Why does a read appear multiple times in the output SAM?

BWA-SW and BWA-MEM perform local alignments. If there is a translocation, a gene fusion or a long deletion, a read bridging the break point may have two hits, occupying two lines in the SAM output. With the default setting of BWA-MEM, one and only one line is primary and is soft clipped; other lines are tagged with 0x800 SAM flag (supplementary alignment) and are hard clipped.

这种情况容易在三代测序数据中出现

5. 同一份BAM文件samtools depth和samtools mpileup的结果不同？

如果你用的是Single-End的数据，那么差异应该比较小，不会太明显，而在Pair-End上差异可能会比较大，之所以会产生这些差异，原因有两点：

（1）mpileup会默认将PE reads中比对情况异常的那些reads（包括双端都比上，但是两条配对reads之间的比对距离明显偏离了插入片段的长度分布，或者一端比对上而另一端没比对上）丢弃，则在计算depth时，这些被丢弃的reads不参与统计，而depth则不会，depth默认不做任何过滤

mpileup中提供了一个-A选项来保留异常比对的reads，如果设置了这个选项，那么在这点上mpileup的depth统计就和samtools depth相同了

（2）从上图的第二个红框可以看出，在默认情况下，mpileup还会过滤掉测序质量值低于13的碱基，depth默认不过滤

从上面列出的两点差异可以看出，mpileup默认输出的是高质量的覆盖深度，这是有历史原因的：当场mpileup功能被开发出来就是为了与bcftools组合，将samtools mpileup的输出作为bcftools的输入用于下游的snp-calling，当然需要保证数据的质量

当然可以通过设置对应的参数使得它的属于结果与depth的一致，但是不推荐这么做

6. 为什么samtools flagstat与hisat2的mapping rate不同？

下面是对同一个样本的paired-end Fastaq文件比对结果（比对使用hisat2），hisat2和samtools分别给出的mapping rate的统计

hisat2：

39928651 reads; of these:
  39928651 (100.00%) were paired; of these:
    8439896 (21.14%) aligned concordantly 0 times
    30158798 (75.53%) aligned concordantly exactly 1 time
    1329957 (3.33%) aligned concordantly >1 times
    ----
    8439896 pairs aligned concordantly 0 times; of these:
      298072 (3.53%) aligned discordantly 1 time
    ----
    8141824 pairs aligned 0 times concordantly or discordantly; of these:
      16283648 mates make up the pairs; of these:
        11658022 (71.59%) aligned 0 times
        4210488 (25.86%) aligned exactly 1 time
        415138 (2.55%) aligned >1 times
85.40% overall alignment rate

samtools flagstat：

85207970 + 0 in total (QC-passed reads + QC-failed reads)
5350668 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
73549948 + 0 mapped (86.32% : N/A)
79857302 + 0 paired in sequencing
39928651 + 0 read1
39928651 + 0 read2
62977510 + 0 properly paired (78.86% : N/A)
63902424 + 0 with itself and mate mapped
4296856 + 0 singletons (5.38% : N/A)
173078 + 0 with mate mapped to a different chr
127302 + 0 with mate mapped to a different chr (mapQ>=5)

从上面可以看出，hisat2给出的mapping rate为85.40%，而samtools给出的为86.32%，两个的统计结果不一样，而且samtools的统计会大一些，what？

介四什么鬼？

我们来简单地分析一下：

hisat2中，

samtools中，

计算没问题，那问题出在哪呢？

有没有注意到上面的两个式子中的分子和分母，计算它们的差值：

分子：

分母：

发现了没有，它们的差值正好都等于samtools flagstat的输出结果的第二行：

5350668 + 0 secondary

所以，hisat2和samtools计算mapping rate的公式实际上分别为：

• hisat2

  

• samtools

  

  其中，recorder number表示SAM文件中除去头信息部分的比对记录数，每一行是一条比对记录

一般来说，我们想得到的是hisat2计算公式所得到的统计结果，hisat2统计结果在比对结束后会以标准错误形式给出，我们可以将标准错误重定向到一个log文件中，但是如果我们忘了保持这个统计结果，怎么办？

最简单的办法就是重新跑一遍hisat2，但是这样太耗费时间和计算资源了，这时我们可以利用samtools flagstat对SAM文件的统计结果，以及它的部分统计值与hisat2计算公式的关系，快速地算出准确的mapping rate：

参考资料：

(1) 【简书】从零开始完整学习全基因组测序数据分析：第5节 理解并操作BAM文件

(2) 【生信技能树】【直播】我的基因组（十五）:提取未比对的测序数据

(3) BWA's README in github

(4) 【简书】黄树嘉《样本量重要，还是测序深度重要? 生物信息工程师可以分为多少种类型? |《解螺旋技术交流圈》精华第3期》

(5) 生信媛《HISAT2的比对率计算结果和SAMTools flagstat不同，你想过原因吗？ 》