Bash
使用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
我想將行提取到一個新的文本文件中,其中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