解析 VCF 文件的 INFO 字段答案

【问题标题】：Parse VCF file's INFO field解析 VCF 文件的 INFO 字段
【发布时间】：2017-10-08 14:47:48
【问题描述】：

请帮我解析一个 VCF 文件。我正在粘贴一个真实的例子。

输入：

1   1014143 rs786201005 C   T   .   .   RS=786201005;RSPOS=1014143;dbSNPBuildID=144;SSR=0;SAO=1;VP=0x050068000605000002110100;GENEINFO=ISG15:9636;WGT=1;VC=SNV;PM;PMC;NSN;REF;ASP;LSD;OM;CLNALLE=1;CLNHGVS=NC_000001.11:g.1014143C>T;CLNSRC=OMIM_Allelic_Variant;CLNORIGIN=1;CLNSRCID=147571.0003;CLNSIG=5;CLNDSDB=MedGen:OMIM:Orphanet;CLNDSDBID=C4015293:616126:ORPHA319563;CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification;CLNREVSTAT=no_criteria;CLNACC=RCV000162196.3
1   1014228 rs1921  G   A,C .   .   RS=1921;RSPOS=1014228;dbSNPBuildID=36;SSR=0;SAO=0;VP=0x050328000a0517053f000100;GENEINFO=ISG15:9636;WGT=1;VC=SNV;PM;PMC;S3D;SLO;NSM;REF;ASP;VLD;G5A;G5;HD;GNO;KGPhase1;KGPhase3;CLNALLE=1;CLNHGVS=NC_000001.11:g.1014228G>A;CLNSRC=.;CLNORIGIN=1;CLNSRCID=.;CLNSIG=2;CLNDSDB=MedGen;CLNDSDBID=CN169374;CLNDBN=not_specified;CLNREVSTAT=single;CLNACC=RCV000455759.1;CAF=0.6611,0.3389,.;COMMON=1
1   1014316 rs672601345 C   CG  .   .   RS=672601345;RSPOS=1014319;dbSNPBuildID=142;SSR=0;SAO=1;VP=0x050068001205000002110200;GENEINFO=ISG15:9636;WGT=1;VC=DIV;PM;PMC;NSF;REF;ASP;LSD;OM;CLNALLE=1;CLNHGVS=NC_000001.11:g.1014319dupG;CLNSRC=OMIM_Allelic_Variant;CLNORIGIN=1;CLNSRCID=147571.0002;CLNSIG=5;CLNDSDB=MedGen:OMIM:Orphanet;CLNDSDBID=C4015293:616126:ORPHA319563;CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification;CLNREVSTAT=no_criteria;CLNACC=RCV000148989.5
1   1014359 rs672601312 G   T   .   .   RS=672601312;RSPOS=1014359;dbSNPBuildID=142;SSR=0;SAO=1;VP=0x050068000605000002110100;GENEINFO=ISG15:9636;WGT=1;VC=SNV;PM;PMC;NSN;REF;ASP;LSD;OM;CLNALLE=1;CLNHGVS=NC_000001.11:g.1014359G>T;CLNSRC=OMIM_Allelic_Variant;CLNORIGIN=1;CLNSRCID=147571.0001;CLNSIG=5;CLNDSDB=MedGen:OMIM:Orphanet;CLNDSDBID=C4015293:616126:ORPHA319563;CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification;CLNREVSTAT=no_criteria;CLNACC=RCV000148988.5
1   1020183 rs539283387 G   C   .   .   RS=539283387;RSPOS=1020183;dbSNPBuildID=142;SSR=0;SAO=0;VP=0x050000000a05040026000100;GENEINFO=AGRN:375790;WGT=1;VC=SNV;NSM;REF;ASP;VLD;KGPhase3;CLNALLE=1;CLNHGVS=NC_000001.11:g.1020183G>C;CLNSRC=.;CLNORIGIN=1;CLNSRCID=.;CLNSIG=3;CLNDSDB=MedGen;CLNDSDBID=CN169374;CLNDBN=not_specified;CLNREVSTAT=single;CLNACC=RCV000424799.1;CAF=0.9904,0.009585;COMMON=1
1   1020216 rs764659938 C   G   .   .   RS=764659938;RSPOS=1020216;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000a05040002000100;GENEINFO=AGRN:375790;WGT=1;VC=SNV;NSM;REF;ASP;VLD;CLNALLE=1;CLNHGVS=NC_000001.11:g.1020216C>G;CLNSRC=.;CLNORIGIN=1;CLNSRCID=.;CLNSIG=0;CLNDSDB=MedGen;CLNDSDBID=CN221809;CLNDBN=cancer;CLNREVSTAT=single;CLNACC=RCV000422793.1

我需要一个输出：

1014143 rs786201005 C   T   CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1014228 rs1921  G   A,C CLNSIG=2    CLNDBN=not_specified
1014316 rs672601345 C   CG  CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1014359 rs672601312 G   T   CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1020183 rs539283387 G   C   CLNSIG=3    CLNDBN=not_specified
1020216 rs764659938 C   G   CLNSIG=0    CLNDBN=not_provided

