SAM/BAM相关的进阶知识

samtools和picard的排序问题

samtoolspicard都有对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文件的头信息:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
> @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的排列顺序:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
> >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文件中染色体的排列顺序为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
> >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这样的顺序进行组织就可以了:

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

SAM文件中FLAG值的理解

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

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

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

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

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

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

序号 简写 说明
1 PAIRED paired-end (or multiple-segment) sequencing technology
2 PROPER_PAIR each segment properly aligned according to the aligner
3 UNMAP segment unmapped
4 MUNMAP next segment in the template unmapped
5 REVERSE SEQ is reverse complemented
6 MREVERSE SEQ of the next segment in the template is reversed
7 READ1 the first segment in the template
8 READ2 the last segment in the template
9 SECONDARY secondary alignment
10 QCFAIL not passing quality controls
11 DUP PCR or optical duplicate
12 SUPPLEMENTARY supplementary alignment

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

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

SAM文件中那些未比对的reads

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

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

  • SAM文件的第3列是*的(如果是PE数据,需要考虑第3,7列)
1
2
> $ samtools view sample.bam | perl -alne '{print if $F[2] eq "*" or $F[6] eq "*" }' sample.unmapped.sam
>
  • 或者SAM文件的flag标签包含0x4的
1
2
3
> # 小写的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
2
3
4
> # 对于取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
2
3
4
> # 对于取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
2
3
> # 对于取1的位点采用提取的策略,用-f参数,值设为4+8=12
> $ samtools view -u -f 12 alignments.bam >pairs_unmap.bam
>

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

为什么一条read会有多条比对记录?

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 将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)中发现了不同:

1
2


在输出中出现了一些不太和谐的结果:有极少部分的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.


参考资料:

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

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