samtools

对于一个生物信息从业人员来说,玩的就是序列,这里面离不开的是序列的比对,作为一个强大的查看序列比对结果的工具,有必要好好认识一下samtools。

 一、samtools的简介

SAMtools: Primer / Tutorial  http://biobits.org/samtools_primer.html


1,samtools的安装
下载samtools工具:http://sourceforge.net/projects/
samtools/files/tar zxvf samtools-0.1.19.tar.bz2
cd samtools-0.1.19/
make                    #为什么他的安装只需要make?
samtoolpath=pwd
PATH=PATH:$samtoolpath
修改bash文件:
export SAMTOOLS_HOME=/Users/ecerami/libraries/samtools-0.1.19

[/crayon]

 二、常用命令

2.1 help

Command: view        SAM<->BAM conversion 两种文件的互换
         sort        sort alignment file
         mpileup     multi-way pileup
         depth       compute the depth
         faidx       index/extract FASTA
         tview       text alignment viewer
         index       index alignment
         idxstats    BAM index stats (r595 or later)
         fixmate     fix mate information
         flagstat    simple stats
         calmd       recalculate MD/NM tags and ‘=’ bases
         merge       merge sorted alignments
         rmdup       remove PCR duplicates
         reheader    replace BAM header
         cat         只支持合并BAM文件
         bedcov      read depth per BED region
         targetcut   cut fosmid regions (for fosmid pool only)
         phase       phase heterozygotes
         bamshuf     shuffle and group alignments by name
2.2 View
将sam文件转换成bam文件;然后对bam文件进行各种操作,比如数据的排序(不属于本命令的功能)和提取(这些操作是对bam文件进行的,因而当输入为sam文件的时候,不能进行该操作);最后将排序或提取得到的数据输出为bam或sam(默认的)格式。
bam文件优点:bam文件为二进制文件,占用的磁盘空间比sam文本文件小;利用bam二进制文件的运算速度快。

samtools 将SAM文件转化为BAM文件 
    samtools view -bS eg2.sam > eg2.bam
也可以输入:

[/crayon]
   -b: indicates that the output is BAM.
   -S: indicates that the input is SAM.
   -o: specifies the name of the output file.
因为bam文件是一个压缩的二进制文件,不能直接看,所以 需要使用view功能

[/crayon]
也可以使用view来仅仅展示那些匹配到specific filtering criteria的reads

[/crayon]
    -f int :仅仅提取匹配specified SAM flag的reads,上面的那个例子中,我们仅仅提取那些flag value 为4的匹配上去的reads

[/crayon]
     -F init: 上面的例子中是删掉匹配上flag value 为4的reads

[/crayon]
     上面不产生结果文件,仅仅甘肃你多少个reads没有匹配到reference genome

[/crayon]
    显示匹配的分数在42分以上的reads数

2.3 sort,index

用samtools排序和然后建立索引

eg2.sorted.bam
[/crayon]
第一步生成eg2.sorted.bam文件,第二步生成sim_reads_aligned.sorted.bam.bai,如果没有排序的话,是没有办法建立索引的。
建立索引后将产生后缀为.bai的文件,用于快速的随机处理。很多情况下需要有bai文件的存在,特别是显示序列比对情况下。比如samtool的tview命令就需要;gbrowse2显示reads的比对图形的时候也需要。

2.4 merge/cat

merge将2个或2个以上的已经sort了的bam文件融合成一个bam文件。融合后的文件不需要则是已经sort过了的。

2.5 faidx

对fasta文件建立索引,生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条(子)序列

2.6  tview

tview能直观的显示出reads比对基因组的情况,和基因组浏览器有点类似。

2.7  flagstat 统计比对结果

2.8 depth 得到每个碱基位点的测序深度

得到每个碱基位点的测序深度,并输出到标准输出。

2.9. reheader   替换bam文件的头

2.10  idxstats  评估map的质量

注意,这一步之前需要经过sort和index。结果会显示:

其中第一列是染色体名称,第二列是序列长度,第三列是mapped reads数,第四列是unmapped reads数,paired reads中有一端能匹配到该scaffold上,而另外一端不匹配到任何scaffolds上的reads数。idxstats 统计一个表格,4列,分别为”序列名,序列长度,比对上的reads数,unmapped reads number”。第4列应该是paired reads中有一端能匹配到该scaffold上,而另外一端不匹配到任何scaffolds上的reads数。

2.11 mpileup  寻找突变通过VCF格式

samtools还有个非常重要的命令mpileup,以前为pileup。该命令用于生成bcf文件,再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件,在samtools的安装文件夹中可以找到。

最常用的参数有2:

-f 来输入有索引文件的fasta参考序列;
-g 输出到bcf格式。用法和最简单的例子如下

mpileup不使用-u或-g参数时,则不生成二进制的bcf文件,而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况;该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如:

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]+’代表在该位点后缺失的碱基;

samtools mpileup -uf $BT2_HOME/example/reference/lambda_virus.fa eg2.sorted.bam | bcftools view -bvcg – > eg2.raw.bcf

samtools mpileup -g -f genomes/NC_008253.fna alignments/sim_reads_aligned.sorted.bam > variants/sim_variants.bcf
    通过工具mpileup统计通过比对基因型相似的reads
