搜索
查看: 9720|回复: 5

使用bcftools、freebayes和varscan2 call SNV

[复制链接]

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-9-18 22:06:44 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-9-18 22:16 编辑

看了直播基因组系列,学习了下这3种不同的方法来call SNV,简单做个笔记

虽然软件来看比较老,但是网上教程和使用的人也不少,说明还是值得学习下的

  • bwa mem算法 比对human_g1k_v37.fasta参考序列
  • 用samtools sort或者picard.jar SortSam对比对后的sam进行默认排序,然后转化为bam文件
  • 用picard MarkDuplicates对bam文件进行mark duplicates,以免PCR重复reads对后续call snp造成影响,没必要去除,可参照https://www.biostars.org/p/3917/
  • 用samtools merge将多个lane合并(如果是多个lane的WGS的话),生成merged_marked.bam文件
  • 使用samtools mpileup + bcftools call SNP

    • 先使用samtools mpileup将bam文件生成bcf文件(二进制文件)
        mpileup是samtools中call snp的工具,可以不使用-g参数,则会生成一个文本格式的文件,我们可以看到参考序列上每个碱基的比对结果:
        总共6列,分别是参考序列名(染色体),位置,参考序列的碱基,比对上的reads数,比对情况,比对上的碱基的质量
        生成bcf文件的主要参数:
      -g Compute genotype likelihoods and output them in the binary call format (BCF)(计算genotype likelihoods并以bcf格式输出)  
      -f The faidx-indexed reference file in the FASTA format(用samtools faidx对参考序列建索引.fai文件,用其他软件建索引也可以的)  
      -u Generate uncompressed VCF/BCF output(输出不压缩的bcf文件)
      -t Comma-separated list of FORMAT and INFO tags to output (case-insensitive)(按照自己的需求,输出vcf格式中常见的各种tags,比如-t DP输出vcf中的DP值;-t SP输出类似GATK call snp 出的VCF的GQ值?)
        
      [AppleScript] 纯文本查看 复制代码
      samtools mpileup -go samtools.bcf -f ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta -t DP -t SP merged_marked.bam
        后来发现,这里-t DP -t SP加不加都行,因为最后的bcftools call snp时,都会有DP和SP的tag        
    • 再用bcftools将bcf文件生成常用的vcf格式文件
        bcftools call snp的主要参数如下:
      -O Output compressed BCF (b), uncompressed BCF (u), compressed VCF (z), uncompressed VCF (v)(一般可以选择z,压缩的vcf格式)  
      —threads Number of output compression threads to use in addition to main thread(主要用于压缩时的多线程,在-O b/z时使用)  
      -m alternative modelfor multiallelic and rare-variant calling designed to overcome known limitations in -c calling model(一种更优的新算法,克服了旧算法的缺陷)  
      -v output variant sites only(只输出突变位点信息)
      -f —format-fields 可以在输出的VCF文件中增加新的tag,例如增加GQ
      -V, —skip-variants 可以跳过输出snps|indels(比如我想只输出snps,则跳过indels)
        下面我就用一个常见的输出:
        
      [AppleScript] 纯文本查看 复制代码
      bcftools call -vmO z -o bcftools_raw.vcf.gz samtools.bcf
        当然上述两步可以合起来运行(但是不能提交后台运行,反正会报错。。。)
        
      [AppleScript] 纯文本查看 复制代码
      samtools mpileup -ugf ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta -t DP -t SP merged_marked.bam |bcftools call -vmO z -o bcftools_raw.vcf.gz
    • 用bcftools过滤不可靠的位点
        bcftools filter的主要参数如下:
      -e  —exclude 主要用表达式方式去除过匹配上的位点,这个参数很关键,很多过滤都依靠这个表达式
      -g, —SnpGap 过滤INDEL附近的snp位点,比如-SnpGap 5则过滤INDEL附近5个碱基距离内的SNP
      -G, —IndelGap 过滤INDEL附近的INDEL位点。。。应该是这个意思
      -o, —output 输出文件的名称
      -O, —output-type 输出的格式,一般z和v都行
      -s, —soft-filter 将过滤掉的位点用字符串注释
      -S, —set-GTs set genotypes of failed samples to missing value (.) or reference allele (0)(将不符合要求的个体基因型改为./.)
        下面就简单的过滤QUAL小于10,DP值小于5,INDEL附近的位点
        
      [AppleScript] 纯文本查看 复制代码
      bcftools filter -O v -o bcftools_filter.vcf -s LOWQUAL -e 'QUAL<10 || FMT/DP <5' --SnpGap 5 --set-GTs . bcftools_raw.vcf.gz
    • 提取过滤后的SNP位点
        
      [AppleScript] 纯文本查看 复制代码
      bcftools view -v snps bcftools_filter.vcf >bcftools_snp_filter.vcf
        或者在vcf文件中INFO列里,如果是INDEL的话,会标注出INDEL,因此提取SNP也可以:
        
      [AppleScript] 纯文本查看 复制代码
      grep -v 'INDEL' bcftools_fiter.vcf >bcftools_snp_filter.vcf


  • 使用Freebayes Call SNP

    • 下载及安装
        git clone —recursive git://github.com/ekg/freebayes.git
        make
        make install
    • 软件介绍及原理
    • 使用freebayes从bam到vcf文件
      [AppleScript] 纯文本查看 复制代码
       freebayes -h
        可看到参数众多,对几个主要参数做个解释
      —min-alternate-count (default: 2):Require at least this count of observations supporting an alternate allele within a single individual in order to evaluate the position(最少的alt碱基数,可以根据你测序深度来定)
      —min-alternate-fraction (default: 0.2):Require at least this fraction of observations supporting an alternate allele within a single individual in the
        in order to evaluate the position(跟—min-alternate-count类似,就是变成比例了,双倍体用默认值0.2即可,如果多倍体(>2)的话,建议使用>0.2的值)
      —genotype-qualities :Calculate the marginal probability of genotypes and report as GQ in each sample field in the VCF output(在输出的VCF文件中有GQ的tag)
      —min-mapping-quality(default: 1):Exclude alignments from analysis if they have a mapping quality less than Q
      —min-base-quality (default: 0):Exclude alleles from analysis if their supporting base quality is less than Q (这个和上面两个参数是对input进行过滤的,但我觉得也可以在最后vcf文件,对结果再过滤也行?)
      -X —no-mnps:Ignore multi-nuceotide polymorphisms, MNPs.
      -u —no-complex:Ignore complex events (composites of other classes)(官网是不建议这么做,其认为减少这些信息会降低检测变异的可行度,其建议在最后文件中再去过滤获得snp or indel)
        从上述可看出,基本上如果单纯的call变异的话,并且前期reads已做质控,那么在freebayes call变异这步几乎可以说使用默认参数即可(对于2倍体生物),或者后续有其他需求可再加入其他参数,所以命令如下:
        
      [AppleScript] 纯文本查看 复制代码
      freebayes -f ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta merged_marked.bam >freebayes_raw.vcf
    • 提取SNP以及INDEL
        freebayes产生的VCF文件中INFO一列有专门的一个tag来注释是snp、ins(插入)、del(缺失)、mnp(连续两个snp位点,如ref为AT, alt为CG)以及complex(composite insetion and substitution events)
        所以只要用下面的shell命令即可提取SNP
        
      [AppleScript] 纯文本查看 复制代码
      grep 'TYPE=snp' freebayes_raw.vcf > freebayes_snp_raw.vcf
        其他类型的变异也是类似
    • 最后就是过滤低质量的SNP位点了
        可以用bcftools过滤,可看之间的文章http://www.bioinfo-scrounger.com/archives/248
        也可以用vcftools,粗略看过,功能比bcftools的filter功能要强大点
        当然如果是简单过滤,直接脚本就行了


  • 使用VarScan2 Call SNP

    • 为什么选用VarScan来call variant,VarScan是这么说的:对于NGS数据,大多数variant caller都是基于贝叶斯算法来识别变异位点,然后通过可信度来评估,但会受extreme read depth,pooled samples,contaminated or impure samples的影响;而VarScan是基于启发式算法来寻找变异位点,可以很好的meet desired thresholds for read depth, base quality, variant allele frequency, and statistical significance
    • 安装很简单,下载jar文件即可使用,https://sourceforge.net/projects/varscan/files/
    • 用samtools mpileup生成mpileup文件作为VarScan的输入文件,这里跟samtools+bcftools流程的是不同的,所以只需要输出mpileup文件即可,那么就不用加-g参数了,还有其他参数倒是可以加上
      -q, -min-MQ 比对的mapping quality
      -d, —max-depth 最大测序深度,过滤掉超深度测序的位点
        samtools mpileup -q 1 -d 30000 -f ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta merged_marked.bam 1>merged_marked.mpileup 2>mpileup.log
    • 既然是call snp,那么肯定是用VarScan2的mpileup2snp命令,具体参数可以查看http://dkoboldt.github.io/varscan/using-varscan.html#v2.3_mpileup2snp,其他命令也类似
        
      [AppleScript] 纯文本查看 复制代码
      USAGE: java -jar VarScan.jar mpileup2snp [mpileup file] OPTIONS mpileup file - The SAMtools mpileup file
        参数不多,全部列出来也可以:
      —min-coverage    Minimum read depth at a position to make a call(位点的测序深度,覆盖度)默认 8
      —min-reads2    Minimum supporting reads at a position to call variants(paired reads,同上)默认 2
      —min-avg-qual    Minimum base quality at a position to count a read(位点上的reads碱基质量)默认 15
      —min-var-freq    Minimum variant allele frequency threshold(不太懂这个阈值是怎么卡的。。。。。。)默认 0.01
      —min-freq-for-hom    Minimum frequency to call homozygote(这个还是不懂)默认 0.75
      —p-value    Default p-value threshold for calling variants(call出来的varians的可行度)默认是0.99
      —strand-filter     Ignore variants with >90% support on one strand 默认1
      —output-vcf    If set to 1, outputs in VCF format
      —variants    Report only variant (SNP/indel) positions (mpileup2cns only)(对于mpileup2snp就无视这个了) [0]
        由上可看出,大部分默认参数即可,特定情况再依情况修改
        
      [AppleScript] 纯文本查看 复制代码
      java -jar VarScan.jar mpileup2snp merged_marked.mpileup --output-vcf 1 >varscan.snp.vcf
        但是VarScan的http://dkoboldt.github.io/varscan/germline-calling.html里面对于input文件(mpileup)又一个卡阈值
        简单的说就是满足minimum base quality=20,minimum mapping quality=1,minimum coverage=20的位点才call variant,怎么跟mpileup2snp的默认参数又有些不一致了呢,不知道有无理解错
    • 过滤位点,VarScan call出来的vcf文件有点跟其他软件不同,最明显的地方就是没有QUAL值,毕竟算法不一样嘛,但是其他常用的参数比如GQ,DP还是都有值的
        varscan也有自己filter命令,内容很简单http://dkoboldt.github.io/varscan/using-varscan.html#v2.3_filter,可以看看
        所以也可以从vcf的第10列的信息来过滤,脚本可以简单的过滤低GQ和DP值的位点

    也参考自 生信技能树、生物信息学与机器学习等公众号
