samtools

在前面，我们向大家介绍了 sambamba — samtools 的高效平替工具，今天让我们为大家带来生信分析最重要的工具之一 samtools。

简介

samtools 是由李恒大神开发的一套用于读/写/编辑/索引/查看 SAM/ BAM/CRAM 格式的工具集。其功能包括：

• view 进行 SAM 和 BAM 之间的文件格式转换或者查看二进制文件；
• sort 用于对文件进行排序，也是我们用得最多的子命令。
• index 生成索引文件，必须是 sort 之后的 bam 文件。
• stats：生成关于比对数据的统计信息，如测序覆盖度、比对质量等；
• faidx 对 fasta 文件建立索引，生成的索引文件以 .fai 后缀结尾。
• tview 查看 reads 比对到基因组的情况，类似基因组浏览器的功能。
• markdup 标记重复序列，但并没有将其去除。

samtools 安装方式

下载地址：http://www.htslib.org/download/

方式一：编译后安装

tar jxvf samtools-1.19.2.tar.bz2
cd samtools-1.19.2    # and similarly for bcftools and htslib
./configure --prefix=/where/to/install # 如果报错，跟着报错信息走
make
make install

典型报错信息：

configure: error: curses development files not found

yum install ncurses-devel

configure: error: libbzip2 development files not found

yum install bzip2-devel

configure: error: liblzma development files not found

yum install xz-devel

安装完成后添加到环境变量：

export PATH=/where/to/install/bin:$PATH    # for sh or bash users

或者创建软链接：

ln -s /where/to/install/bin/samtools /usr/local/bin/samtools

方法二：使用 conda 安装

$ conda install -y samtools

使用方法

按照不同场景，samtools 的多个子命令可以按照以下类别划分：

1. 格式准备

比对完参考基因组之后得到的 SAM 格式文件一般需要经过以下三连操作转换为 BAM 格式用于下游分析（IGV 可视化、variant calling 等）：

1.1 view

SAM 文件转 BAM：

samtools view -bS -o out.bam in.sam

-b 输出 bam 文件；

-S 输入 SAM 文件；

查看 bam 文件头部信息：

samtools view -H out.bam

1.2 sort

对 bam 文件排序

samtools sort -@ 4 -o out.bam in.bam

值得注意的是：

1. samtools1.9 版本以后 sort 算法从 qsort 升级为了 radixsort，性能提升非常大。此外，相对于 samtools+zlib，samtools+libdeflate 编译后的性能速度也得到了大幅提高，但 sambamba 仍然是我们的最佳选择。参考：https://bioinformatics.stackexchange.com/questions/18538/samt...
2. 此外，排序方式可分为两种，默认按 coordinate 进行排序；也可使用 -t 选项按 tag 排序；

1.3 index

对排序后 bam 格式的比对文件构建索引。

samtools index sorted.bam

注意：

1. 必须对 bam 文件进行排序后，才能建立索引；
2. BAI 索引格式支持最长 512 Mbp（2^29 碱基）的单个染色体。如果输入文件可能包含比对到更远位置的 reads，建议使用 sambamba。

2. 数据统计

depth

统计每个碱基位点的测序深度

samtools depth -b test.bed in.bam > test.depth

输出结果：

• 第一列 序列名称
• 第二列 碱基位置
• 第三列 测序深度

flagstat

统计 bam 文件中 read 的比对情况

samtools flagstat in.bam

输出：

in total (QC-passed reads + QC-failed reads),332843724,0
secondary,0,0
supplementary,0,0
duplicates,73213258,0
mapped,328386001:98.66%,0:N/A
paired in sequencing,332843724,0
read1,166421862,0
read2,166421862,0
properly paired,317070676:95.26%,0:N/A
with itself and mate mapped,327763284,0
singletons,622717:0.19%,0:N/A
with mate mapped to a different chr,1720040,0
with mate mapped to a different chr (mapQ>=5),275739,0

结果一共有13 项，根据 QC 结果，每一项都分 PASS 和 FAIL，一般都是 PASS，因为 QC 都在比对前做过了，低质量的大都过滤了。

idxstats

统计 bam 索引文件里的比对信息

samtools idxstats in.bam

idxstats 统计一个表格，4列，分别为 ”序列名，序列长度，比对上的 reads 数量，未比对上的 reads 数量”

3. 文件处理

merge

用于合并多个已排序的比对文件，生成一个包含所有输入记录的单一排序输出文件，同时保持现有的排序顺序。在表观基因组学中用的较多。

samtools merge -o merge_out.bam test_sort.bam test2_sort.bam

4. variant calling

markdup

识别并标记那些在进行基因组坐标排序后被视为重复的比对记录（默认情况下不去除）

samtools markdup in.bam out.markdup.bam

-r 删除重复 read

mpileup

没错，samtools 也可以直接用来进行 calling SNP/InDel。与 GATK 中使用复杂的局部组装和机器学习算法不同，samtools 检测变异的方法其实就是对每个位点做碱基堆积，然后提取与参考基因组不同的位点。但相对于 GATK ，samtools 在速度上并没有优势，所以我们一般不推荐，这里也只简单介绍一下。

samtools mpileup [-EB] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]

简单用法：

samtools  mpileup -q 1 -C 50 -S -D -m 2 -F 0.002 -u -f ref.fa -b bam.list -l split1.txt > test.mpileup

• -B, --no-BAQ 不计算单碱基比对质量分值（BAQ）

  BAQ (Base Alignment Quality)

  BAQ 是读取碱基错配的 Phred-scaled 概率。其有助于识别由错配引起的假阳性 SNP。关于 BAQ 的计算参考：“Improving SNP discovery by base alignment quality”, Heng Li, Bioinformatics, Volume 27, Issue 8 https://doi.org/10.1093/bioinformatics/btr076

• -f：输入有索引文件的 fasta 参考序列
• -A：在检测变异中，不忽略异常的 reads 对
• -C：用于调节比对质量的系数，如果 reads 中含有过多的错配，不能设置为零；
• -D：输出每个样本的 reads 深度
• -l：BED 文件或者包含区域位点的位置列表文件。注意：位置文件包含两列，染色体和位置，从1开始计数。BED文件至少包含3列，染色体、起始和终止位置，开始端从0开始计数。
• -r：在指定区域产生 pileup，需已建立索引的 bam 文件，通常和 -l 参数一起使用
• -o/g/v：输出文件类型（标准格式文件或者 VCF、BCF 文件）
• -t：设置 FORMAT 和 INFO 的列表内容，以逗号分割
• -u：生成未压缩的 VCF 和 BCF 文件
• -I：跳过 INDEL 检测
• -m：候选 INDEL 的最小间隔的 reads
• -F：含有间隔 reads 的最小片段

扫码关注微信公众号【生信F3】获取文章完整内容，分享生物信息学最新知识。

本文由mdnice多平台发布