samtools 是一个用于操作sam和bam文件的工具合集。能够实现二进制查看、格式转换、排序及合并等功能，结合sam格式中的flag、tag等信息，可以完成比对结果的统计汇总。同时利用linux中的grep、awk等操作指令，还可以大大扩展samtools的使用范围与功能。
samtools都可以通过哪些方式实现如上👆功能？
如果系统中已经安装过samtools，直接在命令行中敲击samtools查看对应文档
$ samtools
Program: samtools (Tools for alignments in the SAM format)
Version: 1.9 (using htslib 1.9)
Usage: samtools <command> [options]
Commands:
-- Indexing
dict create a sequence dictionary file
faidx index/extract FASTA
fqidx index/extract FASTQ
index index alignment
-- Editing
calmd recalculate MD/NM tags and '=' bases
fixmate fix mate information
reheader replace BAM header
targetcut cut fosmid regions (for fosmid pool only)
addreplacerg adds or replaces RG tags
markdup mark duplicates
-- File operations
collate shuffle and group alignments by name
cat concatenate BAMs
merge merge sorted alignments
mpileup multi-way pileup
sort sort alignment file
split splits a file by read group
quickcheck quickly check if SAM/BAM/CRAM file appears intact
fastq converts a BAM to a FASTQ
fasta converts a BAM to a FASTA
-- Statistics
bedcov read depth per BED region
depth compute the depth
flagstat simple stats
idxstats BAM index stats
phase phase heterozygotes
stats generate stats (former bamcheck)
-- Viewing
flags explain BAM flags
tview text alignment viewer
view SAM<->BAM<->CRAM conversion
depad convert padded BAM to unpadded BAM
大致可以看出来，分为5类命令块，Indexing,Editing,File Operations,Statistics,Viewing
这里就介绍一些常用的命令，应用时如果忘记参数的使用方式，还是建议命令行中敲击samtools + 命令块来查询 eg： samtools view
1. samtools view (文件转换、数据提取)
1.1 功能总结
view是经常用到的一个命令块，总结它的功能如下：
- bam 文件 和 sam 文件 互换
- 数据的提取
- 将排序或提取到的数据输出为bam或sam格式
Usage: samtools view [options]
主要认识以下几个参数的含义
-b #output BAM (输出bam格式，默认下输出SAM格式文件)
-h # print header for the SAM output (默认下输出的sam格式文件不带header,该参数设定输出sam文件时带header信息)
-H #print SAM header only (no alignments 仅仅输出文件的头文件)
-S #ignored (input format is auto-detected) (默认下输入是BAM文件，若是输入是SAM文件，最好加上这个参数，否则会报错)
-u # uncompressed BAM output (force -b该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间)
-c # print only the count of matching records(仅输出匹配的统计记录)
-L FILE # only include reads overlapping this bed FILE(仅包括和bed文件存在overlap的reads)
-o FILE # output file name[stdout](输出文件的名称)
-F INT # only include reads with none of the FLAGS in INT present(过滤包含flag=INT的reads,仅输出指定FLAG值的序列)
-f INT # only include reads with all of the FLAGs in INT present
-q INT # only include reads with mapping quality >= INT(比对的最低质量值，一般认为20就为unique比对了，可以结合上述-bF参数使用提取特定的比对结果)
-@ Number of additional threads to use #使用的线程数
1.2 实例认识 samtools view
Tips:
其中出现的“提取数字” 4、12是 展示template mapping情况的FLAG，不懂的话可以去详细了解一下bam文件格式 bam&sam文件浅析
1.将sam文件转换成bam文件
samtools view -bS abc.sam > abc.bam
2.BAM转换为SAM
samtools view -h -o out.sam out.bam（默认输入为bam，输出为sam）
3.提取比对到参考序列上的比对结果
samtools view -bF 4 abc.bam > abc.F.bam
4.提取paired reads中两条reads都比对到参考序列上的比对结果，只要把两个4+8的值12作为过滤参数即可
samtools view -bF 12 abc.bam > abc.F12.bam
5.提取没有比对到参考序列上的比对结果
Samtools view -bf 4 abc.bam > abc.f.bam
6.提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式
Samtools view abc.bam scaffold1 > scaffold1.sam
7.提取scaffold1上能比对到30k到100k区域的比对结果
Samtools view abc.bam scaffold1:30000-100000 $gt;scaffold1_30k-100k.sam
8.根据fasta文件，将header加入到sam 或 bam 文件中
Samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam
2.sort（数据排序）
2.1 功能总结
sort一般用来对bam文件进行排序
Usage: samtools sort [option] <in.bam> -o <out.prefix>
-n # sort by read name(设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序即header和序列从左往右的位点顺序)
-m INT # set maximum memory per thread; suffix K/M/G recognized[768M默认]#设置每个线程的最大内存，单位可以是K/M/G，默认是768M。对于处理大数据时，如果内存够用，则设置大点的值，以节约时间
-t TAG # sort by value of TAG. Uses position as secondary index(or read name if -n is set 按照TAG值排序)
-o FILE # write final output to FILE rather than standard output(输出到XX文件中)
2.2 实例认识 samtools sort
实例：tem.bam按照序列位置排序，并将结果输出到tem.sort.bam
Samtools sort -n tmp.bam -o tmp.sort.bam
Samtools view tmp.sort.bam
3.merge 和 cat （数据合并）
3.1 merge
merge将多个已经sort了的bam文件融合成1个bam文件，融合后的文件不需要再sort。
Usage: samtools merge [-nurlf][-h inh.sam][-b <bamlist.fofn>]<out.bam><in1.bam>[<in2.bam>...<inN.bam>]
Options:
-n #input files are sorted by read name
-t TAG #Input files are sorted by TAG value
-r #Attach RG tag (inferred from file names添加RG标签)
-u #uncompressed bam output 输出未压缩的bam
-f #overwrite the output bam if exist(覆盖已经存在的bam)
-1 #compress level 1 1倍压缩
-l INT #compression level, from 0 to 9 指定压缩倍数
-R STR #merge file in the specified region STR [all] 方括号内代表默认
-h FILE #copy the header in FILE to <out.bam>
3.2 cat
cat命令不需要将bam文件进行sort。
$samtools cat
Usage: samtools cat [options] <in1.bam> […<inN.bam>]
samtools cat [options] <in1.bam> […<inN.cram>]
Options:
-b FILE list of input BAM/CRAM file names, one per line
-h FILE copy the header from FILE [default is 1st input file]
-o FILE output BAM/CRAM
4. index 和 faidx (建立索引)
对排序后的序列建立索引，并输出bai文件，用于快速随机处理。在很多情况下，特别是需要显示比对序列的时候（samtools tview 或 IGV），bai文件必不可少。
-
samtools index为bam文件建立索引
Usage: samtools index <in.bam>[out.index]
Samtools index abc.sort.bam
-
samtools faidx为fasta文件建立索引
Usage: samtools faidx <file.falfile.fa.gz>[<reg>[…]]samtools faidx<file.falfile.fa.gz>[<reg>[...]]
5.tview（reads比对可视化）
tview能直观的显示出reads比对基因组的情况，和基因组浏览器有点类似。
Usage:samtools tview <aln.bam>[ref.fasta] #samtools tview bam文件 参考序列文件
当给出参考基因组的时候，会在第一排显示参考基因组的序列，否则，第一排全用N表示。
颜色标注的是碱基质量：30-40 白色；20-30 黄色；10-20 绿色；0-10 蓝色。
使用方法：
- 到达基因组的某一个位点：按下g，然后输入位置eg: 10:1000
- 移动显示界面：键盘上的 H(左)J(上)K(下)L(右)按键就可以
- 切换碱基和点号：按键“.”
- 切换显示read name：按键“r”
- 其他使用说明可以按”?”键来查看
6.flagstat （统计bam文件比对结果）
给出bam文件的比对结果，并输出比对统计结果。
Usage: samtools flagstat
实例：
屏幕快照 2019-11-13 下午9.51.04.png
Final.bam总共有153250612条reads
0条reads secondary
3227750条reads supplementary
44754021条reads duplicates
153250612 条reads mapped到参考基因组，总体上reads匹配率100%
150022862 条reads 是属于paired reads
Reads1 中 reads数 75011431
Reads2 中 reads数 75011431
148020994条reads完美匹配：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
150018176条reads 它和对应的paired reads都比对到参考序列上
2257条reads 单独一条匹配到参考序列上
Paired reads中两条分别比对到两条不同的参考序列有 1758676 条
Paired reads中两条分别比对到两条不同的参考序列有 1331738 条 比对质量大于等于5
7.depth （计算测序深度）
得到每个碱基位点的测序深度，并输出到标准输出。（前提：输入的bam文件必须先做samtools index）
Usage: samtools depth [-r reg][-q baseQthres][-Q mapQthres][-b in.bed]<in1.bam>[…]
-r <chr:from-to> region #后面跟染色体号和位置
-a output all positions (including zero depth) #输入所有位置的序列，包括深度为0的
-q <int> #base quality threshold[0]碱基质量阈值
-Q <int> #mapping quality threshold 比对的质量阈值
实例：
samtools depth tmp.index.bam > tmp.depth.bam（只是将深度信息保存到tmp.depth.bam文件，如果不定向的话会直接打印到界面）
8. mpileup （生成位点匹配信息）
Usage: samtools mpileup [-EBug][-C capQcoef][-r reg][-f in.fa][-l list][-M capMapQ][-Q minBaseQ][-q minMapQ] in.bam [in2.bam [...]]
-f FILE #来输入有索引文件的fasta参考序列
-C INT #adjust mapping quality; recommended:50, disable:0 [默认0]
-r chr$:pos1-pos2 #region in which pileup is generated 输出chr$:pos1-pos2区域的pileup结果
-q, --min-MQ INT skip alignments with mapQ smaller than INT [0]
-Q, --min-BQ INT skip bases with baseQ/BAQ smaller than INT [13]
mpileup生成的结果包含6列：参考序列名称；位置；参考碱基；比对上的reads数；比对情况；比对上的碱基的质量
其中第5列需要理解一下pileup的展现形式
.代表与参考序列正链匹配
,代表与参考序列负链匹配
ATCGN代表在正链上的不匹配
atcgn代表在负链上的不匹配
*碱基缺失
$代表一个read的结束，该符号修饰的是它前面的碱基
8. rmdup（去除dup）
在进行了reads比对后需要将由PCR duplicates获得的reads去掉，并只保留最高比对质量的read。使用rmdup命令完成
Usage: samtools rmdup [-sS]
-s # rmdup for SE reads 默认情况下，只针对single-end reads
-S # treat PE reads as SE in rmdup (force -s) 将paired-end reads作为single-end reads 处理
9.idxstats（统计表格）
Usage: samtools idxstats <aln.bam>
表格内容分4列：序列名，序列长度，比对上的reads数，unmapped reads number（应该是paired reads中有一端能匹配到该scaffold上，而另外一端不匹配到任何scaffolds上的reads数）
samtools idxstats
网友评论