記錄一下我用到的samtools的用法。
samtools的說明文檔：http://samtools.sourceforge.net/samtools.shtml
bam文件優(yōu)點(diǎn)：bam文件為二進(jìn)制文件，占用的磁盤空間比sam文本文件小；利用bam二進(jìn)制文件的運(yùn)算速度快覆获。
首先需要意識(shí)到的是samtools是一個(gè)非常強(qiáng)大的工具，想要熟練的使用它，還需要不斷的摸索小压。
samtools的用法
（1）View
samtools view -bS abc.sam > abc.bam? ? #將sam文件轉(zhuǎn)換為bam文件
?參數(shù)：
-b bam 輸出bam
-S sam 輸入sam
-@ 線程
在比對(duì)完成的sam文件中，包含著mapped reads 和unmapped reads
$ samtools view -bF 4? abc.bam > abc.F.bam? ? ? ?#提取沒有比對(duì)到參考序列上的比對(duì)結(jié)果椰于，步包含標(biāo)簽
$ samtools view -bF 12 abc.bam > abc.F12.bam? ?#提取paired reads中兩條reads都比對(duì)到參考序列上的比對(duì)結(jié)果怠益，只需要把兩個(gè)4+8的值12作為過濾參數(shù)即可
$ samtools view -bf 4 abc.bam > abc.f.bam? ? #提取沒有比對(duì)到參考序列上的比對(duì)結(jié)果，包含標(biāo)簽
$ samtools view abc.bam scaffold1 > scaffold1.sam? ? ?#提取bam文件中比對(duì)到caffold1上的比對(duì)結(jié)果瘾婿，并保存到sam文件格式
$ samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam? ? #提取scaffold1上能比對(duì)到30k到100k區(qū)域的比對(duì)結(jié)果
$ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam? ? #根據(jù)fasta文件蜻牢，將 header 加入到 sam 或 bam 文件中
samtools的view不就可以進(jìn)行格式轉(zhuǎn)換烤咧，還可以進(jìn)行數(shù)據(jù)的提取
例：提取1號(hào)染色體上1234~123456區(qū)域的以對(duì)read
samtools view SRR3589957_sorted.bam chr1:1234-123456| head
?samtools view SRR3589957_sorted.bam chr1:1234-123456 > sub.bam?
使FLAG更具可讀性
samtools view -X sample.sorted.bam | head -n 5
計(jì)算總的比對(duì)數(shù)量
samtools view sample.sorted.bam | wc -l
顯示標(biāo)題，-H選項(xiàng)
samtools view -H sample.sorted.bam
將bam文件轉(zhuǎn)換為sam文件
samtools view -h abc.bam > abc.sam
（2）Sort
samtools sort對(duì)bam文件進(jìn)行排序抢呆，不能對(duì)sam文件進(jìn)行排序髓削。
以leftmost coordinates的方式對(duì)比對(duì)結(jié)果進(jìn)行排序，或者使用-n參數(shù)以read名稱進(jìn)行排序镀娶。將會(huì)添加適當(dāng)?shù)腀HD-SO排序順序標(biāo)頭標(biāo)簽或者如果有必要的話立膛，將會(huì)更新現(xiàn)存的一個(gè)排序順序標(biāo)頭標(biāo)簽。sort命令的輸出默認(rèn)是標(biāo)準(zhǔn)輸出寫入梯码，或者使用-o參數(shù)時(shí)宝泵，指定bam文件輸出名。sort命令還會(huì)在內(nèi)存不足時(shí)創(chuàng)建臨時(shí)文件tmpprefix.%d.bam轩娶。
也就是說：samtools的排序方式有兩種（常用）
默認(rèn)方式儿奶，按照染色體的位置進(jìn)行排序
samtools sort test.bam default
參數(shù)-n則是根據(jù)read名進(jìn)行排序。
samtools sort -n test.bam sort_left
usage: samtools sort [-l?level] [-m?maxMem] [-o?out.bam] [-O?format] [-n] [-T?tmpprefix] [-@?threads] [in.sam|in.bam|in.cram]
例如：samtools sort abc.bam abc.sort
samtools sort -O bam -@ 2 SRR1909070.bam -o SRR1909070.sorted.bam
RNA-seq 的數(shù)據(jù)比對(duì)結(jié)果 BAM 文件使用 samtools 進(jìn)行 sort 之后文件壓縮比例變化會(huì)比DNA-seq 更甚鳄抒。另外闯捎，samtools 對(duì) BAM 文件進(jìn)行排序之后那些沒有比對(duì)上的 reads 會(huì)被放在文件的末尾。
?參數(shù)：
-l INT 設(shè)置輸出文件壓縮等級(jí)许溅。0-9瓤鼻，0是不壓縮，9是壓縮等級(jí)最高贤重。不設(shè)置此參數(shù)時(shí)茬祷，使用默認(rèn)壓縮等級(jí)；
-m INT 設(shè)置每個(gè)線程運(yùn)行時(shí)的內(nèi)存大小并蝗，可以使用K祭犯，M和G表示內(nèi)存大小。
-n 設(shè)定排序方式按short reads的ID排序滚停。默認(rèn)下是按序列在fasta文件中的順序（即header）和序列從左往右的位點(diǎn)排序沃粗。
-o FILE 設(shè)置最終排序后的輸出文件名；
-O FORMAT 設(shè)置最終輸出的文件格式键畴，可以是bam最盅，sam或者cram，默認(rèn)為bam镰吵；
-T PREFIX 設(shè)置臨時(shí)文件的前綴檩禾；
-@ INT 設(shè)置排序和壓縮是的線程數(shù)量，默認(rèn)是單線程疤祭。
（3）index
samtools index 建立索引盼产，在建立索引之前應(yīng)該先對(duì)bam文件進(jìn)行排序。必須對(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í)候也需要些己。
samtools index abc.sort.bam
如果想要建立索引的，具體可以看看比對(duì)的內(nèi)部的算法嘿般，鏈接具體是怎么建立索引的
建立索引的目的應(yīng)該是為了提高比對(duì)的效率
以下兩種命令結(jié)果一樣
$ samtools index abc.sort.bam
$ samtools index abc.sort.bam abc.sort.bam.bai
（4）flagstat
samtools flagstat? 給出BAM文件的比對(duì)結(jié)果
samtools flagstat [options] <in.bam>
-@ 線程
-O FORMAT 設(shè)置最終輸出的文件格式段标，可以是txt，json或者tsv炉奴，默認(rèn)為json逼庞，tsv；
samtools flagstat輸出結(jié)果解釋:
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)
#paired reads中兩條分別比對(duì)到兩條不同的參考序列的reads數(shù),并且其中比對(duì)質(zhì)量>=5的reads的數(shù)量
（5）depth
得到每個(gè)堿基位點(diǎn)的測(cè)序深度,并輸出到標(biāo)準(zhǔn)輸出璧南。
usage: samtools depth [options] in.bam [in.bam ...]
注意：做depth之前必須做samtools index；
示例：
samtools depth in.bam? >? out.depth.txt
注意： in.bam 必須經(jīng)過了排序痹兜。
（6）samtools rmdup
NGS上機(jī)測(cè)序前需要進(jìn)行PCR一步穆咐，使一個(gè)模板擴(kuò)增出一簇，從而在上機(jī)測(cè)序的時(shí)候表現(xiàn)出為1個(gè)點(diǎn)字旭，即一個(gè)reads。若一個(gè)模板擴(kuò)增出了多簇崖叫，結(jié)果得到了多個(gè)reads遗淳，這些reads的坐標(biāo)(coordinates)是相近的。在進(jìn)行了reads比對(duì)后需要將這些由PCRduplicates獲得的reads去掉心傀，并只保留最高比對(duì)質(zhì)量的read屈暗。使用rmdup命令即可完成.
Usage:
samtools rmdup[-sS]
-s對(duì)single-end reads。默認(rèn)情況下脂男，只對(duì)paired-endreads
-S將Paired-endreads作為single-endreads處理养叛。
$samtools?rmdup input.sorted.bam output.bam
（7）mpileup
samtools還有個(gè)非常重要的命令mpileup，以前為pileup宰翅。該命令用于生成bcf文件弃甥，再使用bcftools進(jìn)行SNP和Indel的分析。bcftools是samtool中附帶的軟件汁讼，在samtools的安裝文件夾中可以找到淆攻。
最常用的參數(shù)有2個(gè)：
-f來輸入有索引文件的fasta參考序列阔墩；
-g輸出到bcf格式。用法和最簡(jiǎn)單的例子如下
Usage:samtoolsmpileup[-EBug][-CcapQcoef][-rreg][-fin.fa][-llist][-McapMapQ][-QminBaseQ][-qminMapQ]in.bam[in2.bam[...]]
$samtoolsmpileup-fgenome.fastaabc.bam>abc.txt
$samtoolsmpileup-gSDfgenome.fastaabc.bam>abc.bcf
$samtoolsmpileup-guSDfgenome.fastaabc.bam|\bcftoolsview-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é)果忘苛。比如：
（8）faidx
對(duì)fasta文件建立索引，比如基因組的文件唱较，生成的索引文件以.fai后綴結(jié)尾扎唾。該命令也能依據(jù)索引文件快速提取fasta文件中的某一條（子）序列
Usage: samtools faidx <in.bam> [ [...]]
對(duì)基因組文件建立索引
$ samtools faidx genome.fasta
生成了索引文件genome.fasta.fai,是一個(gè)文本文件，分成了5列绊汹。
第一列是子序列的名稱稽屏；
第二列是子序列的長(zhǎng)度；
第三列是序列所在的位置西乖，因?yàn)樵摂?shù)字從上往下逐漸變大狐榔，最后的數(shù)字是genome.fasta文件的大小获雕；
第4和5列不知是啥意思薄腻。于是通過此文件，可以定
位子序列在fasta文件在磁盤上的存放位置届案，直接快速調(diào)出子序列庵楷。
由于有索引文件，可以使用以下命令很快從基因組中提取到fasta格式的子序列
$ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta
拓展：bcftools軟件
bcftools和samtools類似楣颠，用于處理vcf(variant call format)文件和bcf(binary call format)文件尽纽。前者為文本文件，后者為其二進(jìn)制文件童漩。
bcftools使用簡(jiǎn)單弄贿，最主要的命令是view命令，其次還有index和cat等命令矫膨。index和cat命令和samtools中類似差凹。此處主講使用view命令來進(jìn)行SNP和Indel calling。該命令的使用方法和例子為：
$ 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)隔開；9 FORMAT of the genotype fields, separated by colon (optional)船逮； 10 SAMPLE genotypes and per-sample information (optional)顾腊。
參考鏈接：
原文鏈接：https://blog.csdn.net/u013553061/article/details/53179945
https://www.cnblogs.com/emanlee/p/4316581.html
http://events.jianshu.io/p/794d82bccf6c
http://blog.sina.com.cn/s/blog_13de3725c0102v7rd.html
https://www.cnblogs.com/shuaihe/articles/6802246.html
