一文了解并掌握samtools

一、写在前面

随着高通量测序技术的发展，新的数据格式为数据交互性、紧凑存储和高效数据分析提供了必要条件(毕竟买硬盘要花实打实的银子)。而目前最常用的SAM文件和VCF文件是由制表符分隔的文本文件，使用自定义脚本可以轻松处理这些文件，但解析速度慢，存储效率低下(3GB clip reads约生成34GB的sam文件，而bam仅有3GB)。因此2009年发表的samtools可以将SAM转换为二进制的BAM来应对这一问题，而在2013年3月发布0.1.19版本后，多线程的支持也大大节省了计算时间。时至今日，bam文件依然是测序数据比结果储存的"硬通货"。

1.webp.jpg

如果你对下面的教程比较迷茫，那么你可以先行学习Linux教程：

十小时学会Linux

生信Linux及服务器使用技巧

如果你的计算机不足以支持下面流程的计算，可按需选用适合自己的计算资源：

共享(经济实惠)：有root权限的共享服务器,报我名字立减200￥

独享(省电省心)：生信分析不求人

实体(稳定高效)：为实验室准备一份生物信息学不动产

访问链接：https://biomamba.xiyoucloud.net/

二、工作流程

1、fastq to BAM/CRAN

测序

产生的未比对文件如fastq文件可以存储未比对数据至BAM中，这种方法其实最佳，因为其记录了原始数据和辅助标签信息。但是这个教程是产生sorted并且aligned的BAM数据。

1.1 alignment/mapping to a known reference

步骤一般为：

Map / align

Fix mate-pair issues

Mark duplicates (part 1: preparation)

Sort to positional order

Mark duplicates (part 2: marking)

Convert to final file format

1.1.1 Mapping:

作者以minimap2为例，但无论哪种方式或软件，生成SAM文件即可：

minimap2 -t 8 -a -x sr C.Elegans.fa SRR065390_1.fastq SRR065390_2.fastq -o CE.sam

输出的SAM文件中reads的顺序与原始fastq文件相同（可能由于多线程的原因有一定程度的波动）。建议用-o输出而不是重定向。

1.1.2 Fixing of mate-pair issues

矫正比对时引入的缺陷，保证SAM FLAG, RNEXT, PNEXT and TLEN fields一致，老版本可能需要-S来指定SAM文件，而新版版的软件无需：

samtools fixmate -O bam,level=1 CE.sam fixmate.bam

-O bam,level=1会压缩bam文件，你也可以通过level=0 (或 -u)来取消压缩输出的bam文件

1.1.3 marking duplicates(part1:preparation)

标记重复需要对数据进行一些分析，这些分析按照读名称排序，有些则按照基因组位置顺序进行。而不是使用markdup工具内部进行排序，我们将其分为两个阶段，这样它就可以插入到我们的流程中，而不需要任何额外的排序需求，只要在上面的命令中加一个-m即可。

samtools fixmate -O bam,level=1 -m CE.sam fixmate.bam

1.1.4 sort to positinal order

可以根据基因组的染色体位置来对BAM进行排序：

samtools sort -l 1 -@8 -o pos.srt.bam -T /tmp/example_prefix fixmate.bam

-@8用于开启并行计算，samtools sort中的-m用于提供更多的内存以加快计算速度(每个线程需要的内存)，-l 1 意味着压缩，与-O bam,level=1类似。

1.1.5 marking duplicates(part2:marking)

这时可以按以位置排序后的位置来标记重复了：

samtools markdup -O bam,level=1 pos.srt.bam markdup.bam

1.1.6 输出最终文件：

BAM:

samtools view -@8 markdup.bam -o final.bam

1.2 合并流程

上面的命令可以用管道符传递从而避免临时文件的生成从而节省磁盘空间

minimap2 -t 8 -a -x sr C.Elegans.fa SRR065390_[12].fastq  | /

notes：

samtools只会调用其需要的线程，所以选项参数指定的过大时未必会白白占用计算资源。

alignment是整个流程的限速步骤

建议在shell脚本中加上以下设置，这样任何一个错误都会导致整个脚本的终止

set -o pipefail

上面的步骤没有丢弃任何没有mapping的数据，因此可以通过BAM文件转换回fastq文件，只是顺序会与原来的有所不同：

samtools sort -n -@8 final.cram | /

1.2 De-nove assembly

略

2、WGS/WES Mapping to Variant Calls

DNA测序(如WGS与WES)通常包含Mapping、improvement、variant calling三个步骤。

(1)mapping

这里samtools的作者使用的mapping工具为bwa

bwa index  #构建索引

由于BWA和其它的比对工具在生成sam文件时会记录一些不正常的alignment，因此首先对数据进行清洗与整理：

samtools fixmate -O bam  

从fastq中的name顺序转换为 coordinate order:

samtools sort -O bam -o  -T  

(2)improvement

为了减少数据中插入带来的错误，重新align包含gap的alignment是有帮助的：

java -Xmx2g -jar GenomeAnalysisTK.jar -T RealignerTargetCreator /

BQSR可以减少人工引入的测序效应，第一步测定协变量，第二部根据协变量进行一定的质量分数补偿：

java -Xmx4g -jar GenomeAnalysisTK.jar -T BaseRecalibrator /

把多个文库合并为一个BAM，这样有利于同时标记PCR重复与光学重复

java -Xmx2g -jar MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT /

也可以如此合并：

samtools merge    

如果计算资源与时间充裕，可以把插入缺失也重新align一下：

java -Xmx2g -jar GenomeAnalysisTK.jar -T RealignerTargetCreator /

最后将你的BAM索引化：

samtools index 

(3)variant caling 变异检测

为了将BAM转换为基因组位置，需要将其”堆积“为包含所有位置信息的BCF文件。可以利用这些信息检测基因型与变异位点：

bcftools mpileup -Ou -f     | /

如果希望展示一些未被BCF检测过程中检测到的特定位点，可以把上面的步骤拆成两步

bcftools mpileup -Ob -o  -f    

构建vcf索引：

tabix -p vcf 

下一步将有助于你过滤变量时制备图片或统计：

bcftools stats -F  -s -  > 

可以如此过滤VCF文件：

bcftools filter -O z -o  -s LOWQUAL -i'%QUAL>10' 

@SQ记录的是染色体相关的位置，若没有这一列则在转BAM时：

samtools faidx ref.fa

若SAM中已包含@SQ，则：

samtools view -b aln.sam > aln.bam

参考：

https://academic.oup.com/gigascience/article/10/2/giab008/6137722?login=true

https://www.htslib.org/workflow