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进行突变筛选的示例
摘自公众号生信开发者
#选出位于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突变筛选摘自公众号生信开发者
文章作者 zzx
上次更新 2018-05-24