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.digitsORAMR_AF=digits，每当此正则表达式在线获得匹配项时，它都会设置RSTART并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 从上到下处理每个输入行的脚本，其中每个脚本语句都有

[条件] [{动作}]

如方括号所示，这两个部分都是可选的 - 如果缺少条件，则无条件执行操作;如果缺少操作，则默认为{ 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手册中的更多选项，请参阅。