eve*_*007 9 python dictionary numpy pandas defaultdict
这个问题遵循问题:如何从文件中读取两行并在for循环中创建动态键？
但是,问题的性质已演变为我想要解决的某种复杂性.
下面是按空格分隔的数据结构.
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,而另一个字符串值是TCATGACAA.其中一个字符串来自M-group而另一个来自S-group,但我需要做一些统计分析.但是,我可以直观地说TCATGACAA字符串很可能来自M-group.
程序:
我读了第一行并使用列信息创建了一个唯一键.
然后我读了第二行和第三行以及F1_hybrid中的值,即C | T和G | C.现在,我需要计算在M组与S组之间存在多少GgC(解释为G给定C)与CgT(C给定T).
然后在F1_hybrid中用第4(G | A)行读取第3个(G | C).所以,州是GgG和AgC.同样,我现在认为M对S组中存在许多GcG对抗AgC.
因此,我正在尝试构建一个Markov-model从F1中计算分阶段字符串的状态数并在group Mvs中计算观察计数的数量group S.
我现在解释,如何基于F1_hyrbid计算任何XgY(X给定Y)的数量:
条件01:
该M1样品具有状态为(T/T与C/C)为第二和第三线.由于分隔符是斜杠(/)而不是管道(|),我们无法分辨M1样本所处的确切状态.但是,我们可以创建组合矩阵(对于具有当前状态的先前状态)
T T
C CgT CgT
C CgT CgT
现在,我们可以看出共有4个CgT
如果这个条件满足,我们会继续做同样的矩阵.
条件02
对于来自M组的其他样品的情况也是如此,除了Mg /,其中G/T在A/C之前.所以,矩阵是:
G T
A AgG AgT
C CgG CgT
所以,在这里我们观察到1个CgT计数.
条件03:
但是,如果早期的状态 - 存在状态在两个状态下都被管道定相(如位置16230007处的A | T,样本Sk1处的位置16230011处的C | G),我们可以直接计算该位置处观察状态的相位状态,只有CgA和GgT,因此CgT的计数为0.
条件04: 如果其中一个状态具有管道(|)但是其他状态具有斜杠(/),则该条件将与具有斜杠的两个状态相同.
条件05: 如果previous_state或present_state中的任何一个是句点(.),则对于从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。适应熊猫应该不难。
\n\n自上而下的方法对我来说意味着：如果你不知道如何写它，就直接说出它的名字！
\n\n您有一个文件（或一个字符串）作为输入，并且您想要输出一个文件（或一个字符串）。这看起来很简单，但您想要合并行对来构建每个新行。这个想法是：
\n\n您现在不知道如何编写行生成器。您也不知道如何为每一对构建新行。不要被困难所阻碍，只需说出解决方案即可。想象一下你有一个函数get_rows和一个函数build_new_row。让我们这样写：
def build_new_rows(f):\n """generate the new rows. Output may be redirected to a file"""\n rows = get_rows(f) # get a generator on rows = dictionaries.\n r1 = next(rows) # store the first row\n for r2 in rows: # for every following row\n yield build_new_row(r1, r2) # yield a new row built of the previous stored row and the current row.\n r1 = r2 # store the current row, which becomes the previous row\n现在，检查两个“缺失”的函数：get_rows和build_new_row。\n该函数get_rows非常容易编写。这是主要部分：
header = process_line(next(f))\nfor line in f:\n yield {k:v for k,v in zip(header, process_line(line))}\n其中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):\n """build a row"""\n h0, h1 = pipe_compute(r1["F1_hybrid"], r2["F1_hybrid"])\n\n # initialize the dict whith the pos, H0 and H1\n new_row = {"pos":r2["pos"], "H0":h0, "H1":h1}\n\n for key in r1.keys():\n if key[0] in ("M", "S"):\n new_row[key] = build_count(r1[key], r2[key], h1)\n\n return new_row\n你现在几乎拥有了一切。看一下pipe_compute：这正是你在条件 03 中所写的。
def pipe_compute(v1, v2):\n """build H0 H1 according to condition 03"""\n xs = v1.split("|")\n ys = v2.split("|")\n return [ys[0]+"g"+xs[0], ys[1]+"g"+xs[1]]\n对于buid_count，坚持自上而下的方法：
def build_count(v1, v2, to_count):\n """nothing funny here: just follow the conditions"""\n if is_slash_count(v1, v2): # are conditions 01, 02, 04 true ?\n 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.\n elif "|" in v1 and "|" in v2: # condition 03\n c = pipe_count(v1, v2)[to_count]\n elif "." in v1 or "." in v2: # condition 05\n return \'0\'\n else:\n raise Exception(v1, v2)\n\n return "{}-{}".format(c, to_count) # n-XgY\n我们还在往下走。我们什么时候有is_slash_count？两条斜线（条件 01 和 02）或一条斜线和一根竖线（条件 04）：
def is_slash_count(v1, v2):\n """conditions 01, 02, 04"""\n return "/" in v1 and "/" in v2 or "/" in v1 and "|" in v2 or "|" in v1 and "/" in v2\n该函数slash_count只是条件 01 和 02 的 2 x 2 表：
def slash_count(v1, v2):\n """count according to conditions 01, 02, 04"""\n cnt = collections.Counter()\n for x in re.split("[|/]", v1): # cartesian product\n for y in re.split("[|/]", v2): # cartesian product\n cnt[y+"g"+x] += 1\n return cnt # a dictionary XgY -> count(XgY)\n该函数pipe_count甚至更简单，因为您只需计算 的结果pipe_compute：
def pipe_count(v1, v2):\n """count according to condition 03"""\n return collections.Counter(pipe_compute(v1, v2))\n现在你已经完成了（并且结束了）。我得到这个结果，与您的期望略有不同，但您肯定已经看到了我的错误：
\n\npos M1 M2 Mk Mg1 H0 H1 S1 Sk1 S2 Sj\n16229783 4-CgT 4-CgT 4-CgT 1-CgT GgC CgT 0 1-CgT 1-CgT 1-CgT\n16229992 4-AgC 4-AgC 4-AgC 1-AgC GgG AgC 2-AgC 2-AgC 2-AgC 1-AgC\n16230007 4-TgA 4-TgA 4-TgA 1-TgA AgG TgA 2-TgA 2-TgA 2-TgA 0-TgA\n16230011 4-GgT 4-GgT 4-GgT 2-GgT CgA GgT 1-GgT 1-GgT 1-GgT 1-GgT\n16230049 4-AgG 4-AgG 4-AgG 4-AgG TgC AgG 1-AgG 0 1-AgG 1-AgG\n16230174 0 0 0 4-CgA TgT CgA 1-CgA 0 1-CgA 1-CgA\n16230190 0 0 0 4-AgC TgT AgC 0-AgC 0-AgC 0-AgC 0-AgC\n16230260 4-AgA 4-AgA 4-AgA 4-AgA GgT AgA 0-AgA 0-AgA 0-AgA 0-AgA\n奖励：在线尝试！
\n\n重要的是，除了这个具体问题的解决方案之外，我使用的方法以及在软件开发中广泛使用的方法。代码可能会改进很多。
\n