如何在R中重写基因组期望覆盖深度函数并优化运行效率？

原函数低效的原因 #

• 循环次数过多：1000万次循环在R中性能极差，R的for循环对大规模迭代的处理效率天生不足。
• 冗余计算：每次循环生成整个片段的fragment_length个位置再判断是否等于随机位置，完全没必要——单个read覆盖某固定位置的概率仅为fragment_length / genome，无需生成所有位置验证。

基于统计分布的高效实现 #

单个随机位置的覆盖深度本质服从二项分布：

• 试验次数n = reads（总reads数）
• 单次试验成功概率p = fragment_length / genome（单个read覆盖目标位置的概率）

当n极大且p极小时，二项分布可通过泊松分布近似，泊松参数lambda = n * p = (reads * fragment_length) / genome，计算效率更高且结果几乎一致。

代码实现 #

# 定义参数
genome_size <- 3e9
fragment_len <- 600
reads_num <- 10e6
num_samples <- 100000  # 用于绘制密度曲线的样本数量

# 方法1：泊松分布生成样本（最快）
set.seed(123)  # 设置随机种子保证结果可复现
lambda <- (reads_num * fragment_len) / genome_size
depth_samples_pois <- rpois(n = num_samples, lambda = lambda)

# 方法2：二项分布生成样本（更精确，速度仍极快）
depth_samples_binom <- rbinom(n = num_samples, size = reads_num, prob = fragment_len / genome_size)

# 绘制概率密度曲线
plot(density(depth_samples_pois), 
     main = "随机位置覆盖深度概率密度曲线",
     xlab = "覆盖深度",
     ylab = "概率密度",
     col = "blue",
     lwd = 2)
# 叠加二项分布曲线对比（几乎重合）
lines(density(depth_samples_binom), col = "red", lty = 2, lwd = 2)
legend("topright", legend = c("泊松近似", "二项分布"), col = c("blue", "red"), lty = c(1,2), lwd=2)

若需贴近“模拟read生成”逻辑（非必要） #

如果一定要模拟read起始位置并判断覆盖情况，可用向量化操作替代循环，避免逐次迭代：

depth_of_coverage_fast <- function(genome = 3E9, fragment_length = 600, reads = 10E6, target_pos = sample(1:genome, 1)) {
  # 生成所有read的起始位置
  start_pos <- sample(1:(genome - fragment_length + 1), reads, replace = TRUE)
  # 判断每个read是否覆盖目标位置
  covers <- (start_pos <= target_pos) & (target_pos <= start_pos + fragment_length - 1)
  # 计算覆盖深度
  sum(covers)
}

# 生成多个目标位置的覆盖深度样本
set.seed(123)
depth_samples <- replicate(100000, depth_of_coverage_fast())

这种方式比原函数快很多，但仍不如直接用统计分布高效，因为需要生成1000万个起始位置的向量，内存占用较大。

结果说明 #

给定参数下，lambda = (10e6 * 600) / 3e9 = 2，覆盖深度服从泊松(2)分布，大部分位置的覆盖深度集中在0-4之间，概率密度曲线符合泊松分布形态。

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