如何在 vcf（变体调用格式）文件中进行上下文相关的替换（更新字段）？
How can I do context depended substitution (updating of fields) in the vcf (variant call format) files?
我有一个 vcf 文件，如下所示：
CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 2ms01e 2ms02g 2ms03g 2ms04h
2 15882505 . T A 12134.90 PASS AC=2;AF=0.250;AN=8;BaseQRankSum=-0.021;ClippingRankSum=0.000;DP=695;ExcessHet=3.6798;FS=0.523;MLEAC=2;MLEAF=0.250;MQ=60.00;MQRankSum=0.000;QD=25.18;ReadPosRankSum=1.280;SOR=0.630 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/1:59,89:148:99:3620,0,2177:1|0:.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.:1452:|:0.5 0/1:125,209:334:99:8549,0,4529:.:.:.:.:. 0/0:130,0:130:99:0,400,5809:.:.:.:.:. 0/0:82,0:82:99:0,250,3702:.:.:.:.:.
2 15882583 . G T 1221.33 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=-2.475;ClippingRankSum=0.000;DP=929;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.125;MQ=60.00;MQRankSum=0.000;QD=9.25;ReadPosRankSum=0.299;SOR=0.686 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/0:178,0:178:99:0,539,7601:0/0:.:.:0/0:. 0/0:446,0:446:99:0,1343,16290:.:.:.:.:. 0/0:172,0:172:99:0,517,6205:.:.:.:.:. 0/1:75,57:132:99:1253,0,2863:.:.:.:.:.
第一行是 header（前面还有其他 header 信息，已在此处删除），各列以制表符分隔。
为了方便理解数据结构，我在此 link 中共享文件的 sub-sample（在可以下载的保管箱上）： https://www.dropbox.com/sh/coihujii38t5prd/AABDXv8ACGIYczeMtzKBo0eea?dl=0
请下载文件，只有300Kb左右，可以用文本编辑器打开。这将有助于更好地理解问题。
问题：
我需要进行上下文相关的替换（值更新）。
- 所有 header 信息（在行前用 # 标记）需要保持不变。
不同行的值对应于最后一个 header（即 CHROM POS ID ....）
首先，我们需要查看 FORMAT 列中的 PG（阶段基因型）字段。字段的不同标签以“:”分隔。 SAMPLE 列中的特定字段有一个相应的值（目前为 2ms01e）。因此，对于第一行，样本 (2ms01e) 的 PG 值为 1|0.
- 现在，我们需要查看同一行FORMAT 列中的GT 字段，并将其对应的值更新为与PG 相同的值。即，将 0/1 更改为 1|0。保持 PG 中的顺序很重要（如果是 1|0 或 0|1，则需要准确）。
但是，如果 PG 字段的值为 0/1、0/0、1/0 或任何其他值（即有一个斜杠），则 GT 字段需要不会更改（或更新）。
最终输出：
因此，第一行数据的 GT 值应更改为：
GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/1:59,89:148:99:3620,0,2177:1|0:.....
到
GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 1|0:59,89:148:99:3620,0,2177:1|0:.....
你可以看到这里只有GT字段的值发生了变化，而其他所有字段的值都保持不变。
对于第二行，GT 值保持不变 - 即 0/0 到 0/0，因为该行的 PG 值为 0/0（与 GT 值相同，因此没有变化）。
简单方法：
我认为最好是 PG 字段的值可以 copy-pasted 到 SAMPLE (2ms01e) 列中的 GT 字段值。 GT字段值在第1位，PG字段在第6位，不同字段之间用“:”分隔。因此，我们需要做的就是用第 6 个字段的值更新第一个字段中的值。
这个简单的方法可能会奏效，因为当 PG 有斜杠“/”时，GT 也会有斜杠，而且顺序无关紧要。
但我不确定它是否适用于每一行。但是，这将是一个简单的解决方案，我至少可以检查并确保它是否有效。
硬方法：
如果 easy 方法没有按预期工作，我认为考虑每个上下文变得很重要。
上下文：
是PG字段值用竖线(|)。如果是，则需要更改。
如果 FORMAT 列中没有 PG 字段 - 则跳过它。
FORMAT 列中字段的分隔符是":"。因此，在SAMPLE 列中。因此，计算字段之间的距离很重要。 GT 和 PG 字段是第 1 和第 6 个位置。
任何类型的解决方案都值得赞赏，但如果它 python 更好，那么如果我的上下文发生变化，我可以阅读和操作它。
此外，对给定解决方案的解释会有很大帮助。
提前致谢，抱歉这么挑剔。我的计算机技能很一般，但仍然不会编程。
干杯！ :))
-
下面的答案不会包含任何逻辑的细节，但它会给你一个可能的起点，你可以自己玩：
class VCFProcessor():
def __init__(self):
pass
def load(self, filename):
with open(filename, "rb") as f:
self.load_string(f.read())
def load_string(self, data):
index = 0
for line in data.split("\n"):
# Skip empty rows
if line.strip() == "":
continue
# Assuming there is only header and valid rows
if self.is_valid_row(line):
self.process_row(index, line)
else:
self.process_header(line)
index += 1
def is_valid_row(self, line):
columns = line.split(":")
if len(columns) == 46:
return True
def process_row(self, index, line):
print "Processing line {0} with {1} columns".format(index, len(line.split(":")))
def process_header(self, line):
print "Header has {0} columns".format(len(line.split(":")))
if __name__ == "__main__":
data = """CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 2ms01e 2ms02g 2ms03g 2ms04h
2 15882505 . T A 12134.90 PASS AC=2;AF=0.250;AN=8;BaseQRankSum=-0.021;ClippingRankSum=0.000;DP=695;ExcessHet=3.6798;FS=0.523;MLEAC=2;MLEAF=0.250;MQ=60.00;MQRankSum=0.000;QD=25.18;ReadPosRankSum=1.280;SOR=0.630 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/1:59,89:148:99:3620,0,2177:1|0:.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.:1452:|:0.5 0/1:125,209:334:99:8549,0,4529:.:.:.:.:. 0/0:130,0:130:99:0,400,5809:.:.:.:.:. 0/0:82,0:82:99:0,250,3702:.:.:.:.:.
2 15882583 . G T 1221.33 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=-2.475;ClippingRankSum=0.000;DP=929;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.125;MQ=60.00;MQRankSum=0.000;QD=9.25;ReadPosRankSum=0.299;SOR=0.686 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/0:178,0:178:99:0,539,7601:0/0:.:.:0/0:. 0/0:446,0:446:99:0,1343,16290:.:.:.:.:. 0/0:172,0:172:99:0,517,6205:.:.:.:.:. 0/1:75,57:132:99:1253,0,2863:.:.:.:.:.
"""
v = VCFProcessor()
v.load_string(data)
希望对您有所帮助。
$ cat > another_mess.awk
[=10=]!="" {
n=split(,a,":") # split at ":" to a array
if(substr(a[6],2,1)=="|") { # if "|" in PG
a[1]=a[6] # copy PG to GT
="" # empty
for(i=1;i<=n;i++) # rebuild
=a[i] (i<n?":":"") # use ":" as delimiter
}
print # PRINT TO TEST, CHANGE TO [=10=]
}
$ awk -f another_mess.awk mess.in
1|0:59,89:148:99:3620,0,2177:1|0:.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.:1452:|:0.5
0/0:178,0:178:99:0,539,7601:0/0:.:.:0/0:.
我有一个 vcf 文件，如下所示：
CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 2ms01e 2ms02g 2ms03g 2ms04h
2 15882505 . T A 12134.90 PASS AC=2;AF=0.250;AN=8;BaseQRankSum=-0.021;ClippingRankSum=0.000;DP=695;ExcessHet=3.6798;FS=0.523;MLEAC=2;MLEAF=0.250;MQ=60.00;MQRankSum=0.000;QD=25.18;ReadPosRankSum=1.280;SOR=0.630 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/1:59,89:148:99:3620,0,2177:1|0:.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.:1452:|:0.5 0/1:125,209:334:99:8549,0,4529:.:.:.:.:. 0/0:130,0:130:99:0,400,5809:.:.:.:.:. 0/0:82,0:82:99:0,250,3702:.:.:.:.:.
2 15882583 . G T 1221.33 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=-2.475;ClippingRankSum=0.000;DP=929;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.125;MQ=60.00;MQRankSum=0.000;QD=9.25;ReadPosRankSum=0.299;SOR=0.686 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/0:178,0:178:99:0,539,7601:0/0:.:.:0/0:. 0/0:446,0:446:99:0,1343,16290:.:.:.:.:. 0/0:172,0:172:99:0,517,6205:.:.:.:.:. 0/1:75,57:132:99:1253,0,2863:.:.:.:.:.
第一行是 header（前面还有其他 header 信息，已在此处删除），各列以制表符分隔。
为了方便理解数据结构，我在此 link 中共享文件的 sub-sample（在可以下载的保管箱上）： https://www.dropbox.com/sh/coihujii38t5prd/AABDXv8ACGIYczeMtzKBo0eea?dl=0
请下载文件，只有300Kb左右，可以用文本编辑器打开。这将有助于更好地理解问题。
问题：
我需要进行上下文相关的替换（值更新）。 - 所有 header 信息（在行前用 # 标记）需要保持不变。
不同行的值对应于最后一个 header（即 CHROM POS ID ....）
首先，我们需要查看 FORMAT 列中的 PG（阶段基因型）字段。字段的不同标签以“:”分隔。 SAMPLE 列中的特定字段有一个相应的值（目前为 2ms01e）。因此，对于第一行，样本 (2ms01e) 的 PG 值为 1|0.
- 现在，我们需要查看同一行FORMAT 列中的GT 字段，并将其对应的值更新为与PG 相同的值。即，将 0/1 更改为 1|0。保持 PG 中的顺序很重要（如果是 1|0 或 0|1，则需要准确）。
但是，如果 PG 字段的值为 0/1、0/0、1/0 或任何其他值（即有一个斜杠），则 GT 字段需要不会更改（或更新）。
最终输出：
因此，第一行数据的 GT 值应更改为：
GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/1:59,89:148:99:3620,0,2177:1|0:.....
到
GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 1|0:59,89:148:99:3620,0,2177:1|0:.....
你可以看到这里只有GT字段的值发生了变化，而其他所有字段的值都保持不变。
对于第二行，GT 值保持不变 - 即 0/0 到 0/0，因为该行的 PG 值为 0/0（与 GT 值相同，因此没有变化）。
简单方法： 我认为最好是 PG 字段的值可以 copy-pasted 到 SAMPLE (2ms01e) 列中的 GT 字段值。 GT字段值在第1位，PG字段在第6位，不同字段之间用“:”分隔。因此，我们需要做的就是用第 6 个字段的值更新第一个字段中的值。
这个简单的方法可能会奏效，因为当 PG 有斜杠“/”时，GT 也会有斜杠，而且顺序无关紧要。 但我不确定它是否适用于每一行。但是，这将是一个简单的解决方案，我至少可以检查并确保它是否有效。
硬方法： 如果 easy 方法没有按预期工作，我认为考虑每个上下文变得很重要。
上下文：
是PG字段值用竖线(|)。如果是，则需要更改。
如果 FORMAT 列中没有 PG 字段 - 则跳过它。
FORMAT 列中字段的分隔符是":"。因此，在SAMPLE 列中。因此，计算字段之间的距离很重要。 GT 和 PG 字段是第 1 和第 6 个位置。
任何类型的解决方案都值得赞赏，但如果它 python 更好，那么如果我的上下文发生变化，我可以阅读和操作它。 此外，对给定解决方案的解释会有很大帮助。
提前致谢，抱歉这么挑剔。我的计算机技能很一般，但仍然不会编程。
干杯！ :))
-
-
下面的答案不会包含任何逻辑的细节，但它会给你一个可能的起点，你可以自己玩：
class VCFProcessor():
def __init__(self):
pass
def load(self, filename):
with open(filename, "rb") as f:
self.load_string(f.read())
def load_string(self, data):
index = 0
for line in data.split("\n"):
# Skip empty rows
if line.strip() == "":
continue
# Assuming there is only header and valid rows
if self.is_valid_row(line):
self.process_row(index, line)
else:
self.process_header(line)
index += 1
def is_valid_row(self, line):
columns = line.split(":")
if len(columns) == 46:
return True
def process_row(self, index, line):
print "Processing line {0} with {1} columns".format(index, len(line.split(":")))
def process_header(self, line):
print "Header has {0} columns".format(len(line.split(":")))
if __name__ == "__main__":
data = """CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 2ms01e 2ms02g 2ms03g 2ms04h
2 15882505 . T A 12134.90 PASS AC=2;AF=0.250;AN=8;BaseQRankSum=-0.021;ClippingRankSum=0.000;DP=695;ExcessHet=3.6798;FS=0.523;MLEAC=2;MLEAF=0.250;MQ=60.00;MQRankSum=0.000;QD=25.18;ReadPosRankSum=1.280;SOR=0.630 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/1:59,89:148:99:3620,0,2177:1|0:.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.:1452:|:0.5 0/1:125,209:334:99:8549,0,4529:.:.:.:.:. 0/0:130,0:130:99:0,400,5809:.:.:.:.:. 0/0:82,0:82:99:0,250,3702:.:.:.:.:.
2 15882583 . G T 1221.33 PASS AC=1;AF=0.125;AN=8;BaseQRankSum=-2.475;ClippingRankSum=0.000;DP=929;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.125;MQ=60.00;MQRankSum=0.000;QD=9.25;ReadPosRankSum=0.299;SOR=0.686 GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC 0/0:178,0:178:99:0,539,7601:0/0:.:.:0/0:. 0/0:446,0:446:99:0,1343,16290:.:.:.:.:. 0/0:172,0:172:99:0,517,6205:.:.:.:.:. 0/1:75,57:132:99:1253,0,2863:.:.:.:.:.
"""
v = VCFProcessor()
v.load_string(data)
希望对您有所帮助。
$ cat > another_mess.awk
[=10=]!="" {
n=split(,a,":") # split at ":" to a array
if(substr(a[6],2,1)=="|") { # if "|" in PG
a[1]=a[6] # copy PG to GT
="" # empty
for(i=1;i<=n;i++) # rebuild
=a[i] (i<n?":":"") # use ":" as delimiter
}
print # PRINT TO TEST, CHANGE TO [=10=]
}
$ awk -f another_mess.awk mess.in
1|0:59,89:148:99:3620,0,2177:1|0:.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.:1452:|:0.5
0/0:178,0:178:99:0,539,7601:0/0:.:.:0/0:.