sam格式的简单了解

Posted 2023-05-05

参考技术A

部分转自生信菜鸟团公众号

SAM的全称是sequence alignment/map format。而BAM就是SAM的二进制文件(B源自binary)。
SAM 格式主要包括两大部分:
1 .标头注释部分（header section）
2 .比对结果部分（alignment section）

SAM格式是用来来支持高通量测序数据分析：

(1):快速查找与坐标重叠的比对。例如，选择与染色体2上的坐标323,567,334重叠的比对。

(2):根据read的属性进行选择和过滤。例如，我们希望能够快速选择能过比对到反向链上的read。

(3):有效地存储数据。例如，从SAM格式转化成BAM格式，单个压缩文件包含所有样本的数据，每个样本都以某种方式标记。

标头注释部分
标头信息可有可无，都是以@开头，用不同的tag表示不同的信息

比对结果部分
每一列表示一个read的比对信息，包括11个必须的字段和一个可选字段，字段之间用tag分割。必须的字段有11个，顺序固定。这11个字段包括：

第一列: Query Name (QNAME) ：
这一列代表着比对片段的（template）的编号

第二列：FLAG ：
这是一种常用且高效的保存多个布尔特征值的方法。

举个简单的例子: 在 SAM 格式中，当 flag 为 1，也即对应的二进制为 01 时，表示该 read 有多个测序数据 ， 一般理解为有双端测序数据 (另一条没被过滤掉)， 而 flag 为 2， 也即二进制 10 时， 表示这条 read 的多个片断都有比对结果， 通常理解为双端 reads 都比对上了， 那么就可以推断出 flag 为 3 时， 也即二进制的 11， 表示该 read 有另一端的 read 并且比对成功， 可以看到， 其实就是 01 加 10。
一般flag值不需要自己去算，直接将flag值导入网站即可
http://broadinstitute.github.io/picard/explain-flags.html

所有flag对应值的含义
1 ： 代表这个序列采用的是PE双端测序

2： 代表这个序列和参考序列完全匹配，没有插入缺失

4： 代表这个序列没有mapping到参考序列上

8： 代表这个序列的另一端序列没有比对到参考序列上，比如这条序列是R1,它对应的R2端序列没有比对到参考序列上

16：代表这个序列比对到参考序列的负链上

32 ：代表这个序列对应的另一端序列比对到参考序列的负链上

64 ： 代表这个序列是R1端序列， read1;

128 : 代表这个序列是R2端序列，read2；

256： 代表这个序列不是主要的比对，一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置，其他都是次要的

512： 代表这个序列在QC时失败了，被过滤不掉了（# 这个标签不常用）

1024: 代表这个序列是PCR重复序列（#这个标签不常用）

2048: 代表这个序列是补充的比对（#这个标签具体什么意思，没搞清楚，但是不常用）

第三列： Reference Name (RNAME) ：
reference sequence name，实际上就是比对到参考序列上的染色体号。若是无法比对，则是*

第四列： Position (POS) ：
比对上的位置，注意是从1开始计数，没有比对上，此处为0

第五列：Mapping Quality (MAPQ) ：
比对的质量；比对的质量分数，越高说明该read比对到参考基因组上的位置越准确

第六列：Compact Idiosyncratic Gapped Alignment Representation (CIGAR)
CIGAR 代表着简要比对信息表达式，其以参考序列为基础，使用数字加字母表示比对结果
例如 3S6M1P1I4M
前三个碱基被剪切去除了，然后6个比对上了，然后打开了一 个缺口，有一个碱基插入，最后是4个比对上了。

这里的总长度对应的就是测出来的一条序列的长度，如果是PE100，这里加起来就应该是100，如果是PE150，这里加起来就是150，这里的长度和第10列的长度是一致的

第七列：RNEXT ：
双端测序中下一个reads比对的参考系列的名称。“＊”是完全没有比对上，“＝”代表完全比对

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

第八列：PNEXT ：
如果是双端测序，是指另一端匹配到参考基因组的位置，如果设置为0，那么该列不可用

第九列：TLEN Template的长度
最左边得为正，最右边的为负，中间的不用定义正负，不分区段（single-segment)的比对上，或者不可用时，此处为0

区别于第6列和第10列是对应测出来的序列的长度。这里第9列的长度是对应插入片段的长度，insert size，也就是建库时，将DNA片段打断成的长度。

