解析VCF文件的INFO字段

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

输入：

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。

谢谢你的任何想法。

来源

2017-10-08 Geroge

@EdMorton我只是想复制所有文本以更好地解释输入文件。输入和输出是明确的 - 提取列由制表符2,3,4,5,8分隔，然后从最后一列提取仅CLNSIG + CLNDBN及其值。 – Geroge

我建议使用perl或python等语言的小脚本来做到这一点。然后很容易将最后一个字段拆分为“;”并使用生成的元素来构建一个“词典”（用python词汇表），使您可以选择您感兴趣的信息。 – bli

随着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开始）
见http://perldoc.perl.org/perlrun.html#Command-Switches的详细信息，-lane选项

来源

2017-10-08 15:11:54 Sundeep

纯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 ${[email protected]} ${[email protected]} 
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壳，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

来源

2017-10-08 15:15:31 agc

它可以用awk来完成：

脚本。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中的第一个子表达式（第一个圆括号）放在pos 1，将第二个子表达式放置在pos 2等等。

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

来源

2017-10-08 18:19:20

解析VCF文件的INFO字段

回答

相关问题