☰

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

関連情報