samtools 工具处理SAM文件

image

SAM 是用来存储核苷酸序列比对信息的文件格式。SAM Tools 工具包提供各种工具处理SAM文件。包括功能：sorting, merging, indexing and generating alignments。安装教程见https://www.jianshu.com/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文件转化为bam文件

samtools view -bS in.sam > in.bam

是view的一个应用-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参数。

随机取出bam文件的某一部分出来，必须要有index文件，否则会报错： [bam_index_load] fail to load BAM index. [main_samview] random alignment retrieval only works for indexed BAM files.

关于view更详细的参数说明：

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文件转化为sam文件

samtools view -h R12.merged.bam > test.sam

-h是将bam文件中的header也加入到sam文件中。比如htseq-count老版本只接受sam文件

二、建立索引Indexing

注意

如果你执行命令的地方和参考序列不在同一个目录，参考序列用全路径，最后的index结果和参考序列在同一个目录里面，而不是执行命令的目录。
在fasta文件中，对于某一个序列，除了最后一行，其他行所含碱基数应该一样。

1、对fasta文件建立索引

samtools faidx ref.fasta`

2、对BAM文件建立索引

samtools index in.bam

结果文件名为in.bam.bai

3、知道位置信息查找对应的序列信息

samtools faidx ref.fa Chr1:33667-33667

指查看染色体一上的第33667个碱基。

三、将bam文件进行sort

只能对bam文件进行sort, 不能对sam文件。

samtools sort aln.bam anl.sorted

默认是根据coordinate进行sort，如果输入bam文件为in.bam , 则输出文件名为in.sorted.bam
如果要按照read name进行sort, 需要加-n, 如heseq-count 就要求文件时按照read name 而不是coordinate。

samtools sort -n aln.bam anl.sorted

四、去除bam文件中pcr导致的重复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是某个某样本的三个重复，我们可以将他们合并为一个bam文件。

samtools merge -R chr1 out.bam in1.bam in2.bam in3.bam

如果想对部分合并，如至合并一号染色的上的bam文件合并，chr1必须为序列的名字，一号染色体序列的名字为Chr1，那么就应为-R Chr1

注意：要合并的bam文件，必须有对应的index文件。

samtools index in.bam  #结果文件名为in.bam.bai

六、统计bam文件中的reads情况，有多少reads比对上

命令：

samtools flagstat RF12.merged.bam

结果如下：

66196754 + 0 in total (QC-passed reads + QC-failed reads)
#bam文件中所有的reads数。

0 + 0 duplicates
66196754 + 0 mapped (100.00%:-nan%)
#比对到基因组上的reads数目， 我用的比对软件是tophat， 结果中只保留比对上的reads信息。

66196754 + 0 paired in sequencing
#属于paired reads数目， 我的数据都是双端测序。所以都是paired reads。

33493586 + 0 read1
#来自于read1中的reads数目

32703168 + 0 read2
#来自于read2 中的reads数目
#read1 + read2 = paired reads

45729393 + 0 properly paired (69.08%:-nan%)
#完美匹配的reads数， 即一对reads比对到同一条参考序列，并且这对reads之间的举例符合设置的阈值

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)
#和上面一样，只不过比对的质量大于等于5的reads数目 </pre>