008-处理SAM、BAM你需要Samtools

Jul 5, 2018 7 min read cnposts

关于SAM

Sam: I Want You!

I Want You

sam是短序列比对默认的标准格式，由sanger制定，是以TAB为分割符的文本格式。主要应用于测序序列mapping到基因组上的结果表示，另外也可以表示其他的多重比对结果。一般把测序reads比对到参考基因组以后，通常得到的就是sam文件，全称是sequence alignment/map format，BAM就是SAM的二进制文件(B即：Binary) 官方参考文档http://samtools.github.io/hts-specs/SAMv1.pdf

需要知道的名词

**template：**DNA/RNA序列的一部分，用它进行测序或者由原始数据拼接得到它，作为研究的模版； 【参考序列和比对上的序列共同组成的序列为template】

read： 测序仪下机得到的原始数据，这些数据依照测序时间的先后被打上不同标签；双端测序存在read1，read2

**segment：**比对到read的时候的一段连续序列或者被分开几条子序列。一条read可以包含多条segment；

linear alignment: 一条read比对到参考序列上，可以存在插入(insert)、缺失(delete)、跳跃(skip)、剪切(clip)，但是方向不变（不能是一部分和正链匹配，另一部分又和负链匹配），sam文件中只占用一行记录；

chimeric alignment: “嵌合比对” 的形成是由于一条测序read比对到基因组上时分别比对到两个不同的区域，而这两个区域基本没有overlap。因此它在sam文件中需要占用多行记录显示。只有第一个记录被称作"representative”,其他的都是"supplementary”【Chimeric reads are also called split reads】

嵌合 reads 与嵌合/融合基因（重组由来源与功能不同的基因序列剪接而形成的杂合基因）并不完全一样，RNA-seq中的chimeric read或许可以说明有融合基因存在，但在基因组中一般作为结构变异的证据

【下面👇的图1是原始数据(其中Coor是坐标的简写；ref是参考序列)；图2是使用bwa、bowtie2、hisat2等软件比对的结果sam文件】 【其中r001/1与r001/2是paired end read； r003是嵌合read，有两个记录】

1

2

read alignment： 不管黑猫(linear alignment)、白猫(chimeric alignment)，抓住老鼠(能够完整表示一个read)就是好猫(read alignment)

multiple mapping: 假如有一个短序列，他在比对的时候看到哪哪都有他的影子，这种就是受重复区域影响比较大，所以read越长出现这种的可能性越低。一般指定首先匹配上的为最优匹配结果primary。

坐标系统 ：

  • 1-based coordinate system ：序列的第一个碱基设为数字1，如SAM,VCF,GFF,wiggle格式
  • 0-based coordinate system ：序列的第一个碱基设为数字0，如BAM, BCFv2, BED, PSL格式

SAM要处理好的问题：

  • 非常多的序列（read)，mapping到多个参考基因组（reference）上；
  • 同一条序列，分多段（segment）比对到参考基因组上；
  • 各种结构化信息表示，包括错配、删除、插入等比对信息；

具体的Sam格式：

花花之前有介绍，并且在Linux考试题中也有所提及

sam格式

sam格式的头文件可有可无，主要有@HD，说明符合标准的版本、对比序列的排列顺序；@SQ，参考序列说明；@RG，比对上的序列（read）说明；@PG，使用的程序说明；@CO，任意的说明信息

关于Samtools

一般测序文件都比较大，生成的单个sam文件从几百M到十几G不等，为了学习samtools的使用首先要有sam文件

##准备原始数据
ascp -v -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -T -l500m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra 
./raw_data
## 转换成fq格式
fastq-dump --gzip --split-3 *.sra
##下载参考基因组
mkdir -p genome/hg19  && cd genome/hg19
for i in $(seq 1 22) X Y M; #指定染色体号下载
do echo $i;
# 一般下载UCSC的数据
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M; # 指定染色体号拼成整体
do cat chr${i}.fa >> hg19.fasta
done
rm -fr chr*.fasta
##构建参考基因组比对索引
mkdir -p ~/reference/genome/hg19 && cd ~/reference/genome/hg19
nohup time bwa index -a bwtsw -p hg19 ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1 &
## 有了索引以后，才能用bwa进行比对
bwa mem -t 10 -M ~/reference/index/bwa/hg19/hg19 SRR3589956_1.fastq.gz SRR3589956_2.fastq.gz 1>SRR3589956.sam 2>SRR3589956.bwa.align.log
done

测试数据

1. View

作用：转换sam与bam，对bam进行排序（sort）和提取的操作

  • sam转bam，-S指输入文件格式（不加-S默认输入是bam），-b指定输出文件（默认输出sam）

    samtools view -Sb SRR3589956.sam >SRR3589956.bam

  • 【如果要bam转sam，-h设置输出sam时带上头注释信息】

    samtools view -h SRR3589956.bam > SRR3589956.sam

  • 过滤功能-F：后接flag数字，常用有4（表示序列没比对上）、8（配对的另一条序列，即mate序列没比对上）以及12（两条序列都没比对上）。加上-F就表示过滤掉这些情况

