monocle细胞轨迹追踪原理

作者: 小潤澤 | 来源:发表于2021-08-31 00:47 被阅读0次

我们用seurat选取HVG来进行细胞轨迹推断，并且过滤了一些低质量的细胞

加载相应的工具包

library(igraph)
library(monocle)

本文的数据和预处理来源于：单细胞之轨迹分析-2：monocle2 原理解读+实操

DDRTree树降维

假设采用DDRTree进行降维，DDRTree降维的意义是将高维空间的数据点降至二维，方便可视化

ncenter <- cal_ncenter(ncol(FM))

ddrtree_res <- DDRTree(FM, max_components, verbose = verbose,ncenter = ncenter)

## 其中ncenter代表经过树降维后所选取的细胞数
## ncol(FM)代表单细胞稀疏矩阵所测的细胞数


## cal_ncenter目的是筛选出 ncenter ，即低维空间（Z）中心或者根位置的一小群细胞的数量
cal_ncenter <- function(ncells, ncells_limit = 100){
  round(2 * ncells_limit * log(ncells)/ (log(ncells) + log(ncells_limit)))
}


## 定义：Y represents latent points as the center of Z
## Z代表树降维后，所有细胞在低维空间的新坐标
## Y代表树降维后处于低维空间（Z）中心或者根位置的一小群细胞
reducedDimK(cds) <- ddrtree_res$Y
## reducedDimK(cds) 代表树降维后处于低维空间（Z）中心位置的一小群细胞的点坐标

其中 reducedDimK(cds) <- ddrtree_res$Y 有如下几个属性：

定义：Y represents latent points as the center of Z

Z代表树降维后，所有细胞在低维空间的新坐标

Y代表树降维后处于低维空间（Z）中心位置的一小群细胞

reducedDimK(cds) 代表树降维后处于低维空间（Z）中心的一小群细胞的点坐标，也就是这 ncenter（126）个特征细胞在低维空间（Z）的坐标

其中FM为单细胞稀疏矩阵(count)：

FM稀疏矩阵

这里不清楚为什么要选取 ncenter 个细胞，推测是DDRTree降维后，选取贡献大的 ncenter 个细胞作为代表，该例子一共有126个细胞

计算最小生成树距离

当我们获得上一步的 reducedDimK(cds) ，那么下一步就是计算这些中心细胞（126个）的最小生成树，计算最小生成树的目的是衡量这些细胞之间的距离（以最小生成树的权重来表示）

adjusted_K <- Matrix::t(reducedDimK(cds))
dp <- as.matrix(dist(adjusted_K))

gp <- graph.adjacency(dp, mode = "undirected", weighted = TRUE)
dp_mst <- minimum.spanning.tree(gp)
minSpanningTree(cds) <- dp_mst

上述语句是利用了最小生成树算法，最小生成树是一种图的模型，可以衡量每个点之间的距离

上图每一个点可以理解为不同的细胞，而每条边上的权重可以理解为细胞之间的距离，那么最小生成树的特点是：1.两两细胞之间的权重和是最小的；2. 没有孤立的点（细胞），即所有细胞都是相连接的，任意细胞出发都可以到达所有细胞

下图即为最小生成树两两细胞之间的配对，每个细胞用Y来表示，一共126个细胞

上图表示为低维空间中心附近的细胞，其中一共有126个节点（细胞），125条边，igraph将细胞之间的距离关系用图来表示

软件支持自定义根细胞（根据生物学意义制定根细胞）
如果未定义根细胞，那么软件会自己定义为这126个细胞构成的网络中 网络直径所对应的起点细胞（网络直径在伪时间部分会讲解）为我们的根细胞，在该例子中为 Y_22

# 计算网络直径
diameter <- get.diameter(minSpanningTree(cds))
root_cell = names(diameter[1])

关于分支点

