☰

使用 awk 根据列中的子字符串提取行

tag-icon • tag-icon bash awk
使用 awk 根据列中的子字符串提取行

我有一个制表符分隔的 vcf 文件，格式如下

#CHROM  POS   REF   ALT       INFO
chr1    111    A    TT;C     AC=0;AN=33
chr1    111    A     G;t     AC=0;AN=100
chr1    111    G     A       AC=110;AN=51
chr2    737    T     Q       AC=99;AN=10003
chr2    888    G     G       AC=100;AN=1636

我想将行提取到一个新的文本文件中，其中 AC 在信息列大于 100，因此预期输出变为：

#CHROM  POS   REF   ALT  INFO
chr1    111    G     A   AC=110;AN=51

到目前为止我拥有的 awk 命令是：


awk 'NR==1 || /AC=[0-9][0-9][0-9]+/ && !/AC=100/'  file.vcf > output.txt

不过，我的文件很大，需要很长时间才能完成。有没有办法提取它，我指定 $5 中的 AC（即信息列）应大于 100。见解将不胜感激。

答案1

$ awk -F'[\t=]' 'NR==1 || ($6+0)>100' file
#CHROM  POS     REF     ALT     INFO
chr1    111     G       A       AC=110;AN=51

或者如果您愿意：

$ awk '{split($NF,p,/[=;]/)} NR==1 || p[2]>100' file
#CHROM  POS     REF     ALT     INFO
chr1    111     G       A       AC=110;AN=51

答案2

不要为此使用 awk。我的意思是，你可以，但是有更好的工具。如果这确实是一个有效的 VCF 文件，则如下所示：

##fileformat=VCFv4.3
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr1>
##contig=<ID=chr2>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  foo
chr1    111 .   A   TT,C    100 PASS    AC=0;AN=33  GT  0/1
chr1    111 .   A   G,t 100 PASS    AC=0;AN=100 GT  0/1
chr1    111 .   G   A   100 PASS    AC=110;AN=51    GT  0/1
chr2    737 .   T   Q   100 PASS    AC=99;AN=10003  GT  0/1
chr2    888 .   G   G   100 PASS    AC=100;AN=1636  GT  1/1

然后你可以使用bcftools：

$ bcftools view -i "AC[*]>100" foo.vcf
##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr1>
##contig=<ID=chr2>
##bcftools_viewVersion=1.16+htslib-1.16
##bcftools_viewCommand=view -i AC[*]>100 foo.vcf; Date=Sat Nov  5 12:40:53 2022
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  foo
chr1    111 .   G   A   100 PASS    AC=110;AN=51    GT  0/1

如果它不是真正的 VCF，并且正如您在问题中所示，您可以执行以下操作：

$ perl -ne '/AC=(\d+)/; print if /^#/ || $1 > 100' foo.notVcf
#CHROM  POS   REF   ALT       INFO
chr1    111    G     A       AC=110;AN=51

答案3

使用awk

$ awk '{split($NF,array,";");split(array[1],var,"=")}  NR==1 || var[2]>100'
#CHROM  POS   REF   ALT       INFO
chr1    111    G     A       AC=110;AN=51

相关内容