|
samtools的說明文檔：http://samtools./samtools.shtml 1. viewview命令的主要功能是：將sam文件轉(zhuǎn)換成bam文件,；然后對bam文件進行各種操作,，比如數(shù)據(jù)的排序(不屬于本命令的功能)和提取(這些操作 是對bam文件進行的，因而當(dāng)輸入為sam文件的時候,，不能進行該操作),；最后將排序或提取得到的數(shù)據(jù)輸出為bam或sam（默認的）格式。 bam文件優(yōu)點：bam文件為二進制文件,，占用的磁盤空間比sam文本文件?。焕胋am二進制文件的運算速度快,。 view命令中,，對sam文件頭部的輸入(-t或-T）和輸出(-h)是單獨的一些參數(shù)來控制的。 Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]
默認情況下不加 region,，則是輸出所有的 region.
Options: -b output BAM
默認下輸出是 SAM 格式文件,，該參數(shù)設(shè)置輸出 BAM 格式
-h print header for the SAM output
默認下輸出的 sam 格式文件不帶 header，該參數(shù)設(shè)定輸出sam文件時帶 header 信息
-H print header only (no alignments)
-S input is SAM
默認下輸入是 BAM 文件,，若是輸入是 SAM 文件,，則最好加該參數(shù)，否則有時候會報錯,。
-u uncompressed BAM output (force -b)
該參數(shù)的使用需要有-b參數(shù),，能節(jié)約時間，但是需要更多磁盤空間,。
-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]
數(shù)字4代表該序列沒有比對到參考序列上
數(shù)字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文件轉(zhuǎn)換成bam文件
$ samtools view -bS abc.sam > abc.bam
$ samtools view -b -S abc.sam -o abc.bam
提取比對到參考序列上的比對結(jié)果
$ samtools view -bF 4 abc.bam > abc.F.bam
提取paired reads中兩條reads都比對到參考序列上的比對結(jié)果，只需要把兩個4+8的值12作為過濾參數(shù)即可
$ samtools view -bF 12 abc.bam > abc.F12.bam
提取沒有比對到參考序列上的比對結(jié)果
$ samtools view -bf 4 abc.bam > abc.f.bam
提取bam文件中比對到caffold1上的比對結(jié)果，并保存到sam文件格式
$ samtools view abc.bam scaffold1 > scaffold1.sam
提取scaffold1上能比對到30k到100k區(qū)域的比對結(jié)果
$ samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam
根據(jù)fasta文件,，將 header 加入到 sam 或 bam 文件中
$ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam
2. sortsort對bam文件進行排序,。 Usage: samtools sort [-n] [-m <maxMem>] <in.bam> <out.prefix> -m 參數(shù)默認下是 500,000,000 即500M（不支持K，M,，G等縮寫）,。對于處理大數(shù)據(jù)時，如果內(nèi)存夠用,，則設(shè)置大點的值,，以節(jié)約時間。 -n 設(shè)定排序方式按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個以上的已經(jīng)sort了的bam文件融合成一個bam文件,。融合后的文件不需要則是已經(jīng)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,。否則會報錯,。 建立索引后將產(chǎn)生后綴為.bai的文件，用于快速的隨機處理,。很多情況下需要有bai文件的存在,，特別是顯示序列比對情況下。比如samtool的tview命令就需要,；gbrowse2顯示reads的比對圖形的時候也需要,。 Usage: samtools index <in.bam> [out.index] 例子： 以下兩種命令結(jié)果一樣
$ samtools index abc.sort.bam
$ samtools index abc.sort.bam abc.sort.bam.bai
5. faidx對fasta文件建立索引,生成的索引文件以.fai后綴結(jié)尾。該命令也能依據(jù)索引文件快速提取fasta文件中的某一條（子）序列 Usage: samtools faidx <in.bam> [ [...]] 對基因組文件建立索引 $ samtools faidx genome.fasta 生成了索引文件genome.fasta.fai,是一個文本文件,，分成了5列,。第一列是子序列的名稱； 第二列是子序列的長度,；個人認為“第三列是序列所在的位置”,，因為該數(shù)字從上往下逐漸變大， 最后的數(shù)字是genome.fasta文件的大??；第4和5列不知是啥意思。于是通過此文件,，可以定 位子序列在fasta文件在磁盤上的存放位置,，直接快速調(diào)出子序列。 由于有索引文件,，可以使用以下命令很快從基因組中提取到fasta格式的子序列 $ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta 6. tviewtview能直觀的顯示出reads比對基因組的情況,，和基因組瀏覽器有點類似,。 Usage: samtools tview <aln.bam> [ref.fasta] 當(dāng)給出參考基因組的時候，會在第一排顯示參考基因組的序列,，否則,，第一排全用N表示。 按下 g ,，則提示輸入要到達基因組的某一個位點,。例子“scaffold_10:1000"表示到達第 10號scaffold的第1000個堿基位點處。 使用H(左）J（上）K（下）L（右）移動顯示界面,。大寫字母移動快,，小寫字母移動慢。 使用空格建向左快速移動（和 L 類似）,，使用Backspace鍵向左快速移動（和 H 類似）,。 Ctrl+H 向左移動1kb堿基距離； Ctrl+L 向右移動1kb堿基距離 可以用顏色標(biāo)注比對質(zhì)量,，堿基質(zhì)量,，核苷酸等。30～40的堿基質(zhì)量或比對質(zhì)量使用白色表示,； 20～30黃色,；10～20綠色；0～10藍色,。 使用點號'.'切換顯示堿基和點號,；使用r切換顯示read name等 還有很多其它的使用說明，具體按 ,？ 鍵來查看,。 7. flagstat給出BAM文件的比對結(jié)果 Usage: samtools flagstat <in.bam> $ samtools flagstat example.bam 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ù)：比對到同一條參考序列，并且兩條reads之間的距離符合設(shè)置的閾值 6899708 + 0 with itself and mate mapped #paired reads中兩條都比對到參考序列上的reads數(shù) 636656 + 0 singletons (5.33%:-nan%) #單獨一條匹配到參考序列上的reads數(shù),，和上一個相加,，則是總的匹配上的reads數(shù)。 469868 + 0 with mate mapped to a different chr #paired reads中兩條分別比對到兩條不同的參考序列的reads數(shù) 243047 + 0 with mate mapped to a different chr (mapQ>=5) #同上一個,，只是其中比對質(zhì)量>=5的reads的數(shù)量 7. depth得到每個堿基位點的測序深度,并輸出到標(biāo)準(zhǔn)輸出,。 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 統(tǒng)計一個表格，4列,，分別為”序列名,，序列長度，比對上的reads數(shù),，unmapped reads number”,。第4列應(yīng)該是paired reads中有一端能匹配到該scaffold上,，而另外一端不匹配到任何scaffolds上的reads數(shù)。 $ samtools idxstats <aln.bam> 9. 將bam文件轉(zhuǎn)換為fastq文件有時候,，我們需要提取出比對到一段參考序列的reads,，進行小范圍的分析,，以利于debug等,。這時需要將bam或sam文件轉(zhuǎn)換為fastq格式。 $ wget http://www./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. mpileupsamtools還有個非常重要的命令mpileup,，以前為pileup,。該命令用于生成bcf文件，再使用bcftools進行SNP和Indel的分析,。bcftools是samtool中附帶的軟件,，在samtools的安裝文件夾中可以找到。 最常用的參數(shù)有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參數(shù)時，則不生成二進制的bcf文件,，而生成一個文本文件(輸出到標(biāo)準(zhǔn)輸出),。該文本文件統(tǒng)計了參考序列中每個堿基位點的比對情況；該文件每一行代表了參考序列中某一個堿基位點的比對結(jié)果,。比如： 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生成的結(jié)果包含6行：參考序列名,；位置；參考堿基,；比對上的reads數(shù),；比對情況；比對上的堿基的質(zhì)量,。其中第5列比較復(fù)雜,解釋如下： 1 ‘.’代表與參考序列正鏈匹配,。 pileup具體的參數(shù)如下： 輸入?yún)?shù)
-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]
輸出參數(shù)
-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. 使用bcftoolsbcftools和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
生成的結(jié)果文件為vcf格式，有10列,，分別是：1 參考序列名,；2 varianti所在的left-most位置；3 variant的ID（默認未設(shè)置,，用’.'表示),；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 的具體參數(shù)如下： 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]
當(dāng)有多個sample用于variants calling時,，比如多個轉(zhuǎn)錄組數(shù)據(jù)或多個重測序
數(shù)據(jù)需要比對到參考基因組上，設(shè)置該值,，表明至少有該<float 0~1>比例的
samples在該位點都有覆蓋才計算入variant.所以對于只有一個sample的情況
下,，該值設(shè)置在0～1之間沒有意義，大于1則得不到任何結(jié)果,。
-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結(jié)果后,。需要對結(jié)果再次進行過濾。主要依據(jù)比對結(jié)果中第8列信息,。其中的 DP4 一行尤為重要,，提供了4個數(shù)據(jù)：1 比對結(jié)果和正鏈一致的reads數(shù)、2 比對結(jié)果和負鏈一致的reads數(shù),、3 比對結(jié)果在正鏈的variant上的reads數(shù),、4 比對結(jié)果在負鏈的variant上的reads數(shù)?？梢栽O(shè)定 （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 rmdupNGS上機測序前需要進行PCR一步，使一個模板擴增出一簇,，從而在上機測序的時候表現(xiàn)出為1個點,，即一個reads。若一個模板擴增出了多簇,，結(jié) 果得到了多個reads,，這些reads的坐標(biāo)(coordinates)是相近的。在進行了reads比對后需要將這些由PCR duplicates獲得的reads去掉,，并只保留最高比對質(zhì)量的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 |
|
|
