PCR STR扩增模型技术咨询：批量返回多基因座数据及高效存储方案

一、批量生成21个基因座的返回数据 #

手动编写基因座列的方式效率极低，可通过定义基因座列表+动态列生成实现批量处理：

1. 先列出所有21个目标基因座（示例补充常见STR基因座，可按需替换）
2. 为每个基因座生成对应的7类数据列（Allele1/2、Stutter Allele1/2、Stutter2 Allele1/2、PCR效率）
3. 用循环批量映射数据到列

核心代码修改（返回部分） #

# 定义21个STR基因座列表（按需补充完整）
loci <- c("D3S1358", "vWA", "FGA", "TH01", "TPOX", "CSF1PO", "D5S818", 
          "D13S317", "D7S820", "D16S539", "D2S1338", "D19S433", "D18S51", 
          "D8S1179", "D21S11", "D10S1248", "D1S1656", "D12S391", "D22S1045", 
          "SE33", "PentaE")

# 批量构建所有基因座的列数据
result_df <- do.call(cbind, lapply(loci, function(locus) {
  data.frame(
    `Allele 1` = amplicon_count_allele1,
    `Allele 2` = amplicon_count_allele2,
    `Stutter Allele 1` = stutter_count_allele1,
    `Stutter Allele 2` = stutter_count_allele2,
    `Stutter2 Allele 1` = stutter2_count_allele1,
    `Stutter2 Allele 2` = stutter2_count_allele2,
    `PCR Efficiency` = efficiency_allele
  ) %>% 
    setNames(paste(locus, names(.), sep = " "))
}))

return(result_df)

二、更高效的数据存储与调用方式 #

原数据框结构不利于按基因座快速提取变量，推荐使用命名嵌套列表结构：每个基因座作为独立子列表，包含其所有扩增数据，调用时直接通过基因座名称+变量名访问，便捷性大幅提升。

优化后的完整函数 #

# 优化后的PCR扩增模型函数
pcr <- function(initial_efficiency, num_cycles) {
  prob_stutter <- 0.006
  
  # 初始化所有跟踪变量（预分配内存，提升运行效率）
  amplicon_count_allele1 <- numeric(num_cycles)
  amplicon_count_allele2 <- numeric(num_cycles)
  stutter_count_allele1 <- numeric(num_cycles)
  stutter_count_allele2 <- numeric(num_cycles)
  stutter2_count_allele1 <- numeric(num_cycles)
  stutter2_count_allele2 <- numeric(num_cycles)
  efficiency_allele <- numeric(num_cycles)
  
  efficiency_allele[1] <- initial_efficiency
  amplicon_count_allele1[1] <- 1
  amplicon_count_allele2[1] <- 1
  stutter_count_allele1[1] <- 0
  stutter_count_allele2[1] <- 0
  stutter2_count_allele1[1] <- 0
  stutter2_count_allele2[1] <- 0
  
  # PCR循环扩增逻辑
  for (i in 2:num_cycles) {
    # 等位基因扩增（修正原代码bug：使用上一轮的效率值）
    amplicon_count_allele1[i] <- amplicon_count_allele1[i-1] + 
      sum(rbinom(amplicon_count_allele1[i-1], 1, efficiency_allele[i-1]))
    amplicon_count_allele2[i] <- amplicon_count_allele2[i-1] + 
      sum(rbinom(amplicon_count_allele2[i-1], 1, efficiency_allele[i-1]))
    
    # 一级拖峰生成
    stutter_count_allele1[i] <- stutter_count_allele1[i-1] + 
      sum(rbinom(amplicon_count_allele1[i-1], 1, prob_stutter)) + 
      sum(rbinom(stutter_count_allele1[i-1], 1, prob_stutter))
    stutter_count_allele2[i] <- stutter_count_allele2[i-1] + 
      sum(rbinom(amplicon_count_allele2[i-1], 1, prob_stutter)) + 
      sum(rbinom(stutter_count_allele2[i-1], 1, prob_stutter))
    
    # 二级拖峰生成
    stutter2_count_allele1[i] <- stutter2_count_allele1[i-1] + 
      sum(rbinom(stutter_count_allele1[i-1], 1, prob_stutter)) +
      sum(rbinom(stutter2_count_allele1[i-1],1,prob_stutter))
    stutter2_count_allele2[i] <- stutter2_count_allele2[i-1] + 
      sum(rbinom(stutter_count_allele2[i-1], 1, prob_stutter)) +
      sum(rbinom(stutter2_count_allele2[i-1],1, prob_stutter))
    
    # 更新PCR效率
    total_amplicon <- amplicon_count_allele1[i] + amplicon_count_allele2[i] +
      stutter_count_allele1[i] + stutter_count_allele2[i] +
      stutter2_count_allele1[i] + stutter2_count_allele2[i]
    efficiency_allele[i] <- efficiency_allele[i-1] * exp(-1.12e-12 * total_amplicon)
  }
  
  # 定义21个STR基因座列表（按需补充完整）
  loci <- c("D3S1358", "vWA", "FGA", "TH01", "TPOX", "CSF1PO", "D5S818", 
            "D13S317", "D7S820", "D16S539", "D2S1338", "D19S433", "D18S51", 
            "D8S1179", "D21S11", "D10S1248", "D1S1656", "D12S391", "D22S1045", 
            "SE33", "PentaE")
  
  # 构建嵌套列表结构：每个基因座为一个子列表
  nested_result <- lapply(loci, function(locus) {
    list(
      Allele1 = amplicon_count_allele1,
      Allele2 = amplicon_count_allele2,
      Stutter_Allele1 = stutter_count_allele1,
      Stutter_Allele2 = stutter_count_allele2,
      Stutter2_Allele1 = stutter2_count_allele1,
      Stutter2_Allele2 = stutter2_count_allele2,
      PCR_Efficiency = efficiency_allele
    )
  })
  names(nested_result) <- loci
  
  # 同时返回嵌套列表（便于调用）和平铺数据框（便于查看）
  return(list(
    nested_data = nested_result,
    flat_data = do.call(cbind, lapply(loci, function(locus) {
      data.frame(
        `Allele 1` = amplicon_count_allele1,
        `Allele 2` = amplicon_count_allele2,
        `Stutter Allele 1` = stutter_count_allele1,
        `Stutter Allele 2` = stutter_count_allele2,
        `Stutter2 Allele 1` = stutter2_count_allele1,
        `Stutter2 Allele 2` = stutter2_count_allele2,
        `PCR Efficiency` = efficiency_allele
      ) %>% setNames(paste(locus, names(.), sep = " "))
    }))
  ))
}

调用示例 #

# 运行模型
model_result <- pcr(initial_efficiency = 0.9, num_cycles = 30)

# 提取D3S1358基因座的Allele1最终循环值
d3_allele1_final <- model_result$nested_data$D3S1358$Allele1[30]

# 提取vWA基因座的全程PCR效率数据
vwa_efficiency <- model_result$nested_data$vWA$PCR_Efficiency

内容的提问来源于stack exchange，提问作者Advent