翻译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中有特殊含义。比如

1
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突变筛选摘自公众号生信开发者