根据其他列的内容从VCF中提取可变位置

Question

我有一个vcf文件，我试图从这些列中提取信息：

#CHROM  POS   REF     ALT

但是，只有当示例-1列包含字符串DeNovo (而不是DeNovoSV)，而该示例-1、SAMPLE-2和SAMPLE-3都包含PASS时，我才想提取它们。

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SAMPLE-1  SAMPLE-2  SAMPLE-3
chr1    10230   .       AC      A       186.90  .       AC=4;AF=0.667;AN=6;DP=77;FS=0.000;MQ=26.38;MQRankSum=0.436;QD=3.89;ReadPosRankSum=0.000;SOR=1.162       GT:AD:AF:DP:GQ:FT:F1R2:F2R1:PL:GP:PP
:DN 0/1:6,12:0.667:18:32:PASS:2,3:4,9:69,0,30:3.9669e+01,2.7888e-03,3.2724e+01:295,0,215:Inherited  0/1:5,15:0.750:20:11:PASS:3,6:2,9:60,0,8:3.0340e+01,3.6694e-01,1.0964e+01:172,0,137:.   1/1:0,10
:1.000:10:26:PASS:0,2:0,8:93,29,0:6.1212e+01,2.6342e+01,1.0101e-02:378,0,183:.
chr1    61871   .       C       CT      60.27   .       AC=3;AF=0.500;AN=6;DP=29;FS=11.290;MQ=33.00;MQRankSum=-0.423;QD=2.51;ReadPosRankSum=0.705;SOR=0.478     GT:AD:AF:DP:GQ:FT:F1R2:F2R1:PL:GP:PP:DPL:DN:DQ  0/0:5,0:0.000:5:15:PASS:.:.:0,15,182:.:0,7,93:0,15,100:DeNovo:2.9227e-07        1/1:0,2:1.000:2:5:PASS:0,1:0,1:42,6,0:2.4787e+01,4.7870e+00,1.7754e+00:29,0,9:40,6,0:.  0/1:15,7:0.318:22:26:PASS:6,3:9,4:43,0,41:2.6538e+01,9.8206e-03,4.4010e+01:65,0,250:74,0,234:.
chr1    66369   .       TA      T       116.77  .       AC=2;AF=0.500;AN=4;DP=56;FS=10.138;MQ=173.59;MQRankSum=1.468;QD=4.32;ReadPosRankSum=0.929;SOR=0.367     GT:AD:AF:DP:GQ:FT:F1R2:F2R1:PL:GP:PP    ./.:11,5:0.312:16:0:LowGQ:.:.   0/1:8,4:0.333:12:40:PASS:3,2:5,2:71,0,43:4.1762e+01,4.0824e-04,4.5625e+01:71,0,43       0/1:8,7:0.467:15:45:PASS:4,4:4,3:77,0,47:4.7400e+01,1.2244e-04,5.0000e+01:77,0,47
chr1    934273  .       G       C       8.67    .       AC=1;AF=0.167;AN=6;DP=26;FS=0.000;MQ=19.17;MQRankSum=-1.179;QD=0.96;ReadPosRankSum=1.666;SOR=0.223      GT:AD:AF:DP:GQ:FT:F1R2:F2R1:PL:GP:PP
:DPL:DN:DQ  0/1:7,2:0.222:9:11:PASS:5,1:2,1:45,0,32:1.0868e+01,3.7242e-01,3.5372e+01:45,0,60:45,0,32:DeNovoSV:4.3945e-09    0/0:6,3:0.333:9:0:LowGQ:.:.:0,0,140:.:55,0,191:46,0,186:.       0/0:8,0:
0.000:8:23:PASS:.:.:0,23,190:.:0,25,195:0,23,190:.
chr1    934274  .       G       C       8.68    .       AC=1;AF=0.167;AN=6;DP=26;FS=0.000;MQ=19.17;MQRankSum=-1.179;QD=0.96;ReadPosRankSum=1.666;SOR=0.223      GT:AD:AF:DP:GQ:FT:F1R2:F2R1:PL:GP:PP
:DPL:DN:DQ  0/1:7,2:0.222:9:11:PASS:5,1:2,1:45,0,32:1.0868e+01,3.7242e-01,3.5372e+01:45,0,60:45,0,32:DeNovoSV:4.3945e-09    0/0:6,3:0.333:9:0:PASS:.:.:0,0,140:.:55,0,191:46,0,186:.       0/0:8,0:
0.000:8:23:PASS:.:.:0,23,190:.:0,25,195:0,23,190:.

我试过使用bcftools，见下文。

bcftools query -f '%CHROM %POS %REF %ALT\n' file.vcf | head -3
chr1 10230 AC A
chr1 61871 C CT
chr1 66369 TA T

有没有一种方法可以使用bcftools，或者将bcftools与awk结合，以获得我正在寻找的vcf文件格式的输出？

非常感谢

Answer 1

FWIW，我不熟悉bcftools，但如果目的是将bcftools输出导入awk，那么我们就可以在awk中完成全部工作。

假设：

没有一个字段(在文件中)包含空白，ALT

• these

• 总是要为列( #CHROM / POS / REF /4)提取数据，列总是位于相同的字段/位置(#CHROM = $1 / POS = $2 / REF = $4 / ALT = $5)
• DeNovo总是出现在两个冒号之间(即，我们正在寻找字段中的$10)
• the，三个测试字段总是出现在相同的字段/位置(SAMPLE_1 = $10 / SAMPLE_2 = $11 / SAMPLE_3 = $12)
• PASS总是出现在两个冒号之间(即，我们在$10 / $11 / $12)

字段中寻找:PASS: )

awk的一个想法是：

awk '
$10 ~ ":DeNovo:" &&
$10 ~ ":PASS:"   &&
$11 ~ ":PASS:"   &&
$12 ~ ":PASS:"      { print $1,$2,$4,$5 }
' file.vcf

这就产生了：

chr1 61871 C CT

注：显而易见(?)这种方法的缺点是我们已经硬编码了列引用；我们当然可以修改awk代码以提供一个更动态的接口(复制bcftools功能？)但我不确定增加的复杂性是否值得努力，ymmv。

假设OP需要bcftools的功能(例如，指定一组变量列)，一种想法是修改当前的bcftools调用以包含SAMPLE_X列，然后将输出输送到awk；类似于：

bcftools query -f '%CHROM %POS %REF %ALT %SAMPLE_1 %SAMPLE_2 %SAMPLE_3\n' file.vcf | awk '
$5 ~ ":DeNovo:" &&
$5 ~ ":PASS:"   &&
$6 ~ ":PASS:"   &&
$7 ~ ":PASS:"      { print $1,$2,$3,$4 }'

备注：

• 的进一步改进将包括参数化DeNovo和PASS bcftools，返回到消除bcftools的想法，并在awk中完成全部工作。可能需要两个列列表(要显示的列，要测试的列)以及字符串列表，以便将每个“test”列与.可行，但现在我们正在进行相当大的编码工作(相对于我们迄今所做的)