上面讨论了最小生成树的原理，那么最小生成树的特点是大部分细胞（点）只有一条边，仅有少数细胞（点）会有两条边，因此monocle的分支点选取的就是大于2条边的细胞（点）

mst_branch_nodes <- V(minSpanningTree(cds))[which(degree(minSpanningTree(cds)) > 2)]$name

## 选择最小生成树中大于2条边的细胞（点）

分支点
如上图所示的1，2，3即为分支点，该例子中对应的分支点细胞为：

轨迹线的画法

那么上图黑色的轨迹线是怎么画出来的呢？经过DDRTree降维后：

# reducedDimK(cds)代表这 126 个细胞在低维空间（二维）的坐标
reducedDimK(cds) <- ddrtree_res$Y
reduced_dim_coords <- reducedDimK(cds)

reducedDimK(cds)
reducedDimK(cds)代表reducedDimK(cds)代表这 126 个细胞在低维空间（二维）的坐标，Y_1，Y_2.....代表不同的细胞，一共126个，每个细胞对应下面的值代表二维坐标的横纵坐标

ica_space_df <- Matrix::t(reduced_dim_coords) %>%
    as.data.frame() %>%
    select_(prin_graph_dim_1 = x, prin_graph_dim_2 = y) %>%
    mutate(sample_name = rownames(.), sample_state = rownames(.))

ica_space_df
ica_space_df 的每一行表示不同的细胞对应的在低维（二维）空间的坐标，整理为一个大的文件

igraph::as_data_frame(dp_mst)

上图这张表交代了在最小生成树中两两细胞之间的联系（图的边连接两个细胞）

edge_df <- dp_mst %>%
    igraph::as_data_frame() %>%
    select_(source = "from", target = "to") %>%
    left_join(ica_space_df %>% select_(source="sample_name", 
                                       source_prin_graph_dim_1="prin_graph_dim_1", 
                                       source_prin_graph_dim_2="prin_graph_dim_2"), by = "source") %>%
    left_join(ica_space_df %>% select_(target="sample_name", target_prin_graph_dim_1="prin_graph_dim_1", 
                                       target_prin_graph_dim_2="prin_graph_dim_2"), by = "target")

edge_df

edge_df 表示的意义如下：

第一列表示不同的细胞，与第二列的细胞是用最小生成树的边相连接的，我们简称为最小生成树的细胞对

第二列表示不同的细胞，与第一列的细胞是用最小生成树的边相连接的，我们简称为最小生成树的细胞对

第三列和第四列表示的是第一列细胞在低维（二维）空间里的坐标

第五列和第六列表示对的第二列细胞在低维（二维）空间里的坐标

每一行的两个细胞代表最小生成树的细胞对，这两个细胞在二维空间的连线构成了其中的一条轨迹图

用这126个细胞画轨迹，推测是DDRTree降维后，选取贡献大的这 126 个细胞作为代表

事实上，黑色轨迹线即为 edge_df 中第一列和第二列细胞（对应最小生成树的细胞对）在低维（二维）空间两个细胞坐标之间的连线

当我们选取前30个细胞画时：

当我们选取前80个细胞画时：

当我们选取全部细胞画时：

画轨迹的代码：

g <- ggplot(data=data_df, aes(x=data_dim_1, y=data_dim_2)) 
g <- g + geom_segment(aes_string(x="source_prin_graph_dim_1", y="source_prin_graph_dim_2", 
                           xend="target_prin_graph_dim_1", yend="target_prin_graph_dim_2"), 
                           size=cell_link_size, linetype="solid", na.rm=TRUE, data=edge_df)

综上，轨迹线即为 edge_df 中第一列和第二列细胞（对应最小生成树的细胞对）在低维（二维）空间两个细胞坐标之间的连线

计算伪时间（Pseudotime）

