R语言计算多态性信息含量(PIC)的代码公式正确性问询

• 你对PIC计算的核心逻辑理解没有偏差：双等位SNP场景下，仅需代入参考等位、第一替代等位的频率即可完成计算，暂不涉及第二替代等位时不需要额外扩展逻辑。
• 你当前的R代码实现存在算术错误，和标准PIC计算公式不符。

错误原因说明 #

双等位基因场景下，设p为替代等位基因频率（即你代码中的a2），q为参考等位基因频率，即q = 1 - p，代入标准PIC公式可化简为：

PIC = 1 - p² - q² - 2p²q²
将q=1-p替换后，公式为：
PIC = 1 - p² - (1-p)² - 2p²(1-p)²

你现有代码的错误出现在最后一项的计算：你写的是2*(var_freq$a2^2)*(1-(var_freq$a2^2))，错误将(1 - a2)²写成了1 - a2²，二者计算结果差异很大。

修正后的代码 #

你可以直接替换为以下正确的计算代码：

# 写法1：直接对应你的原有变量
var_freq$PIC <- 1 - var_freq$a2^2 - (1 - var_freq$a2)^2 - 2 * var_freq$a2^2 * (1 - var_freq$a2)^2

# 写法2：简化重复计算，可读性更高
p <- var_freq$a2
q <- 1 - p
var_freq$PIC <- 1 - p^2 - q^2 - 2 * p^2 * q^2

验证示例 #

举个简单的测试场景：当替代等位频率a2=0.5时，正确PIC计算结果为0.375，使用你原有错误代码计算得到的结果为0.125，差异明显，可验证修正后代码的正确性。

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