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
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
將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.
必須對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
對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
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等 還有不少其它的使用說明,具體按 ? 鍵來查看。
給出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值來標識出其對應的全部含義。
獲得每一個鹼基位點的測序深度,並輸出到標準輸出。
Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] <in1.bam> [...]
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>
有時候,咱們須要提取出比對到一段參考序列的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>
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]
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
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