拆分并计算 AWK答案

【问题标题】：Split and count AWK拆分并计算 AWK
【发布时间】：2020-10-02 11:52:04
【问题描述】：

我有像这样的文件（VCF）

##fileformat=VCFv4.0
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS     ID        REF ALT    QUAL FILTER INFO    FORMAT      NA00001
Chr02   259 .   A   .   20  .   .   GT:DP:A:C:G:T:PP:GQ 0/0:1:0,1:0,0:0,0:0,0:0,26,23,75,33,33,33,47,52,49:23
Chr02   260 .   C   .   13  .   .   GT:DP:A:C:G:T:PP:GQ 0/0:1:0,0:0,1:0,0:0,0:24,0,70,17,25,49,43,25,25,44:16
Chr02   261 .   C   .   13  .   .   GT:DP:A:C:G:T:PP:GQ 0/0:1:0,0:0,1:0,0:0,0:24,0,194,18,25,49,44,25,25,45:16
Chr02   262 .   C   A   21  .   .   GT:DP:A:C:G:T:PP:GQ 0/1:1:0,0:0,1:0,0:0,0:387,0,342,348,25,368,376,25,25,368:25
Chr02   263 .   C   .   24  .   .   GT:DP:A:C:G:T:PP:GQ 0/0:2:0,0:1,1:0,0:0,0:541,0,529,495,29,556,508,29,29,499:29
Chr02   264 .   A   .   31  .   .   GT:DP:A:C:G:T:PP:GQ 0/0:2:1,1:0,0:0,0:0,0:0,280,192,317,36,36,36,178,302,219:36
Chr02   265 .   G   C   25  .   .   GT:DP:A:C:G:T:PP:GQ 0/1:2:0,0:0,0:1,1:0,0:255,414,0,328,284,29,284,29,351,29:29
Chr02   266 .   A   .   31  .   .   GT:DP:A:C:G:T:PP:GQ 0/0:2:1,1:0,0:0,0:0,0:0,281,323,440,36,36,36,209,309,315:36
Chr02   267 .   C   .   24  .   .   GT:DP:A:C:G:T:PP:GQ 0/0:2:0,0:1,1:0,0:0,0:595,0,541,481,28,567,512,28,28,512:

我需要根据行数（1000）拆分此文件，并从每个文件（拆分文件）中计算“0/1”。为此，我使用此命令拆分文件

head -10 Chromosome02.vcf| tee header subset.1.VCF >/dev/null ; awk -v header="`cat header`" -v count=1 '( (NR>10) && !( (NR-1) % 2) ) { count++ ; print header >"subset." count ".VCF";} {print $0 >>"subset." count ".VCF";}' Chromosome02.vcf | grep “0/1”

但它不起作用。有没有办法在不生成拆分文件的情况下做到这一点？

预期输出

Chromosome02
00
00
01
00
01
00

【问题讨论】：

对不起，您拆分文件的逻辑在这里不清楚。你想创建一个新的输出文件是每 1000 行的数量吗？或者是否有任何其他逻辑可以从 Input_file 中看到任何值？由于您显示的输出不是连续显示的，请在您的问题中添加清晰的细节，以便我们更好地理解它，干杯。
Without splitting the files 你的意思是每 1000 行或每 990 行获取计数，就好像你在每个块前面添加了标题一样？
如果我根据上面的脚本拆分文件，将给出 1010 行包含标题。我只想从每 1010 行中获取计数（0/1）
所以它是每 1000 行的计数，让我们忘记标题，事实是您不想实际拆分。

标签： awk vcf-variant-call-format

【解决方案1】：

tail -n +11 file | 
    awk -v n=1000 '/0\/1/{c++} NR%n==0{print c; c=0} END {if (NR%n!=0) print c}'

tail 命令将从文件中排除标题。然后我们计算包含该模式的行数，但每 N 行我们打印一次并重置计数器。

【讨论】：

【解决方案2】：

您能否根据您显示的示例尝试以下操作，并假设我们需要在每 1011 行中计算 0/1 并打印它们。

使用tail + awk`：

tail -n +11 Input_file | 
awk '
FNR%1000==0{
  if(++count==1){ print "Chromosome02" }
  print total
  total=""
}
{
  total+=gsub(/0\/1/,"&")
}
END{
  if(++count==1){ print "Chromosome02" }
  if(total){ print total }
}'

只有awk：

awk '
FNR==11{
  start=1
}
start && ++line && line%1000==0{
  if(++count==1){ print "Chromosome02" }
  print total
  total=""
}
{
  total+=gsub(/0\/1/,"&")
}
END{
  if(++count==1){ print "Chromosome02" }
  if(total){ print total }
}' Input_file

注意： 这认为 OP 的 Input_file 需要从整行中计数 0/1，以防您要检查特定字段，然后可以针对特定字段更改上述替换。

【讨论】：

您的回答很棒。但是对于像 vcf fasta sam bed gff 这样的生物输出文件；您还可以探索 bioawk（awk 的生物版本），它可以让您更好地控制。就像 Qs 一样，有 11 个评论行；但它也可以不同。此评论仅供您参考。