samtools使用總結

1. view

view命令的主要功能是:將sam文件轉換成bam文件;而後對bam文件進行各類操做,好比數據的排序(不屬於本命令的功能)和提取(這些操做是對bam文件進行的,於是當輸入爲sam文件的時候,不能進行該操做);最後將排序或提取獲得的數據輸出爲bam或sam(默認的)格式。瀏覽器

bam文件優勢:bam文件爲二進制文件,佔用的磁盤空間比sam文本文件小;利用bam二進制文件的運算速度快。app

標準sam/bam 文件less


1.read名稱,一般包括測序平臺等信息
2.SAM標記(Flag),沒有mapping的標記爲「 * 」
3.chromosome
4.比對上的位置,注意是從1開始計數。
5.MAPQ(mapping quality,描述比對的質量,數字越大,特異性越高,說明該read比對到參考基因組上的位置越惟一)
6.CIGAR字串,記錄插入,刪除,錯配以及splice junctions(後剪切拼接的接頭)
7.mate名稱,記錄mate pair信息
8.mate的位置
9.模板的長度
10.read序列
11.read質量
12.程序用標記ide

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)和序列從左往右的位點排序。

例子:ui

$ samtools sort abc.bam abc.sort
$ samtools view abc.sort.bam | less -S

3.merge

將2個或2個以上的已經sort了的bam文件融合成一個bam文件。融合後的文件不須要則是已經sort過了的。this

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。不然會報錯。spa

創建索引後將產生後綴爲.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的數量

 

 

Flag 標籤說明
在sam/bam文件中輸出的比對falg信息,是這條reads符合含義的十進制數值之和,網上也有許多在線工具能夠經過輸入結果flag值來標識出其對應的全部含義。

 

 

7. depth

獲得每一個鹼基位點的測序深度,並輸出到標準輸出。

Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] <in1.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
相關文章
相關標籤/搜索