当我们好不容易得到了还凑合的cell cluster，又做了pseudotime trajectory，然后如何为我们的细胞找到一个合适的路径呢？
比如在scanpy的trajectory inference教程中，run PAGA之后，得到了如下结果：
图片.png
并且后面出现了这样的path：
paths = [('erythrocytes', [16, 12, 7, 13, 18, 6, 5, 10]),
('neutrophils', [16, 0, 4, 2, 14, 19]),
('monocytes', [16, 0, 4, 11, 1, 9, 24])]
我很好奇，为什么erythrocyte的路径为上图，为啥不能是[16, 12, 7, 13, 18, 6,15, 3,17,10]等等其他组合呢？
通过询问小伙伴，他们告诉我，在得到paga的结果后，需要根据自己的background knowledge来确定path。作为一个半文盲的外科医学生，我并没有什么胚胎发育的background knowledge，此时我的心情是：
图片.png
然后，在PI大人的帮助下，我们开始作弊。（期间我还在没有任何基础的情况下被Gary忽悠硬去看了最大流最小割）
我们来看一下paga的结果：
数据来自scanpy tutorial
https://scanpy-tutorials.readthedocs.io/en/latest/paga-paul15.html
adata.uns['paga']
{'connectivities': <25x25 sparse matrix of type '<class 'numpy.float64'>'
with 70 stored elements in Compressed Sparse Row format>,
'connectivities_tree': <25x25 sparse matrix of type '<class 'numpy.float64'>'
with 22 stored elements in Compressed Sparse Row format>,
'groups': 'louvain_anno',
'pos': array([[ -828.751990457 , 5207.8173336207],
[ 7739.1243479772, 501.297186083 ],
[ 4710.602631487 , 6958.3421748608],
[ -1227.3342416429, -8336.2604694421],
[ 4014.7802194193, 4006.1770106812],
[ 238.0243621847, -5954.9137182842],
[ -3813.7119249991, -5807.4355041023],
[ -7147.3500300809, -1873.8699750282],
[ -7073.1481397142, 3800.1887626503],
[ 9185.0028776231, -1083.1673649154],
[ 2207.4910493818, -7830.277016848 ],
[ 7035.6931062009, 1687.6925420044],
[ -7058.0693477259, 2230.5216592432],
[ -6837.8460264989, -3101.0364176838],
[ 4175.5076840799, 9565.9115138965],
[ -4117.15456722 , -8148.2464387789],
[ -3097.2839569304, 3831.915714829 ],
[ 1850.347505207 , -9555.3736777067],
[ -6354.2341646177, -4841.4294204418],
[ 2992.4039702936, 11111.1153167889],
[ -9663.5270537571, 5103.458104299 ],
[ 5296.8206062395, -13759.453603083 ],
[ -732.970275698 , 2272.8366737242],
[ -9438.3334273383, 12743.9548016194],
[ 10113.6409562296, -3324.471455179 ]])}
发现它有两个martix，paga和paga tree 的connectivities matrix，大小为25x25，与cluster的数目对应，类型都是sparse matrix。
我们把第一个matrix取出来，然后取出其中的非0值的坐标，可以得到两个array
connection_matrix = adata.uns['paga']['connectivities'].todense()
ind = connection_matrix.nonzero()
ind
(array([ 0, 0, 0, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4,
4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 8, 8, 8, 8,
9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14,
14, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 18, 18, 18, 18, 19, 20,
22, 24]),
array([ 4, 11, 16, 4, 9, 11, 4, 14, 5, 6, 10, 15, 17, 0, 1, 2, 11,
16, 3, 6, 10, 13, 3, 5, 13, 15, 18, 12, 13, 18, 12, 16, 18, 20,
1, 11, 24, 3, 5, 17, 0, 1, 4, 9, 7, 8, 16, 5, 6, 7, 2,
19, 3, 6, 18, 0, 4, 8, 12, 22, 3, 10, 6, 7, 8, 15, 14, 8,
16, 9]))
其中，array1 为paga graph中的每一个node，array 为nodes之间的edges（这里的edge是双向的，即node0到node1的edge和node1到node0的edge是不同的）
我们把node和edge之间的对应关系取出，放到一个字典里
list_1 = [chr(x+97).upper() for x in ind[0]]
list_2 = [chr(x+97).upper() for x in ind[1]]
graph = {}
# traversal all the position and key
for pos,key in enumerate(list_1):
# wheter key is in the graph
if key in graph:
graph[key].append(list_2[pos])
else:
graph[key] = [list_2[pos]]
graph
{'A': ['E', 'L', 'Q'],
'B': ['E', 'J', 'L'],
'C': ['E', 'O'],
'D': ['F', 'G', 'K', 'P', 'R'],
'E': ['A', 'B', 'C', 'L', 'Q'],
'F': ['D', 'G', 'K', 'N'],
'G': ['D', 'F', 'N', 'P', 'S'],
'H': ['M', 'N', 'S'],
'I': ['M', 'Q', 'S', 'U'],
'J': ['B', 'L', 'Y'],
'K': ['D', 'F', 'R'],
'L': ['A', 'B', 'E', 'J'],
'M': ['H', 'I', 'Q'],
'N': ['F', 'G', 'H'],
'O': ['C', 'T'],
'P': ['D', 'G', 'S'],
'Q': ['A', 'E', 'I', 'M', 'W'],
'R': ['D', 'K'],
'S': ['G', 'H', 'I', 'P'],
'T': ['O'],
'U': ['I'],
'W': ['Q'],
'Y': ['J']}
(ps：这段优美简洁的代码是我跑去洲更老师的星球提问洲更老师写的，相比之下我自己写的简直没眼看，我还很无耻的发邮件催回答，然后发现大神不仅知识渊博，脾气还特别好，人美心善！)
我们给他起点和终点，寻找连接起点和终点的所有可能路径
def findAllPath(graph,start,end,path=[]):
path = path +[start]
if start == end:
return [path]
paths = [] #store all paths
for node in graph[start]:
if node not in path:
newpaths = findAllPath(graph,node,end,path)
for newpath in newpaths:
paths.append(newpath)
return paths
allpath = findAllPath(graph,'Q','K')
len(allpath)
178
我们找到了178条path
（PS2： findAllPath这段代码是我从网上搜到的，看了两遍没看懂，在纸上画了一边大概懂了，后来Gary告诉我这种写法叫做recursion，又去星球提问洲更老师，大神指点说这种嵌套的本质是图的遍历，又去看了BFS和DFS，好像又明白了一些些，忽然觉得临床医学还是很简单的）
然后，我们就可以想办法来从这些path中来条一条合眼缘的啦。
设想从node（x）跳到node（x+1），因为总要往下跳嘛，所以跳到所有可能的x+1的node的概率之和为1，根据这个想法，我们可以把connectivity matrix 转换成transition matrix
keylist = [chr(i) for i in range(97, 122)]
keylist = [item.upper() for item in keylist]
connection_dataframe=pd.DataFrame(connection_matrix.A,index = keylist,columns = keylist)
row_sum = connection_dataframe.sum(axis = 1)
propobility_matrix = connection_dataframe.div(row_sum, axis='rows')
propobility_matrix.sum(axis = 1)
然后计算每条路径的概率
result = [ ]
for j in range(len(allpath)):
times = 1
for i in range(len(allpath[j])-1):
times = times * propobility_matrix[allpath[j][i]][allpath[j][i+1]]
result.append(times)
print(result)
终于有两段是自己写的了。。我也就只能写写这种简单的小loop
挑一挑概率最大的那条，瞅一眼
result_sorted = sorted(result,reverse= True)
allpath[result.index(result_sorted[0])]
['Q', 'M', 'H', 'N', 'G', 'F', 'K']
tutorial中给的路径是
['Q', 'M', 'H', 'N', 'S', 'G', 'F', 'K']
我们的路径中少了S，查看一下graph中，S和N之间并没有边，所以这就是所谓的靠强大的background knowledge确定path吧。
把这种方法用到我自己的数据里，然后对应着去搜一搜文献，发现还挺make sense的。最后心得就是，虽然觉得make sense，还是得有背景知识做支撑，否则特别心虚。
idea来自PI大人，code来自网友，whatever， I learned a lot。
写在最后的话：
今天挺特别的，第一件事，今天早晨我的好朋友离开达拉斯去开始新的旅途啦，为他感到高兴，同时又很不舍。他是第一个开始带我看文章，做实验的人，告诉我应该怎么做事情，怎么规划自己的生活，在地图上帮我标好要去哪里玩，开车带我去踩点party，告诉我最重要的是要明白自己想要什么，随波逐流也可以，但是有可能会错过真正感兴趣的东西。萍水相逢一场，却在一开始就站在我的立场上帮我想问题。希望他在新的旅途上，能够看到更好的风景。
第二件事是下面的图：
2dc5bc41a84de5e8ef9f228dc69b5ed.jpg
我完成了第100个小时啦，一直觉得自己啥也没学会，很焦虑，现在想想，相比2个月前，我连python是什么都不知道，现在好像也还可以啊。（定100个小时也是因为Russell说其他的他可以教我，基础的100个小时是没有任何人能帮我完成的。现在他不能教我了，sad。但还是为他感到高兴）
米国进入混乱模式2.0，彻底打消了我road trip的心，踏实宅吧。
第二个100个小时结束后，希望我更喜欢我自己！
图片.png
以一张达拉斯的夜空结束今天的熬夜，不知道我在这里看了多少晚霞，一抬头就是满眼的绚丽或温柔。
网友评论