从模拟网格的历史中找到细胞的聚结时间
Find coalescent times of cells from the history of a simulated grid
我正在编写一个由 n x n 单元格组成的模拟。在模拟的不同时间，随机抽取一个单元格到'divide'。当一个细胞分裂时，它会死亡并产生两个子细胞。一个女儿替换原来的单元格，另一个女儿随机替换网格上 8 个邻居中的一个。
网格由数据帧编码，开头有 n^2 行，每个单元格一行（每个单元格有 birth_time=0，death_time=50 和 parent=0 开头）。随着模拟的进行，为每个分裂事件添加代表子细胞的两行，并更新 parent（和相邻前体）的死亡时间。女儿们被分配 birth_time!=0、death_time=50 和 parent（示例见下文）。
在模拟 运行 一段时间后（在下面的示例中为 50），我对所有具有相同 x-coordinate 的单元格进行了采样。对于这些细胞，我想使用 grid-dataframe 中编码的历史信息来查找它们的聚结时间，即最终样本中两个或多个细胞的祖先的所有细胞的死亡时间。我正在寻找一个函数来在 R 中完成此操作（或帮助构建一个我可以自己在 R 中编写代码的算法）。
下面是三个示例，我希望它们能使我的要求更加明确：
测试 1：
> grid1
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 1 1 1 1 0 0 50
2 2 2 1 1 0 0 50
3 3 3 1 1 0 0 2
4 4 4 1 1 0 0 50
5 5 5 1 1 0 0 50
6 6 1 2 1 0 0 50
7 7 2 2 0 0 0 50
8 8 3 2 0 0 0 2
9 9 4 2 0 0 0 50
10 10 5 2 1 0 0 50
11 11 1 3 1 0 0 50
12 12 2 3 0 0 0 50
13 13 3 3 0 0 0 12
14 14 4 3 0 0 0 50
15 15 5 3 1 0 0 50
16 16 1 4 1 0 0 50
17 17 2 4 0 0 0 50
18 18 3 4 0 0 0 21
19 19 4 4 0 0 0 50
20 20 5 4 1 0 0 50
21 21 1 5 1 0 0 50
22 22 2 5 1 0 0 50
23 23 3 5 1 0 0 50
24 24 4 5 1 0 0 50
25 25 5 5 1 0 0 50
26 26 3 2 0 8 2 12
27 27 3 1 1 8 2 50
28 28 3 2 0 26 12 33
29 29 3 3 0 26 12 21
30 30 3 3 0 29 21 33
31 31 3 4 0 29 21 45
32 32 3 3 0 30 33 45
33 33 3 2 0 30 33 50
34 34 3 4 0 31 45 50
35 35 3 3 0 31 45 50
我对结束时间 (50) 存在的 crypts 进行了采样，并且 x-coordinate=3。请注意，虽然我在此测试用例中对所有 5 个 crypt 进行了采样，但在实际模拟中将对一个子集进行采样。
> sample1
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
23 23 3 5 1 0 0 50
27 27 3 1 1 8 2 50
33 33 3 2 0 30 33 50
34 34 3 4 0 31 45 50
35 35 3 3 0 31 45 50
在此示例中，位于 (3,5) 的单元格与其他单元格无关（除了所有单元格 (0) 的 pseudo-parent 节点。其他四个单元格均相关，并且有 3为系统发育提供信息的分裂事件。它们是：
> res1
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 8 3 2 0 0 0 2
3 29 3 3 0 26 12 21
5 31 3 4 0 29 21 45
下面的树显示了我试图捕获的关系
这里还有另外两个例子：
测试 2：
> grid2
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 1 1 1 1 0 0 50
2 2 2 1 1 0 0 2
3 3 3 1 1 0 0 50
4 4 4 1 1 0 0 45
5 5 5 1 1 0 0 50
6 6 1 2 1 0 0 50
7 7 2 2 0 0 0 2
8 8 3 2 0 0 0 45
9 9 4 2 0 0 0 21
10 10 5 2 1 0 0 21
11 11 1 3 1 0 0 50
12 12 2 3 0 0 0 50
13 13 3 3 0 0 0 33
14 14 4 3 0 0 0 50
15 15 5 3 1 0 0 50
16 16 1 4 1 0 0 50
17 17 2 4 0 0 0 33
18 18 3 4 0 0 0 12
19 19 4 4 0 0 0 50
20 20 5 4 1 0 0 50
21 21 1 5 1 0 0 50
22 22 2 5 1 0 0 50
23 23 3 5 1 0 0 50
24 24 4 5 1 0 0 12
25 25 5 5 1 0 0 50
26 26 2 2 0 7 2 50
27 27 2 1 1 7 2 50
28 28 3 4 0 18 12 50
29 29 4 5 1 18 12 50
30 30 4 2 0 9 21 50
31 31 5 2 1 9 21 50
32 32 2 4 0 17 33 50
33 33 3 3 0 17 33 50
34 34 3 2 0 8 45 50
35 35 4 1 1 8 45 50
> sample2
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
3 3 3 1 1 0 0 50
23 23 3 5 1 0 0 50
28 28 3 4 0 18 12 50
33 33 3 3 0 17 33 50
34 34 3 2 0 8 45 50
sample2中的细胞完全是un-related（它们的most-recent共同祖先是0pseudo-node）。该函数应该 return 什么都没有（或者只是时间 0）。
测试 3：
> grid3
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 1 1 1 1 0 0 50
2 2 2 1 1 0 0 50
3 3 3 1 1 0 0 50
4 4 4 1 1 0 0 50
5 5 5 1 1 0 0 50
6 6 1 2 1 0 0 50
7 7 2 2 0 0 0 31
8 8 3 2 0 0 0 34
9 9 4 2 0 0 0 37
10 10 5 2 1 0 0 50
11 11 1 3 1 0 0 50
12 12 2 3 0 0 0 22
13 13 3 3 0 0 0 8
14 14 4 3 0 0 0 8
15 15 5 3 1 0 0 6
16 16 1 4 1 0 0 50
17 17 2 4 0 0 0 2
18 18 3 4 0 0 0 2
19 19 4 4 0 0 0 3
20 20 5 4 1 0 0 50
21 21 1 5 1 0 0 50
22 22 2 5 1 0 0 50
23 23 3 5 1 0 0 50
24 24 4 5 1 0 0 50
25 25 5 5 1 0 0 50
26 26 2 4 0 17 2 50
27 27 3 4 0 17 2 3
28 28 3 4 0 27 3 45
29 29 4 4 0 27 3 6
30 30 4 4 0 29 6 50
31 31 5 3 1 29 6 50
32 32 4 3 0 14 8 50
33 33 3 3 0 14 8 22
34 34 3 3 0 33 22 45
35 35 2 3 0 33 22 31
36 36 2 3 0 35 31 50
37 37 2 2 0 35 31 34
38 38 2 2 0 37 34 50
39 39 3 2 0 37 34 37
40 40 3 2 0 39 37 49
41 41 4 2 0 39 37 50
42 42 3 3 0 34 45 49
43 43 3 4 0 34 45 50
44 44 3 3 0 42 49 50
45 45 3 2 0 42 49 50
> sample3 <- subset(grid3, x_coordinate==3 & death_time==50)
> sample3
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
3 3 3 1 1 0 0 50
23 23 3 5 1 0 0 50
43 43 3 4 0 34 45 50
44 44 3 3 0 42 49 50
45 45 3 2 0 42 49 50
此网格有许多事件与 x-coordinate 3 重叠，但只有两个具有信息性：
> res3
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 42 3 3 0 34 45 49
2 34 3 3 0 33 22 45
如果有人觉得它有帮助，这是我 semi-crude 每个时间点每个网格状态的绘图（忽略前两行）：
非常感谢您的帮助！
你的问题很难理解，我没有完全理解你需要什么以及为什么选择每个数据行作为结果。我的职能是检查在当地社区中幸存下来的每一代人的祖先以及 returns 他们的信息。也许这会为您解决问题提供指导。
find.elders = function(x, dead, dat){
locals = dat[dat$x_coordinate == x & dat$death_time != dead,]
survivors = dat[dat$x_coordinate == x & dat$death_time == dead,]
anc = survivors$parent
res = NULL
while(any(anc != 0)){
anc = anc[anc > 0]
cat("Ancestors:", anc, "\n")
res = c(res, which(locals$parent %in% anc))
survivors = locals[locals$cellID %in% anc,]
anc = survivors$parent
}
#res = c(res, which(locals$parent %in% anc))
locals[res,]
}
find.elders(3, 50, grid1)
对于未来的读者，我意识到这个问题相当复杂且难以理解。我为此道歉。
我最终通过在表格 0/1、0/1/27 等上向网格数据帧添加一个属性 'pathString' 来解决我的问题，其中对于每个单元格，该单元格也存储它的所有祖先作为 'itself'。
然后我可以使用 R 中 data.tree 包中的 as.Node() 函数将我的网格转换为树对象，随后可以使用 as.phylo() 函数将其转换为 phylo 对象在猿。
一旦采样的单元格存储为树，ape 和 ggtree 中的现有函数使剩下的事情变得容易。
请参阅 data.tree 包的 data(acme) 和此处的 #Tree 示例：https://rdrr.io/cran/data.tree/man/as.Node.data.frame.html
我正在编写一个由 n x n 单元格组成的模拟。在模拟的不同时间，随机抽取一个单元格到'divide'。当一个细胞分裂时，它会死亡并产生两个子细胞。一个女儿替换原来的单元格，另一个女儿随机替换网格上 8 个邻居中的一个。
网格由数据帧编码，开头有 n^2 行，每个单元格一行（每个单元格有 birth_time=0，death_time=50 和 parent=0 开头）。随着模拟的进行，为每个分裂事件添加代表子细胞的两行，并更新 parent（和相邻前体）的死亡时间。女儿们被分配 birth_time!=0、death_time=50 和 parent（示例见下文）。
在模拟 运行 一段时间后（在下面的示例中为 50），我对所有具有相同 x-coordinate 的单元格进行了采样。对于这些细胞，我想使用 grid-dataframe 中编码的历史信息来查找它们的聚结时间，即最终样本中两个或多个细胞的祖先的所有细胞的死亡时间。我正在寻找一个函数来在 R 中完成此操作（或帮助构建一个我可以自己在 R 中编写代码的算法）。
下面是三个示例，我希望它们能使我的要求更加明确：
测试 1：
> grid1
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 1 1 1 1 0 0 50
2 2 2 1 1 0 0 50
3 3 3 1 1 0 0 2
4 4 4 1 1 0 0 50
5 5 5 1 1 0 0 50
6 6 1 2 1 0 0 50
7 7 2 2 0 0 0 50
8 8 3 2 0 0 0 2
9 9 4 2 0 0 0 50
10 10 5 2 1 0 0 50
11 11 1 3 1 0 0 50
12 12 2 3 0 0 0 50
13 13 3 3 0 0 0 12
14 14 4 3 0 0 0 50
15 15 5 3 1 0 0 50
16 16 1 4 1 0 0 50
17 17 2 4 0 0 0 50
18 18 3 4 0 0 0 21
19 19 4 4 0 0 0 50
20 20 5 4 1 0 0 50
21 21 1 5 1 0 0 50
22 22 2 5 1 0 0 50
23 23 3 5 1 0 0 50
24 24 4 5 1 0 0 50
25 25 5 5 1 0 0 50
26 26 3 2 0 8 2 12
27 27 3 1 1 8 2 50
28 28 3 2 0 26 12 33
29 29 3 3 0 26 12 21
30 30 3 3 0 29 21 33
31 31 3 4 0 29 21 45
32 32 3 3 0 30 33 45
33 33 3 2 0 30 33 50
34 34 3 4 0 31 45 50
35 35 3 3 0 31 45 50
我对结束时间 (50) 存在的 crypts 进行了采样，并且 x-coordinate=3。请注意，虽然我在此测试用例中对所有 5 个 crypt 进行了采样，但在实际模拟中将对一个子集进行采样。
> sample1
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
23 23 3 5 1 0 0 50
27 27 3 1 1 8 2 50
33 33 3 2 0 30 33 50
34 34 3 4 0 31 45 50
35 35 3 3 0 31 45 50
在此示例中，位于 (3,5) 的单元格与其他单元格无关（除了所有单元格 (0) 的 pseudo-parent 节点。其他四个单元格均相关，并且有 3为系统发育提供信息的分裂事件。它们是：
> res1
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 8 3 2 0 0 0 2
3 29 3 3 0 26 12 21
5 31 3 4 0 29 21 45
下面的树显示了我试图捕获的关系
这里还有另外两个例子： 测试 2：
> grid2
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 1 1 1 1 0 0 50
2 2 2 1 1 0 0 2
3 3 3 1 1 0 0 50
4 4 4 1 1 0 0 45
5 5 5 1 1 0 0 50
6 6 1 2 1 0 0 50
7 7 2 2 0 0 0 2
8 8 3 2 0 0 0 45
9 9 4 2 0 0 0 21
10 10 5 2 1 0 0 21
11 11 1 3 1 0 0 50
12 12 2 3 0 0 0 50
13 13 3 3 0 0 0 33
14 14 4 3 0 0 0 50
15 15 5 3 1 0 0 50
16 16 1 4 1 0 0 50
17 17 2 4 0 0 0 33
18 18 3 4 0 0 0 12
19 19 4 4 0 0 0 50
20 20 5 4 1 0 0 50
21 21 1 5 1 0 0 50
22 22 2 5 1 0 0 50
23 23 3 5 1 0 0 50
24 24 4 5 1 0 0 12
25 25 5 5 1 0 0 50
26 26 2 2 0 7 2 50
27 27 2 1 1 7 2 50
28 28 3 4 0 18 12 50
29 29 4 5 1 18 12 50
30 30 4 2 0 9 21 50
31 31 5 2 1 9 21 50
32 32 2 4 0 17 33 50
33 33 3 3 0 17 33 50
34 34 3 2 0 8 45 50
35 35 4 1 1 8 45 50
> sample2
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
3 3 3 1 1 0 0 50
23 23 3 5 1 0 0 50
28 28 3 4 0 18 12 50
33 33 3 3 0 17 33 50
34 34 3 2 0 8 45 50
sample2中的细胞完全是un-related（它们的most-recent共同祖先是0pseudo-node）。该函数应该 return 什么都没有（或者只是时间 0）。
测试 3：
> grid3
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 1 1 1 1 0 0 50
2 2 2 1 1 0 0 50
3 3 3 1 1 0 0 50
4 4 4 1 1 0 0 50
5 5 5 1 1 0 0 50
6 6 1 2 1 0 0 50
7 7 2 2 0 0 0 31
8 8 3 2 0 0 0 34
9 9 4 2 0 0 0 37
10 10 5 2 1 0 0 50
11 11 1 3 1 0 0 50
12 12 2 3 0 0 0 22
13 13 3 3 0 0 0 8
14 14 4 3 0 0 0 8
15 15 5 3 1 0 0 6
16 16 1 4 1 0 0 50
17 17 2 4 0 0 0 2
18 18 3 4 0 0 0 2
19 19 4 4 0 0 0 3
20 20 5 4 1 0 0 50
21 21 1 5 1 0 0 50
22 22 2 5 1 0 0 50
23 23 3 5 1 0 0 50
24 24 4 5 1 0 0 50
25 25 5 5 1 0 0 50
26 26 2 4 0 17 2 50
27 27 3 4 0 17 2 3
28 28 3 4 0 27 3 45
29 29 4 4 0 27 3 6
30 30 4 4 0 29 6 50
31 31 5 3 1 29 6 50
32 32 4 3 0 14 8 50
33 33 3 3 0 14 8 22
34 34 3 3 0 33 22 45
35 35 2 3 0 33 22 31
36 36 2 3 0 35 31 50
37 37 2 2 0 35 31 34
38 38 2 2 0 37 34 50
39 39 3 2 0 37 34 37
40 40 3 2 0 39 37 49
41 41 4 2 0 39 37 50
42 42 3 3 0 34 45 49
43 43 3 4 0 34 45 50
44 44 3 3 0 42 49 50
45 45 3 2 0 42 49 50
> sample3 <- subset(grid3, x_coordinate==3 & death_time==50)
> sample3
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
3 3 3 1 1 0 0 50
23 23 3 5 1 0 0 50
43 43 3 4 0 34 45 50
44 44 3 3 0 42 49 50
45 45 3 2 0 42 49 50
此网格有许多事件与 x-coordinate 3 重叠，但只有两个具有信息性：
> res3
cellID x_coordinate y_coordinate onEdge parent birth_time death_time
1 42 3 3 0 34 45 49
2 34 3 3 0 33 22 45
如果有人觉得它有帮助，这是我 semi-crude 每个时间点每个网格状态的绘图（忽略前两行）：
非常感谢您的帮助！
你的问题很难理解，我没有完全理解你需要什么以及为什么选择每个数据行作为结果。我的职能是检查在当地社区中幸存下来的每一代人的祖先以及 returns 他们的信息。也许这会为您解决问题提供指导。
find.elders = function(x, dead, dat){
locals = dat[dat$x_coordinate == x & dat$death_time != dead,]
survivors = dat[dat$x_coordinate == x & dat$death_time == dead,]
anc = survivors$parent
res = NULL
while(any(anc != 0)){
anc = anc[anc > 0]
cat("Ancestors:", anc, "\n")
res = c(res, which(locals$parent %in% anc))
survivors = locals[locals$cellID %in% anc,]
anc = survivors$parent
}
#res = c(res, which(locals$parent %in% anc))
locals[res,]
}
find.elders(3, 50, grid1)
对于未来的读者，我意识到这个问题相当复杂且难以理解。我为此道歉。
我最终通过在表格 0/1、0/1/27 等上向网格数据帧添加一个属性 'pathString' 来解决我的问题，其中对于每个单元格，该单元格也存储它的所有祖先作为 'itself'。
然后我可以使用 R 中 data.tree 包中的 as.Node() 函数将我的网格转换为树对象，随后可以使用 as.phylo() 函数将其转换为 phylo 对象在猿。
一旦采样的单元格存储为树，ape 和 ggtree 中的现有函数使剩下的事情变得容易。
请参阅 data.tree 包的 data(acme) 和此处的 #Tree 示例：https://rdrr.io/cran/data.tree/man/as.Node.data.frame.html