samtools的說明文檔：http://samtools.sourceforge.net/samtools.shtml
samtools是一個用於操作sam和bam文件的工具合集。包含有許多命令。以下是常用命令的介紹
1. view
view命令的主要功能是：將sam文件轉換成bam文件；然后對bam文件進行各種操作，比如數據的排序(不屬於本命令的功能)和提取(這些操作 是對bam文件進行的，因而當輸入為sam文件的時候，不能進行該操作)；最后將排序或提取得到的數據輸出為bam或sam（默認的）格式。
bam文件優點：bam文件為二進制文件，占用的磁盤空間比sam文本文件小；利用bam二進制文件的運算速度快。
view命令中，對sam文件頭部的輸入(-t或-T）和輸出(-h)是單獨的一些參數來控制的。
Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]
默認情況下不加 region，則是輸出所有的 region.
Options: -b output BAM
默認下輸出是 SAM 格式文件，該參數設置輸出 BAM 格式
-h print header for the SAM output
默認下輸出的 sam 格式文件不帶 header，該參數設定輸出sam文件時帶 header 信息
-H print header only (no alignments)
-S input is SAM
默認下輸入是 BAM 文件，若是輸入是 SAM 文件，則最好加該參數，否則有時候會報錯。
-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)
-X output FLAG in string (samtools-C specific)
-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]
使用一個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]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
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]
-r STR only output reads in read group STR [null]
-s FLOAT fraction of templates to subsample; integer part as seed [-1]
-? longer help
例子：
將sam文件轉換成bam文件
$ samtools view -bS abc.sam > abc.bam
$ samtools view -b -S abc.sam -o abc.bam
提取比對到參考序列上的比對結果
$ samtools view -bF 4 abc.bam > abc.F.bam 提取paired reads中兩條reads都比對到參考序列上的比對結果，只需要把兩個4+8的值12作為過濾參數即可 $ samtools view -bF 12 abc.bam > abc.F12.bam 提取沒有比對到參考序列上的比對結果 $ samtools view -bf 4 abc.bam > abc.f.bam 提取bam文件中比對到caffold1上的比對結果，並保存到sam文件格式 $ samtools view abc.bam scaffold1 > scaffold1.sam 提取scaffold1上能比對到30k到100k區域的比對結果 $ samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam 根據fasta文件，將 header 加入到 sam 或 bam 文件中 $ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam
2. sort
sort對bam文件進行排序。
Usage: samtools sort [-n] [-m <maxMem>] <in.bam> <out.prefix> -m 參數默認下是 500,000,000 即500M（不支持K，M，G等縮寫）。對於處理大數據時，如果內存夠用，則設置大點的值，以節約時間。 -n 設定排序方式按short reads的ID排序。默認下是按序列在fasta文件中的順序（即header）和序列從左往右的位點排序。
例子：
$ samtools sort abc.bam abc.sort ###注意 abc.sort 是輸出文件的前綴，實際輸出是 abc.sort.bam
$ samtools view abc.sort.bam | less -S
3.merge
將2個或2個以上的已經sort了的bam文件融合成一個bam文件。融合后的文件不需要則是已經sort過了的。
Usage: samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam>[...]
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 <out.bam> [in1.bam]
Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users
must provide the correct header with -h, or uses Picard which properly maintains
the header dictionary in merging.
4.index
必須對bam文件進行默認情況下的排序后，才能進行index。否則會報錯。
建立索引后將產生后綴為.bai的文件，用於快速的隨機處理。很多情況下需要有bai文件的存在，特別是顯示序列比對情況下。比如samtool的tview命令就需要；gbrowse2顯示reads的比對圖形的時候也需要。
Usage: samtools index <in.bam> [out.index]
例子：
以下兩種命令結果一樣
$ samtools index abc.sort.bam
$ samtools index abc.sort.bam abc.sort.bam.bai
5. faidx
對fasta文件建立索引,生成的索引文件以.fai后綴結尾。該命令也能依據索引文件快速提取fasta文件中的某一條（子）序列
Usage: samtools faidx <in.bam> [ [...]] 對基因組文件建立索引 $ samtools faidx genome.fasta 生成了索引文件genome.fasta.fai,是一個文本文件，分成了5列。第一列是子序列的名稱； 第二列是子序列的長度；個人認為“第三列是序列所在的位置”，因為該數字從上往下逐漸變大， 最后的數字是genome.fasta文件的大小；第4和5列不知是啥意思。於是通過此文件，可以定 位子序列在fasta文件在磁盤上的存放位置，直接快速調出子序列。 由於有索引文件，可以使用以下命令很快從基因組中提取到fasta格式的子序列 $ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta
6. tview
tview能直觀的顯示出reads比對基因組的情況，和基因組瀏覽器有點類似。
Usage: samtools tview <aln.bam> [ref.fasta] 當給出參考基因組的時候，會在第一排顯示參考基因組的序列，否則，第一排全用N表示。 按下 g ，則提示輸入要到達基因組的某一個位點。例子“scaffold_10:1000"表示到達第 10號scaffold的第1000個鹼基位點處。 使用H(左）J（上）K（下）L（右）移動顯示界面。大寫字母移動快，小寫字母移動慢。 使用空格建向左快速移動（和 L 類似），使用Backspace鍵向左快速移動（和 H 類似）。 Ctrl+H 向左移動1kb鹼基距離； Ctrl+L 向右移動1kb鹼基距離 可以用顏色標注比對質量，鹼基質量，核苷酸等。30～40的鹼基質量或比對質量使用白色表示； 20～30黃色；10～20綠色；0～10藍色。 使用點號'.'切換顯示鹼基和點號；使用r切換顯示read name等 還有很多其它的使用說明，具體按 ？ 鍵來查看。
7. flagstat
給出BAM文件的比對結果
Usage: samtools flagstat <in.bam> $ 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的數量
7. depth
得到每個鹼基位點的測序深度,並輸出到標准輸出。
Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] <in1.bam> [...]
-r 后面跟染色體號（region）
注意：做depth之前必須做samtools index；
示例：
samtools depth in.bam > out.depth.txt
注意： in.bam 必須經過了排序。
8. 其它有用的命令
reheader 替換bam文件的頭
$ samtools reheader <in.header.sam> <in.bam>
cat 連接多個bam文件，適用於非sorted的bam文件
$ samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [ ... ]
idxstats 統計一個表格，4列，分別為”序列名，序列長度，比對上的reads數，unmapped reads number”。第4列應該是paired reads中有一端能匹配到該scaffold上，而另外一端不匹配到任何scaffolds上的reads數。
$ samtools idxstats <aln.bam>
9. 將bam文件轉換為fastq文件
有時候，我們需要提取出比對到一段參考序列的reads，進行小范圍的分析，以利於debug等。這時需要將bam或sam文件轉換為fastq格式。
該網站提供了一個bam轉換為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 <in.bam>
10. mpileup
samtools還有個非常重要的命令mpileup，以前為pileup。該命令用於生成bcf文件，再使用bcftools進行SNP和Indel的分析。bcftools是samtool中附帶的軟件，在samtools的安裝文件夾中可以找到。
最常用的參數有2： -f 來輸入有索引文件的fasta參考序列； -g 輸出到bcf格式。用法和最簡單的例子如下
Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] 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參數時，則不生成二進制的bcf文件，而生成一個文本文件(輸出到標准輸出)。該文本文件統計了參考序列中每個鹼基位點的比對情況；該文件每一行代表了參考序列中某一個鹼基位點的比對結果。比如：
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
mpileup生成的結果包含6行：參考序列名；位置；參考鹼基；比對上的reads數；比對情況；比對上的鹼基的質量。其中第5列比較復雜,解釋如下：
1 ‘.’代表與參考序列正鏈匹配。
2 ‘,’代表與參考序列負鏈匹配。
3 ‘ATCGN’代表在正鏈上的不匹配。
4 ‘atcgn’代表在負鏈上的不匹配。
5 ‘*’代表模糊鹼基
6 ‘^’代表匹配的鹼基是一個read的開始；’^'后面緊跟的ascii碼減去33代表比對質量；這兩個符號修飾的是后面的鹼基，其后緊跟的鹼基(.,ATCGatcgNn)代表該read的第一個鹼基。
7 ‘$’代表一個read的結束，該符號修飾的是其前面的鹼基。
8 正則式’\+[0-9]+[ACGTNacgtn]+’代表在該位點后插入的鹼基；比如上例中在scaffold_1的2847后插入了9個長度的鹼基acggtgaag。表明此處極可能是indel。
9 正則式’-[0-9]+[ACGTNacgtn]+’代表在該位點后缺失的鹼基；
pileup具體的參數如下：
輸入參數
-6 Assume the quality is in the Illumina 1.3+ encoding. -A Do not skip anomalous read pairs in variant calling.
-B 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.
-b FILE List of input BAM files, one file per line [null]
-C INT Coefficient for downgrading mapping quality for reads containing excessive mismatches. Given a read with a phred-scaled probability q of being generated from the mapped position, the new mapping quality is about sqrt((INT-q)/INT)*INT. A zero value disables this functionality; if enabled, the recommended value for BWA is 50. [0]
-d INT At a position, read maximally INT reads per input BAM. [250]
-E Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt specificity a little bit.
-f FILE The faidx-indexed reference file in the FASTA format. The file can be optionally compressed by razip. [null]
-l FILE BED or position list file containing a list of regions or sites where pileup or BCF should be generated [null]
-M INT cap mapping quality at INT [60]
-q INT Minimum mapping quality for an alignment to be used [0]
-Q INT Minimum base quality for a base to be considered [13]
-r STR Only generate pileup in region STR [all sites]
輸出參數
-D Output per-sample read depth (require -g/-u)
-g Compute genotype likelihoods and output them in the binary call format (BCF). -S Output per-sample Phred-scaled strand bias P-value (require -g/-u) -u Similar to -g except that the output is uncompressed BCF, which is preferred for piping. Options for Genotype Likelihood Computation (for -g or -u): -e INT Phred-scaled gap extension sequencing error probability. Reducing INT leads to longer indels. [20] -h INT Coefficient for modeling homopolymer errors. Given an l-long homopolymer run, the sequencing error of an indel of size s is modeled as INT*s/l. [100] -I Do not perform INDEL calling -L INT Skip INDEL calling if the average per-sample depth is above INT. [250] -o INT Phred-scaled gap open sequencing error probability. Reducing INT leads to more indel calls. [40] -P STR Comma dilimited 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]
11. 使用bcftools
bcftools和samtools類似，用於處理vcf(variant call format)文件和bcf(binary call format)文件。前者為文本文件，后者為其二進制文件。
bcftools使用簡單，最主要的命令是view命令，其次還有index和cat等命令。index和cat命令和samtools中類似。此處主講使用view命令來進行SNP和Indel calling。該命令的使用方法和例子為：
$ bcftools view [-AbFGNQSucgv] [-D seqDict] [-l listLoci] [-s listSample]
[-i gapSNPratio] [-t mutRate] [-p varThres] [-P prior]
[-1 nGroup1] [-d minFrac] [-U nPerm] [-X permThres]
[-T trioType] in.bcf [region]
$ bcftools view -cvNg abc.bcf > snp_indel.vcf
生成的結果文件為vcf格式，有10列，分別是：1 參考序列名；2 varianti所在的left-most位置；3 variant的ID（默認未設置，用’.'表示)；4 參考序列的allele；5 variant的allele(有多個alleles，則用’,'分隔);6 variant/reference QUALity;7 FILTers applied;8 variant的信息，使用分號隔開；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列中顯示了對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 view 的具體參數如下：
Input/Output Options:
-A Retain all possible alternate alleles at variant sites. By default, the view command discards unlikely alleles.
-b Output in the BCF format. The default is VCF.
-D FILE Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
-F Indicate PL is generated by r921 or before (ordering is different).
-G Suppress all individual genotype information.
-l FILE List of sites at which information are outputted [all sites]
-N Skip sites where the REF field is not A/C/G/T
-Q Output the QCALL likelihood format
-s FILE List of samples to use. The first column in the input gives the sample names and the second gives the ploidy, which can only be 1 or 2. When the 2nd column is absent, the sample ploidy is assumed to be 2. In the output, the ordering of samples will be identical to the one in FILE. [null]
-S The input is VCF instead of BCF.
-u Uncompressed BCF output (force -b).
Consensus/Variant Calling Options:
-c Call variants using Bayesian inference. This option automatically invokes option -e. -d FLOAT When -v is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0] 當有多個sample用於variants calling時，比如多個轉錄組數據或多個重測序 數據需要比對到參考基因組上，設置該值，表明至少有該<float 0~1>比例的 samples在該位點都有覆蓋才計算入variant.所以對於只有一個sample的情況 下，該值設置在0～1之間沒有意義，大於1則得不到任何結果。 -e Perform max-likelihood inference only, including estimating the site allele frequency, testing Hardy-Weinberg equlibrium and testing associations with LRT. -g Call per-sample genotypes at variant sites (force -c) -i FLOAT Ratio of INDEL-to-SNP mutation rate [0.15] -p FLOAT A site is considered to be a variant if P(ref|D) -t FLOAT Scaled muttion rate for variant calling [0.001] -T STR Enable pair/trio calling. For trio calling, option -s is usually needed to be applied to configure the trio members and their ordering. In the file supplied to the option -s, the first sample must be the child, the second the father and the third the mother. The valid values of STR are ‘pair’, ‘trioauto’, ‘trioxd’ and ‘trioxs’, where ‘pair’ calls differences between two input samples, and ‘trioxd’ (‘trioxs’) specifies that the input is from the X chromosome non-PAR regions and the child is a female (male). [null] -v Output variant sites only (force -c) Contrast Calling and Association Test Options: -1 INT Number of group-1 samples. This option is used for dividing the samples into two groups for contrast SNP calling or association test. When this option is in use, the following VCF INFO will be outputted: PC2, PCHI2 and QCHI2. [0] -U INT Number of permutations for association test (effective only with -1) [0] -X FLOAT Only perform permutations for P(chi^2)
使用bcftools得到variant calling結果后。需要對結果再次進行過濾。主要依據比對結果中第8列信息。其中的 DP4 一行尤為重要，提供了4個數據：1 比對結果和正鏈一致的reads數、2 比對結果和負鏈一致的reads數、3 比對結果在正鏈的variant上的reads數、4 比對結果在負鏈的variant上的reads數。可以設定 （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
12. samtools rmdup
NGS上機測序前需要進行PCR一步，使一個模板擴增出一簇，從而在上機測序的時候表現出為1個點，即一個reads。若一個模板擴增出了多簇，結 果得到了多個reads，這些reads的坐標(coordinates)是相近的。在進行了reads比對后需要將這些由PCR duplicates獲得的reads去掉，並只保留最高比對質量的read。使用rmdup命令即可完成.
Usage: samtools rmdup [-sS] -s 對single-end reads。默認情況下，只對paired-end reads -S 將Paired-end reads作為single-end reads處理。 $ samtools input.sorted.bam output.bam
sam 和 bam 格式轉換
BAM轉換為SAM
samtools view -h -o out.sam out.bam
SAM轉換為BAM
samtools view -bS out.sam >out.bam
-b 意思使輸出使BAM format
-S 意思使輸入使SAM，如果@SQ 缺剩， 要寫-t
所以如果沒有@SQ
samtools faidx ref.fa
samtools view -bt ref.fa.fai out.sam > out.bam
REF:
http://www.chenlianfu.com/?p=1399
http://en.wikipedia.org/wiki/SAMtools
http://www.htslib.org/doc/samtools-1.1.html
http://sourceforge.net/projects/samtools/files/
http://samtools.sourceforge.net/samtools.shtml