伪时间在下面这幅图表示为颜色的深浅，蓝色越深代表伪时间越小，蓝色越浅代表伪时间越大，而上面的点表示所有的细胞（一共2638多个细胞）；轨迹线是用 edge_df 中第一列和第二列细胞（一共126个细胞，中心细胞或称为特征细胞）即最小生成树的细胞对在低维（二维）空间两个细胞坐标之间的连线画的

我们以这126个细胞为例子：

dp <- cellPairwiseDistances(cds)
dp_mst <- minSpanningTree(cds)

dp
其中，dp 代表的是两个细胞间的距离，值越大代表两个细胞距离比较大，而细胞间的距离推测使用基因表达来计算的

step_1：定义网络的直径

diameter <- get.diameter(minSpanningTree(cds))
# minSpanningTree(cds) 指的是最小生成树的细胞对

diameter
其中diameter所代表的轨迹即为该最小生成树所构成的图的直径，而 get.diameter 表示的是计算该图网络中最长的路径（该路径的权重之和应该是在树里面各个路径中是最小的）
网络直径:一个网络的直径被定义为网络中最长或最短路径
那么上图diameter的路径为由Y_22到Y_120（途经89个细胞，一共91个细胞），软件把这91个细胞（从Y_22开始到Y_120）所连成的轨迹定义为该网络的直径：
由 Y_22 出发到 Y_38（其中一个分支点）；由 Y_38 出发经过 Y_13 到 Y_86（其中一个分支点）；由 Y_86 出发经过 Y_73 到 Y_41（其中一个分支点）；由 Y_41 出发经过 Y_72 到 Y_120
以上路径被称为该最小生成树的直径，轨迹中最长的路径：

step_2 初始化表格，并定义父节点

curr_state <- 1

res <- list(subtree = dp_mst, root = root_cell)

states = rep(1, ncol(dp))
names(states) <- V(dp_mst)$name

# 伪时间的初始值为 0
pseudotimes = rep(0, ncol(dp))
names(pseudotimes) <- V(dp_mst)$name

parents = rep(NA, ncol(dp))
names(parents) <- V(dp_mst)$name

# 该例子中 root_cell 为 Y_22
mst_traversal <- graph.dfs(dp_mst,
                           root=root_cell,
                           neimode = "all",
                           unreachable=FALSE,
                           father=TRUE)
mst_traversal$father <- as.numeric(mst_traversal$father)
curr_state <- 1

其中值得注意的是 mst_traversal <- graph.dfs(dp_mst,root=root_cell,neimode = "all",unreachable=FALSE,father=TRUE) 是利用深度优先搜索来定义这126个细胞的父节点细胞，通过计算遍历这126个细胞会得到如下的顺序：

其中：

第一轨迹是从Y_22 到 Y_124 序号从 1到70的顺序轨迹（紫色路径）

第二条轨迹是从 Y_86—>Y_73—>Y_14.....Y120（省略号的位置代表序号 from 72 to 112；绿色路径）

第三条轨迹是从Y_41—>Y_111—>Y_63 序号从 79—>113—>114，这条轨迹比较短，在图中就没有用颜色进行区分

第四条轨迹是从Y_38—>Y_39—>Y_68.....Y_51（省略号的位置代表序号 from 116 to 126；橙色路径）

回顾之前的分支点 Y_38，Y_41，Y_86（大于两条边的细胞），第二条轨迹的 Y_86 是 Y_73 的父节点细胞，Y_73 是 Y_14 的父节点细胞；第三条轨迹的 Y_41 是 Y_111 的父节点细胞，Y_111 是 Y_63 的父节点细胞；第四条轨迹的 Y_38 是 Y_39 的父节点细胞，Y_39 是 Y_68 的父节点细胞

例如，Y_41多分支的情况，其中Y_41即是Y_72 的父节点细胞，也是Y_111 的父节点细胞

