sam文件header全面解析

发表于 | 分类于

众说周知，最常用的比对软件 bwa 产生的结果文件是 sam 格式，具体格式可以参考 sam 。

其中文件的主体便是比对的结果，而头文件往往是我们最容易忽视的地方，尤其踩过太多的坑，特来解释和记录其 header。

header 中最常见的有3个：@SQ @RG @PG ，@SQ 用来记录参考序列的信息，@RG 用来记录样本的信息，@PG 用来记录 bwa 程序信息。

@SQ 和 @PG 在 bwa 比对之后就会在 sam 文件中存在。而 @RG 则需要自己在 bwa mem 比对的命令中使用 -R 参数来添加。@RG 这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要 。

如果在 bwa 比对期间没有选择 -R 参数，可以 picard 中 AddOrReplaceReadGroups 对 sam 或bam 文件进行加头处理。 

1
2
3
4
5
6
7
8
9
10
-R '@RG\tID:sample\tLB:sample\tSM:sample\tPL:ILLUMINA'   #bwa mem @RG头格式

java -jar picard.jar AddOrReplaceReadGroups \   #sam或bam文件（picard ）
      I = input.bam \ 
      O = output.bam \ 
      RGID = 4 \ 
      RGLB = lib1 \ 
      RGPL = illumina \ 
      RGPU = unit1 \ 
      RGSM = 20
1
2
3
4
5
6
7
8
9
10
11
$ bwa mem -t 4 genome.fa ./samples/A.fastq >A.sam   #非常普通的SE数据bwa比对
$ head -n 3 A.sam
@SQ     SN:I    LN:230218
@PG     ID:bwa  PN:bwa  VN:0.7.17-r1188 CL:bwa mem -t 4 genome.fa -R 
SRR800764.1     4       *       0       0       *       *       0       0       CATCTTTGGAGTAACTATTATTTCGCCCCTTTTGTTTGCTGCATATCGCCCCGCTCTCTGCATACACGATTGGATAATGACCAAAGCAAGGTTTAATACGC…………   #关键看头文件，后面比对信息省略
$ samtools view -Sb A.sam > A.bam   #转为bam文件
$ samtools sort  A.bam > A_sort.bam   #排序
$ samtools view -h A_sort.bam > A_sort.sam   #再转回来sam文件，一定要加-h参数，默认是不输出header的。如果没有-h参数，是无法对sam文件进行转换和排序的。如下：
[E::sam_parse1] missing SAM header
[W::sam_read1] parse error at line 1
[main_samview] truncated file.

上面提到在缺少头文件的前提下是无法将 sam 转为 bam 的，更加无进行 bam 文件排序。但是如果我们使用 picard 这个问题就可以解决了。

1
java -Xmx5g -Djava.io.tmpdir=./tmp -jar ../picard-tools-1.117/SortSam.jar INPUT=A_sort.sam OUTPUT=A_sort.bam SORT_ORDER=coordinate VALIDATION_STRINGENCY=SILENT

picard 的 SortSam.jar 可以直接对没有 header 的 sam 文件排序并转为 bam 文件。

另外，samtools 可以为sam或者bam文件加上 header 中的 @SQ：

1
samtools view -Sb -T genome.fa sample_no_header.sam >sample_with_header.bam  #使用 -T 参考序列加header（@SQ）