使用 AWK 将两个不同的 gsub 结果输出到两列

Question

我正在尝试优化 awk 代码以在我的 VCF 文件中生成新列。 这是输入（最小工作示例）：

#CHROM  POS GENE    REF ALT QUAL    FILTER  INFO    FORMAT  1   2   3   4   5   6
chr14   33925904    EGLN3   A   T   .   .   CSQ=T|  GT:AQ:DP:AD:AB:GQ:PL    0/1:.:0:.:0,0:.:.   1/1:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.
chr14   33925905    EGLN3   G   C   .   .   CSQ=C|  GT:AQ:DP:AD:AB:GQ:PL    ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   0/0:.:0:.:0,0:.:.   1/0:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.
chr14   33925907    EGLN3   G   A   .   .   CSQ=A|  GT:AQ:DP:AD:AB:GQ:PL    ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   0/0:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.
chr14   33925918    EGLN3   T   G   .   .   CSQ=G|  GT:AQ:DP:AD:AB:GQ:PL    ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   1/0:.:0:.:0,0:.:.
chr14   33927014    EGLN3   A   G   .   .   CSQ=G|  GT:AQ:DP:AD:AB:GQ:PL    ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   1/0:.:0:.:0,0:.:.
chr14   33927025    EGLN3   AT  A   .   .   CSQ=-|  GT:AQ:DP:AD:AB:GQ:PL    ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.
chr14   33929079    EGLN3   G   A   .   .   CSQ=A|  GT:AQ:DP:AD:AB:GQ:PL    ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   ./.:.:0:.:0,0:.:.   0/0:.:0:.:0,0:.:.

我使用这个代码：

paste <(bcftools view $MYVCF | awk -F"\t" 'BEGIN {print "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"} !/^#/ {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8}') \ #prints first columns from vcf - Important if file has a vcf header
  <(bcftools query -f '[\t%SAMPLE=%GT]\n' "$MYVCF" |\
    awk 'BEGIN {print "nHet"} {print gsub(/0\|1|1\|0|0\/1|1\/0/, "")}') \ #counts samples (1-6) that have 0/1 or 1/0
    \
  <(bcftools query -f '[\t%SAMPLE=%GT]\n' "$MYVCF" |\
    awk 'BEGIN {print "nHomAlt"} {print gsub(/1\|1|1\/1/, "")}') \
    \ #counts samples with 1/1
  <(bcftools query -f '[\t%SAMPLE=%GT]\n' "$MYVCF" |\
    awk 'BEGIN {print "nHomRef"} {print gsub(/0\|0|0\/0/, "")}')   >> example3.txt

  Which does the following:

首先，它打印 CHR POS ID REF ALT 和 INFO 列，然后在样本 (1-6) 中搜索具有特定 0/0 的 GT 值，并在名为 nHomRef 的新列中打印具有此特定 GT 的样本计数。对 1/0（和组合）执行相同操作并输出到名为 nHet 的列，然后最后对 1/1 组合执行相同操作并输出到名为 hHomRef

的新列

这是输出：

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    nHet    nHomAlt nHomRef
chr14   33925904        EGLN3   A       T       .       .       CSQ=T|  1       1       0
chr14   33925905        EGLN3   G       C       .       .       CSQ=C|  1       0       1
chr14   33925907        EGLN3   G       A       .       .       CSQ=A|  0       0       1
chr14   33925918        EGLN3   T       G       .       .       CSQ=G|  1       0       0
chr14   33927014        EGLN3   A       G       .       .       CSQ=G|  1       0       0
chr14   33927025        EGLN3   AT      A       .       .       CSQ=-|  0       0       0
chr14   33929079        EGLN3   G       A       .       .       CSQ=A|  0       0       1

我正在尝试优化代码，因此不会创建很多子流程。 为了简化我只使用了两列

paste <(bcftools view "$MYVCF" | awk -F"\t" 'BEGIN {print "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"} !/^#/ {print $1"\t"$2"\t"$3"\t"$4"\t"$5}') <(bcftools query -f '[\t%SAMPLE=%GT]\n' "$MYVCF" | awk 'BEGIN {print "nHomAlt\tnHet"} {print gsub(/1\|1|1\/1/, "")} {print gsub(/0\|1|1\|0|0\/1|1\/0/, "")}') >> example3_test.txt

如何将 awk 内的 gsub 输出到正确的列？ 所以，这个片段

awk 'BEGIN {print "nHomAlt\tnHet"} {print gsub(/1\|1|1\/1/, "")} {print gsub(/0\|1|1\|0|0\/1|1\/0/, "")}'

输出正确的列？ 我得到：

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    nHomAlt nHet
chr14   33925904    EGLN3   A   T   .   .   CSQ=T|  1
chr14   33925905    EGLN3   G   C   .   .   CSQ=C|  1
chr14   33925907    EGLN3   G   A   .   .   CSQ=A|  0
chr14   33925918    EGLN3   T   G   .   .   CSQ=G|  1
chr14   33927014    EGLN3   A   G   .   .   CSQ=G|  0
chr14   33927025    EGLN3   AT  A   .   .   CSQ=-|  0
chr14   33929079    EGLN3   G   A   .   .   CSQ=A|  0
    1
    0
    1
    0
    0
    0
    0

Answer 1

首先，任何时候你的目标包括“优化”，你的首要努力应该是消除重复的任何工作......例如重复

bcftools query -f '[\t%SAMPLE=%GT]\n' "$MYVCF" | awk 

您的部分问题是当您已经在

paste

中拥有数据时，尝试将数据推送到

awk

。我认为你真正需要做的是记住

print

将附加一个输出记录分隔符或 ORS，默认情况下它是一个换行符。您可以使用

printf

输出每个字段，但我认为将字段作为参数传递会更容易。您还可以设置 OFS（输出字段分隔符）以默认使用制表符分隔它们。

根据我存储在文件中的输入示例，您可能只需要管理数据收集和格式设置。这是一个可能的例子：

$: awk 'BEGIN {OFS="\t"; print "nHomAlt\tnHet"} {
 nHet=gsub(/0\|1|1\|0|0\/1|1\/0/, "");
 nHomAlt=gsub(/1\|1|1\/1/, "");
 print nHomAlt,nHet;
}' input
nHomAlt nHet
1       3
0       3
3       2
0       1
0       3
0       0
0       7
0       0

你说你想要计数，这就是

gsub

返回的内容，所以我认为这就是你使用它的原因。

你想要的可能看起来像这样：

awk 'BEGIN { OFS="\t"; print "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tnHet\tnHomAlt\tnHomRef" } {
 nHet=gsub(/0[|/]1|1[|/]0/, "");
 nHomAlt=gsub(/1[|/]1/, "");
 nHomRef=gsub(/0[|/]0/, "");
 print $1,$2,$3,$4,$5,$6,$7,$8,nHet,nHomAlt,nHomRef;
}' input

...您显然已经知道如何用命令替换输入。希望这能让你足够接近。