PNAS：模拟微生物群落互作及生命游戏在R中的实现

Listenlii-生物信息知识分享

发布于 2022-03-31 21:44:24

7910

发布于 2022-03-31 21:44:24

文章被收录于专栏：Listenlii的生物信息笔记Listenlii的生物信息笔记

Journal: PNAS

Publish: January 4,2022

Link:https://www.pnas.org/content/119/1/e2020956119

今天看到一篇文章，通过对空间参数和微生物生长状况进行参数设置，计算模拟了多个物种之间的相互作用关系及微生物多样性的形成机制。

参数设置：

图2 群落形成的动态模拟。

代码见：https://github.com/levifussell/MicroEvo

不过这篇文章不是本文的重点，而是其模拟方法让我瞬间想到了生命游戏（game of life）。

生命游戏是英国数学家约翰·何顿·康威在1970年发明的。

简单来说，对于一个网格状的空间，其中一些点可以有细胞存在。在下一时刻，细胞是否存在只依赖于其周围8个格子是否存在细胞。如果周围存在2-3个细胞，则下一时刻该细胞存在。如果周围细胞小于2个，则该细胞由于太过孤独而死亡。如周围细胞大于3个，则该细胞由于太过拥挤或者细胞间的竞争而死亡。

给定一个初始状态，生命就会相互交织纠缠，使得仅仅如此简单的规律即可产生让人叹为观止的生命现象。生命游戏可以产生很多有趣的图形，具体可自行百度。。。

而这篇PNAS似乎是生命游戏在微生物群落中的推广。通过给定参数，模拟群落在时间轴上的多样性变化。并利用随机森林考察了不同的参数对群落多样性的影响程度。

我在网上搜了一下还真搜到了R语言实现生命游戏的代码。看了之后发现思路并不难，有点后悔没有自己先思考一下就直接搜索了。

依据别人代码的思路，我也在R中实现了简单的生命游戏：

# Game of Life
# Refer to: https://zhuanlan.zhihu.com/p/136727731

### 构造初始状态：
set.seed(2022-2-21)
size = 20  # 矩阵的行和列数
d = round(runif(size*size,0,0.6)) # 最大值低一些，保证初始有值的少一些。
start = matrix(data=d,ncol=size,nrow=size)

### 求下一个状态时格子周围值的和：
fun = function(m,x,y,size){ # m为当前状态的矩阵；x和y为坐标；size为矩阵大小
  fun.sum = 0
  for (i in c(x-1,x,x+1)){ #依次遍历一个格子周围3x3的邻居格子
    for (j in c(y-1,y,y+1)){
      #如果格子在角落或者边，则邻居的值直接为0
      if (i  > 0 & i <= size & j > 0 & j <= size)  fun.sum = fun.sum + m[i,j]  # 把9个格子先求和
    }
  }
  fun.sum = fun.sum - m[x,y]  # 减去中间格子的值，即为周围8个值的和
}

### 设置运行次数
time = 100
life = list()
life[[1]] = start
for (k in 2:time){ # k = 3
  life.next = matrix(data = 0,ncol=size,nrow = size)
  for (i in 1:size){
    for (j in 1:size){
      fun.sum = fun(life[[k-1]],i,j,size) # 运行上述函数
      
      # 判断下个状态时当前位置是否有值存在。判断依据来自 http://nonoas.gitee.io/webproj/LifeGame/
      
      #孤单死亡：如果细胞的邻居小于等于1个，则该细胞在下一次状态将死亡；
      
      #拥挤死亡：如果细胞的邻居在4个及以上，则该细胞在下一次状态将死亡；
      
      #稳定：如果细胞的邻居为2个或3个，则下一次状态为稳定存活；
      
      #复活：如果某位置原无细胞存活，而该位置的邻居为2个或3个，则该位置将复活一个细胞。

      life.next[i,j] =  ifelse( fun.sum == 2|fun.sum == 3, 1, 0)
    }
    
  }
  life[[k]] = life.next
}


#输出gif：
# setwd("E:/桌面/")

install.packages("magick")
install.packages("animation")
library(animation)
library(pheatmap)   #加载pheatmap
saveGIF( 
  (for (k in 1:time){ 
      pheatmap(life[[k]],cluster_rows=F,cluster_cols=F,display_numbers=F,
           legend = F,cellwidth = 20, cellheight = 20,
           color = c("white","black"),main = "Life Game")
  }), movie.name = "LifeGame.gif"
)