如何从文件中读取两行并在 for 循环中创建动态键，后续

这个问题紧接着所讨论的问题：如何从文件中读取两行并在 for 循环中创建动态键？ https://stackoverflow.com/q/41929351/868546

但是，问题的本质已经发展到我想要解决的某种复杂性。

下面是我的数据结构，用空格分隔。

chr pos         M1  M2  Mk  Mg1  F1_hybrid     F1_PG    F1_block    S1  Sk1   S2    Sj
2   16229767    T/T T/T T/T G/T C|T 1|0 726  .  T/C T/C T/C
2   16229783    C/C C/C C/C A/C G|C 0|1 726 G/C G/C G/C C|G
2   16229992    A/A A/A A/A G/A G|A 1|0 726 A/A A/A A/A A|G
2   16230007    T/T T/T T/T A/T A|T 1|0 726 A|T A|T A|T A|T
2   16230011    G/G G/G G/G G/G C|G 1|0 726 G/C C|G C|G G/C
2   16230049    A/A A/A A/A A/A T|A 1|0 726 A|T .   A/T A/T
2   16230174    .   .   .   C/C T|C 1|0 726 C|T T|C T|C C|T
2   16230190    A/A A/A A/A A/A T|A 1|0 726 T|G G|T T|G T|G
2   16230260    A/A A/A A/A A/A G|A 1|0 726 G/G G/G G/G G/G

解释：

• 上述文件中有两大类数据。数据来自Group M样本名称以M，并且类似地group S有几个以以下开头的列名称S.

• 还有一个混合列（表示为F1_混合动力).

• 数据是沿着位置线的字符串。 F1_hybrid 通过管道 (|) 区分两个字母。因此，F1 中的两个字符串值是C-G-G-A-C-T-T-T-G，而另一个字符串值是 T-C-A-T-G-A-C-A-A。该字符串之一来自M-group而另一个来自S-group but 为此，我需要进行一些统计分析。然而，我可以从视觉上看出 T-C-A-T-G-A-C-A-A 字符串很可能来自 M 组。

程序：

• 我阅读第一行并使用列信息创建唯一键。

• 然后我读取第二行和第三行以及 F1_hybrid 中的值，即 C|T 和 G|C。现在，我需要计算 M 组与 S 组之间存在多少个 GgC（解释为 G 给定 C）与 CgT（C 给定 T）。

• 然后读取 F1_hybrid 中的第三行 (G|C) 和第四行 (G|A)。因此，状态是 GgG 和 AgC。相似地，我现在统计M vs. S组中存在很多GcG vs. AgC。

因此，我正在尝试建立一个Markov-model它计算 F1 中的相控串的状态数，并将观察到的计数取入group M vs group S.

我现在解释一下，如何根据F1_hyrbid计算任意XgY（X给定Y）的数量：

• 在进行计数之前，请务必注意条件。现有条件可以是分阶段的（通过有管道表示）与不分阶段的（如果两条线至少有一个斜线（/）。

条件01：

The M1样品有状态为（T/T 和 C/C）对于第二行和第三行。因为分隔符是斜杠 (/)并不是pipe (|)我们无法判断 M1-sample 处于哪个确切状态。但是，我们可以创建组合矩阵（对于先前状态与当前状态）

    T     T
C  CgT   CgT
C  CgT   CgT

现在，我们可以说有4 总 CgT

如果满足这个条件，我们就继续做同样的矩阵。

条件02

M 组的其他样品也是如此，但 Mg1 除外，其中 G/T 位于 A/C 之前。所以，矩阵是：

    G     T
A  AgG   AgT
C  CgG   CgT

因此，我们在这里观察到 1 个 CgT 计数。

条件03：

但是，如果早期状态 - 当前状态在两种状态下都通过管道进行分阶段（对于样本 Sk1，A|T 位于位置 16230007，C|G 位于位置 16230011）我们可以直接计数该位置处观测态的相态，只有CgA和GgT，所以CgT的计数为0。

条件04：如果其中一个状态具有竖线 (|)，而另一个状态具有斜杠 (/)，则条件将与两个具有斜杠的状态相同。

条件05：如果 previous_state 或 present_state 中的任何一个是 period(.)，则对于 F1_hybrid 预期的状态，观察计数自动为零 (0)。

因此，预期的输出应该是这样的：

pos     M1  M2  Mk  Mg1 H0  H1  S1  Sk1 S2  Sj
16..9783    4-CgT   4-CgT   4-CgT   1-CgT   GgC CgT 0   1-CgT   1-CgT   1-CgT
16..9992    4-AgC   4-AgC   4-AgC   2-AgC   GgG AgC 1-AgC   1-AgC   1-AgC   1-AgC,1-GgG
16..0007    4-TgA   4-TgA   4-TgA   1-AgG,1-TgA AgG TgA 2-TgA   2-TgA   2-TgA1  1-TgA
..................contd

或者，每列的字典格式的值同样有效。就像是['4-CgT','4-CgT','4-CgT','4-CgT']对于位置 16..9783 处的第一个 M1，其他位置相同。

