61.《Bioinformatics Data Skills》之
2021-08-30 本文已影响0人
DataScience
我们已经知道 BED 与 GTF 文件用于存储基因组区域与基因组特征注释信息。还有一种文件用于存储高通量read比对数据，叫做SAM文件（Sequene Alignment Mapping，其压缩的二进制形式叫BAM）。SAM 文件应用非常广泛，有必要了解一下其格式。
SAM 头部
数据地址：https://github.com/vsbuffalo/bds-files/tree/master/chapter-11-alignment
这里以一个小 SAM 文件为例子（使用线虫的参考基因组进行模拟read对比创建），以 head 命令查看文件的前10行内容：
image-20210830205603602
SAM 文件的头部以@开头，紧接的两个字母代表信息的类别，接着以 tab 分隔的 key:value 对展示具体信息，以下是每行的详细解释：
- @SQ 代表参考基因组序列信息，所有分离的染色体都对应一行这样的信息。SN 代表序列名字，比如 SN:I 为 1 号染色体，LN 代表染色体长度。
- @RG 代表 read 组(read group)，这里的ID是唯一的字符串，代表 read 的来源。通常为产生read的通道，这样一些程序可能以此判断序列数据是否有批次效应。ID 是唯一必要的信息，其它可选信息包括SM代表的样本名，以及PL代表的测序技术如ILLUMINA, PACBIO等
- @PG 代表生成SAM文件的程序，ID 也是唯一的，VN 代表程序版本，一般是程序自动生成。
- 这行没有以@开头，所以是SAM文件的内容第一行，这部分的说明请参考下一节的内容。
对于 SAM/BAM 文件，标准的交互方式是采用 samtools 工具，通过这个工具的子命令 view 结合 -H 参数 查看 SAM 文件头部内容，如下：
$ samtools view -H celegans.sam
@SQ SN:I LN:15072434
@SQ SN:II LN:15279421
@SQ SN:III LN:13783801
@SQ SN:IV LN:17493829
@SQ SN:MtDNA LN:13794
@SQ SN:V LN:20924180
@SQ SN:X LN:17718942
@RG ID:VB00023_L001 SM:celegans-01
@PG ID:bwa PN:bwa VN:0.7.10-r789 CL:bwa mem -R @RG\tID:VB00023_L001\tSM:celegans-01 Caenorhabditis_elegans.WBcel235.dna.toplevel.fa celegans-1.fq celegans-2.fq
其中-H代表头部，此命令也可以直接用于bam文件：
$ samtools view -H celegans.bam
@SQ SN:I LN:15072434
@SQ SN:II LN:15279421
...
可以与Linux命令组合，例如查找@RG信息：
$ samtools view -H celegans.bam | grep "^@RG"
@RG ID:VB00023_L001 SM:celegans-01
如果不加任何参数的话则返回不包括头部的内容：
$ samtools view celegans.sam | head -n 1
I_2011868_2012306_0:0:0_0:0:0_2489 83 I 2012257 40 50M [...]