第十列：Sequence ：
序列片段的序列信息，如果不存储此类信息，此处为’*‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度；就是read的碱基序列，如果是比对到互补链上则是reverse completed。
就是测序的reads序列信息

第十一列：ASCII ：
read质量值

其实很容易发现，如果将第1，10，11列提取出来的话，就能还原成我们常见的fastq格式信息。

第十二列：Optional fields ：
可选的区域
格式如：TAG:TYPE:VALUE，其中TAG有两个大写字母组成，每个TAG代表一类信息，每一行一个TAG只能出现一次，TYPE表示TAG对应值的类型，可以是字符串、整数、字节、数组等。

备注 ：
看一下KPGP-00001这个韩国人的测序reads比对到hg38的其中一个lane的sam格式部分信息：

可以看出这个是用的PE90测序，测序read长度为90bp，建库打断成约490bp，这个read名称是B80KJTABXX:4:1:1404:2065#CTAGTTAT，flag值是163，代表着

reads是比对到7号染色体，比对的位置是50962731，比对的质量值是60，"90M"意味着90个碱基都match（当然可能是mismatch），“=”意味着双端测序的另一条read也比对上，并且是比对到同一个片段，另一条read比对的位置是
50963137 ，这条read的序列信息是“AGAAAATTATTTAAATGACCCGAGCCTCACCTTCAACATGAGGAACATCAT
ATTCCCTTTGATAAAATGTTGCTGGTGTAAGTGCTCCAT ”
对应质量值分ASCII值为“gggfgfggeggdgggadegggdegegeeggeegedggegggeggegedgggedgggfggeceeggaedgcgggggecgaQ_`X``BBBBB ”

以上。

SAM/BAM文件处理

当测序得到的fastq文件map到基因组之后，我们通常会得到一个sam或者bam为扩展名的文件。SAM的全称是sequence alignment/map format。而BAM就是SAM的二进制文件(B取自binary)。 那么SAM文件的格式是什么样子的呢？如果你想真实地了解SAM文件，可以查看它的说明文档。SAM由头文件和map结果组成。头文件由一行行以@起始的注释构成。而map结果是类似下面的东西：

HWI-ST1001:137:C12FPACXX:7:1115:14131:66670 0 chr1 12805 1 42M4I5M * 0 0 TTGGATGCCCCTCCACACCCTCTTGATCTTCCCTGTGATGTCACCAATATG CCCFFFFFHHGHHJJJJJHJJJJJJJJJJJJJJJJIJJJJJJJJJJJJIJJ AS:i:-28 XN:i:0 XM:i:2 XO:i:1XG:i:4 NM:i:6 MD:Z:2C41C2 YT:Z:UU NH:i:3 CC:Z:chr15 CP:i:102518319 XS:A:+ HI:i:0
HWI-ST1001:137:C12FPACXX:7:2313:17391:30032 272 chr1 13494 1 51M * 0 0 ACTGCCTGGCGCTGTGCCCTTCCTTTGCTCTGCCCGCTGGAGACAGTGTTT [email protected] AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0NM:i:1 MD:Z:44G6 YT:Z:UU XS:A:+ NH:i:3 CC:Z:chr15 CP:i:102517626 HI:i:0
HWI-ST1001:137:C12FPACXX:7:1109:17518:53305 16 chr1 13528 1 51M * 0 0 CGCTGGAGCCGGTGTTTGTCATGGGCCTGGGCTGCAGGGATCCTGCTACAA #############AB=?:*B?;A?<2+233++;A+A2+<[email protected],A<A<=> AS:i:-5 XN:i:0 XM:i:2 XO:i:0 XG:i:0NM:i:2 MD:Z:8A21T20 YT:Z:UU XS:A:+ NH:i:4 CC:Z:chr15 CP:i:102517592 HI:i:0

看上去很类似fastq文件，它也有read名称，序列，质量等信息，但是又不完全一样。首先，每个read只占一行，只是它被tab分成了很多列，一共有12列，分别记录了：

1. read名称

2. SAM标记

3. chromosome

4. 5′端起始位置

5. MAPQ（mapping quality，描述比对的质量，数字越大，特异性越高）

6. CIGAR字串，记录插入，删除，错配以及splice junctions（后剪切拼接的接头）

7. mate名称，记录mate pair信息

8. mate的位置

9. 模板的长度

10. read序列

11. read质量

12. 程序用标记

显然，其中chromosome至CIGAR的信息都是非常重要的。但是这些对我们不重要，我们只需要了解SAM/BAM文件是什么，就可以了。重要的是如果进行下游的操作。 要操作SAM/BAM文件，首先需要安装samtools。它的安装过程和所有的linux/unix程序一样，都是经过make之后生成可执行程序，然后把它的路径告知系统，或者放在系统可以找到的位置就可以了。 比如：

 tar zxvf samtools-0.1.18.tar.bz2
 cd samtools-0.1.18/
 make
 samtoolpath=`pwd`
 PATH=PATH:$samtoolpath

然后就可以按照samtools主页上介绍的工具进行各种操作了。我们最常见的几步操作比如 0. SAM，BAM转换

 samtools view -h file.bam > file.sam
 samtools view -b -S file.sam > file.bam

1. sorting BAM文件。大多数下游程序都要求BAM文件是被排过序的。

samtools sort file.bam outputPrefix

2. 创建BAM index。这也是被大多数下游程序所要求。

samtools index sorted.bam

3. index模板基因组。这也是被大多数下游程序所要求。

samtools faidx Homo_sapiens_assembly19.fasta

在很多时候，我们还会看到一种扩展名为BED的mapping文件。其具体格式也是几经变化，但是现在以UCSC的描述为准。从BAM文件转换成BED文件，我们需要安装BEDtools。下载安装就不多说了。示例一个如何从BAM文件转换成BED文件的命令：

bamToBed -i reads.bam > reads.bed

更多的具体内容可以参见其说明文档。 当然，还有很多种格式来记录mapping的结果，大多数都收录在UCSC的帮助文档中。比如上次有人问及的.bw是什么文件（bigWig文件）之类的，都可以在那里找到答案。 上次谈及fastq文件时，有讲过其质量评估的问题，那么在mapping之后，如何对mapping的结果进行评估呢？ 最简单的，就是通过samtools来评估mapping质量了。

samtools idxstats aln.sorted.bam

注意，这一步之前需要经过sort和index。结果会显示：

 chr1 195471971 6112404 0
 chr10 130694993 3933316 0
 chr11 122082543 6550325 0
 chr12 120129022 3876527 0
 chr13 120421639 5511799 0
 chr14 124902244 3949332 0
 chr15 104043685 3872649 0
 chr16 98207768 6038669 0
 chr17 94987271 13544866 0
 chr18 90702639 4739331 0
 chr19 61431566 2706779 0
 chr2 182113224 8517357 0
 chr3 160039680 5647950 0
 chr4 156508116 4880584 0
 chr5 151834684 6134814 0
 chr6 149736546 7955095 0
 chr7 145441459 5463859 0
 chr8 129401213 5216734 0
 chr9 124595110 7122219 0
 chrM 16299 1091260 0
 chrX 171031299 3248378 0
 chrY 91744698 259078 0
 * 0 0 0

其中第一列是染色体名称，第二列是序列长度，第三列是mapped reads数，第四列是unmapped reads数。 如果是RNAseq，我们可以使用broad institute的RNA-SeQC来得到更加完整的报告。下载到文件之后，也许需要安装BWA来获取更精准的结果，但是如果不安装的话，也可以进行分析。一般来说，这一步不需要特别精准的结果，所以我很少使用BWA选项。下载的文件如果是.zip结尾的，直接把它改写成.jar就可以运行了。 在它的主页上下载所需要的Example RNA-seq Data。下载结束之后，该解压的解压缩。接下来运行：

 samtools index example/ThousandReads.bam
 samtools faidx example/Homo_sapiens_assembly19.fasta
 java -Xmx2048m -jar RNA-SeQC_v1.1.7.jar -n 1000 -s "TestId|example/ThousandReads.bam|TestDesc" -t example/gencode.v7.annotation_goodContig.gtf -r example/Homo_sapiens_assembly19.fasta -o ./testReport/ -start gc -gc example/gencode.v7.gc.txt

以上的参数只有一个与其说明文档不一样的地方就是使用了-Xmx2048m来指定java虚拟机的内存大小为2G。如果遇到java.lang.OutOfMemoryError，还可以指定得再大些。

当然如果是自己的文件的话，还需要多两步：

1.BAM，reference及GTF文件的基因组名称必须一致。

2.需要使用picard工具包中的CreateSequenceDictionary来构建一个dictionary文件。

原文来自：http://pgfe.umassmed.edu/ou/archives/3050

生物信息学交流论坛 http://bbs.bbioo.com/forum-76-1.html