这个问题有点老了，但是很有趣，因为您有一个非常清晰的规范，并且需要帮助来编写代码。我将采用自上而下的方法公开一个解决方案，这是一种众所周知的方法，使用普通的旧 python。适应熊猫应该不难。

自上而下的方法对我来说意味着：如果你不知道怎么写，就直接命名吧！

您有一个文件（或一个字符串）作为输入，并且您想要输出一个文件（或一个字符串）。这看起来很简单，但您想要合并行对来构建每个新行。这个想法是：

1. 获取输入的行，作为字典
2. 把他们分成两个
3. 为每对建立一个新行
4. 输出结果

您现在不知道如何编写行生成器。您也不知道如何为每对构建新行。不要被困难所阻碍，只提出解决方案。假设你有一个函数get_rows和一个函数build_new_row。我们这样写：

def build_new_rows(f):
    """generate the new rows. Output may be redirected to a file"""
    rows = get_rows(f) # get a generator on rows = dictionaries.
    r1 = next(rows) # store the first row
    for r2 in rows: # for every following row
        yield build_new_row(r1, r2) # yield a new row built of the previous stored row and the current row.
        r1 = r2 # store the current row, which becomes the previous row

现在，检查两个“缺失”的函数：get_rows and build_new_row。 功能get_rows很容易写。这是主要部分：

header = process_line(next(f))
for line in f:
    yield {k:v for k,v in zip(header, process_line(line))}

where process_line只是在空间上分割线，例如与一个re.split("\s+", line.strip()).

第二部分是build_new_row。仍然是自上而下的方法：您需要从预期的表中构建 H0 和 H1，然后根据您暴露的条件为每个 M 和 S 构建 H1 的计数。假装你有一个pipe_compute计算 H0 和 H1 的函数，以及build_count为每个 M 和 S 构建 H1 计数的函数：

def build_new_row(r1, r2):
    """build a row"""
    h0, h1 = pipe_compute(r1["F1_hybrid"], r2["F1_hybrid"])

    # initialize the dict whith the pos, H0 and H1
    new_row = {"pos":r2["pos"], "H0":h0, "H1":h1}

    for key in r1.keys():
        if key[0] in ("M", "S"):
            new_row[key] = build_count(r1[key], r2[key], h1)

    return new_row

你现在几乎拥有了一切。看一眼pipe_compute: 这正是你在条件03中所写的。

def pipe_compute(v1, v2):
    """build H0 H1 according to condition 03"""
    xs = v1.split("|")
    ys = v2.split("|")
    return [ys[0]+"g"+xs[0], ys[1]+"g"+xs[1]]

And for buid_count，坚持自上而下的方法：

def build_count(v1, v2, to_count):
    """nothing funny here: just follow the conditions"""
    if is_slash_count(v1, v2): # are conditions 01, 02, 04 true ?
        c = slash_count(v1, v2)[to_count] # count how many "to_count" we find in the 2 x 2 table of condtions 01 or 02.
    elif "|" in v1 and "|" in v2: # condition 03
        c = pipe_count(v1, v2)[to_count]
    elif "." in v1 or "." in v2: # condition 05
        return '0'
    else:
        raise Exception(v1, v2)

    return "{}-{}".format(c, to_count) # n-XgY

我们还在往下走。我们什么时候有is_slash_count？两条斜线（条件 01 和 02）或一条斜线和一根竖线（条件 04）：

def is_slash_count(v1, v2):
    """conditions 01, 02, 04"""
    return "/" in v1 and "/" in v2 or "/" in v1 and "|" in v2 or "|" in v1 and "/" in v2

功能slash_count就是条件 01 和 02 的 2 x 2 表：

def slash_count(v1, v2):
    """count according to conditions 01, 02, 04"""
    cnt = collections.Counter()
    for x in re.split("[|/]", v1): # cartesian product
        for y in re.split("[|/]", v2): # cartesian product
            cnt[y+"g"+x] += 1
    return cnt # a dictionary XgY -> count(XgY)

功能pipe_count甚至更简单，因为你只需要计算结果pipe_compute:

def pipe_count(v1, v2):
    """count according to condition 03"""
    return collections.Counter(pipe_compute(v1, v2))

现在你已经完成了（并且结束了）。我得到这个结果，与您的期望略有不同，但您肯定已经看到了我的错误：

pos M1  M2  Mk  Mg1 H0  H1  S1  Sk1 S2  Sj
16229783    4-CgT   4-CgT   4-CgT   1-CgT   GgC CgT 0   1-CgT   1-CgT   1-CgT
16229992    4-AgC   4-AgC   4-AgC   1-AgC   GgG AgC 2-AgC   2-AgC   2-AgC   1-AgC
16230007    4-TgA   4-TgA   4-TgA   1-TgA   AgG TgA 2-TgA   2-TgA   2-TgA   0-TgA
16230011    4-GgT   4-GgT   4-GgT   2-GgT   CgA GgT 1-GgT   1-GgT   1-GgT   1-GgT
16230049    4-AgG   4-AgG   4-AgG   4-AgG   TgC AgG 1-AgG   0   1-AgG   1-AgG
16230174    0   0   0   4-CgA   TgT CgA 1-CgA   0   1-CgA   1-CgA
16230190    0   0   0   4-AgC   TgT AgC 0-AgC   0-AgC   0-AgC   0-AgC
16230260    4-AgA   4-AgA   4-AgA   4-AgA   GgT AgA 0-AgA   0-AgA   0-AgA   0-AgA

