BCFTOOLS支持的表达式

翻译https://samtools.github.io/bcftools/bcftools-man.html#expressions

有效的表达式可能包含

  • 数字常量,字符常量,文件名
  • 1, 1.0, 1e-4
    "String"
    @file_name
  • 算数运算符:
  • +,*,-,/
  • 比较运算符:
  • == (同=), >, >=, <=, <, !=
  • 正则运算符"~"及其非值"!~"。大小写敏感,加"/i"这大小写不敏感:
  • INFO/HAYSTACK ~ "needle"
    INFO/HAYSTACK ~ "NEEDless/i"
  • 括号:
  • (, )
  • 逻辑运算符:
  • && (同&), ||,  |
  • VCF/BCF中的INFO tags, FORMAT tags, column names:
  • INFO tags: INFO/DP or DP
    FORMAT tags: FORMAT/DV, FMT/DV, or DV
    column names:  FILTER, QUAL, ID, POS, REF, ALT[0]
  • 1 (或0) 检测flag的存在 (或不存在) :
  • FlagA=1 && FlagB=0

  • "." 检测缺失值:
  • DP=".", DP!=".", ALT="."
  • 不论分型phase或者倍体ploid(".|.", "./.", ".") ,缺失的基因型genotype可以通过以下表达式匹配:
  • GT~"\.", GT!~"\."
  • 包括分型phase或者倍体ploid(".|.", "./.", ".") ,可以通过以下表达式匹配:
  • GT=".|.", GT="./.", GT="."
  • 样本基因型genotype:参考
  • 参考基因型(单倍体或二倍体):  GT="ref"
    突变基因型 (纯合或杂合,单倍体或二倍体): GT="alt"
    缺失基因型:GT="mis"
    纯合:GT="hom"
    杂合:GT="het"
    单倍体: GT="hap"
    纯合参考:GT="RR"
    纯合突变: GT="AA"
    ref-alt杂合: GT="RA" or GT="AR"
    alt-alt杂合:GT="Aa" or GT="aA"
    haploid ref: GT="R"
    haploid alt:GT="A"
    大小写不敏感
  • 在REF和ALT列中的突变类型:(indel,snp,mnp,ref,bnd(断点breakend?),other),"\~" 表示至少一个等位基因是这种类型,"=" 要求所有等位基因是这个类型:
  • TYPE="snp"
    TYPE~"snp"
    TYPE!="snp"
    TYPE!~"snp"
  • 数组下标(0开始), "*" 任意元素, "-" 表示范围。在检索FORMAT向量时,冒号":" 可以用来选择样本和元素,如下:
  • INFO/AF[0] > 0.3             ..   第一个AF值大于0.3
    FORMAT/AD[0:0] > 30          ..   第一个样本的第一个AD值大于30
    FORMAT/AD[0:1]               ..   第一个样本的第二个AD值
    FORMAT/AD[1:0]               ..   第二个样本的第一个AD值
    DP4[*] == 0                  ..   任意DP4的值
    FORMAT/DP[0]   > 30         ..   第一个样本的DP值大于30
    FORMAT/DP[1-3] > 10         ..  第二到四的样本
    FORMAT/DP[1-]  < 7          ..    除了第一个样本的其他所有样本
    FORMAT/DP[0,2-4] > 20       ..   第一,三到五个样本
    FORMAT/AD[0:1]              ..   第一个样本的第二个AD field 
    FORMAT/AD[0:*], AD[0:] or AD[0] ..  第一个样本得到所有AD field
    FORMAT/AD[*:1] or AD[:1]        ..   所有样本的AD field
    (DP4[0]+DP4[1])/(DP4[2]+DP4[3]) > 0.3
    CSQ[*] ~ "missense_variant.*deleterious"
  • 基于FORMAT tags和INFO tags的函数:
  • MAX, MIN, AVG, SUM, STRLEN, ABS, COUNT
  • 变量可以快速计算得到:
  • alternate allele的数目:N_ALT
    样本数目:N_SAMPLES
    alternate allele的的计数:AC
    次等位基因计数,与AC相似,但统计的是小于0.5的等位基因):MAC,
    alternate等位基因的频率(AF=AC/AN):AF
    次等位基因频率(MAF=MAC/AN): MAF
    在called genotypes中等位基因的数目: AN
    有missing genotype的样本数目:N_MISSING
    有missing genotype的样本比例:F_MISSING
  • 配置perl过滤。注意这个命令默认是没有编译的,在安装的时候可以查看INSTALL文件的Optional Compilation with Perl部分,misc/demo-flt.pl 是一个示例,这个示例中定义了一个名为severity的perl子程序,并被以下命令行调用:
  • perl:path/to/script.pl; perl.severity(INFO/CSQ) > 3

    注意

    字符和正则表达式都是不区分大小写的,除了tag名外,变量和函数名也是不区分大小写的。例如"qual" 等同 "QUAL","strlen()"等同"STRLEN()" ,但是"dp" 不等同 "DP"。
    当检索多个值时,所有的元素都会被测试,通过或逻辑(OR)得到结果。比如检索"TAG=1,2,3,4"时:

    -i 'TAG[*]=1'   .. true, 真,该条记录会被打印
    -i 'TAG[*]!=1'  .. true
    -e 'TAG[*]=1'   .. false,否,该条打印会被丢弃
    -e 'TAG[*]!=1'  .. false
    -i 'TAG[0]=1'   .. true
    -i 'TAG[0]!=1'  .. false
    -e 'TAG[0]=1'   .. false
    -e 'TAG[0]!=1'  .. true

    示例

    MIN(DV)>5
    MIN(DV/DP)>0.3
    MIN(DP)>10 & MIN(DV)>3
    FMT/DP>10  & FMT/GQ>10 ..   同一个样本需满足两个条件 
    FMT/DP>10 && FMT/GQ>10 ..  允许不同样本间可以满足这两个条件
    QUAL>10 |  FMT/GQ>10   ..   对QUAL>10 的位点或者 FMT/GQ>10的样本为真,但只选择FMT/GQ>10的样本
    QUAL>10 || FMT/GQ>10   .. 对QUAL>10 的位点或者 FMT/GQ>10的样本为真,选择该位点的所有样本
    TYPE="snp" && QUAL>=10 && (DP4[2]+DP4[3] > 2)
    COUNT(GT="hom")=0
    MIN(DP)>35 && AVG(GQ)>50
    ID=@file       ..   选择出现在文件file中的ID
    ID!=@~/file    ..   忽略出现在文件file中的ID
    MAF[0]<0.05    ..   以5%为阈值,选择rare variant
    POS>=100   ..   限制检索范围,比如20:100-200内的突变

    Shell 扩展

    表达式需要用单引号引起来,因为一些字符可能在shell中有特殊含义。比如

    bcftools view -i '%ID!="." & MAF[0]<0.01' 

    用bcftools进行突变筛选的示例

    摘自公众号生信开发者
    #选出位于bed文件中的所有区域的突变位点

    bcftools view -R target.bed sample.filter.vcf.bgz 

    #选出INFO中SAO=1的所有位点

    bcftools view -i 'SAO=1' sample.filter..b37.vcf 

    #选出INFO中AC>2的所有位点

    bcftools view -i "AC>=2" ExAC.r1.sites.vep.vcf.gz 

    #选出遗传方式是AR或XR的位点(需INFO字段中已有INHERITANCE注释)

    bcftools view -i 'INHERITANCE[*] = "AR" || INHERITANCE[*] = "XR"' sample.filter.annovar.hg19_multianno.Fun.ar.vcf 

    #选出位于区域X:31136335-33358725的所有位点

    bcftools view sample.filter.vcf.bgz X:31136335-33358725

    #选出除了INFO字段CLINSIG匹配"Benign"以外的所有

    bcftools view -e 'CLINSIG~"Benign"' sample.filter.hg19_multianno.Fun.vcf 

    #选出最小的depth>500而且,肿瘤样品的VAF>0.05的所有位点

    bcftools view -i " MIN(FMT/DP)>500 && FORMAT/AD[1:1]/FORMAT/DP[1]>0.05 " annovar.exdbsnp.exCOSMIC.exTSG.cds.vcf | sed '/^#/d' 

    参考

    表达式手册翻译自https://samtools.github.io/bcftools/bcftools-man.html#expressions
    bcftools突变筛选摘自公众号生信开发者

    发表评论

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

    This site uses Akismet to reduce spam. Learn how your comment data is processed.