Bash

使用awk根據列中的子字元串提取行

  • November 5, 2022

我有一個製表符分隔的 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

我想將行提取到一個新的文本文件中,其中INFO列中的 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 的地方提取它。將不勝感激。

不要為此使用 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

引用自:https://unix.stackexchange.com/questions/723752