-g: 指导SAMtools产生基因型相似性的二进制访问格式(binary call format)(BCF).
-f: 指导SAMtools使用指定的相关基因组,相关基因组必须被指定,上面的例子指定的是E. coli
    mpileup命令能够自动浏览咩个位置的比对的reads,计算所有基因型的可能性,然后统计在我们样品中这些基因型的每一个的可能性。

2.12  bcftools突变位点

bcftools和samtools类似,用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件,后者为其二进制文件。
bcftools使用简单,最主要的命令是view命令,其次还有index和cat等命令。index和cat命令和samtools中类似。此处主讲使用view命令来进行SNP和Indel calling。该命令的使用方法和例子为:

生成的结果文件为vcf格式,有10列,分别是:
1 参考序列名;
2 variant所在的left-most位置;
3 variant的ID(默认未设置,用’.’表示);
4 参考序列的allele;
5 variant的allele(有多个alleles,则用’,’分隔);
6 variant/reference QUALity;
7 FILTers applied;
8 variant的信息,使用分号隔开;
9 FORMAT of the genotype fields, separated by colon (optional);
10 SAMPLE genotypes and per-sample information (optional)。

例如:

第8列中显示了对variants的信息描述,比较重要,其中的 Tag 的描述如下:

bcftools view 的具体参数如下:

使用bcftools得到variant calling结果后。需要对结果再次进行过滤。主要依据比对结果中第8列信息。其中的 DP4 一行尤为重要,提供了4个数据:1 比对结果和正链一致的reads数、2 比对结果和负链一致的reads数、3 比对结果在正链的variant上的reads数、4 比对结果在负链的variant上的reads数。可以设定 (value3 + value4)大于某一阈值,才算是variant。比如:

    bcftools view eg2.raw.bcf

[/crayon]

-c: directs bcftools to call SNPs.
-v: directs bcftools to only output potential variants

bcftools view命令用之前产生的基因型相似性来寻找SNP和节点,结果保存在vcf中

2.13 rmdup

NGS上机测序前需要进行PCR一步,使一个模板扩增出一簇,从而在上机测序的时候表现出为1个点,即一个reads。若一个模板扩增出了多簇,结果得到了多个reads,这些reads的坐标(coordinates)是相近的。在进行了reads比对后需要将这些由PCR duplicates获得的reads去掉,并只保留最高比对质量的read。使用rmdup命令即可完成.

三、讨论

3.1 depth和mpileup的区别

3.2 samtools和igv看到的depth不一致

这里面有两个原因:序列的质量以及duplicates

samtools在计算depth的时候,会过滤掉duplicates的序列,解决办法:

3.3  将bam文件转换为fastq文件

有时候,我们需要提取出比对到一段参考序列的reads,进行小范围的分析,以利于debug等。这时需要将bam或sam文件转换为fastq格式。

该网站提供了一个bam转换为fastq的程序:http://www.hudsonalpha.org/gsl/information/software/bam2fastq

3.4 报错: hts_open_format

问题原因:

解决办法:

1.去1000g数据库中下载对应的bai文件(。。。。)
2.修改系统容许打开的文件,提高命令中threads的值(未这么做)

四、我的案例

3.1求得coverage
cd /sam/uncltured/contig7/index  #建立索引
/sam/bowtie2/bowtie2-build  ../contig.fasta scaffold  #这里的contig.fasta是我变形以后的scaffold,
先用paired 去比对
/sam/bowtie2/bowtie2 -p 12 -x scaffold -1 /sam/idba_ud/myproject/reads2/output_forward_paired.fq -2 /sam/idba_ud/myproject/reads2/output_reverse_paired.fq -S eg2.sam
219079809 reads; of these:
  219079809 (100.00%) were paired; of these:
    9992084 (4.56%) aligned concordantly 0 times
    190538854 (86.97%) aligned concordantly exactly 1 time
    18548871 (8.47%) aligned concordantly >1 times
    —-
    9992084 pairs aligned concordantly 0 times; of these:
      985097 (9.86%) aligned discordantly 1 time
    —-
    9006987 pairs aligned 0 times concordantly or discordantly; of these:
      18013974 mates make up the pairs; of these:
        9341217 (51.86%) aligned 0 times
        4121203 (22.88%) aligned exactly 1 time
        4551554 (25.27%) aligned >1 times
97.87% overall alignment rate
再用unpair去比对
/sam/bowtie2/bowtie2 -p 12 -x scaffold -U /sam/idba_ud/myproject/reads2/unpair.fq -S eg1.sam
9910673 reads; of these:
  9910673 (100.00%) were unpaired; of these:
    303363 (3.06%) aligned 0 times
    8085141 (81.58%) aligned exactly 1 time
    1522169 (15.36%) aligned >1 times
96.94% overall alignment rate
samtools view -bS eg1.sam > eg1.bam
[samopen] SAM header is present: 144759 sequences.
samtools view -bS eg2.sam > eg2.bam            [11:23PM]
[samopen] SAM header is present: 144759 sequences
samtools cat eg1.bam eg2.bam >eg.bam
samtools sort eg.bam eg.sorted
[bam_sort_core] merging from 184 files…

参考资料:
(推荐)SAMtools: Primer / Tutorial  http://biobits.org/samtools_primer.html
菜菜的一天的博客 http://blog.sina.com.cn/s/blog_670445240101l30k.html
糗世界的博客
https://www.biostars.org/p/57968/
https://www.plob.org/article/7112.html

发表评论

电子邮件地址不会被公开。 必填项已用*标注