生物信息学

10.18 bcftools下载使用 samtools语句

2018-10-18  本文已影响0人  KK_f2d5

师兄让下载bcftools,学习一下。

1.下载

电脑太不给力了,wget都没有,用brew下载一个

太慢了,在github直接下载。

https://github.com/samtools/bcftools/releases

安装看readme

共享目录出错,

这个时候就要在/etc/ld.so.conf中加入xxx.so所在的目录。

一般而言,有很多so档会在/usr/local/lib这个目录下

在/etc/ld.so.conf中加入/usr/local/lib这一行后,执行 在etc中执行ldconfig 命令来更新生效。

./bftools 使用


2.samtools语句

介绍samtools和bcftools的语句与用法 :https://www.jianshu.com/p/da32e3c3a168

samtools基本语句和功能:

1.view 查看bam,sam文件的内容

2.sort用来对bam文件进行排序

3.merge将2个或2个以上的已经sort的bam文件合并成一个bam文件,合并后的文件不需要再次sort,是已经sort过了的。

4.index对bam文件建立索引,生成后缀为.bai的文件,用于快速检索reads。需要对bam文件先进行排序,否则会报错。很多时候都需要索引文件的存在,特别是显示序列比对的情况下。比如samtools tview,gbrowse2等。

5.faidx对fasta序列建立索引,生成后缀为.fai的文件。该命令也能依据索引文件快速提取fasta文件中的某一条序列。

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

7.flagstat给出bam文件的比对结果的summary。

8.depth得到每个位点的测序深度,并输出到标准输出

samtools depth a.bam | less

9.rmdup去除PCR重复和光学重复,重复的reads仅保留一对。

10.reheader替换bam文件的头文件

11.cat链接多个bam文件,适用于非sorted的bam文件

12.mpileup(现在也属于bcftools)

使用samtools的子命令mpileup分析参考序列上的每个碱基位点的比对结果,并生成VCF/BCF格式文件,BCF是VCF的二进制文件。再使用Bcftools对VCF/BCF格式文件进行SNP/Indel calling。其中,Bcftools是附属于samtools的程序。

后面接Bcftools进行variation calling

bcftools进行SNP calling:  www.jianshu.com/p/af93ca26ba44

time $bcf mpileup ${precallbam} --fasta-ref ${reference} | $bcf call -mv -o ./out/bcfcall.vcf && echo " ** bcf snp calling done **"

在进行SNP calling 时,必须选择一种算法,有两种calling算法可供选择,分别对应-c和-m参数。-c参数对应consensus-caller算法,-m参数对应multiallelic-caller算法,后者更适合多种allel和罕见变异的calling。

-v参数也是常用参数,作用是只输出变异位点的信息,如果一个位点不是snp/indel, 不会输出。


3.bam文件格式

samtools view -H SRR1770413.sorted.markdup.bam

详细讲解bam文件:https://www.jianshu.com/p/364e640d3c9f



4.bftools合并vcf文件

师兄用bftools合并vcf文件,这里的vcf文件不是gatk得到的,没有combine过。

https://www.jianshu.com/p/59728f20d38c


5 vcf 文件格式

https://www.jianshu.com/p/38f734ae47f5

https://www.jianshu.com/p/957efb50108f


刚刚看到一个实战:

https://www.jianshu.com/p/6f3198b7a070

发现gatk在github上有workflow: https://github.com/gatk-workflows

其中vcf文件生成步骤:

1.HaplotyperCaller      2.MergeVcfs     3.GenomicsDBImport     4.GenotypeGVCFs

在github上选择workflow打开,查看wdl文件中的command了解流程

上一篇下一篇

猜你喜欢

热点阅读