溫馨提示×
您好，登錄后才能下訂單哦！
點擊 登錄注冊 即表示同意《億速云用戶服務條款》
您好，登錄后才能下訂單哦！
這篇“samtools怎么使用”文章的知識點大部分人都不太理解，所以小編給大家總結了以下內容，內容詳細，步驟清晰，具有一定的借鑒價值，希望大家閱讀完這篇文章能有所收獲，下面我們一起來看看這篇“samtools怎么使用”文章吧。
$samtoolsProgram: samtools (Tools for alignments in the SAM format) Version: 1.9 (using htslib 1.9) Usage: samtools <command> [options] Commands: -- Indexing dict create a sequence dictionary file faidx index/extract FASTA fqidx index/extract FASTQ index index alignment -- Editing calmd recalculate MD/NM tags and '=' bases fixmate fix mate information reheader replace BAM header targetcut cut fosmid regions (for fosmid pool only) addreplacerg adds or replaces RG tags markdup mark duplicates -- File operations collate shuffle and group alignments by name cat concatenate BAMs merge merge sorted alignments mpileup multi-way pileup sort sort alignment file split splits a file by read group quickcheck quickly check if SAM/BAM/CRAM file appears intact fastq converts a BAM to a FASTQ fasta converts a BAM to a FASTA -- Statistics bedcov read depth per BED region depth compute the depth flagstat simple stats idxstats BAM index stats phase phase heterozygotes stats generate stats (former bamcheck) -- Viewing flags explain BAM flags tview text alignment viewer view SAM<->BAM<->CRAM conversion depad convert padded BAM to unpadded BAM
從上面我們可以看到，大致我5類命令塊：Indexing,Editing,File operations,Statistics,Viewing，下面我們來看看幾個常用的命令
view命令的主要功能是：將sam文件與bam文件互換；然后對bam文件進行各種操作，比如數據的排序(sort)和提取(這些操作是對bam文件進行的，因而當輸入為sam文件的時候，不能進行該操作)；最后將排序或提取得到的數據輸出為bam或sam（默認的）格式。
bam文件優點：bam文件為二進制文件，占用的磁盤空間比sam文本文件小；利用bam二進制文件的運算速度快。
view命令中，對sam文件頭部（序列ID）的輸入(-t或-T）和輸出(-h)是單獨的一些參數來控制的。
Usage: samtools view [options] | [region1 [...]]
下面的view命令的部分參數
默認情況下不加 region，則是輸出所有的 region.options:
-b output BAM # 該參數設置輸出 BAM 格式，默認下輸出是 SAM 格式文件-h print header for the SAM output # 默認下輸出的 sam 格式文件不帶 header，該參數設定輸出sam文件時帶 header 信息 -H print SAM header only (no alignments) # 僅僅輸出文件的頭文件-S input is SAM # 默認下輸入是 BAM 文件，若是輸入是 SAM 文件，則最好加該參數，否則有時候會報錯。 -u uncompressed BAM output (force -b) # 該參數的使用需要有-b參數，能節約時間，但是需要更多磁盤空間。 -c print only the count of matching records# 僅輸出匹配的統計記錄 -L FILE only include reads overlapping this BED FILE [null] # 僅包括和bed文件存在overlap的reads-o FILE output file name [stdout] # 輸出文件的名稱-F INT only include reads with none of the FLAGS in INT present [0] # 過濾flag，僅輸出指定FLAG值的序列-q INT only include reads with mapping quality >= INT [0] # 比對的最低質量值，一般認為20就為unique比對了，可以結合上述-bF參數使用使用提取特定的比對結果-@ Number of additional threads to use [0]# 指使用的線程數
下面來看幾個例子，如果想要比較直觀的結果，可以用軟件自帶的測試數據，在example文件夾中
# 將sam文件轉換成bam文件samtools view -bS abc.sam > abc.bam# BAM轉換為SAMsamtools view -h -o out.sam out.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 [option] <in.bam> -o <out.prefix> -n Sort by read name#設定排序方式按short reads的ID排序。默認下是按序列在fasta文件中的順序（即header）和序列從左往右的位點排序。-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]# 設置每個線程的最大內存，單位可以是K/M/G，默認是 768M。對于處理大數據時，如果內存夠用，則設置大點的值，以節約時間。-t TAG Sort by value of TAG. Uses position as secondary index (or read name if -n is set)# 按照TAG值排序-o FILE Write final output to FILE rather than standard output # 輸出到文件中，加文件名
例子：
# tmp.bam 按照序列位置排序，并將結果輸出到tmp.sort.bam samtools sort -n tmp.bam -o tmp.sort.bam samtools view tmp.sort.bam
merge將多個已經sort了的bam文件融合成一個bam文件。融合后的文件不需要則是已經sort過了的。而cat命令不需要將bam文件進行sort。
Usage: samtools merge [-nurlf] [-h inh.sam] [-b <bamlist.fofn>] <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>] Options: -n Input files are sorted by read name# 輸入文件是經過sort -n的 -t TAG Input files are sorted by TAG value# 輸入文件是經過sort -t的 -r Attach RG tag (inferred from file names)# 添加上RG標簽 -u Uncompressed BAM output# 輸出未壓縮的bam -f Overwrite the output BAM if exist# 覆蓋已經存在的bam -1 Compress level 1# 1倍壓縮 -l INT Compression level, from 0 to 9 [-1]# 指定壓縮倍數 -R STR Merge file in the specified region STR [all] -h FILE Copy the header in FILE to <out.bam> [in1.bam]
$samtools catUsage: samtools cat [options] <in1.bam> [... <inN.bam>] samtools cat [options] <in1.cram> [... <inN.cram>] Options: -b FILE list of input BAM/CRAM file names, one per line -h FILE copy the header from FILE [default is 1st input file] -o FILE output BAM/CRAM
對排序后的序列建立索引，并輸出為bai文件，用于快速隨機處理。在很多情況下，特別是需要顯示比對序列的時候，bai文件是必不可少的，例如之后的tview命令。
Usage: samtools index <in.bam> [out.index]samtools index abc.sort.bam
對fasta文件建立索引,生成的索引文件以.fai后綴結尾。該命令也能依據索引文件快速提取fasta文件中的某一條（子）序列
Usage: samtools faidx <file.fa|file.fa.gz> [<reg> [...]]samtools faidx <file.fa|file.fa.gz> [<reg> [...]]
如對基因組文件建立索引
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比對基因組的情況，和基因組瀏覽器有點類似。可視化一般用IGV比較好，不建議用tview
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 [options] <in.bam> samtools flagstat tmp.bam 20000 + 0 in total (QC-passed reads + QC-failed reads)# 總共的reads數0 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 18995 + 0 mapped (94.98% : N/A)# 總體上reads的匹配率20000 + 0 paired in sequencing# 有多少reads是屬于paired reads10000 + 0 read1# reads1中的reads數10000 + 0 read2# reads2中的reads數18332 + 0 properly paired (91.66% : N/A)# 完美匹配的reads數和比例：比對到同一條參考序列，并且兩條reads之間的距離符合設置的閾值18416 + 0 with itself and mate mapped# paired reads中兩條都比對到參考序列上的reads數579 + 0 singletons (2.90% : N/A)# 單獨一條匹配到參考序列上的reads數，和上一個相加，則是總的匹配上的reads數。0 + 0 with mate mapped to a different chr# paired reads中兩條分別比對到兩條不同的參考序列的reads數0 + 0 with mate mapped to a different chr (mapQ>=5)# 同上一個，只是其中比對質量>=5的reads的數量
得到每個堿基位點的測序深度，并輸出到標準輸出。輸入的bam文件必須先做samtools index
Usage: samtools depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] <in1.bam> [...] -r <chr:from-to> region# 后面跟染色體號（region）-a output all positions (including zero depth)# 輸入所有位置的序列，包括測序深度為0的-q <int> base quality threshold [0]# 堿基質量閾值-Q <int> mapping quality threshold [0]# 比對的質量閾值
舉例：
samtools depth tmp.index.bam > tmp.depth.bam
reheader 替換bam文件的頭
$ samtools reheader <in.header.sam> <in.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.tgztar zxf bam2fastq-1.1.0.tgz cd bam2fastq-1.1.0make ./bam2fastq <in.bam>
samtools還有個非常重要的命令mpileup，以前為pileup。該命令用于生成bcf文件，再使用bcftools進行SNP和Indel的分析。bcftools是samtool中附帶的軟件，在samtools的安裝文件夾中可以找到。
用法：
Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
最常用的參數有2：
-f 來輸入有索引文件的fasta參考序列；
-g 輸出到bcf格式。用法和最簡單的例子如下
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$,...,..... !0000000000scaffold_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]+’代表在該位點后缺失的堿基；
NGS上機測序前需要進行PCR一步，使一個模板擴增出一簇，從而在上機測序的時候表現出為1個點，即一個reads。若一個模板擴增出了多簇，結 果得到了多個reads，這些reads的坐標(coordinates)是相近的。在進行了reads比對后需要將這些由PCR duplicates獲得的reads去掉，并只保留最高比對質量的read。使用rmdup命令即可完成.
Usage: samtools rmdup [-sS] -s rmdup for SE reads# 對single-end reads。默認情況下，只對paired-end reads-S treat PE reads as SE in rmdup (force -s)# 將Paired-end reads作為single-end reads處理。
以上就是關于“samtools怎么使用”這篇文章的內容，相信大家都有了一定的了解，希望小編分享的內容對大家有幫助，若想了解更多相關的知識內容，請關注億速云行業資訊頻道。
免責聲明：本站發布的內容（圖片、視頻和文字）以原創、轉載和分享為主，文章觀點不代表本網站立場，如果涉及侵權請聯系站長郵箱：is@yisu.com進行舉報，并提供相關證據，一經查實，將立刻刪除涉嫌侵權內容。
