0. samtools背景
（1）如果集群已经安装过conda，安装samtools使用命令：conda install -c bioconda samtools
（2）如果没有conda，这里参考👉conda安装链接
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.pngFinal.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 （生成位点匹配信息）
mpileup命令是samtools中一个很重要命令，也是samtools中call snp的工具。它的主要功能是，生成一个或多个bam文件的BCF、VCF文件或者pileup文本（染色体上每个碱基的比对情况汇总文件）。以 @RG中的样本名 作为区分标识符，对不同的样本记录。如果样本标识符缺失，那么每一个输入文件则视为一个样本。
pileup格式：每一行记录一个位点的测序信息，每一列的信息分别对应为染色体号（chromosome name）；物理位置（1-based coordinate）；参考基因组对应碱基类型（reference base）；该位置的测序深度（the number of reads covering the site）；reads读到的碱基类型（read bases）；reads对应碱基的测序质量值（base qualities）；比对质量值（alignment mapping qualities）。
其中，比对信息（如mismatch, indel等）都体现在read bases这一列：
.代表与参考序列正链匹配
,代表与参考序列负链匹配
ATCGN代表在正链上的不匹配
atcgn代表在负链上的不匹配
*碱基缺失
$代表一个read的结束，该符号修饰的是它前面的碱基
^代表一个read的开头
参考文件：http://www.htslib.org/doc/1.2/samtools.html
Usage: samtools mpileup [-EBug][-C capQcoef][-r reg][-f in.fa][-l list][-M capMapQ][-Q minBaseQ][-q minMapQ] in.bam [in2.bam [...]]
1. 输入相关参数:
-A --count-orphans ##在检测变异时，不忽略异常的reads对
-b --bam-list FILE #存储bam文件名称的文本（一行一个bam）
-B --no-BAQ #Disable probabilistic realignment for the computation of base alignment quality (BAQ). BAQ is the Phred-scaled probability of a read base being misaligned. Applying this option greatly helps to reduce false SNPs caused by misalignments.
-C --adjust-MQ INT #指定 mapping quality; BWA推荐为50（默认为0）
-d --max-depth INT #指定单碱基的最大测序深度（感觉没什么用，默认250）
-E --redo-BAQ #Recalculate BAQ on the fly, ignore existing BQ tags
-f FILE #输入有索引文件的fasta参考序列
-G --exclude-RG FILE #除去bam文件中指定readgroups的reads(在FILE文件中一行一个@RG-ID)
-l --positions FILE #接bed文件或包含区域位点的位置列表文件（BED or position list file containing a list of regions or sites where pileup or BCF should be generated. If BED, positions are 0-based half-open）
-q --min-MQ INT # 忽略比对值小于 INT的reads（默认为0）
-Q --min-BQ INT # 忽略碱基质量值小于INT的碱基（默认为13）
-r chr$:pos1-pos2 #在指定区域产生pileup结果（输出chr$:pos1-pos2区域的pileup结果，不加该参数默认all sites）
-R --ignore-RG #忽略 RG 标签，可以将Bam文件中所有的reads视为一个样本
--rf --incl-flags STR|INT #设置必备标签，忽略不包含“必备标签”的reads
--ff --excl-flags STR|INT #设置排除标签，忽略包含“排除标签”的reads（默认标签UNMAP,SECONDARY,QCFAIL,DUP）
-x --ignore-overlaps #忽略存在overlap的read-pair(存在overlap说明插入片段过短)
2. 输出文件的相关参数：
-o –output FILE #生成pileup格式文件或者VCF、BCF文件而不是默认的标准输出
-g --BCF FILE # 计算基因型的似然值和输出bcf形式文件，As of v1.0, this is BCF2 which is incompatible with the BCF1 format produced by previous (0.1.x) versions of samtools
-v --VCF FILE # 计算基因型的似然值和输出vcf形式文件，Output is bgzip-compressed VCF unless -u option is set
3. 设置mpileup输出格式(输出参数中没有 -g or -v)：
-O --output-BP #输出碱基在reads中的位置
-s --output-MQ #输出mapping quality
4. 设定VCF/BCF格式的参数(with -g or -v)：
-D #输出每个样本的测序深度 [DEPRECATED - use -t DP instead]
-S #Output per-sample Phred-scaled strand bias P-value [DEPRECATED - use -t SP instead]
-t --output-tags LIST #Comma-separated list of FORMAT and INFO tags to output (case-insensitive): DP (Number of high-quality bases, FORMAT), DV (Number of high-quality non-reference bases, FORMAT), DPR (Number of high-quality bases for each observed allele, FORMAT), INFO/DPR (Number of high-quality bases for each observed allele, INFO), DP4 (Number of high-quality ref-forward, ref-reverse, alt-forward and alt-reverse bases, FORMAT), SP (Phred-scaled strand bias P-value, FORMAT) [null]
-u --uncompressed #生成不压缩的VCF/BCF
-V #输出每个样本中未比对到参考基因组的reads数目 [DEPRECATED - use -t DV instead]
5. 计算 SNP/INDEL突变类型所使用参数(for -g or -v)：
-e --ext-prob INT #设定gap extension的测序错误概率（降低INT值：更容易call到gap较长的indel突变，默认值为20）#Phred-scaled gap extension sequencing error probability. Reducing INT leads to longer indels
-F --gap-frac FLOAT #设置gap reads最低频率（默认0.002）
-h --tandem-qual INT #用来矫正连续相同碱基区（homopolymer）的indel错误（默认值为100）：对于一段由l个碱基组成的homopolymer序列，其检测到s bp 插入或缺失，它的测序错误率建模为INT*s/l
-I --skip-indels #不检测indel
-L --max-idepth INT #如果平均每个样本的测序深度高于INT，不再检测indel（默认250）
-m --min-ireads INT #设定 支持indel的reads数目最小值（默认为1）
-o --open-prob INT #设定gap open的测序错误概率（降低INT值：更容易call到indel突变，默认值为40）#(The same short option is used for both --open-prob and --output. When -o's argument contains only an optional + or - sign followed by the digits 0 to 9, it is interpreted as --open-prob.)
-p --per-sample-mF #Apply -m and -F thresholds per sample to increase sensitivity of calling. By default both options are applied to reads pooled from all samples.
-P --platforms STR #Comma-delimited list of platforms (determined by @RG-PL) from which indel candidates are obtained. It is recommended to collect indel candidates from sequencing technologies that have low indel error rate such as ILLUMINA. [all]
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数）
PS：
纯属个人见解，如有错误欢迎批评指正～
