awk 或其他生物信息学工具来过滤 vcf

Question

我正在尝试过滤 vcf 文件中的一些行，以下是行示例：

1   10505   rs548419688 A   T   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9632;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;E
UR_AF=0;SAS_AF=0;AA=.|||;VT=SNP
1   10506   rs568405545 C   G   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9676;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;E
UR_AF=0;SAS_AF=0;AA=.|||;VT=SNP
1   10511   rs534229142 G   A   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9869;EAS_AF=0;AMR_AF=0.0014;AFR_AF=0;E
UR_AF=0;SAS_AF=0;AA=.|||;VT=SNP
1   10539   rs537182016 C   A   100 PASS    AC=3;AF=0.000599042;AN=5008;NS=2504;DP=9203;EAS_AF=0;AMR_AF=0.0014;AFR_AF=0;E
UR_AF=0.001;SAS_AF=0.001;AA=.|||;VT=SNP
1   10542   rs572818783 C   T   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9007;EAS_AF=0.001;AMR_AF=0;AFR_AF=0;EU
R_AF=0;SAS_AF=0;AA=.|||;VT=SNP

假设我想提取 AMR_AF 大于 0.5 的行，但不知道如何使用 Awk 正则表达式来完成这项工作。试过vcftools，但没用。

Answer 1

请您尝试关注一下。

awk 'match($0,/AMR_AF=[0-9]+\.[0-9]+|AMR_AF=[0-9]+/) && substr($0,RSTART+7,RLENGTH-7)>0.5'  Input_file

解释： 使用 awk 的 match 函数来匹配正则表达式 AMR_AF= digits.digits 或 AMR_AF=digits 并且每当此正则表达式在线匹配时，它就会设置 @ 987654326@ 和 RLENGTH 变量。 &&（AND 条件）检查 RSTART+7 的子字符串值是否到 RLENGTH-7 值大于 0.5，然后打印该行。

Answer 2

您可以在您选择的字段上拆分线，并检查拆分后元素的数值是否大于您的阈值。

更详细地说，在,bar= 上拆分输入yes,foo=2,bar=0.23,baz=1 将产生一个包含yes,foo=2 和0.23,baz=1 的数组。在 Awk 中，如果您将第二个元素与 0.2 进行比较，它会简单地将值的开头尽可能多地转换为数字，然后执行数字比较。

这样

awk '{ split($0, x, /[\t;]AMR_AF=/) } x[2]>0.5' file.vcf

应该做你想做的。我们将该行拆分为x，并检查x[2] 的数值。

正则表达式中的[\t;] 允许在字段名称前使用制表符或分号；完全笼统地说，也许您甚至应该使用 (^|[\t;]) 来允许匹配发生在行首。

如果你想参数化这个，也许试试

awk -v field="AMR_AF" -v thres=0.5 '{ split($0, x, "(^|[\t;])" field "=")) } x[2]>thres' file.vcf

回想一下，Awk 从上到下处理每个输入行的脚本，其中每个脚本语句都有格式

[ 条件 ] [ { 动作 } ]

如方括号所示，这两个部分都是可选的——如果缺少条件，则动作将无条件执行；如果缺少 action，则默认为 { print $0 }。所以我们的脚本会先无条件地分割行，如果x[2]大于阈值则有条件地打印出来。

GNU Awk 可以在多字符字段分隔符上进行拆分，因此您也可以使用 -F '[\t;]AMR_AF='。

awk -F '[\t;]AMR_AF=' '$2>0.5' file.vcf

Answer 3

使用bcftools：

bcftools view -i 'INFO/AMR_AF > 0.5' myFile.vcf

更多选项请参阅bcftools manuals。