在R中计算含3种状态的网络可行路径数并绘制路径的方法

我写了一个R函数get.fit.land，用来绘制两种类型的基因型网络：

• 2状态网络：每个位点仅能取0或1
• 3状态网络：每个位点能取0、1或2

现有函数代码 #

par(mfrow=c(1,1))
get.fit.land <- function(nb.g = 3, # Number of loci
                         alleles.state = FALSE, 
                         main = "Fitness landscape",
                         draw.network = TRUE, 
                         add.geno.txt = TRUE) {
  # Empty list 
  l.gen = list()
  
  # 计算基因型状态总数
  nb.states = ifelse(alleles.state, 2^nb.g, 3^nb.g)
  
  if (nb.g > 10) {
    stop("Too many combinations for computation!")
  }
  
  # 生成每个位点的状态列表
  if (alleles.state) {
    for (i in 1:(nb.g)) {
      l.gen[[i]] <- 0:1
    }
    
  } else {
    for (i in 1:(nb.g)) {
      l.gen[[i]] <- 0:2
    }
    
  }  
  # 生成所有可能的基因型组合
  comb.gen = expand.grid(l.gen)
  comb.gen$sort = 1:nrow(comb.gen)
  comb.gen$sq.space = apply(comb.gen[,-which(names(comb.gen)=="sort")], 
                            MARGIN = 1, 
                            FUN = function(x) paste(x, collapse = ""))
  
  # 计算每个基因型的总突变数（所有位点数值之和）
  comb.gen$nb.mut = apply(X = comb.gen[,-which(names(comb.gen) %in% c("sort","sq.space"))], 
                          MARGIN = 1, 
                          FUN = sum)
  
  # 按突变数排序
  comb.gen.od = comb.gen[order(comb.gen$nb.mut),]
  
  # 初始化空图
  plot(0, type = "n",
       ylab = "",
       xlab = "",
       yaxt="n",
       xaxt="n",
       main = main, 
       xlim = c(-.50, max(comb.gen.od$nb.mut)+0.5), 
       ylim = c(-(max(table(comb.gen.od$nb.mut))-1)/1.5, 
                (max(table(comb.gen.od$nb.mut))-1)/1.5));abline(h = 0, lty = 2)
  
  # 为每个基因型分配绘图坐标
  pos_df = data.frame()
  unik_vals = unique(comb.gen.od$nb.mut)
  for (j in seq_along(unik_vals)) {
    current_mut = unik_vals[j]
    genos = comb.gen.od[comb.gen.od$nb.mut == current_mut, ]
    y_pos = seq(-(nrow(genos)-1)/2, (nrow(genos)-1)/2, length.out = nrow(genos))
    pos_df = rbind(pos_df, 
                   data.frame(sq.space = genos$sq.space,
                              x = current_mut,
                              y = y_pos,
                              geno_vec = I(lapply(1:nrow(genos), function(k) unlist(genos[k,1:nb.g])))))
  }
  
  # 绘制基因型节点
  points(pos_df$x, pos_df$y, pch = 19, cex = 1, col = scales::alpha("black", .5))
  
  # 添加基因型标签（如果开启）
  if(add.geno.txt){
    text(pos_df$x, pos_df$y, pos = 2, labels = pos_df$sq.space)
  }
  
  # 绘制合法路径：仅允许单个位点+1操作
  if (draw.network) {
    for (i in 1:nrow(pos_df)) {
      current_vec = pos_df$geno_vec[[i]]
      current_x = pos_df$x[i]
      current_y = pos_df$y[i]
      
      # 寻找下一个突变数层级的基因型
      target_x = current_x + 1
      target_genos = pos_df[pos_df$x == target_x, ]
      
      for (k in 1:nrow(target_genos)) {
        target_vec = target_genos$geno_vec[[k]]
        diff = target_vec - current_vec
        # 验证：仅一个位点+1，其余位点不变，且无数值减小
        if (all(diff >= 0) && sum(diff) == 1 && sum(diff == 1) == 1) {
          segments(x0 = current_x, y0 = current_y,
                   x1 = target_genos$x[k], y1 = target_genos$y[k],
                   col = scales::alpha("red", .4))
        }
      }
    }
  }
  
} # end function 

# 测试调用
get.fit.land(nb.g = 3, alleles.state = TRUE, draw.network = TRUE)  
get.fit.land(nb.g = 3, alleles.state = FALSE, draw.network = TRUE)  

现有网络图示 #

• 2状态网络（仅0/1）：
  
• 3状态网络（0/1/2）：
  

当前需求与问题 #

我需要绘制仅允许对任意位点执行+1操作的可行路径（即Hamming图，不能出现数值减小的情况）：

• 已知2状态下，3个位点的合法路径数是3! = 6条（每个位点依次从0→1的全排列）
• 不清楚3状态下的路径数计算方式
• 原有函数会绘制所有相邻节点的连接，不符合“仅+1”的规则，需要修改

解决方案 #

1. 合法路径数计算

• 2状态网络（每个位点仅0→1一次突变）：n个位点的合法路径数为n!，本质是n个位点的全排列，每个路径对应依次将不同位点从0变为1的顺序。
• 3状态网络（每个位点需0→1→2两次突变）：n个位点需要完成2n次+1操作，每个位点恰好被选中2次，属于多重排列问题，路径数公式为：

  (2n)! / (2!^n)
  

  以n=3为例，路径数为6!/(2!×2!×2!) = 720/8 = 90条。

2. 函数修改说明

修改后的函数核心优化点：

• 为每个基因型存储对应的数值向量，方便验证突变合法性
• 替换原有的全连接逻辑，仅保留满足“单个位点+1、其余位点不变”的节点连接
• 确保所有突变操作都是非递减的，避免出现数值减小的情况

修改后的函数会精准绘制符合规则的红色路径，同时保留原有节点展示和标签功能。

内容的提问来源于stack exchange，提问作者M. Beausoleil