如何用awk提取VCF文件的ID列与INFO列中的AN值

输入VCF文件内容 #

##bcftools_annotateVersion=1.3.1+htslib-1.3.1
##bcftools_annotateCommand=annotate 
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HG005
chr1    817186  rs3094315   G   A   50  PASS    platforms=2;platformnames=Illumina,CG;datasets=3;datasetnames=HiSeq250x250,CGnormal,HiSeqMatePair;callsets=5;callsetnames=HiSeq250x250Sentieon,CGnormal,HiSeq250x250freebayes,HiSeqMatePairSentieon,HiSeqMatePairfreebayes;datasetsmissingcall=IonExome,SolidSE75bp;callable=CS_HiSeq250x250Sentieon_callable,CS_CGnormal_callable,CS_HiSeq250x250freebayes_callable;AN=2;AF=1;AC=2 GT:PS:DP:ADALL:AD:GQ    1/1:.:809:0,363:78,428:237
chr1    817341  rs3131972   A   G   50  PASS    platforms=3;platformnames=Illumina,CG,Solid;datasets=4;datasetnames=HiSeq250x250,CGnormal,HiSeqMatePair,SolidSE75bp;callsets=6;callsetnames=HiSeq250x250Sentieon,CGnormal,HiSeq250x250freebayes,HiSeqMatePairSentieon,HiSeqMatePairfreebayes,SolidSE75GATKHC;datasetsmissingcall=IonExome;callable=CS_HiSeq250x250Sentieon_callable,CS_CGnormal_callable,CS_HiSeq250x250freebayes_callable;AN=2;AF=1;AC=2   GT:PS:DP:ADALL:AD:GQ    1/1:.:732:1,330:99,391:302

期望输出 #

ID         INFO
rs3094315   2
rs3131972   2

问题分析 #

你之前尝试的awk '/^[^#]/ { print $3, gsub(/^[^AN=])/,"",$8)}' file.vcf命令存在两处核心问题：

• 正则表达式语法错误，/^[^AN=])/多了冗余的右括号，且逻辑无法匹配AN=目标内容
• gsub是替换操作，无法直接提取目标值，且未正确定位AN=后的数字

正确解决方案 #

方案一：兼容多数awk版本 #

awk 'BEGIN{printf "%-10s %s\n", "ID", "INFO"} /^[^#]/{split($8, info_arr, /;/); for(i in info_arr){if(info_arr[i] ~ /^AN=/){sub(/^AN=/,"",info_arr[i]); printf "%-10s %s\n", $3, info_arr[i]}}}' file.vcf

命令说明：

1. BEGIN块：提前打印格式化的表头，保证输出对齐
2. 匹配非注释行：/^[^#]/跳过所有以#开头的注释行
3. 拆分INFO字段：用split($8, info_arr, /;/)把INFO字段按分号拆分成数组
4. 提取AN值：遍历数组找到以AN=开头的元素，用sub去掉前缀后，格式化输出ID和对应值

方案二：简洁正则捕获（需awk支持数组捕获） #

awk 'BEGIN{print "ID         INFO"} /^[^#]/{match($8, /AN=([0-9]+)/, an_val); printf "%-10s %s\n", $3, an_val[1]}' file.vcf

命令说明：

用match函数直接匹配AN=后的数字并捕获到数组，代码更简洁，适合较新版本的awk。

内容的提问来源于stack exchange，提问作者Anna