轉(zhuǎn)載：https://biozx.top/samtools.html
image
SAM 是用來存儲核苷酸序列比對信息的文件格式。SAM Tools 工具包提供各種工具處理SAM文件。包括功能：sorting, merging, indexing and generating alignments。安裝教程見http://www.itdecent.cn/p/53de170927a7
安裝過程中有許多依賴的庫需要安裝的，可能每個人缺的庫都不盡相同，不懂的就百度一下吧：
sudo apt-get update
sudo apt install gcc
sudo apt install git
sudo apt install make
sudo apt-get install libncurses5-dev ##bam_tview_curses.c:41:20: fatal error: curses.h: No such file or directory
sudo yum install bzip2-devel ##cram/cram_io.c:57:19: fatal error: bzlib.h: No such file or directory
sudo apt-get install liblzma-dev ##cram/cram_io.c:60:18: fatal error: lzma.h: No such file or directory
裝好之后有如下這么多命令，下面我們只介紹samtools：
image
rqq@ubuntu:~/bio$ samtools
Program: samtools (Tools for alignments in the SAM format)Version: 1.5 (using htslib 1.5)
Usage: samtools <command> [options]
Commands: -- Indexing
dict create a sequence dictionary file
faidx index/extract FASTA
index index alignment
-- Editing
calmd recalculate MD/NM tags and '=' bases
fixmate fix mate information
reheader replace BAM header
rmdup remove PCR duplicates
targetcut cut fosmid regions (for fosmid pool only)
addreplacerg adds or replaces RG tags
-- 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
一、samtools view功能
1、將sam文件轉(zhuǎn)化為bam文件
samtools view -bS in.sam > in.bam
是view的一個應(yīng)用-b指定輸出文件為bam, -S 指定輸入文件為sam
2、查看bam文件的header信息
samtools view -H in.bam
3、取出bam文件的一部分
samtools view -b your.bam Chr1 > Chr1.bam
- -b 指定是輸出文件為bam,
- Chr1指定你要看的是哪一部分， 這里指看Chr1那一部分，然后重定向到一個新的bam文件，注意，這個bam文件是沒有header的，如果想要包括header 可以使用-h參數(shù)。
隨機取出bam文件的某一部分出來， 必須要有index文件，否則會報錯： [bam_index_load] fail to load BAM index. [main_samview] random alignment retrieval only works for indexed BAM files.
關(guān)于view更詳細(xì)的參數(shù)說明：
Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]
Options:
-b output BAM
-h print header for the SAM output
-H print header only (no alignments)
-S input is SAM
-u uncompressed BAM output (force -b)
-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
-B collapse the backward CIGAR operation
-@ INT number of BAM compression threads [0]
-L FILE output alignments overlapping the input BED FILE [null]
-t FILE list of reference names and lengths (force -S) [null]
-T FILE reference sequence file (force -S) [null]
-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]
-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
4、bam文件 轉(zhuǎn)化為sam文件
samtools view -h R12.merged.bam > test.sam
-h是將bam文件中的header也加入到sam文件中。 比如htseq-count老版本只接受sam文件
二、建立索引Indexing
注意
如果你執(zhí)行命令的地方和參考序列不在同一個目錄，參考序列用全路徑，最后的
index結(jié)果和參考序列在同一個目錄里面，而不是執(zhí)行命令的目錄。在fasta文件中，對于某一個序列，除了最后一行，其他行所含堿基數(shù)應(yīng)該一樣。
1、對fasta文件建立索引
samtools faidx ref.fasta`
2、對BAM文件建立索引
samtools index in.bam
- 結(jié)果文件名為in.bam.bai
3、知道位置信息查找對應(yīng)的序列信息
samtools faidx ref.fa Chr1:33667-33667
指查看染色體一上的第33667個堿基。
三、將bam文件進(jìn)行sort
只能對bam文件進(jìn)行sort, 不能對sam文件。
samtools sort aln.bam anl.sorted
- 默認(rèn)是根據(jù)coordinate進(jìn)行sort， 如果輸入bam文件為in.bam , 則輸出文件名為in.sorted.bam
- 如果要按照read name進(jìn)行sort, 需要加
-n, 如heseq-count 就要求文件時按照read name 而不是coordinate。
samtools sort -n aln.bam anl.sorted
四、去除bam文件中pcr導(dǎo)致的重復(fù)reads信息
samtools rmdup in.bam in.rmp.bam
五、合并bam文件
samtools merge out.bam in1.bam in2.bam in3.bam
假如in1.bam, in2.bam, in3.bam是某個某樣本的三個重復(fù)，我們可以將他們合并為一個bam文件。
samtools merge -R chr1 out.bam in1.bam in2.bam in3.bam
如果想對部分合并，如至合并一號染色的上的bam文件合并，chr1必須為序列的名字，一號染色體序列的名字為Chr1，那么就應(yīng)為-R Chr1
注意：要合并的bam文件，必須有對應(yīng)的index文件。
samtools index in.bam #結(jié)果文件名為in.bam.bai
六、統(tǒng)計bam文件中的reads情況，有多少reads比對上
命令：
samtools flagstat RF12.merged.bam
結(jié)果如下：
66196754 + 0 in total (QC-passed reads + QC-failed reads)
#bam文件中所有的reads數(shù)。
0 + 0 duplicates
66196754 + 0 mapped (100.00%:-nan%)
#比對到基因組上的reads數(shù)目， 我用的比對軟件是tophat， 結(jié)果中只保留比對上的reads信息。
66196754 + 0 paired in sequencing
#屬于paired reads數(shù)目， 我的數(shù)據(jù)都是雙端測序。所以都是paired reads。
33493586 + 0 read1
#來自于read1中的reads數(shù)目
32703168 + 0 read2
#來自于read2 中的reads數(shù)目
#read1 + read2 = paired reads
45729393 + 0 properly paired (69.08%:-nan%)
#完美匹配的reads數(shù)， 即一對reads比對到同一條參考序列，并且這對reads之間的舉例符合設(shè)置的閾值
61118410 + 0 with itself and mate mapped
#一對reads 都比對到了參考序列上，但是不一定比對到同一條染色體
5078344 + 0 singletons (7.67%:-nan%)
#一對reads中，只有一條匹配到了參考序列上。
#61118410+5078344=66196754
703397 + 0 with mate mapped to a different chr
#一對reads比對到了不同的染色體上。針對的是with itself and mate mapped
346693 + 0 with mate mapped to a different chr (mapQ>=5)
#和上面一樣，只不過比對的質(zhì)量大于等于5的reads數(shù)目 </pre>
