Samtools安裝及常用命令詳解

Samtools是一個用于操作sam和bam文件的工具合集。包含有許多命令。以下是常用命令的介紹

下載安裝包：

http://www.htslib.org/download/

安裝依賴：

yum install bzip2-devel ncurses-libs ncurses-devel xz-devel zlib-devel

編譯安裝Samtools：

tar xvf samtools-1.9.tar.bz2

cd samtools-1.9

./configure --prefix=/opt/samtools1.9

make

make install

配置環(huán)境變量：

gedit ~/.bashrc

#Samtools1.9

export PATH=/opt/samtools1.9/bin:$PATH

source ~/.bashrc

運行

samtools

samtools常用命令詳解

1. view

view命令的主要功能是：將sam文件轉(zhuǎn)換成bam文件；然后對bam文件進行各種操作，比如數(shù)據(jù)的排序(不屬于本命令的功能)和提取(這些操作是對bam文件進行的，因而當(dāng)輸入為sam文件的時候，不能進行該操作)；最后將排序或提取得到的數(shù)據(jù)輸出為bam或sam（默認(rèn)的）格式。

bam文件優(yōu)點：bam文件為二進制文件，占用的磁盤空間比sam文本文件小；利用bam二進制文件的運算速度快。

view命令中，對sam文件頭部的輸入(-t或-T）和輸出(-h)是單獨的一些參數(shù)來控制的。

Usage: samtools view [options] | [region1 [...]]

默認(rèn)情況下不加 region，則是輸出所有的 region.

Options: -b output BAM

默認(rèn)下輸出是 SAM 格式文件，該參數(shù)設(shè)置輸出 BAM 格式

-h print header for the SAM output

默認(rèn)下輸出的 sam 格式文件不帶 header，該參數(shù)設(shè)定輸出sam文件時帶 header 信息

-H print header only (no alignments)

-S input is SAM

默認(rèn)下輸入是 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. sort

sort對bam文件進行排序。

Usage: samtools sort [-n] [-m ]

-m 參數(shù)默認(rèn)下是 500,000,000 即500M（不支持K，M，G等縮寫）。對于處理大數(shù)據(jù)時，如果內(nèi)存夠用，則設(shè)置大點的值，以節(jié)約時間。

-n 設(shè)定排序方式按short reads的ID排序。默認(rèn)下是按序列在fasta文件中的順序（即header）和序列從左往右的位點排序。

3.merge

將2個或2個以上的已經(jīng)sort了的bam文件融合成一個bam文件。融合后的文件不需要則是已經(jīng)sort過了的。

Usage: samtools merge [-nr] [-h inh.sam] [...]

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 [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文件進行默認(rèn)情況下的排序后，才能進行index。否則會報錯。

建立索引后將產(chǎn)生后綴為.bai的文件，用于快速的隨機處理。很多情況下需要有bai文件的存在，特別是顯示序列比對情況下。比如samtool的tview命令就需要；gbrowse2顯示reads的比對圖形的時候也需要。

例子：

以下兩種命令結(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 [ [...]]

對基因組文件建立索引

$ samtools faidx genome.fasta

生成了索引文件genome.fasta.fai,是一個文本文件，分成了5列。第一列是子序列的名稱；

第二列是子序列的長度；個人認(rèn)為“第三列是序列所在的位置”，因為該數(shù)字從上往下逐漸變大，

最后的數(shù)字是genome.fasta文件的大小；第4和5列不知是啥意思。于是通過此文件，可以定

位子序列在fasta文件在磁盤上的存放位置，直接快速調(diào)出子序列。

由于有索引文件，可以使用以下命令很快從基因組中提取到fasta格式的子序列

$ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta

6. tview

tview能直觀的顯示出reads比對基因組的情況，和基因組瀏覽器有點類似。

Usage: samtools tview [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. 將bam文件轉(zhuǎn)換為fastq文件

有時候，我們需要提取出比對到一段參考序列的reads，進行小范圍的分析，以利于debug等。這時需要將bam或sam文件轉(zhuǎn)換為fastq格式。

8. 使用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

參考：

http://www.chenlianfu.com/?p=1399

http://www.bio-info-trainee.com/518.html

版權(quán)聲明：本文內(nèi)容由網(wǎng)絡(luò)用戶投稿，版權(quán)歸原作者所有，本站不擁有其著作權(quán)，亦不承擔(dān)相應(yīng)法律責(zé)任。如果您發(fā)現(xiàn)本站中有涉嫌抄襲或描述失實的內(nèi)容，請聯(lián)系我們jiasou666@gmail.com 處理，核實后本網(wǎng)站將在24小時內(nèi)刪除侵權(quán)內(nèi)容。

版權(quán)聲明：本文內(nèi)容由網(wǎng)絡(luò)用戶投稿，版權(quán)歸原作者所有，本站不擁有其著作權(quán)，亦不承擔(dān)相應(yīng)法律責(zé)任。如果您發(fā)現(xiàn)本站中有涉嫌抄襲或描述失實的內(nèi)容，請聯(lián)系我們jiasou666@gmail.com 處理，核實后本網(wǎng)站將在24小時內(nèi)刪除侵權(quán)內(nèi)容。