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

Question

这个问题跟在问题之后:How to read two lines from a file and create dynamics keys in a for-loop?

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

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

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_hybrid表示)。

• 数据是位置线上的字符串。 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，但我需要做一些统计分析才能做到这一点。但是，我可以从视觉上看出 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 中的第 3 (G|C) 和第 4 (G|A) 行。所以，状态是 GgG 和 AgC。同样，我现在数出在 M 与 S 组中存在许多 GcG 与 AgC。

因此，我正在尝试构建一个 Markov 模型，该模型计算 F1 阶段字符串的状态数，并在 group M 与 S组。

我现在解释一下，如何根据F1_hyrbid计算任意XgY(X给定Y)的个数:

• 在计数之前注意条件很重要。现有条件可能是分阶段的(用竖线表示)与非分阶段的(如果两行至少有一个斜杠 (/))。

条件 01:

M1 示例的第 2 行和第 3 行的状态为(T/T with C/C)。因为分隔符是斜杠(/)而不是竖线(|)，所以我们无法判断 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:

但是，如果较早的状态 - 当前状态在两种状态中都由管道分阶段(例如 A|T 在位置 16230007 和 C|G 在位置 16230011 对于样本 Sk1) 我们可以直接做该位置观测状态的相态计数，只有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'] 第一个 M1 在位置 16..9783 和其他相同。

Answer 1

这个问题有点老，但很有趣，因为你有一个非常明确的规范，你需要帮助来编写代码。我将按照自上而下的方法公开一个解决方案，这是一种众所周知的方法，使用普通的旧 python。适应 pandas 应该不难。

自上而下的方法对我来说意味着:如果您不知道如何编写它，只需说出它的名字即可!

你有一个文件(或字符串)作为输入，你想输出一个文件(或字符串)。看起来很简单，但您想合并成对的行来构建每个新行。想法是:

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 和 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))}

process_line 只是在空间上分割线，例如使用 re.split("\s+", line.strip())。

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

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]]

对于 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

奖金:Try it online!

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