samtools使用详解 - ricket-sjtu/bioinformatics GitHub Wiki
samtools是由Heng Li开发的针对序列比对结果标准格式sam及其二进制格式bam的分析处理工具包：
-
samtools view：将sam与bam之间进行相互转换； -
samtools sort：对bam文件进行排序，可以按照染色体坐标，也可按照基因名等序列标识进行排序； -
samtools merge：合并多个排序后的bam文件，结果也是排序好的； -
samtools index：为排序后的bam文件建立索引，索引文件的后缀为.bai； -
samtools faidx：为FASTA文件建立索引，索引文件的后缀为.fai； -
samtools tview：将比对结果reads与基因组显示在同一张图中，类似于基因组浏览器； -
samtools flagstat：给出bam的比对结果； -
samtools depth：输出每个位点的测序深度，并输出到stdout； -
samtools mpileup：输出突变位点的信息，保存为vcf或bcf
$ samtools flags
About: Convert between textual and numeric flag representation
Usage: samtools flags INT|STR[,...]
Flags:
0x1 PAIRED .. paired-end (or multiple-segment) sequencing technology
表示采用了PE测序技术，否则为SE技术
0x2 PROPER_PAIR .. each segment properly aligned according to the aligner
根据比对结果，表明序列片段有正常的pairing
0x4 UNMAP .. segment unmapped
比对结果说明片段没有比对到参考序列上
0x8 MUNMAP .. next segment in the template unmapped
模板对应的Pair片段没有比对到参考序列上
0x10 REVERSE .. SEQ is reverse complemented
序列与目标参考序列反向互补，说明是另一条链
0x20 MREVERSE .. SEQ of the next segment in the template is reversed
配对序列是反向互补
0x40 READ1 .. the first segment in the template
是模板中的第一个片段
0x80 READ2 .. the last segment in the template
是模板中的最后一个片段
0x100 SECONDARY .. secondary alignment
0x200 QCFAIL .. not passing quality controls
比对没有达到mapping质量的要求
0x400 DUP .. PCR or optical duplicate
是PCR复制的结果
0x800 SUPPLEMENTARY .. supplementary alignment
嵌合（chimeric）比对
Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]
默认情况下不加 region，则是输出所有的 region.
Options: -b output BAM
默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式
-S input is SAM
默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
-h print header for the SAM output
默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息
-H print header only (no alignments)
只输出header部分，不输出比对结果
-u uncompressed BAM output (force -b)
该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。
-c Instead of printing the alignments, only count them and print the
total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ ,
are taken into account.
只输出比对的数量。如果用了过滤选项，只计算符合条件的比对结果
-1 fast compression (force -b)
快速压缩
-x output FLAG in HEX (samtools-C specific)
输出的FLAG用十六进制表示
-X output FLAG in string (samtools-C specific)
输出的FLAG用字符串表示
-c print only the count of matching records
只输出匹配的记录数
-L FILE output alignments overlapping the input BED FILE [null]
输出与输入BED文件具有overlapping的比对结果
-t FILE list of reference names and lengths (force -S) [null]
使用一个list文件来作为header的输入
-T FILE reference sequence file (force -S) [null]
使用序列fasta文件作为header的输入
-o FILE output file name [stdout]
输出文件名
-R FILE list of read groups to be outputted [null]
指定需要输出的RG
-f INT required flag, 0 for unset [0]
只输出符合FLAG的比对结果
-F INT filtering flag, 0 for unset [0]
过滤掉符合FLAG的结果
Skip alignments with bits present in INT [0]
数字4代表该序列没有比对到参考序列上
数字8代表该序列的mate序列没有比对到参考序列上
-q INT minimum mapping quality [0]
最小比对质量分数
-l STR only output reads in library STR [null]
只输出某个库中的reads
-r STR only output reads in read group STR [null]
只输出某个RG中的reads
-s FLOAT fraction of templates to subsample; integer part as seed [-1]
用于采样的模板的比例，整数部分作为seed，如1234.60表示的是采样比例是0.60，而seed=1234
-? longer help
1. 将sam文件转换成bam文件
$ samtools view -bS abc.sam > abc.bam
$ samtools view -b -S abc.sam -o abc.bam
2. 提取比对到参考序列上的比对结果
$ samtools view -bF 4 abc.bam > abc.F.bam
3. 提取paired-end reads中两条reads都比对到reference的比对结果，只需要把两个4+8的值12作为过滤参数即可
$ samtools view -bF 12 abc.bam > abc.F12.bam
4. 提取没有比对到reference的比对结果
$ samtools view -bf 4 abc.bam > abc.f.bam
5. 提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式
$ samtools view abc.bam scaffold1 > scaffold1.sam
6. 提取scaffold1上能比对到30k到100k区域的比对结果
$ samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam
7. 根据fasta文件模板，将 header 加入到 sam 或 bam 文件中
$ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam
Usage: samtools sort [options...] [in.bam]
Options:
-l INT Set compression level, from 0 (uncompressed) to 9 (best)
设置压缩比，从0（不压缩）到9（最大压缩）
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]
设置每个线程能使用的最大内存（单位为K/M/G）
-n Sort by read name
按照read名进行排序，默认按照参考序列的名字顺序和坐标排序
-o FILE Write final output to FILE rather than standard output
将结果写到文件FILE中
-T PREFIX Write temporary files to PREFIX.nnnn.bam
指定生成临时文件的前缀
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
指定每个输入文件的格式
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
指定输出格式（SAM/BAM/CRAM）
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
指定每个输出文件的格式
--reference FILE
Reference sequence FASTA FILE [null]
参考序列FASTA文件
-@, --threads INT
Number of additional threads to use [0]
指定线程数
$ samtools sort -O BAM -o abc.sort.bam abc.bam
Usage: samtools merge [-nurlf] [-h inh.sam] [-b <bamlist.fofn>] <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]
将in1.bam等合并到out.bam文件中。
Options:
-n Input files are sorted by read name
输入文件按照read name排序
-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
最小压缩
-l INT Compression level, from 0 to 9 [-1]
设置压缩水平
-R STR Merge file in the specified region STR [all]
合并在特定region的文件
-h FILE Copy the header in FILE to <out.bam> [in1.bam]
将FILE中的header复制到out.bam中
-c Combine @RG headers with colliding IDs [alter IDs to be distinct]
合并@RG header
-p Combine @PG headers with colliding IDs [alter IDs to be distinct]
合并@PG header
-s VALUE Override random seed
重载随机seed
-b FILE List of input BAM filenames, one per line [null]
输入BAM文件的列表，每行一个
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
指定输入文件类型
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
指定输出文件类型
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
指定单个输出文件类型
--reference FILE
Reference sequence FASTA FILE [null]
参考序列FASTA文件
-@, --threads INT
Number of additional threads to use [0]
设置线程数
Usage: samtools index [-bc] [-m INT] <in.bam> [out.index]
Options:
-b Generate BAI-format index for BAM files [default]
为BAM文件创建bai索引
-c Generate CSI-format index for BAM files
为BAM文件创建csi索引
-m INT Set minimum interval size for CSI indices to 2^INT [14]
设置CSI索引的最小区间大小
-@ INT Sets the number of threads [none]
设置线程数
Usage: samtools faidx <file.fa|file.fa.gz> [<reg> [...]]
Usage: samtools flagstat [options] <in.bam>
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
指定输入文件格式
-@, --threads INT
Number of additional threads to use [0]
设置线程数
$ samtools flagstat example.bam
11945742 + 0 in total (QC-passed reads + QC-failed reads)
#总共的reads数
0 + 0 duplicates
7536364 + 0 mapped (63.09%:-nan%)
#总体上reads的匹配率
11945742 + 0 paired in sequencing
#有多少reads是属于paired reads
5972871 + 0 read1
#reads1中的reads数
5972871 + 0 read2
#reads2中的reads数
6412042 + 0 properly paired (53.68%:-nan%)
#完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
6899708 + 0 with itself and mate mapped
#paired reads中两条都比对到参考序列上的reads数
636656 + 0 singletons (5.33%:-nan%)
#单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。
469868 + 0 with mate mapped to a different chr
#paired reads中两条分别比对到两条不同的参考序列的reads数
243047 + 0 with mate mapped to a different chr (mapQ>=5)
#同上一个，只是其中比对质量>=5的reads的数量
samtools depth
Usage: samtools depth [options] in1.bam [in2.bam [...]]
Options:
-a output all positions (including zero depth)
输出包括测序深度为0的所有位点
-a -a (or -aa) output absolutely all positions, including unused ref. sequences
输出包含一直未用的参考序列
-b <bed> list of positions or regions
包含位点或区域的列表文件
-f <list> list of input BAM filenames, one per line [null]
输入BAM文件列表，每行一个
-l <int> read length threshold (ignore reads shorter than <int>) [0]
读长的阈值（忽略短于该长度的read）
-d/-m <int> maximum coverage depth [8000]
最大深度
-q <int> base quality threshold [0]
设置纳入计算的最小测序质量分值阈值
-Q <int> mapping quality threshold [0]
比对质量分值的阈值
-r <chr:from-to> region
设置需要计算的区域
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
输入文件格式
--reference FILE
Reference sequence FASTA FILE [null]
参考序列FASTA文件
The output is a simple tab-separated table with three columns: reference name,
position, and coverage depth. Note that positions with zero coverage may be
omitted by default; see the -a option
samtools mpileup
Usage: samtools mpileup [options] in1.bam [in2.bam [...]]
Input options (输入选项):
-6, --illumina1.3+ quality is in the Illumina-1.3+ encoding
质量分值采用Illumina-1.3+编码
-A, --count-orphans do not discard anomalous read pairs
保留没有正常匹配的read对
-b, --bam-list FILE list of input BAM filenames, one per line
输入BAM文件的列表，每行一个
-B, --no-BAQ disable BAQ (per-Base Alignment Quality)
不计算单碱基比对质量分值
-C, --adjust-MQ INT adjust mapping quality; recommended:50, disable:0 [0]
校正比对质量
-d, --max-depth INT max per-file depth; avoids excessive memory usage [250]
每个文件的最大深度
-E, --redo-BAQ recalculate BAQ on the fly, ignore existing BQs
重新计算BAQ
-f, --fasta-ref FILE faidx indexed reference sequence file
用samtools faidx建立过索引的参考序列文件
-G, --exclude-RG FILE exclude read groups listed in FILE
排除所有包含在FILE中的RG
-l, --positions FILE skip unlisted positions (chr pos) or regions (BED)
排除没有包含在FILE中的位置（chr pos）和区域（BED文件）
-q, --min-MQ INT skip alignments with mapQ smaller than INT [0]
忽略mapQ小于INT的比对结果
-Q, --min-BQ INT skip bases with baseQ/BAQ smaller than INT [13]
忽略BAQ小于该阈值的碱基
-r, --region REG region in which pileup is generated
只分析该区域
-R, --ignore-RG ignore RG tags (one BAM = one sample)
忽略RG标记
--rf, --incl-flags STR|INT required flags: skip reads with mask bits unset []
只包含FLAG已经标记的结果
--ff, --excl-flags STR|INT filter flags: skip reads with mask bits set
[UNMAP,SECONDARY,QCFAIL,DUP]
忽略FLAG被标记的结果
-x, --ignore-overlaps disable read-pair overlap detection
不检查read pair的overlap
Output options（输出选型）:
-o, --output FILE write output to FILE [standard output]
输出文件
-g, --BCF generate genotype likelihoods in BCF format
输出文件格式为BCF
-v, --VCF generate genotype likelihoods in VCF format
输出文件格式为VCF
Output options for mpileup format (without -g/-v):
-O, --output-BP output base positions on reads
输出reads中的碱基位置信息
-s, --output-MQ output mapping quality
输出mapping质量分值
-a output all positions (including zero depth)
输出所有位点
-a -a (or -aa) output absolutely all positions, including unused ref. sequences
输出所有位点，包含未比对的参考序列
Output options for genotype likelihoods (when -g/-v is used):
-t, --output-tags LIST optional tags to output:
DP,AD,ADF,ADR,SP,INFO/AD,INFO/ADF,INFO/ADR []
输出哪些可选标签
-u, --uncompressed generate uncompressed VCF/BCF output
结果不压缩
SNP/INDEL genotype likelihoods options (effective with -g/-v):
-e, --ext-prob INT Phred-scaled gap extension seq error probability [20]
空位延伸的测序错误概率用Phred分数表示，不能大于20
-F, --gap-frac FLOAT minimum fraction of gapped reads [0.002]
空位reads的最小比例
-h, --tandem-qual INT coefficient for homopolymer errors [100]
-I, --skip-indels do not perform indel calling
不call INDEL
-L, --max-idepth INT maximum per-file depth for INDEL calling [250]
进行INDEL calling的单文件最大深度
-m, --min-ireads INT minimum number gapped reads for indel candidates [1]
最少需要多少条gapped reads才能作为候选的INDEL
-o, --open-prob INT Phred-scaled gap open seq error probability [40]
可作为gap open候选的测序错误分值的最大分值
-p, --per-sample-mF apply -m and -F per-sample for increased sensitivity
同时设置-m和-F以提高结果的敏感性
-P, --platforms STR comma separated list of platforms for indels [all]
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPT=VALUE
--reference FILE
Reference sequence FASTA FILE [null]
Notes: Assuming diploid individuals.假定是二倍体
scaffold_1 2841 A 11 ,,,...,.... BHIGDGIJ?FF
scaffold_1 2842 C 12 ,$,,...,....^I. CFGEGEGGCFF+
scaffold_1 2843 G 11 ,,...,..... FDDDDCD?DD+
scaffold_1 2844 G 11 ,,...,..... FA?AAAA<AA+
scaffold_1 2845 G 11 ,,...,..... F656666166*
scaffold_1 2846 A 11 ,,...,..... (1.1111)11*
scaffold_1 2847 A 11 ,,+9acggtgaag.+9ACGGTGAAT.+9ACGGTGAAG.+9ACGGTGAAG,+9acggtgaag.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG %.+....-..)
scaffold_1 2848 N 11 agGGGgGGGGG !!$!!!!!!!!
scaffold_1 2849 A 11 c$,...,..... !0000000000
scaffold_1 2850 A 10 ,...,..... 353333333
文件一般包含6列，每一列的意义是：
- 参考序列
- 参考序列中的位置（坐标）
- 参考序列对应位置的核苷酸
- 比对上的reads的数量
- 比对的具体情况
- ","表示与参考序列的正链匹配
- "."表示与参考序列的负链匹配
- "[ATCGN]"表示在该位置上的碱基与参考序列的正链不匹配（mismatch）
- "[atcgn]"表示在该位置上的碱基与参考序列的负链不匹配（mismatch）
- "+[ATCGNatcgn]"表示在该位置插入了序列，大小写分别表示正链或者负链
- "-[ATCGNagcgn]"表示在该位置发生了DEL，大小写分别表示对应正链还是负链
- "^"表示一个read的开始，其后跟的ASCII码对应的是这个碱基的mappling质量的PHRED分值，再后面的才是这个位置的碱基
- "$"表示一个read的结束，修饰的是其前面的碱基
- "*"表示模糊（ambiguous）碱基
- 比对上的reads在该位置的mapping打分