这意味着打印列 2,3,4,5，然后解析最后一列并仅打印 CLNSIG 和 CLNDBN。问题是，这些值并不总是在同一个位置。

我的尝试是：

awk -v OFS="\t"'{print $2,$3,$4,$5,$8}' input

...然后我不知道如何获得 CLNSIG 和 CLNDBN。

感谢您的任何想法。

【问题讨论】：

@EdMorton 我只是想复制所有文本以更好地解释输入文件。输入和输出是明确的 - 提取由制表符 2,3,4,5,8 分隔的列，然后从最后 8 列中仅提取 CLNSIG+CLNDBN 及其值。
我建议使用 perl 或 python 等语言的小脚本来执行此操作。然后很容易拆分“;”上的最后一个字段并使用生成的元素构建一个“字典”（在 python 词汇表中），允许您选择您感兴趣的信息。

标签： bash awk sed bioinformatics vcf-variant-call-format

【解决方案1】：

纯bash，使用bash 解析剩余变量在$h 中，带有parameter tranformation 输出：

while read a b c d e f g h ; do 
    declare ${h//;/ }
    printf "%s\t%-10s\t%s\t%s\t%s\t%s\n" $b $c $d $e ${CLNSIG@A} ${CLNDBN@A}
done < input

输出：

1014143 rs786201005 C   T   CLNSIG='5'  CLNDBN='Immunodeficiency_38_with_basal_ganglia_calcification'
1014228 rs1921      G   A,C CLNSIG='2'  CLNDBN='not_specified'
1014316 rs672601345 C   CG  CLNSIG='5'  CLNDBN='Immunodeficiency_38_with_basal_ganglia_calcification'
1014359 rs672601312 G   T   CLNSIG='5'  CLNDBN='Immunodeficiency_38_with_basal_ganglia_calcification'
1020183 rs539283387 G   C   CLNSIG='3'  CLNDBN='not_specified'
1020216 rs764659938 C   G   CLNSIG='0'  CLNDBN='cancer'

POSIX shell，grep 和 printf 方法：

while read a b c d e f g h ; do
    printf "%s\t%-10s\t%s\t%s\t%s\t%s\n" $b $c $d $e \
          $( echo "$h" | grep -o 'CLN\(SIG\|DBN\)=[^;]*' ) ; 
done < input

输出：

1014143 rs786201005 C   T   CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1014228 rs1921      G   A,C CLNSIG=2    CLNDBN=not_specified
1014316 rs672601345 C   CG  CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1014359 rs672601312 G   T   CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1020183 rs539283387 G   C   CLNSIG=3    CLNDBN=not_specified
1020216 rs764659938 C   G   CLNSIG=0    CLNDBN=cancer

【讨论】：

【解决方案2】：

与perl

$ perl -lane 'print join "\t",(@F[1..4], /(?:CLNSIG|CLNDBN)=[^;]+/g)' ip.txt
1014143 rs786201005 C   T   CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1014228 rs1921  G   A,C CLNSIG=2    CLNDBN=not_specified
1014316 rs672601345 C   CG  CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1014359 rs672601312 G   T   CLNSIG=5    CLNDBN=Immunodeficiency_38_with_basal_ganglia_calcification
1020183 rs539283387 G   C   CLNSIG=3    CLNDBN=not_specified
1020216 rs764659938 C   G   CLNSIG=0    CLNDBN=cancer

-a 用于在空白处拆分输入的选项，保存在 @F 数组中
/(?:CLNSIG|CLNDBN)=[^;]+/g 将返回 CLNSIG 和 CLNDBN 字段
@F[1..4] 给出第 2 到第 5 个字段（索引从 0 开始）
有关-lane 选项的详细信息，请参阅http://perldoc.perl.org/perlrun.html#Command-Switches

【讨论】：

【解决方案3】：

可以使用 awk 来完成：

script.awk

BEGIN { OFS="\t" }
      { clnsig = clndbn = ""
        if( match( $8, /CLNSIG=[^;]+/ ) ) {
          clnsig = substr( $8, RSTART, RLENGTH )
        }
        if( match( $8, /CLNDBN=[^;]+/ ) ) {
          clndbn = substr( $8, RSTART, RLENGTH )
        }
        print $2, $3, $4, $5, clnsig, clndbn
      }

或者更紧凑，如果CLNDBN总是在CLNSIG之后：

script.awk

BEGIN { OFS="\t" }
    { match($8,/(CLNSIG=[^;]+).*(CLNDBN=[^;]+)/, tmp) 
      print $2,$3,$4,$5, tmp[1], tmp[2]
    }

函数match 匹配一个正则表达式。第一种形式设置变量RSTART 和RLENGTH，以便您可以提取带有substring 的文本。

第二种形式将第一个子表达式（第一个括号）放入数组tmp的位置1，第二个子表达式放在位置2，依此类推。

正则表达式 CLNSIG=[^;]+ 匹配文字 CLNSIG= 后跟一个子字符串，直到（但不包括）;。

【讨论】：