samtools 的使用
主要介绍一下 SAMtools 的用法。它被誉为 NGS 里面的 瑞士军刀,功能很是强大,所以值得好好学习一番。
文档参考
推荐: blog SAMtools: Primer/Tutorial
安装
由于 samtools 与 htslib , bcftools 经常结合使用,因此,最好一次就把这三个一起安装好,届时使用更方便。
下载,安装
以 samtools 为例
下载,解压到 samtools/ 目录下, 里面有一个 INSTALL 文件,可以看看说明
$ cd samtools-1.2 # and similarly for bcftools and htslib
$ make
很有可能会出现如下的错误信息:
bam_tview_curses.c:30:20: fatal error: curses.h: No such file or directory
compilation terminated.
make: *** [bam_tview_curses.o] Error 1
这说明缺少curses 文件,google 后,找到解决方法:
解决方法
$ sudo apt-get update
$ sudo apt-get install ncurses-dev
这样,问题就解决啦
$ make prefix=/home/user/samtools-1.2/ install
此时,可运行的软件已经安装到了 /home/user/samtools-1.2/bin/ 下,所以下一步把它加入到 PATH 中即可:
$ export PATH=/home/user/samtools-1.2/bin:$PATH
这样就基本算安装,配置完了,可以测试一下:
$ samtools
应该会输出各种信息,包括版本,语法,用法等等,这样就算成功了。
剩下的 htslib 与 bcftools 的安装与这个一样。
这三个的安装比较简单,一切妥当后, 就可以进行使用了。
使用samtools
参考推荐: blog SAMtools: Primer/Tutorial
前提是有已经通过比对软件 mapping 后的初始 sam 文件
方法在之前的博客中已经说过了
比对程序 BWA 的简单用法
比对程序 Bowtie 的简单使用
比对程序 Bowtie2 的简单使用
将 sam 转换成 bam 文件
$ samtools view # 查看 view 的用法
$ samtools view -b -S -o /media/文档/aln-pe.bam /media/文档/aln-pe.sam
# -b: indicates that the output is BAM.
# -S: indicates that the input is SAM.
# -o: specifies the name of the output file
这样就把 aln-pe.sam 文件转换成 aln-pe.bam 文件
sorting
$ samtools sort # 查看 sort 的用法
$ samtools sort /media/文档/aln-pe.bam /media/文档/aln-pe.sorted
这样就把 上一步的 aln-pe.bam 变为 sorted file : aln-pe.sorted.bam
index
$ samtools index # 查看 index 的用法
$ samtools index /media/文档/aln-pe.sorted.bam
This enables tools, including SAMtools itself, and other genomic viewers to perform efficient random access on the BAM file, resulting in greatly improved performance
This reads in the sorted BAM file, and creates a BAM index file with the .bai extension: aln-pe.sorted.bam.bai
至此,现在就可以利用软件 IGV 去查看 mapping 结果了。
下面就可以鉴定基因组的变异位点了,也就是所谓的 SNP 以及 INDEL 等。
samtools mpileup
$ samtools mpileup # 查看 mpileup 的用法
$ samtools mpileup -g -f /media/文档/tigr7.fa /media/文档/samtools-test/aln-pe.sorted.bam > /media/文档/samtools-test/sim_variants.bcf
-g: directs SAMtools to output genotype likelihoods in the binary call format (BCF). This is a compressed binary format.
-f: directs SAMtools to use the specified reference genome
这一步会生成一个 sim_variants.bcf 文件,很大...
The mpileup command automatically scans every position supported by an aligned read, computes all the possible genotypes supported by these reads, and then computes the probability that each of these genotypes is truly present in our sample
bcftools call
$ bcftools call # 查看 call 的用法
$ bcftools call -c -v /media/文档/samtools-test/sim_variants.bcf > /media/文档/samtools-test/sim_variants.vcf
# -c: directs bcftools to call SNPs.
# -v: directs bcftools to only output potential variants
这样就生成了一个 sim_variants.vcf 文件。里面包含了鉴定出的 SNP 以及 INDEL 的详细信息,可以用文本编辑器来查看
Visualizing Reads
除了用软件如 IGV 来查看 mapping 的reads 的情况外,也可以在命令行下使用 samtools 去查看。
$ samtools tview # 查看 tview 的用法
$ samtools tview /media/文档/samtools-test/aln-pe.sorted.bam /media/文档/tigr7.fa
终于算是把 samtools , bcftools 的最简单的用法过了一遍,操作起来并不难,但是很多必要的原理知识还没有去看,去阅读体会,所以这仅仅算是一个入门而已。不过到这一步也算不错的。。