samtools这几年的版本迭代节奏虽然不像普通软件那么频繁但每次更新都会带来实打实的变化。做测序分析的同学应该都清楚只要你的流程里碰过比对文件samtools基本就是绕不开的那个工具——BAM文件的查看、排序、索引、统计、转换绝大多数环节都要跟它打交道。而且它不光是命令好用底层对SAM/BAM/CRAM格式的支持直接影响你后续变异检测、可视化、数据合并能不能顺利跑通。今天这篇笔记我就把samtools安装与使用方法从头到尾捋一遍从源码编译到conda安装从view/sort这类高频命令到mpileup/tview这些进阶操作再把我这些年踩过的坑一并整理出来。不管你是刚入门的学生还是要搭建生产环境pipeline的工程师照着这篇基本都能落地。1. 版本选择与安装前准备别拿到什么就装什么1.1 为什么版本选择很重要很多人习惯直接apt install samtools或者conda install -c bioconda samtools装完发现命令能用就完事了。但等到你要跑一些新参数、新格式或者发现某些命令的行为和文档对不上时才会意识到版本差异有多坑。举例来说samtools 1.17之后对CRAM格式的参考基因组处理方式做过调整旧版本生成的CRAM在1.20以上版本里解码时如果参考序列缺失报错提示反而更严格。再比如markdup命令在1.10之前叫rmdup参数和输出统计都不一样你拿老的脚本在新版本里跑大概率会报错或者得到不同的重复率统计结果。所以我给的建议是在你的项目里固定一个版本或者至少在README里写明用的是哪个版本。如果是在全新环境里部署直接装当前稳定版写这篇笔记时最新的1.20/1.21系列就好老版本除非有特殊兼容需求不然没必要刻意保留。1.2 安装前先检查环境依赖决定成败samtools虽然是一个C写的工具但并不是编译完就完事的。它依赖几个基础库缺一个就会在编译或者运行阶段出幺蛾子zlib处理gz压缩的BAM文件必需这个基本所有Linux发行版都有libbz2bzip2压缩格式支持部分BAM/CRAM文件会用到liblzma用于CRAM和部分索引文件的解压缺少时编译能过但运行可能报错libcurl从远程URL直接读取文件、访问参考序列时用到1.17以后对网络相关功能依赖更明显ncurses支持tview这种交互式文本界面如果不用tview可跳过但建议还是安上。在开始编译之前我习惯先跑一遍系统检查。Ubuntu/Debian环境下通常是sudo apt update sudo apt install -y gcc make libncurses-dev libbz2-dev liblzma-dev libcurl4-openssl-dev zlib1g-devCentOS/RHEL/Fedora类的环境对应的是sudo yum install -y gcc make ncurses-devel bzip2-devel xz-devel libcurl-devel zlib-devel这一步做完后面编译samtools基本不会遇到缺头文件的尴尬。2. 三条安装路线对比与实操步骤2.1 conda/bioconda方式省心省力适合日常开发和教学如果你已经在用conda管理生物信息软件环境那samtools装起来非常简单而且conda会自动处理依赖关系尤其是对老版本的回退支持很好。比如你想在某个项目里固定samtools 1.10直接指定版本号就行。先确认自己有conda或者mamba。推荐用mamba解析依赖速度比conda快很多遇到大环境时体感差别非常大。mamba create -n samtools_env -c bioconda -c conda-forge samtools1.21 conda activate samtools_env samtools --version如果不指定版本号默认会装到当前bioconda通道里的最新版本。多环境隔离的好处是你跑一个老流程可以用Python 3.6 samtools 1.9的环境跑新流程再另建一个环境两边互不干扰。这个优势在做多项目并行时真的能救命。需要注意一点bioconda和conda-forge的通道顺序要写对优先-c bioconda再-c conda-forge否则可能出现依赖解析错误或者拿到非官方构建版本的情况。2.2 源码编译安装生产环境更稳妥的选择源码编译看起来麻烦但在真正的高性能计算集群上往往比conda更可控。因为你可以通过--prefix把samtools装到自己的目录下不用借助任何包管理器不污染系统环境。而且自己编译能针对CPU指令集做一些优化虽然samtools这种IO密集工具在编译优化上提升有限但至少可控。先去samtools官网或者GitHub发布页下载源码包wget https://github.com/samtools/samtools/releases/download/1.21/samtools-1.21.tar.bz2 tar -xjf samtools-1.21.tar.bz2 cd samtools-1.21接下来是常规三步走./configure --prefix/opt/samtools/1.21 make -j 8 make install这里我强调几个细节--prefix指定安装路径后续通过绝对路径调用或者把这个路径加到PATH里make -j 8中的8是并行编译的任务数如果服务器负载高或者内存吃紧可以用-j 4避免编译过程中内存溢出编译完先别急着跑用make test跑一遍自带的功能测试确认所有用例都通过再正式使用。源码编译对依赖库的位置比较敏感如果你之前把某些库装到了非标准路径可能需要在configure前设置CPPFLAGS和LDFLAGS。比如CPPFLAGS-I/usr/local/include LDFLAGS-L/usr/local/lib ./configure --prefix/opt/samtools/1.21这种情况在集群环境里不少见因为管理员不一定把所有库都放到系统默认路径。2.3 预编译包与系统包管理器临时应急可以长期不推荐Ubuntu/Debian下sudo apt install samtools确实能装但版本往往偏旧。以Ubuntu 22.04为例默认仓库里大概率还是1.14或者1.15而官方已经出到1.21。不是说旧版本不能用而是很多新特性你也体验不到。macOS下用Homebrew装brew install samtools这个版本更新相对及时但要留意Homebrew维护者更新有延迟可能不是最新版。Windows用户则建议用WSL或者官方Windows子系统直接原生跑samtools会遇到一些路径兼容和编译困难的问题不太推荐折磨自己。预编译包适合什么场景比如你在对方机器上临时帮忙看一个BAM文件不打算长期使用那直接apt install然后走人就好。如果是要经常跑pipeline我建议还是花5分钟用conda或者源码装一个自己可控的版本。3. 理解SAM/BAM格式命令背后的数据模型3.1 SAM那11列每一列都代表什么samtools所有操作都是围绕SAM/BAM格式展开的不搞懂这个格式你连报错都看不懂。SAM全称是Sequence Alignment/Map格式是纯文本存储比对结果的格式BAM则是它经过BGZF压缩后的二进制形式体积小、读取随机访问快。两者内容一样只是存储方式不同。每一行比对记录有11个核心字段缺一不可字段名称含义1QNAME比对序列的名称也就是read ID2FLAG比对特征标志位用整数表示多个属性3RNAME参考序列名称比对到哪条染色体或contig4POS比对起始位置1-based坐标5MAPQ比对质量值越高质量越高6CIGAR比对操作的描述串如10M1I5M7RNEXT配对read中另一条比对到的参考序列名8PNEXT配对read中另一条比对到的位置9TLEN插入片段长度template length10SEQread的碱基序列11QUALread的碱基质量值ASCII编码的Phred分数view命令输出的SAM就是这11列理解每一列的含义后你基本上能看懂samtools任何输出。比如flagstat里面的“properly paired”统计实际就是根据FLAG和TLEN计算的sort按坐标排序实际上就是按RNAME和POS排序。我平时排查数据问题第一步就是head一个BAM文件看SAM格式的列结构往往能快速发现问题源头比如序列名混乱、FLAG异常、POS越界等。3.2 FLAG值、CIGAR串和MAPQ是读懂比对结果的三把钥匙这三个字段是SAM格式里最容易让人犯迷糊的部分。FLAG值不是简单的0/1而是一个按位编码的十进制整数。它把多个布尔属性打包在一起具体对照关系是1read成对存在paired2paired read中两条都能正确比对proper pair4read自身未比对unmapped8paired read中另一条未比对mate unmapped16read比对的链为负链reverse strand32paired read中另一条比对到负链mate reverse strand64这是read1first in pair128这是read2second in pair256次级比对secondary alignment512QC失败1024PCR或光学重复2048补充比对supplementary alignment比如FLAG99分解一下是643221代表read1、proper pair、另一条负链、成对read一看就知道是正常的双端比对。而FLAG147则是1281621代表read2、本身负链、proper pair、成对。你用samtools view -f 2筛选proper pair底层就是把FLAG值和2做按位与运算结果非零就保留。CIGAR串则是描述比对操作的字符串从左边开始读M匹配或错配比对到参考序列的碱基Iread上有插入参考序列上无对应Dread上有缺失参考序列上有碱基但read没有覆盖N跳过区域比如RNA-seq的intron区域S软剪切read两端未比对的碱基SEQ中保留H硬剪切read两端未比对且SEQ中不保留Ppadding用于表示gap一条CIGAR为76M1I23M的read表示先比对76个碱基中间插入1个碱基再继续比对23个碱基。你算read总长直接76123100就行。MAPQ字段则告诉我们对这个比对位置的信心有多高。不同比对软件输出范围不一样BWA MEM的范围通常到60STAR的MAPQ可能到255。一般过滤阈值取MAPQ20或者30就够用具体取决于你怎么定义“可靠比对”。4. 高频命令的实操要点与参数取舍4.1 viewBAM/SAM/CRAM互转与过滤view是samtools里出场频率最高的命令负责读文件、写文件、格式转换和过滤。最基础的用法# SAM转BAM samtools view -bS aln.sam aln.bam # BAM转SAM加head samtools view -h aln.bam aln.sam # 只输出head samtools view -H aln.bam-b表示输出BAM格式-S表示输入是SAM格式新版samtools即使省略-S也能自动识别但老脚本里会有尊重一下历史习惯就行。过滤是view的重头戏# 只保留比对质量不低于30的read samtools view -q 30 aln.bam # 只保留proper pair samtools view -f 2 aln.bam # 过滤掉unmapped不保留未比对的read samtools view -F 4 aln.bam # 过滤掉secondary和supplementary只留primary比对 samtools view -F 0x900 aln.bam-f和-F的区别要记牢-f是只要FLAG满足某个位就保留-F是只要FLAG满足某个位就过滤掉。很多人最开始搞反结果筛选结果完全不对。另外view还支持直接操作CRAM格式。CRAM相比BAM压缩率更高特别适合存储大样本数据# BAM转CRAM需要参考基因组 samtools view -C -T reference.fa aln.bam aln.cram # CRAM转BAM samtools view -b -T reference.fa aln.cram aln.bamCRAM必须要参考基因组才能解码没有参考序列文件CRAM基本等于一堆乱码。所以分发CRAM文件时一定要记得附带参考基因组的版本信息不然后续协作会非常麻烦。4.2 sort与index排序和索引不只是“跑一下”samtools sort用于按坐标排序BAM文件。排序是对比分析必不可少的一步因为后续很多操作比如变异检测、深度的计算、索引构建都要求输入文件是按坐标排好的。基础用法samtools sort -o aln.sorted.bam aln.bam如果文件特别大一定要记得调整内存和线程参数samtools sort -m 4G - 8 -o aln.sorted.bam aln.bam-m是每个线程使用的内存上限-是线程数。这两者配合起来就决定了排序过程的内存占用。比如-m 4G - 8理论峰值内存是32G如果服务器实际内存只有16G很容易被OOM kill掉。我一般建议按总内存的50%~70%去折算留出一部分给其他进程。索引是构建排序BAM的坐标索引文件samtools index aln.sorted.bam生成的文件默认是aln.sorted.bam.bai。当参考基因组特别长或者酶切位点跨越不同染色体间的大片段时普通BAI索引在超过512M区域时会有问题可以用CSI索引samtools index -c aln.sorted.bamCSI索引没有512M的限制在人类全基因组、泛基因组等大型参考序列场景下是更好的选择。有个非常容易被忽略的细节sort之后的BAM如果接下来想用markdup标记重复建议在排序前先用fixmate补齐配对信息否则markdup的结果会大打折扣。这个我们在后面进阶部分详细说。4.3 flagstat、idxstats、depth用一套命令快速掌握数据全貌拿到一个陌生BAM文件我习惯先跑三个统计命令samtools flagstat aln.sorted.bam samtools idxstats aln.sorted.bam samtools depth -a aln.sorted.bam | headflagstat输出的是比对整体统计包括总条数、proper pair数、unmapped数、secondary数等。读到这个输出你能快速判断比对质量好不好。如果proper pair占比明显低于90%可能是参考基因组不对、测序数据有污染、或者比对软件参数不合适。idxstats则按染色体/contig分别统计每个区域的比对数量和碱基总数。这在判断是否有明显偏倚时很有帮助比如某条染色体比对的read数异常高可能是参考序列包含污染序列、或者样本本身拷贝数异常。depth计算每个位点的测序深度是评估覆盖均匀性的基础。-a参数表示输出所有位点包括深度为0的位置不加-a的话会跳过0深度位点统计得到的结果会偏乐观。我更推荐用depth配合awk做简单汇总samtools depth -a aln.sorted.bam | awk {sum$3; if($30) cnt} END {print mean depth:, sum/NR; print covered bases:, cnt}就能快速估算全基因组的平均深度和覆盖比例。当然做正式报告时可以用mosdepth这类专攻深度的工具但日常摸鱼式抽查samtools一个命令就够了。5. 进阶场景变异检测、去重复与可视化5.1 mpileup配合bcftools做变异检测samtools mpileup生成每个位点的比对信息汇总再交给bcftools call来调用变异。这一套是很多老流程的标配虽然现在GATK HaplotypeCaller在大型项目里更常用但mpileupbcfcall的轻量性和速度优势仍然很明显尤其是做细菌、病毒这类小型基因组时几秒钟就能出来结果。基本流程# 建参考基因组索引 samtools faidx reference.fa # 生成pileup并调用变异 samtools mpileup -uf reference.fa aln.sorted.bam | bcftools call -mv -o variants.vcf这里的关键参数是-u产生未压缩的BCF输出这是给bcftools用的格式不是文本VCF-f指定参考基因组FASTA-g产生包含基因型信息的BCF在部分版本中与-u联用很常见等价于--BCF-A保留所有位点的等位基因信息用于计数等场景-Q过滤比对质量低于指定阈值的read-q过滤碱基质量低于指定阈值的read有些版本是-q注意看帮助。实际使用时我会先做一轮QC过滤掉未正确配对、质量分低的read避免它们对变异位点判断造成干扰samtools view -b -f 2 -q 30 aln.bam aln.filtered.bam samtools sort -o aln.filtered.sorted.bam aln.filtered.bam samtools mpileup -uf reference.fa aln.filtered.sorted.bam | bcftools call -mv -o variants.vcf另外提一句mpileup的输出数据量很大尤其是人类全基因组级别的文件结果直接用管道传给bcftools call处理即可不要先转成中间文件省内存也省磁盘。5.2 markdup标记PCR重复记得先fixmate标记PCR重复是变异检测前很重要的一步因为PCR扩增过程中同一个DNA片段可能被复制多份测序时会产出多条完全相同的read如果不标记掉变异频率会被过高估计。samtools从1.10开始推荐工作流是# 第一步按坐标排序必须 samtools sort -n -o aln.namesorted.bam aln.bam # 第二步fixmate补齐配对信息 samtools fixmate -m aln.namesorted.bam aln.fixmate.bam # 第三步再按坐标排序 samtools sort -o aln.fixmate.sorted.bam aln.fixmate.bam # 第四步标记重复 samtools markdup -r aln.fixmate.sorted.bam aln.markdup.bam第一步用-n是按read名称排序因为fixmate需要在同一对read名称相邻的情况下才能正确工作如果先按坐标排序配对的两条read可能隔得很远fixmate就无法整合mate信息。第三步重新按坐标排序是为了满足后续变异检测输入的要求。第四步的-r表示直接移除重复read如果你只想标记而不删除去掉-r即可。这个流程的排序次数看似多但每一步都有明确目的。实际操作中如果你的输入数据比对质量很高而且确定不需要fixmate的额外信息比如只统计比对、不做变异检测可以跳过第二步。但正规的变异检测流程我还是建议完整走一遍。另外注意markdup的输入文件里不要包含未配对的read或者至少要做好过滤否则会有很多无法妥善处理的单端read进入重复标记流程输出统计会偏差。5.3 tview交互式查看比对结果samtools tview可以在终端里直接展示比对情况这在快速检查单个位点的比对模式时很有用不用打开IGV。用法是samtools tview aln.sorted.bam --reference reference.fa进入交互界面后可以通过方向键移动光标按?查看操作说明按q退出。支持颜色高亮、CIGAR串显示、回文序列标记等。虽然无法和IGV这种图形化工具媲美但胜在轻快便利服务器上没有图形界面时也能快速排查问题。我能想到最实用的tview场景是当变异检测软件在某个位点报出可疑变异时你在屏幕前跑一下tview光标移动到该位点看看到底是同源染色体间的杂合信号还是比对假象。这种即时反馈比反复导数据到本地再加载IGV高效太多。除了tviewsamtools还提供了view -H、reheader这类操作头部区域的命令以及merge合并、bedcov按区间计算覆盖度等实用功能。硬要追求全面功能可以去翻官方文档但上面这些已经是日常pipeline里约90%会用到的高频命令。6. 常见报错与排查经验附速查表6.1 文件损坏、内存爆掉、进度卡死这类问题怎么定位我见过太多人在samtools上栽跟头其实问题原因就那几类提前知道能省不少力气。BAM文件报错“EOF marker is absent”这个提示出现在读取BAM文件尾部时说明文件不完整很可能是中断的下载或者比较早异常退出的流程产生的。遇到这个情况回到上游重新生成文件即可想硬修是不太现实的。报错“Malformed BAM header”通常意味着文件开头信息丢失或者文件根本不是BAM格式比如把txt文件改成了.bam后缀。用file aln.bam看下文件真实类型可避免很多低级错误。报错“Failed to open file ... No such file or directory”看起来像是文件不存在但在集群上更常见的是路径权限问题检查一下目录和文件的读权限即可。另一种可能是指定参考基因组路径时用了相对路径而程序的工作目录和当前目录不一致。运行过程中被“Killed”这是最让人头疼的问题通常意味着内存不够被系统OOM kill。解决方法在sort那一节已经提过调低-m和-给系统留足余量。“Failed to parse reference sequence”多半是参考基因组的FASTA文件有问题比如缺少索引需要faidx生成、FASTA序列行中有非法字符。先用samtools faidx reference.fa生成索引再检查一下FASTA首行是否有开头的有效头。我整理了一个速查表方便你快速定位问题现象可能原因快速排查命令EOF marker is absent文件不完整ls -l 文件大小对比预期大小Malformed BAM后缀名错误或文件损坏file 文件无法打开文件路径或权限ls -l、id被Killed内存不足free -g降低-m/-参考基因组报错缺少faidx或FASTA非法samtools faidx ref.fa版本相关参数失效版本差异打印samtools版本对比文档6.2 性能调优经验线程、内存、压缩级别的权衡samtools很多命令是CPU和IO混合型任务。盲目加线程不一定能加速反而可能因为IO带宽成为瓶颈而更慢。以sort为例设置- 8 -m 4G通常能发挥较高性能但如果你把-调到32甚至64可能效果有限还会把内存吃满。大BAM文件转换时可以尝试换用CRAM格式并将压缩级别调到9samtools view -C -T reference.fa -O CRAM,level9 -o aln.cram aln.bam在保证解码正常的前提下压缩率提升明显特别是人类全基因组数据CRAM体积大约只有BAM的50%~70%。并行同时跑多个samtools命令时也要注意给其他任务留资源。我所在的集群通常一个节点有48核同时跑8个samtools sort任务每个任务占用4~6个线程内存控制在2~4G整个调度就非常平滑不会互相抢资源。还有个小技巧对临时文件、中间文件使用本地高速磁盘如NVMe存放能够显著减少IO等待时间。如果你在共享文件系统上跑大型任务可以试试把TMPDIR指向本地磁盘export TMPDIR/local_disk/tmp samtools sort -o out.bam in.bam这个做法在文件特别大时效果明显排序速度能提升不少。6.3 几个容易被忽略的细节最后再提醒几个细节。第一samtools index只能在文件按坐标排序后才能操作否则生成索引必然报错。第二samtools merge合并多个BAM时要求所有输入文件事先按相同方式排序如果混用sort和sort -n的结果输出会错乱。第三CRAM文件不携带参考序列使用前必须提供和生成时一致的reference.fa。我个人的习惯是每次跑完一步关键操作马上跑一个验证命令确认结果不空、统计正常。比如flagstat检查总数变化、quickcheck检查文件完整性samtools quickcheck aln.sorted.bam echo OK这个quickcheck跑得很快适合放在pipeline的关键节点做自动校验能帮你更早发现文件是否损坏、流程是否有bug。说实话samtools这套命令本身不复杂但它和众多上游下游工具的组合使用才是真正的价值所在。掌握好安装方法和常用命令后续学习bcftools、mosdepth、igv等工具都会顺畅很多。如果你在集群上还需要处理几十甚至上百个样本的批量操作建议再封装一个shell脚本或者用snakemake/nextflow这类流程管理工具把samtools这些命令串起来才算是真正把效率提上来了。