Bonus: 在线尝试一下！ https://tio.run/##lVZtb@M2DP6eX0G4H5qgsWO7XXvNUAxFsKU3oLgB7YcNucBwbDlW41iZrLxd09/ekZLtOD3fhgWgI4lvj0iK0mqvUpFfvr@fQSRins@HsFaJ/alzBqlSq2I4GBQqjBZiw2SSia0TieXg7zUrFBd5MbjyXc@/uXQHqdjaStiShbGttsLOeM4KO5FiaYd2wjNmh3lsR8hXzI73ebjkUWEv2L6weU4iQtqZECt7XSAIe9Xp8OVKSAVcVCPJqlEksoxFGkG1VOzrYVRsOh21U3AHlmVFqYSVKKD6PXpIPtICaY6T37wg3c8kjzUbZ3@My8EsE9GCxk8o9rTwaOTr@UuH/r1r37@9ub6hpefBc01jpNHhGbyDCzf@NYBD/FFFDd1Pl6Q7wvWK7pHGhxG4B0/rjmle0ugwPure3moo94P7msZEh/vab5N3f9S9dF33O8z3RA3MNG5QQ9ejOKCv8Qkhtlq3xKppfNwv6l7dfsRM9NzEjP4cLaMxHXW9myvSdRpE8XrGWFW6FHOaE42amL1b97/8PhNWrT8mOur61626zTh/iAUWXacTswRma57FQc62gRTbopv0hh0yhfw5y5nEkwAqZYACQAIOfFmr1VrBMtzDjGG9x1ximbMYlIAQ6BSRbbJB8ljfc6Yq22bZw8Wc7VSXVs0aniyQPvBcKxkI9NtzlsWnGLvS66Nsr5bR9qRvttNwVm9kXG1EA6LjrmFCxhcMxOwF4TsV5jMoVhlXIHII8z1Owog5mrOSImJFEVDTQIdZuJzFIdBsiFFwtFrX@lpcWH296hRK8lW31@to9RSbDpOo2LTT1WFIescgaOsYhuRjDF4Xw42WWPQ3JPANbRub/VOT9On13lqyW0WuDozmYtKQNwTPARViQAqxZLAJM2yfJlbS/xl8B/vZEvNuqsGwq5ilbh9SSsKKr1hQyqGziVW3LWtKnk8Wyric4Wa44mHGvxnbMY8UbFOuUj3FttiHBxezEcODp1XK7aDDVwvZ1pAs0wCdWA@uNSRA1oOHA@@tU4cW@7guMM@hjt7tHSPME2JO3Cnxu9YjZtB6shoCDa8TlJyiaxPZSKxzRVulVb1FM0i9cnuSqbXMK2WTlJMwbTAnm@9ykvFCQRhFQtJ1R0crEnnM6TYB97IK/I6O18arau9gmTra62X/43IJZbIvcKcX1ty62NGoj/ITr17wptNm6ZgNGpB9xGEWjnBzoVJCmKxzPC4pk3gaXtYIPsHrD3NEKayh1xWDAedFUGRhkZ546A2xHkLZVAEXOS76dq9AyTWDX@B3so@XMx4U5dQ5inDTLRYnFeYp0OsBB4DPAGxfCNeqeBZsGfYELDEsAILsww5JhTNsEyLRcCo0gKXk4nn4HgXLcF8YbjKy8XTFVjN/qL0fU3gCuyyIVtQNy05lGRFUk4@Gfzq2RZPuc/e8NFGwY0HLkBcMft1FbEV6lduTmrVe3@zXN8vBs7MMVTdqpB995vaf879MpfwglVWJtOayviNKX4OTqJUzX@900BpQzTq0ahlU/w6J6qD9fLWhjHJ6qDVedM6ILDDZPXbunW4u9UUwOQym2Ec2XqONkNi@Xcz/0G3Q42RvzuQULu7Aa0YLmSed5P9usdFCKostW2trUz3jdxnyvGqgCUaGC@cJb7t8/vlLFx@1VB/66qAHbv0saN63p1ei6d7YdB89/fX1d6G/c72EXd20dOrM5rswf1r26cUyB4Wu/9MLOh7CJJ4spuYKoOAbt1NT61vJlcaAD3LHTLr4TMfLOxZr1YcYj96SVu/Ov6rzXkOnlMZr1RhsPGTIyw/eVW36FeiupPjiSQ8C6itBAHd3cB4EFO0gODcGTOg77@//AA

重要的是，除了这个具体问题的解决方案之外，我使用的方法以及在软件开发中广泛使用的方法。代码可能会改进很多。