10.18 bcftools下载使用 samtools语句
师兄让下载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了解流程