生信人必會之 samtools
慶祝香港回歸20周年
在這個重要的時刻，讓我為你唱一首歌~~
啊~~努力做好生信人，
為了祖國的統一事業做出大貢獻。
怎麼做？
Jimmy大哥說了學生信人必會之samtools,
於是乎，
魚同志響應號召，給我分享了一些資料。
下面是那首動聽的歌~~~~~~~~~~
samtools
samtools 常用命令詳解
希望自己有朝一日，技能上這樣虐群主嗎？~想來群主也是心甚慰！好好學習唄！
就是瞅！ 將sam文件瞅成bam文件；這是對sam文件最早的操作，然後對bam文件進行各種操作，比如數據的排序(view不行，是別的命令，光瞅能幹啥！)和提取(這些是對bam文件的操作，sam不行，就是不願意)；
瞅你咋的？
最後將排序或提取得到的數據輸出為bam或sam（默認的）格式。
bam文件好在它小啊，小就佔地兒少啊，小就靈活，速度就快啊。
view命令中，對sam文件頭部的輸入(-t或-T）和輸出(-h)是單獨的一些參數來控制的。
慄子，慄子，小寶慄子
#將sam文件轉換成bam文件$ samtools view -bS abc.sam > abc.bam$ samtools view -b -S abc.sam -o abc.bam
Usage: samtools sort [-n] [-m <maxMem>] <in.bam> <out.prefix>
-m 參數默認下是 500,000,000 即500M（不支持K，M，G等縮寫）。對於處理大數據時，如果內存夠用，則設置大點的值，以節約時間。
排隊的時候你得設置人數吧，排到大門外你就不知道有多少人了。所以得設置上限，所以出現了-m 先生，他就是上限（任性，旁邊標他多重，）默認下默認下是 500,000,000 即500M（不支持K，M，G等縮寫）。真重！如果地方夠大，多整點唄！則設置大點的值，以節約時間。
那你問按啥排？-n 先生說按short reads的ID排序唄（你帥你說了算），然後全村的人都知道默認下是按序列在fasta文件中的順序（即header）和序列從左往右的位點排序。
下面是深夜食堂！
兩個慄子（減肥少吃點）
$ samtools sort abc.bam abc.sort$ samtools view abc.sort.bam | less -S
上面排好序，站好隊了，哎（這兩瘦的一樣高，合成一個唄，省地兒），這就有了merge(純屬胡說八道！)！
官方如右：將2個或2個以上的已經sort了的bam文件融合成一個bam文件。融合後的文件不需要則是已經sort過了的。
01Usage: samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam>[...]03Options: -n sort by read names04 -r attach RG tag (inferred from file names)05 -u uncompressed BAM output06 -f overwrite the output BAM if exist07 -1 compress level 108 -R STR merge file in the specified region STR [all]09 -h FILE copy the header in FILE to <out.bam> [in1.bam]11Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users12 must provide the correct header with -h, or uses Picard which properly maintains13 the header dictionary in merging.排好隊後，需要給大家一個編號，就是index,
建立索引後將產生後綴為.bai的文件，用於快速的隨機處理。很多情況下需要有bai文件的存在，特別是顯示序列比對情況下。比如samtool的tview命令就需要；gbrowse2顯示reads的比對圖形的時候也需要。
Usage: samtools index <in.bam> [out.index]
下面是慄子
（媽媽說多吃慄子身體好）
#以下兩種命令結果一樣
$ samtools index abc.sort.bam</span>
看名字就知道是對fasta文件建立索引,生成的索引文件以.fai後綴結尾。該命令也能依據索引文件快速提取fasta文件中的某一條（子）序列。
這就是傳說中的花裡胡哨本人！
tview能直觀的顯示出reads比對基因組的情況，和基因組瀏覽器(IGV,IGB)有點類似。
當給出參考基因組的時候，會在第一排顯示參考基因組的序列，否則，第一排全用N表示。
10號scaffold的第1000個鹼基位點處。
使用H(左）J（上）K（下）L（右）移動顯示界面。大寫字母移動快，小寫字母移動慢。
使用空格建向左快速移動（和 L 類似），使用Backspace鍵向左快速移動（和 H 類似）。
Ctrl+H 向左移動1kb鹼基距離； Ctrl+L 向右移動1kb鹼基距離
20～30黃色；10～20綠色；0～10藍色。
使用點號'.'切換顯示鹼基和點號；使用r切換顯示read name等
PS:建議不要學這個命令了，還是直接用IGV看吧！
扔給你BAM比對結果！！
243047 + 0 with mate mapped to a different chr (mapQ>=5)
#同上一個，只是其中比對質量>=5的reads的數量
PS: 此處應該強勢推薦一個兄弟命令，stats，不僅僅各種各樣的bam統計結果，還可以出圖！！！
一看就是測序深度，就是得到每個鹼基位點的測序深度,並輸出到標準輸出。
PS:這個命令非常重要，尤其是需要做各種各樣的個性化統計的時候。
reheader 替換bam文件的頭（好可怕，換腦袋）
查看原始碼列印幫助
1$ samtools reheader <in.header.sam> <in.bam>cat 連接多個bam文件，適用於非sorted的bam文件
1$ samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [ ... ]
idxstats 統計一個表格，4列，分別為」序列名，序列長度，比對上的reads數，unmapped reads number」。第4列應該是paired reads中有一端能匹配到該scaffold上，而另外一端不匹配到任何scaffolds上的reads數。
存在即合理嘛~
王菲唱：有時候，有時候~我們需要提取出比對到一段參考序列的reads，進行小範圍的分析，以利於debug等。這時需要將bam或sam文件轉換為fastq格式。
該網站提供了一個bam轉換為fastq的程序：http://www.hudsonalpha.org/gsl/information/software/bam2fastq
mpileup，嗯！這是個重要的命令。（我在Jimmy的基因組直播重看到過~~很有趣，詳情關注《生信技能樹》公眾號哦，真的很有趣，看完自己都想做一個測序。快去看12那節，好奇寶寶表示想看哪個出軌基因的結果（壞笑~），求群主原諒）
該命令用於生成bcf文件，再使用bcftools進行SNP和Indel的分析。bcftools是samtool中附帶的軟體，在samtools的安裝文件夾中可以找到。
最常用的參數有2（就不想說兩）：
-f 來輸入有索引文件的fasta參考序列；
-g 輸出到bcf格式。用法和最簡單的例子如下
1Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]3$ samtools mpileup -f genome.fasta abc.bam > abc.txt4$ samtools mpileup -gSDf genome.fasta abc.bam > abc.bcf5$ samtools mpileup -guSDf genome.fasta abc.bam | \6 bcftools view -cvNg - > abc.vcfmpileup不使用-u或-g參數時，則不生成二進位的bcf文件，而生成一個文本文件(輸出到標準輸出)。
該文本文件統計了參考序列中每個鹼基位點的比對情況；該文件每一行代表了參考序列中某一個鹼基位點的比對結果。比如：
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具體的參數如下：
查看原始碼列印幫助
02-6 Assume the quality is in the Illumina 1.3+ encoding. -A Do not skip anomalous read pairs in variant calling.03-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.04-b FILE List of input BAM files, one file per line [null]05-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]06-d INT At a position, read maximally INT reads per input BAM. [250]07-E Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt specificity a little bit.08-f FILE The faidx-indexed reference file in the FASTA format. The file can be optionally compressed by razip. [null]09-l FILE BED or position list file containing a list of regions or sites where pileup or BCF should be generated [null]10-M INT cap mapping quality at INT [60]11-q INT Minimum mapping quality for an alignment to be used [0]12-Q INT Minimum base quality for a base to be considered [13]13-r STR Only generate pileup in region STR [all sites]16-D Output per-sample read depth (require -g/-u)17-g Compute genotype likelihoods and output them in the binary call format (BCF).18-S Output per-sample Phred-scaled strand bias P-value (require -g/-u)19-u Similar to -g except that the output is uncompressed BCF, which is preferred for piping.21Options for Genotype Likelihood Computation (for -g or -u):22-e INT Phred-scaled gap extension sequencing error probability. Reducing INT leads to longer indels. [20]23-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]24-I Do not perform INDEL calling25-L INT Skip INDEL calling if the average per-sample depth is above INT. [250]26-o INT Phred-scaled gap open sequencing error probability. Reducing INT leads to more indel calls. [40]27-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]這種幹不拉幾的，吃完會噎死！！
bcftools和samtools類似，用於處理vcf(variant call format)文件和bcf(binary call format)文件。前者為文本文件，後者為其二進位文件。
bcftools使用簡單，最主要的命令是view命令，其次還有index和cat等命令。index和cat命令和samtools中類似。此處主講使用view命令來進行SNP和Indel calling。該命令的使用方法和例子為：
生成的結果文件為vcf格式，有10列，分別是：
1 參考序列名；
2 variant所在的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)。
例如：
1scaffold_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:632scaffold_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:693scaffold_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:214scaffold_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:455scaffold_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:996scaffold_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:37scaffold_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:168scaffold_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 的描述如下：
bcftools view 的具體參數如下：
02-A Retain all possible alternate alleles at variant sites. By default, the view command discards unlikely alleles.03-b Output in the BCF format. The default is VCF.04-D FILE Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]05-F Indicate PL is generated by r921 or before (ordering is different).06-G Suppress all individual genotype information.07-l FILE List of sites at which information are outputted [all sites]08-N Skip sites where the REF field is not A/C/G/T09-Q Output the QCALL likelihood format10-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]11-S The input is VCF instead of BCF.12-u Uncompressed BCF output (force -b).14Consensus/Variant Calling Options:15-c Call variants using Bayesian inference. This option automatically invokes option -e.16-d FLOAT When -v is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]17 當有多個sample用於variants calling時，比如多個轉錄組數據或多個重測序18 數據需要比對到參考基因組上，設置該值，表明至少有該<float 0~1>比例的19 samples在該位點都有覆蓋才計算入variant.所以對於只有一個sample的情況20 下，該值設置在0～1之間沒有意義，大於1則得不到任何結果。21-e Perform max-likelihood inference only, including estimating the site allele frequency, testing Hardy-Weinberg equlibrium and testing associations with LRT.22-g Call per-sample genotypes at variant sites (force -c)23-i FLOAT Ratio of INDEL-to-SNP mutation rate [0.15]24-p FLOAT A site is considered to be a variant if P(ref|D)25-t FLOAT Scaled muttion rate for variant calling [0.001]26-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]27-v Output variant sites only (force -c)29Contrast Calling and Association Test Options:30-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]31-U INT Number of permutations for association test (effective only with -1) [0]32-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。比如：
NGS上機測序前需要進行PCR一步，使一個模板擴增出一簇，從而在上機測序的時候表現出為1個點，即一個reads。若一個模板擴增出了多簇，結果得到了多個reads，這些reads的坐標(coordinates)是相近的。在進行了reads比對後需要將這些由PCR duplicates獲得的reads去掉，並只保留最高比對質量的read。使用rmdup命令即可完成.
後續還有精彩~~~~