轉(zhuǎn)自：https://blog.csdn.net/sinat_38163598/article/details/72910115
samtools是一個(gè)用于操作sam和bam文件（通常是短序列比對(duì)工具 bwa,bowtie2,hisat2,tophat2等等產(chǎn)生的术陶，具體格式可以在消息框輸入“SAM”查看）的工具合集紊馏，包含有許多命令。以下是常用命令的介紹棉姐。
1.View
view命令的主要功能是：將sam文件與bam文件互換届案；然后對(duì)bam文件進(jìn)行各種操作庵楷，比如數(shù)據(jù)的排序(sort)和提取(這些操作
是對(duì)bam文件進(jìn)行的，因而當(dāng)輸入為sam文件的時(shí)候萝玷，不能進(jìn)行該操作)嫁乘；最后將排序或提取得到的數(shù)據(jù)輸出為bam或sam（默認(rèn)的）格式。
bam文件優(yōu)點(diǎn)：bam文件為二進(jìn)制文件球碉，占用的磁盤(pán)空間比sam文本文件序迅；利用bam二進(jìn)制文件的運(yùn)算速度快睁冬。
view命令中挎春，對(duì)sam文件頭部（序列ID）的輸入(-t或-T）和輸出(-h)是單獨(dú)的一些參數(shù)來(lái)控制的。
Usage: samtools view [options] | [region1 [...]]
默認(rèn)情況下不加 region豆拨，則是輸出所有的 region.
options:
-b?output?BAM
??默認(rèn)下輸出是?SAM?格式文件直奋，該參數(shù)設(shè)置輸出?BAM?格式
?-h?print?header?for?the?SAM?output
??默認(rèn)下輸出的?sam?格式文件不帶?header，該參數(shù)設(shè)定輸出sam文件時(shí)帶?header?信息
?-H?print?header?only?(no?alignments)
??僅僅輸出文件的頭
?-S?input?is?SAM
??默認(rèn)下輸入是?BAM?文件施禾，若是輸入是?SAM?文件脚线，則最好加該參數(shù)，否則有時(shí)候會(huì)報(bào)錯(cuò)弥搞。
?-u?uncompressed?BAM?output?(force?-b)
??該參數(shù)的使用需要有-b參數(shù)邮绿，能節(jié)約時(shí)間渠旁，但是需要更多磁盤(pán)空間。
?-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.
??過(guò)濾功統(tǒng)計(jì)功能
?-c?print?only?the?count?of?matching?records
?-L?FILE??output?alignments?overlapping?the?input?BED?FILE?[null]
?-t?FILE??list?of?reference?names?and?lengths?(force?-S)?[null]
??使用一個(gè)list文件來(lái)作為header的輸入
?-T?FILE??reference?sequence?file?(force?-S)?[null]
??使用序列fasta文件作為header的輸入
?-o?FILE??output?file?name?[stdout]
?-F?INT???filtering?flag,?0?for?unset?[0]?
??Skip?alignments?with?bits?present?in?INT?[0]
??數(shù)字4代表該序列沒(méi)有比對(duì)到參考序列上
??數(shù)字8代表該序列的mate序列沒(méi)有比對(duì)到參考序列上
??過(guò)濾功能船逮。如F12過(guò)濾只有雙端map的
?-q?INT???minimum?mapping?quality?[0]
????比對(duì)的最低質(zhì)量值顾腊，一般認(rèn)為20就為unique比對(duì)了，可以結(jié)合上述-bF參數(shù)使用使用提取特定的比對(duì)結(jié)果
例子：
將sam文件轉(zhuǎn)換成bam文件
samtools view -bS abc.sam > abc.bam
BAM轉(zhuǎn)換為SAM
samtools view -h -o out.sam out.bam
提取比對(duì)到參考序列上的比對(duì)結(jié)果
samtools view -bF 4 abc.bam > abc.F.bam
提取paired reads中兩條reads都比對(duì)到參考序列上的比對(duì)結(jié)果挖胃，只需要把兩個(gè)4+8的值12作為過(guò)濾參數(shù)即可
samtools view -bF 12 abc.bam > abc.F12.bam
提取沒(méi)有比對(duì)到參考序列上的比對(duì)結(jié)果
samtools view -bf 4 abc.bam > abc.f.bam
提取bam文件中比對(duì)到caffold1上的比對(duì)結(jié)果杂靶，并保存到sam文件格式
samtools view abc.bam scaffold1 > scaffold1.sam
提取scaffold1上能比對(duì)到30k到100k區(qū)域的比對(duì)結(jié)果
samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam
根據(jù)fasta文件，將 header 加入到 sam 或 bam 文件中
samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam
2.Sort
sort對(duì)bam文件進(jìn)行排序酱鸭。一些軟件需要sort的bam或者sam文件吗垮，如stringtie，所以必須要sort使用凹髓；求depth時(shí)抱既，也必須要sort；
Usage: samtools sort [-n] [-m ] ?
-m 內(nèi)存參數(shù)默認(rèn)下是 500,000,000 即500M（不支持K扁誓，M，G等縮寫(xiě)）蚀之。對(duì)于處理大數(shù)據(jù)時(shí)蝗敢，如果內(nèi)存夠用，則設(shè)置大點(diǎn)的值足删，以節(jié)約時(shí)間寿谴。
-n 設(shè)定排序方式按short reads的ID排序。默認(rèn)下是按序列在fasta文件中的順序（即header）和序列從左往右的位點(diǎn)排序失受。
例子：
samtools sort accept.bam accept.sort最終產(chǎn)生accept.sort.bam
3.merge
將2個(gè)或2個(gè)以上的已經(jīng)sort了的bam文件融合成一個(gè)bam文件讶泰。融合后的文件已經(jīng)sort過(guò)了的。
Usage:? samtools merge [-nr] [-h inh.sam] [...]Options: -n? ? ? sort by read names? ? ? ? -r? ? ? attach RG tag (inferred from file names)? ? ? ? -u? ? ? uncompressed BAM output? ? ? ? -f? ? ? overwrite the output BAM if exist? ? ? ? -1? ? ? compress level 1? ? ? ? -R STR? merge file in the specified region STR [all]? ? ? ? -h FILE? copy the header in FILE to [in1.bam]
例子：
4.index
必須對(duì)bam文件進(jìn)行默認(rèn)情況下的排序后拂到，才能進(jìn)行index痪署。否則會(huì)報(bào)錯(cuò)。
建立索引后將產(chǎn)生后綴為.bai的文件兄旬，用于快速的隨機(jī)處理狼犯。很多情況下需要有bai文件的存在，特別是顯示序列比對(duì)情況下领铐。比如samtool的tview命令就需要悯森；gbrowse2顯示reads的比對(duì)圖形的時(shí)候也需要。IGV顯示比對(duì)情況也需要绪撵。
Usage: samtools index [out.index]
例子：
以下兩種命令結(jié)果一樣
$ samtools index abc.sort.bam$ samtools index abc.sort.bam abc.sort.bam.bai
5.faidx
對(duì)fasta文件建立索引,生成的索引文件以.fai后綴結(jié)尾瓢姻。該命令也能依據(jù)索引文件快速提取fasta文件中的某一條（子）序列
Usage: samtools faidx [ [...]]對(duì)基因組文件建立索引，方便提取序列音诈。
例子：$ samtools faidx genome.fasta
由于有索引文件幻碱，可以使用以下命令很快從基因組中提取到fasta格式的子序列
$ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta
6.tview
tview能直觀的顯示出reads比對(duì)基因組的情況绎狭，和基因組瀏覽器有點(diǎn)類(lèi)似。
需要事先利用利用上面講的sort和建index命令執(zhí)行完后收班，用下述命令坟岔。
Usage: samtools tview [ref.fasta]
出參考基因組的時(shí)候，會(huì)在第一排顯示參考基因組的序列摔桦，否則社付，第一排全用N表示。按下 g 邻耕，則提示輸入要到達(dá)基因組的某一個(gè)位點(diǎn)鸥咖。例子“scaffold_10:1000"表示到達(dá)第10號(hào)scaffold的第1000個(gè)堿基位點(diǎn)處。使用H(左）J（上）K（下）L（右）移動(dòng)顯示界面兄世。大寫(xiě)字母移動(dòng)快啼辣，小寫(xiě)字母移動(dòng)慢。使用空格建向左快速移動(dòng)（和 L 類(lèi)似）御滩，使用Backspace鍵向左快速移動(dòng)（和 H 類(lèi)似）鸥拧。Ctrl+H 向左移動(dòng)1kb堿基距離； Ctrl+L 向右移動(dòng)1kb堿基距離可以用顏色標(biāo)注比對(duì)質(zhì)量削解，堿基質(zhì)量富弦，核苷酸等。30～40的堿基質(zhì)量或比對(duì)質(zhì)量使用白色表示氛驮；20～30黃色腕柜；10～20綠色；0～10藍(lán)色矫废。使用點(diǎn)號(hào)'.'切換顯示堿基和點(diǎn)號(hào)盏缤；使用r切換顯示read name等還有很多其它的使用說(shuō)明，具體按 蓖扑？ 鍵來(lái)查看唉铜。
7.flagstat
給出BAM文件的比對(duì)結(jié)果
Usage: samtools flagstat
$ samtools flagstat example.bam
11945742 + 0 in total (QC-passed reads + QC-failed reads)
#總共的reads數(shù)
0 + 0 duplicates
7536364 + 0 mapped (63.09%:-nan%)
#總體上reads的匹配率
11945742 + 0 paired in sequencing
#有多少reads是屬于paired reads
5972871 + 0 read1
#reads1中的reads數(shù)
5972871 + 0 read2
#reads2中的reads數(shù)
6412042 + 0 properly paired (53.68%:-nan%)
#完美匹配的reads數(shù)：比對(duì)到同一條參考序列，并且兩條reads之間的距離符合設(shè)置的閾值
6899708 + 0 with itself and mate mapped
#paired reads中兩條都比對(duì)到參考序列上的reads數(shù)
636656 + 0 singletons (5.33%:-nan%)
#單獨(dú)一條匹配到參考序列上的reads數(shù)赵誓，和上一個(gè)相加打毛，則是總的匹配上的reads數(shù)。
469868 + 0 with mate mapped to a different chr
#paired reads中兩條分別比對(duì)到兩條不同的參考序列的reads數(shù)
243047 + 0 with mate mapped to a different chr (mapQ>=5)
#同上一個(gè)俩功，只是其中比對(duì)質(zhì)量>=5的reads的數(shù)量
8.depth
得到每個(gè)堿基位點(diǎn)的測(cè)序深度,并輸出到標(biāo)準(zhǔn)輸出幻枉，所以要用大于號(hào)追加到一個(gè)文件。
Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] [...]
-r 后面跟染色體號(hào)（region）
-q ：計(jì)算深度時(shí)要求測(cè)序堿基質(zhì)量最低質(zhì)量值
-Q：計(jì)算深度時(shí)要求比對(duì)的最低質(zhì)量值
注意：做depth之前必須做samtools index诡蜓；
例子
samtools depth accept.bam >depth
9.其他命令
reheader：替換bam文件的頭
$ samtools reheader
idxstats ：統(tǒng)計(jì)一個(gè)表格熬甫，4列，分別為”序列名蔓罚，序列長(zhǎng)度椿肩，比對(duì)上的reads數(shù)瞻颂，unmapped reads
number。第4列應(yīng)該是paired
reads中有一端能匹配到該scaffold上郑象，而另外一端不匹配到任何scaffolds上的reads數(shù)贡这。
$ samtools idxstats
rmdup：將由PCR duplicates獲得的reads去掉，并只保留最高比對(duì)質(zhì)量的read厂榛。
Usage: ?samtools rmdup [-sS]
-s 對(duì)single-end reads盖矫。默認(rèn)情況下，只對(duì)paired-end reads-S 將Paired-end reads作為single-end reads處理击奶。
10. 將bam文件轉(zhuǎn)換為fastq文件
有時(shí)候辈双，我們需要提取出比對(duì)到一段參考序列的reads，進(jìn)行小范圍的分析柜砾，以利于debug等湃望。這時(shí)需要將bam或sam文件轉(zhuǎn)換為fastq格式。
該網(wǎng)站提供了一個(gè)bam轉(zhuǎn)換為fastq的程序：http://www.hudsonalpha.org/gsl/information/software/bam2fastq
$ wget http://www.hudsonalpha.org/gsl/static/software/bam2fastq-1.1.0.tgz$ tar zxf bam2fastq-1.1.0.tgz$ cd bam2fastq-1.1.0$ make$ ./bam2fastq
11. mpileup
samtools還有個(gè)非常重要的命令mpileup痰驱，以前為pileup证芭。該命令用于生成bcf文件，再使用bcftools進(jìn)行SNP和Indel的分析担映。bcftools是samtool中附帶的軟件檩帐，在samtools的安裝文件夾中可以找到。
最常用的參數(shù)有2：
?-f 來(lái)輸入有索引文件的fasta參考序列另萤； -g 輸出到bcf格式。用法和最簡(jiǎn)單的例子如下
Usage: samtools mpileup [-EBug] [-CcapQcoef] [-rreg] [-fin.fa] [-llist] [-McapMapQ] [-QminBaseQ] [-qminMapQ]in.bam[in2.bam[...]]$ samtools mpileup -f genome.fasta abc.bam > abc.txt
$ samtools mpileup -gSDf genome.fasta abc.bam > abc.bcf
$ samtools mpileup -guSDf genome.fasta abc.bam | \
? ? ? ? ? bcftools view -cvNg - > abc.vcf
mpileup不使用-u或-g參數(shù)時(shí)诅挑，則不生成二進(jìn)制的bcf文件四敞，而生成一個(gè)文本文件(輸出到標(biāo)準(zhǔn)輸出)。該文本文件統(tǒng)計(jì)了參考序列中每個(gè)堿基位點(diǎn)的比對(duì)情況拔妥；該文件每一行代表了參考序列中某一個(gè)堿基位點(diǎn)的比對(duì)結(jié)果忿危。比如：
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
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
mpileup生成的結(jié)果包含6行：參考序列名；位置没龙；參考?jí)A基铺厨；比對(duì)上的reads數(shù)；比對(duì)情況硬纤；比對(duì)上的堿基的質(zhì)量解滓。其中第5列比較復(fù)雜,解釋如下：
1 ‘.’代表與參考序列正鏈匹配。
2 ‘,’代表與參考序列負(fù)鏈匹配筝家。
3 ‘ATCGN’代表在正鏈上的不匹配洼裤。
4 ‘a(chǎn)tcgn’代表在負(fù)鏈上的不匹配。
5 ‘*’代表模糊堿基
6 ‘^’代表匹配的堿基是一個(gè)read的開(kāi)始溪王；’^'后面緊跟的ascii碼減去33代表比對(duì)質(zhì)量腮鞍；這兩個(gè)符號(hào)修飾的是后面的堿基值骇，其后緊跟的堿基(.,ATCGatcgNn)代表該read的第一個(gè)堿基。
7 ‘$’代表一個(gè)read的結(jié)束移国，該符號(hào)修飾的是其前面的堿基吱瘩。
8 正則式’\+[0-9]+[ACGTNacgtn]+’代表在該位點(diǎn)后插入的堿基；比如上例中在scaffold_1的2847后插入了9個(gè)長(zhǎng)度的堿基acggtgaag迹缀。表明此處極可能是indel使碾。
9 正則式’-[0-9]+[ACGTNacgtn]+’代表在該位點(diǎn)后缺失的堿基；
12. 使用bcftools
bcftools和samtools類(lèi)似裹芝，用于處理vcf(variant call format)文件和bcf(binary call format)文件部逮。前者為文本文件，后者為其二進(jìn)制文件嫂易。
bcftools使用簡(jiǎn)單兄朋，最主要的命令是view命令，其次還有index和cat等命令怜械。index和cat命令和samtools中類(lèi)似颅和。此處主講使用view命令來(lái)進(jìn)行SNP和Indel calling。該命令的使用方法和例子為：
$ bcftools view [-AbFGNQSucgv] [-DseqDict] [-llistLoci] [-slistSample] ? ? ? ? ? [-igapSNPratio] [-tmutRate] [-pvarThres] [-Pprior] ? ? ? ? ? [-1nGroup1] [-dminFrac] [-UnPerm] [-XpermThres] ? ? ? ? ? [-TtrioType]in.bcf[region]
$ bcftools view -cvNg abc.bcf > snp_indel.vcf
生成的結(jié)果文件為vcf格式缕允，有10列峡扩，分別是：1 參考序列名；2 varianti所在的left-most位置障本；3 variant的ID（默認(rèn)未設(shè)置教届，用’.'表示)；4 參考序列的allele驾霜；5 variant的allele(有多個(gè)alleles案训，則用’,'分隔);6 variant/reference QUALity;7 FILTers applied;8 variant的信息，使用分號(hào)隔開(kāi)粪糙；9 FORMAT of the genotype fields, separated by colon (optional)强霎； 10 SAMPLE genotypes and per-sample information (optional)。
例如：
scaffold_1? ? ? 2847? ? .? ? ? A? ? ? AACGGTGAAG? ? ? 194? ? .? ? ? INDEL;DP=11;VDB=0.0401;AF1=1;AC1=2;DP4=0,0,8,3;MQ=35;FQ=-67.5? GT:PL:GQ? ? ? ? 1/1:235,33,0:63
scaffold_1? ? ? 3908? ? .? ? ? G? ? ? A? ? ? 111? ? .? ? ? DP=13;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,5,7;MQ=42;FQ=-63? GT:PL:GQ? ? ? ? 1/1:144,36,0:69
scaffold_1? ? ? 4500? ? .? ? ? A? ? ? G? ? ? 31.5? ? .? ? ? DP=8;VDB=0.0034;AF1=1;AC1=2;DP4=0,0,1,3;MQ=42;FQ=-39? ? GT:PL:GQ? ? ? ? 1/1:64,12,0:21
scaffold_1? ? ? 4581? ? .? ? ? TGGNGG? TGG? ? 145? ? .? ? ? INDEL;DP=8;VDB=0.0308;AF1=1;AC1=2;DP4=0,0,0,8;MQ=42;FQ=-58.5? ? GT:PL:GQ? ? ? ? 1/1:186,24,0:45
scaffold_1? ? ? 4644? ? .? ? ? G? ? ? A? ? ? 195? ? .? ? ? DP=21;VDB=0.0198;AF1=1;AC1=2;DP4=0,0,10,10;MQ=42;FQ=-87 GT:PL:GQ? ? ? ? 1/1:228,60,0:99
scaffold_1? ? ? 4827? ? .? ? ? NACAAAGA? ? ? ? NA? ? ? 4.42? ? .? ? ? INDEL;DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=40;FQ=-37.5? ? ? GT:PL:GQ? ? ? ? 0/1:40,3,0:3
scaffold_1? ? ? 4854? ? .? ? ? A? ? ? G? ? ? 48? ? ? .? ? ? DP=6;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,2,1;MQ=41;FQ=-36? ? GT:PL:GQ? ? ? ? 1/1:80,9,0:16
scaffold_1? ? ? 5120? ? .? ? ? A? ? ? G? ? ? 85? ? ? .? ? ? DP=8;VDB=0.0355;AF1=1;AC1=2;DP4=0,0,5,3;MQ=42;FQ=-51? ? GT:PL:GQ? ? ? ? 1/1:118,24,0:45
第8列中顯示了對(duì)variants的信息描述蓉冈，比較重要城舞，其中的 Tag 的描述如下：
Tag Format Description
AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
DP int Raw read depth (without quality filtering)
DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
MQ int Root-Mean-Square mapping quality of covering reads
PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
QCHI2 int Phred-scaled PCHI2
RP int # permutations yielding a smaller PCHI2
CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
UGT string Most probable genotype configuration without the trio constraint
CGT string Most probable configuration with the trio constraint
使用bcftools得到variant calling結(jié)果后。需要對(duì)結(jié)果再次進(jìn)行過(guò)濾寞酿。主要依據(jù)比對(duì)結(jié)果中第8列信息家夺。其中的 DP4
一行尤為重要，提供了4個(gè)數(shù)據(jù)：1 比對(duì)結(jié)果和正鏈一致的reads數(shù)伐弹、2 比對(duì)結(jié)果和負(fù)鏈一致的reads數(shù)秦踪、3
比對(duì)結(jié)果在正鏈的variant上的reads數(shù)、4 比對(duì)結(jié)果在負(fù)鏈的variant上的reads數(shù)∫蔚耍可以設(shè)定 （value3 +
value4）大于某一閾值柠逞，才算是variant。比如：
$ perl -ne 'print $_ if /DP4=(\d+),(\d+),(\d+),(\d+)/ && ($3+$4)>=10 && ($3+$4)/($1+$2+$3+$4)>=0.8' snp_indel.vcf > snp_indel.final.vcf