如果自定义根细胞，那么软件利用深度优先搜索出发的起点将会是用户自定义的根细胞
由于 root_cell = Y_22 所以搜索的起点为 Y_22 ，直到搜索完所有的细胞退出搜索，那么搜索的路径相连接后，就可以定义每个细胞的父节点细胞了，即相邻的两个细胞之间，前面的细胞是后面细胞的父节点细胞；比方说Y_22是Y_109的父节点细胞，Y_108是Y_105的父节点细胞，以此类推

深度优先搜索的意义是，从根细胞出发，沿着图的边先寻找与根细胞最近（权重最小的）细胞A，然后从该细胞A出发沿着图的边去寻找与细胞A最近（权重最小的）细胞B，以此类推，直到遍历完所有的细胞；那么它的生物学意义是，父节点细胞与子细胞的距离比较近，那么子细胞可能是由父节点细胞发育而来的，因此会存在一个发育的时间差，我们将这个时间差称之为伪时间

下面是多分枝的情况：

其中Y_41即是Y_72的父节点细胞，也是Y_111的父节点细胞

step_3 计算伪时间

为方便理解，先用126个中心细胞做伪时间计算的例子

for (i in 1:length(mst_traversal$order)){
# mst_traversal$order表示的是 diameter 中的的细胞（91个细胞）
  curr_node = mst_traversal$order[i]
# curr_node 即这91个细胞其中的一个
  curr_node_name = V(dp_mst)[curr_node]$name
  
  if (is.na(mst_traversal$father[curr_node]) == FALSE){
    # mst_traversal$father 为所有的126个细胞，即父节点细胞
    # 其中 curr_node 为 diameter 中的91个细胞之一
    # 如果 mst_traversal$father[curr_node] 存在，那么从父节点细胞中取出对应的 curr_node 细胞，并定义为 parent_node
    # curr_node 与 parent_node 是一个细胞
    parent_node = mst_traversal$father[curr_node]
    parent_node_name = V(dp_mst)[parent_node]$name
   
    # 从126个细胞的 pseudotimes 中取值，定义为parent_node_pseudotime，这 126 个细胞的初始 pseudotimes 都为 0
    parent_node_pseudotime = pseudotimes[parent_node_name]
    parent_node_state = states[parent_node_name]

   # 那么每个细胞的伪时间等于 parent_node_pseudotime 加上 curr_node_name 和 parent_node_name 细胞之间的距离
   # 比方说 Y_109 的伪时间等于 Y_22的伪时间（0）加上 Y_109与Y_22的距离（0.2273778）
    curr_node_pseudotime = parent_node_pseudotime + dp[curr_node_name, parent_node_name]
    if (degree(dp_mst, v=parent_node_name) > 2){
      curr_state <- curr_state + 1
    }
  }else{
    parent_node = NA
    parent_node_name = NA

   # 如果 mst_traversal$father[curr_node] 不存在，那么伪时间将定义为 0 
    curr_node_pseudotime = 0
  }
  
  curr_node_state = curr_state
  pseudotimes[curr_node_name] <- curr_node_pseudotime
  states[curr_node_name] <- curr_node_state
  parents[curr_node_name] <- parent_node_name
}

ordering_df <- data.frame(sample_name = names(states),
                          cell_state = factor(states),
                          pseudo_time = as.vector(pseudotimes),
                          parent = parents)
row.names(ordering_df) <- ordering_df$sample_name

比方说 Y_109 的伪时间等于 Y_22的伪时间（0）加上 Y_109与Y_22的距离（0.2273778）
又比方说 Y_108 的伪时间等于 Y_109的伪时间（0.2273778）加上 Y_108与Y_109的距离（0.2281914），距离的表示即为上述的dp矩阵

其中父节点细胞关系如下：

以此类推就可以计算出这126个中心细胞的伪时间了

对于全部细胞（2638多个细胞）伪时间的计算如下：
而全部细胞（2638多个细胞，上图蓝色的点）轨迹追踪计算方法和中心细胞的一样

