细胞状态转换图

library(gRbase)
library(igraph)
library(Rgraphviz)
#source("https://bioconductor.org/biocLite.R")
#biocLite(c("Rgraphviz", "RBGL"))
dag0 <- dag(~c0, ~ca*c0,~c0*ca,~c0*cx,~cx*c0,~ca*cx,~cx*ca,~cd*c0,~cd*cx,~cd*ca,~ca*ca,~c0*c0,~cd*cd)
dag0
## A graphNEL graph with directed edges
## Number of Nodes = 4 
## Number of Edges = 12
plot(dag0)

ps：用代码画得实在太丑了，不如手动！

细胞各个状态模拟

设定初始概率

  • 状态依次是c0、ca、cx、cd
  • 最后一个是死亡状态

假定开始细胞停留在c0状态：

init_prob <- c(1,0,0,0)

下面设定各个状态之间的转移概率，其中第四个cd是吸收状态：

tran_matrix <-matrix(c(0.1,0.4,0.4,0.1,0.2,0.2,0.2,0.4,0.3,0.3,0.2,0.2,0,0,0,1),byrow = T,nrow=4,ncol=4)
tran_matrix
##      [,1] [,2] [,3] [,4]
## [1,]  0.1  0.4  0.4  0.1
## [2,]  0.2  0.2  0.2  0.4
## [3,]  0.3  0.3  0.2  0.2
## [4,]  0.0  0.0  0.0  1.0

下面用一个\(400\times 4\)的矩阵显示细胞400次转移在各个状态的概率：

先初始化，全为0

stat_matrix <- matrix(c(0),nrow=400,ncol=4)

第一个是初始状态，位于c0

stat_matrix[1,] <- t(init_prob)

for (i in 2:400)
{
  stat_matrix[i,] <- stat_matrix[i-1,]%*% tran_matrix
}

细胞状态四百次转移位于各个状态的概率：

前二十次的状态

stat_matrix[1:20,]
##              [,1]        [,2]        [,3]      [,4]
##  [1,] 1.000000000 0.000000000 0.000000000 0.0000000
##  [2,] 0.100000000 0.400000000 0.400000000 0.1000000
##  [3,] 0.210000000 0.240000000 0.200000000 0.3500000
##  [4,] 0.129000000 0.192000000 0.172000000 0.5070000
##  [5,] 0.102900000 0.141600000 0.124400000 0.6311000
##  [6,] 0.075930000 0.106800000 0.094360000 0.7229100
##  [7,] 0.057261000 0.080040000 0.070604000 0.7920950
##  [8,] 0.042915300 0.060093600 0.053033200 0.8439579
##  [9,] 0.032220210 0.045094800 0.039791480 0.8828935
## [10,] 0.024178425 0.033844488 0.029865340 0.9121117
## [11,] 0.018146342 0.025399870 0.022413336 0.9340405
## [12,] 0.013618609 0.019062511 0.016821178 0.9504977
## [13,] 0.010220717 0.014306299 0.012624181 0.9628488
## [14,] 0.007670586 0.010736801 0.009474383 0.9721182
## [15,] 0.005756734 0.008057909 0.007110471 0.9790749
## [16,] 0.004320397 0.006047417 0.005336370 0.9842958
## [17,] 0.003242434 0.004538553 0.004004916 0.9882141
## [18,] 0.002433429 0.003406159 0.003005667 0.9911547
## [19,] 0.001826275 0.002556303 0.002255737 0.9933617
## [20,] 0.001370609 0.001918492 0.001692918 0.9950180
# stat_matrix
# stat_matrix[,1]
# plot(stat_matrix[,1],xlim=(0,1))

后二十次的状态：

stat_matrix[380:400,]
##               [,1]         [,2]         [,3] [,4]
##  [1,] 1.828165e-48 2.558949e-48 2.258071e-48    1
##  [2,] 1.372028e-48 1.920477e-48 1.694670e-48    1
##  [3,] 1.029699e-48 1.441308e-48 1.271841e-48    1
##  [4,] 7.727836e-49 1.081693e-48 9.545093e-49    1
##  [5,] 5.799698e-49 8.118049e-49 7.163540e-49    1
##  [6,] 4.352642e-49 6.092551e-49 5.376197e-49    1
##  [7,] 3.266634e-49 4.572426e-49 4.034806e-49    1
##  [8,] 2.451590e-49 3.431581e-49 3.028100e-49    1
##  [9,] 1.839905e-49 2.575382e-49 2.272572e-49    1
## [10,] 1.380839e-49 1.932810e-49 1.705553e-49    1
## [11,] 1.036312e-49 1.450563e-49 1.280008e-49    1
## [12,] 7.777463e-50 1.088640e-49 9.606390e-50    1
## [13,] 5.836943e-50 8.170182e-50 7.209543e-50    1
## [14,] 4.380593e-50 6.131676e-50 5.410722e-50    1
## [15,] 3.287611e-50 4.601789e-50 4.060717e-50    1
## [16,] 2.467334e-50 3.453617e-50 3.047546e-50    1
## [17,] 1.851721e-50 2.591921e-50 2.287166e-50    1
## [18,] 1.389706e-50 1.945222e-50 1.716506e-50    1
## [19,] 1.042967e-50 1.459879e-50 1.288228e-50    1
## [20,] 7.827408e-51 1.095631e-50 9.668080e-51    1
## [21,] 5.874427e-51 8.222649e-51 7.255841e-51    1