R中求和的正确公式，计算多态信息含量（PIC）

Question

（基因组学数据行话警告！）我正在从数千个单核苷酸多态性 (SNP) 的数据集中计算 R 中的多态性信息内容 (PIC)。进行计算所需的数据是等位基因频率。每个观察值都有两种（很少三种）等位基因类型，或者是参考基因类型，或者是第一个替代基因类型。我有一小部分带有第二个替代等位基因。这是我试图在 R 中编写的公式：original publication

其中 Pi 和 Pj 是所选 SNP 标记的第 i 个和第 j 个等位基因的频率。这是我目前在 R

中的表述

var_freq$PIC <- (1-(var_freq$a2^2)-(1-var_freq$a2)^2)-(2*(var_freq$a2^2)*(1-(var_freq$a2^2))))

其中变量a2 是替代等位基因频率。

测试数据集：a1是参考等位基因，a2是替代（仍然忽略a3第二替代等位基因）

library(tidyverse)
set.seed(123)
testdata <- data.frame(a1=rnorm(n=10000, mean = .99, sd=0.15)) %>% 
  filter(., a1<1&a1>0) %>% 
  mutate(., a2=1-a1)

这是 PIC 值的正确 R 代码公式吗，即我对公式的解释正确吗？

Answer 1

如果问题是如何计算给定向量 p 的 PIC 则

1) 从 1 中减去 p 的元素的平方和，并减去 o 矩阵的上三角形中的值之和的两倍。

p <- c(0.1, 0.2, 0.3, 0.4)  # test data

o <- outer(p^2, p^2)
PIC <- 1 - sum(p^2) - 2 * sum(o[upper.tri(o)]); PIC
## [1] 0.6454

1a) 因为 o 是对称的，所以上三角形和下三角形是相同的，所以上三角形的两倍等于除对角线之外的所有元素之和，我们减去它（即添加它，因为它已经否定）所以一个变体是：

o <- outer(p^2, p^2)
PIC <- 1 - sum(p^2) - sum(o) + sum(diag(o)); PIC
## [1] 0.6454

2) 交替使用 listcompr 包，该包有助于编写非常接近问题中所示公式的代码。

library(listcompr)
n <- length(p)
1 - sum(p^2) - 2 * sum(gen.vector(p[i]^2 * p[j]^2, i = 1:(n-1), j = (i+1):n))
## [1] 0.6454

双重检查

# double check
1 - (.1^2+.2^2+.3^2+.4^2) - 2*(.1^2*(.2^2+.3^2+.4^2) + .2^2*(.3^2+.4^2) + .3^2*.4^2)
## [1] 0.6454

Answer 2

在数学上，这相当于：

p <- c(0.1, 0.2, 0.3, 0.4)  # test data

1 - sum((s <- p%o%p) * `diag<-`(s, 1))
[1] 0.6454