[size=0em]​




上一篇:我问大家一个很弱的问题
下一篇:叶绿体基因组上传NCBI,求助大神,急急急!!!
回复

使用道具 举报

4

主题

56

帖子

555

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
555
发表于 2017-9-19 14:19:35 | 显示全部楼层
虽然你给了这3个软件的用法指导,但是目前主流用的是GATK
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
 楼主| 发表于 2017-9-19 15:34:14 | 显示全部楼层
bioinfotrainee 发表于 2017-9-19 14:19
虽然你给了这3个软件的用法指导,但是目前主流用的是GATK

我最先学的也是GATK,似乎现在都是以GATK为金标准?其他的也顺便看看嘛
回复 支持 反对

使用道具 举报

0

主题

8

帖子

86

积分

注册会员

Rank: 2

积分
86
发表于 2017-10-25 15:58:48 | 显示全部楼层
贴主,能分享一下你搜集的关于call SNV的资料不?
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
 楼主| 发表于 2017-10-26 00:00:04 | 显示全部楼层
泡泡 发表于 2017-10-25 15:58
贴主,能分享一下你搜集的关于call SNV的资料不?

我都是看软件说明书,如果看不懂的地方就网上搜一下,所以没资料额。。因为我用过call SNV的软件也不多。。不好意思了。。
回复 支持 反对

使用道具 举报

13

主题

44

帖子

522

积分

高级会员

Rank: 4

积分
522
发表于 2018-1-22 10:08:53 | 显示全部楼层
想知道现在比较主流的一些call SNV的软件的核心算法是不是一样的,像GATK、strelka之类的
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 14:15 , Processed in 0.036418 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.