# 计算全部细胞（2638多个细胞）的伪时间，这里的cds存储了全部细胞（2638多个细胞）的信息
extract_ddrtree_ordering <- function(cds, root_cell, verbose=T)
{
  
  # 计算全部细胞（2638多个细胞）细胞间的距离
  dp <- cellPairwiseDistances(cds)
  # 计算全部细胞（2638多个细胞）的最小生成树路径
  dp_mst <- minSpanningTree(cds)
  

 # 下面的代码是计算伪时间的，与126个中心细胞计算伪时间的方法一样
  curr_state <- 1
  res <- list(subtree = dp_mst, root = root_cell)
  states = rep(1, ncol(dp))
  names(states) <- V(dp_mst)$name
  
  pseudotimes = rep(0, ncol(dp))
  names(pseudotimes) <- V(dp_mst)$name
  
  parents = rep(NA, ncol(dp))
  names(parents) <- V(dp_mst)$name
  
  mst_traversal <- graph.dfs(dp_mst,
                             root=root_cell,
                             neimode = "all",
                             unreachable=FALSE,
                             father=TRUE)
  mst_traversal$father <- as.numeric(mst_traversal$father)
  curr_state <- 1
  
  for (i in 1:length(mst_traversal$order)){
    curr_node = mst_traversal$order[i]
    curr_node_name = V(dp_mst)[curr_node]$name
    
    if (is.na(mst_traversal$father[curr_node]) == FALSE){
      parent_node = mst_traversal$father[curr_node]
      parent_node_name = V(dp_mst)[parent_node]$name
      parent_node_pseudotime = pseudotimes[parent_node_name]
      parent_node_state = states[parent_node_name]
      curr_node_pseudotime = parent_node_pseudotime + dp[curr_node_name, parent_node_name]
      if (degree(dp_mst, v=parent_node_name) > 2){
        curr_state <- curr_state + 1
      }
    }else{
      parent_node = NA
      parent_node_name = NA
      curr_node_pseudotime = 0
    }
    
    curr_node_state = curr_state
    pseudotimes[curr_node_name] <- curr_node_pseudotime
    states[curr_node_name] <- curr_node_state
    parents[curr_node_name] <- parent_node_name
  }
  
  ordering_df <- data.frame(sample_name = names(states),
                            cell_state = factor(states),
                            pseudo_time = as.vector(pseudotimes),
                            parent = parents)
  row.names(ordering_df) <- ordering_df$sample_name
  # ordering_df <- plyr::arrange(ordering_df, pseudo_time)
  return(ordering_df)
}

值得注意的是
# 计算全部细胞（2638多个细胞）细胞间的距离
 dp <- cellPairwiseDistances(cds)
 # 计算全部细胞（2638多个细胞）的最小生成树路径
 dp_mst <- minSpanningTree(cds)  
计算的是所有细胞（2638多个细胞）的距离，这里的 dp 代表所有细胞（2638多个细胞）两两之间的距离
dp 每个值代表细胞间的距离，距离是基于基因表达量进行计算的
而dp_mst 是基于图并利用所有细胞（2638多个细胞）进行计算的最小生成树的轨迹
dp_mst ，之后类似于中心细胞（126个细胞）计算伪时间方法计算所有细胞（2638多个细胞）的伪时间

总结一下，黑色的轨迹线是利用126个中心细胞绘制的，而线旁边的点是所有的细胞（2638多个细胞）绘制的，而所有的细胞（2638多个细胞）点坐标的计算方式和126个中心细胞坐标计算方式一样，参考126个中心细胞的 edge_df

参考：关于图论中的最小生成树(Minimum Spanning Tree)详解

monocle细胞轨迹追踪原理

DDRTree树降维

计算最小生成树距离

关于分支点

轨迹线的画法

计算伪时间（Pseudotime）

step_1：定义网络的直径

step_2 初始化表格，并定义父节点

step_3 计算伪时间

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

机器学习算法

单细胞

单细胞测序