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:.:.:.:.:.

第一行是标题（它前面有其他标题信息，在此处删除），列是制表符分隔的。

为了方便理解数据结构，我在此链接中共享文件的子样本（在可下载的保管箱上）： https ://www.dropbox.com/sh/coihujii38t5prd/AABDXv8ACGIYczeMtzKBo0eea?dl=0

请下载只有大约 300 Kb 的文件，可以通过文本编辑器打开。这将有助于更好地理解问题。

问题：

我需要做一个依赖于上下文的替换（值更新）。 - 所有标题信息（在行前用# 标记）需要保持不变。

  • 不同行的值对应于最后一个标题（即 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 字段的值复制粘贴到 SAMPLE (2ms01e) 列中的 GT 字段值。GT 字段值为第 1 位，PG 字段值为第 6 位，不同字段用“:”分隔。因此，我们需要做的就是使用第 6 个字段中的值更新第一个字段中的值。

这种简单的方法可能会奏效，因为当 PG 有斜线“/”时，GT 也会有斜线，顺序无关紧要。但我不确定它是否适用于每一行。但是，这将是一个简单的解决方案，我至少可以检查并确保它是否有效。

硬方法： 如果简单的方法不能按预期工作，我认为考虑每个上下文变得很重要。

语境：

是带有管道 (|) 的 PG 字段值。如果是，则需要更改。

如果 FORMAT 列中没有 PG 字段 - 则跳过它。

FORMAT 列中字段的分隔符是“：”。因此，在 SAMPLE 列中。因此，计算场之间的距离很重要。GT 和 PG 字段分别位于第 1 和第 6 位。

任何一种解决方案都值得赞赏，但如果它的 python 更好，那么如果我的上下文发生变化，我可以阅读和操作它。此外，对给定解决方案的解释将有很大帮助。

在此先感谢，很抱歉如此挑剔。我有非常中等的计算机技能，但仍然没有编程。

干杯！:))

    -
4

2 回答 2

1 
$ cat > another_mess.awk  
$0!="" {
    n=split($10,a,":")               # split $10 at ":" to a array
    if(substr(a[6],2,1)=="|") {      # if "|" in PG
        a[1]=a[6]                    # copy PG to GT
        $10=""                       # empty $10
        for(i=1;i<=n;i++)            # rebuild $10
            $10=$10a[i] (i<n?":":"") # use ":" as delimiter
    }
    print $10           # PRINT $10 TO TEST, CHANGE TO $0                       
}

$ 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:.
于 2016-08-21T20:00:02.147 回答
0

下面的答案不包含任何逻辑的任何细节，但它会给你一个可能的起点，所以你可以自己玩：

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)

希望能帮助到你。

于 2016-08-21T14:31:09.503 回答