## 提取一条reads比对到参考序列上的序列结果
samtools view -bF4 SRR3589956.bam>SRR3589956.F4.bam
## 提取两条reads都比对到参考序列上的序列结果
samtools view -bF12 SRR3589956.bam>SRR3589956.F12.bam

如果是-f，就是提取指定flag的序列

## 提取没有比对到参考序列上的序列结果
samtools view -bf4 SRR3589956.bam>SRR3589956.f4.bam

转换后的bam文件大小是2.8G，可以看到，相比于sam文件，节省了四分之三的空间；

另外得到的SRR3589956.F4.bam是2.7G，说明本来bam文件中比对上的reads就很多了，才过滤掉了不到四十分之一的reads

2. Sort

作用：对bam进行排序，一些软件需要使用排序后的bam，拿到bam先按照比对位置的顺序排一下，百利无一害

samtools sort -@ 10 -o SRR3589956.sorted.bam SRR3589956.bam

sort后的文件SRR3589956.sorted.bam大小是1.5G

3. merge

作用：将两个及以上的sort过的bam文件融合成一个bam文件

这里就一个示例文件就不演示了

samtools merge xxx.merged.bam xxx.sorted.bam

4. index

作用：对bam文件构建索引，产生.bai文件，方便以后的快速处理 bam文件进行排序sort后，才能进行index，否则报错； 要显示比对结果时，比如用IGV导入bam，就需要有.bai的存在

samtools index SRR3589956.sorted.bam

插播一条

当有多个bam文件时，一般思路就是对每一个bam进行sort、index后，再merge成一个整体merged.bam，然后对merged.bam再进行sort、index，才算能用了，得到最终结果应该是是sorted.merge.bam

5. faidx

作用：对fasta文件建立提取索引，索引文件后缀是.fai。利用索引文件可以快速提取fasta文件中的某些序列

samtools faidx hg19.fasta

6. tview

作用：直观显示reads比对到基因组的情况，与IGV类似 需要先sort和index

samtools tview SRR3589956.sorted.bam hg19.fasta

samtools tview

7. flagstat

作用：统计bam文件的比对结果

samtools flagstat SRR3589956.sorted.bam > SRR3589956.sorted.flagstat.txt

samtools flagstat

8. depth

作用：统计每个碱基位点的测序深度 需要使用重定向定义输出文件；要是用构建过索引的bam

-r：（region）加染色体号； -q：要求测序碱基质量最低值； -Q：要求比对的质量最低值

samtools depth SRR3589956.sorted.bam >SRR3589956.depth

9. mpileup

作用：用于生成bcf文件，之后结合bcftools进行SNP与InDel的分析 安装samtools时，包含了bcftools

-f：输入有索引的参考基因组fasta； -g：输出到二进制的bcf格式【不使用-g，就不生成bcf格式，而是一个文本文件，统计了参考序列中每个碱基位点的比对情况；每一行代表参考序列中某一个碱基位点的比对结果】

samtools mpileup -f hg19.fa SRR3589956.sorted.bam > SRR3589956.mpileup.txt

结果文件大小 7.1G Aug 3 21:48 SRR3589956.mpileup.txt

结果包含6列：参考序列名、匹配位置、参考碱基、比对上的reads数、比对的情况、比对的碱基质量

在第5列比对具体情况中: . 表示与参考序列正链匹配； , 表示与参考序列负链匹配； ATCGN 表示在正链不匹配； atcgn 表示在负链不匹配； * 模糊碱基； ^ 匹配的碱基是一个read的开始，后面的ASCII码-33表示比对质量，再向后修饰的(.,ATCGNatcgn) 表示该read的第一个碱基； $ 表示一个read结束，修饰前面碱基; 正则表达式+[0-9][ATCGNatcgn] 表示在该位点后面插入的碱基； 正则表达式-[0-9][ATCGNatcgn] 表示该位点后面缺失的碱基

samtools mpileup

10. bam转fastq

作用：方便提取出一段比对到参考序列的reads进行分析

利用软件：http://www.hudsonalpha.org/gsl/information/software/bam2fastq

11. rmdup

作用：将测序数据中由于PCR duplicate得到的reads去掉，只保留比对质量最高的reads

samtools rmdup 默认只对PE数据进行处理

12. idxstats

作用：输出一个表格，包含“序列名、序列长度、比对上的reads数、没有比对上的reads数”；其中第四列指PE reads中的一条read能匹配到参考基因组的染色体A，另一条read不能匹配到A上

samtools idxstats

13. reheader

作用：替换bam文件的头文件

samtools reheader

samtools tools
Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related