问题描述：Chargaff规则

核苷酸(nucleotide)是核酸的基本组成单位，它以一个含氮碱基为核心，加上一个五碳糖和一个或者多个磷酸基团组成，下面的图是我从维基百科扒来的，感觉非常清晰。含氮碱基有五种，分别是腺嘌呤（A）、鸟嘌呤（G）、胞嘧啶（C）、胸腺嘧啶（T）和尿嘧啶（U）。五碳糖为脱氧核糖的称为脱氧核糖核苷酸，是DNA的单体基本组成单位；五碳糖为核糖的称为核糖核苷酸，是RNA的基本组成单位。DNA中可以有的碱基是ATCG，RNA中可以有的碱基是AUCG。

核苷酸分布频率的规则是由Elson与Chargaff在1952年发现的（Elson, D, and E Chargaff. 1952. “On the Desoxyribonucleic Acid Content of Sea Urchin Gametes.” Experientia 8 (4). Springer: 143–45.）。下面是Chargaff的一些实验数据：

##                   A    T    C    G
## Human-Thymus   30.9 29.4 19.9 19.8
## Mycobac.Tuber  15.1 14.6 34.9 35.4
## Chicken-Eryth. 28.8 29.2 20.5 21.5
## Sheep-liver    29.3 29.3 20.5 20.7
## Sea Urchin     32.8 32.1 17.7 17.3
## Wheat          27.3 27.1 22.7 22.8
## Yeast          31.3 32.9 18.7 17.1
## E.coli         24.7 23.6 26.0 25.7

第一列表示某种生物的某个部位，每一行的四个数字表示这个部位的四种核苷酸的比例，下面是这些数据的柱状图：

Chargaff根据这些实验数据得到了一个结论：A的含量与T相同，C的含量与G相同，这个结论被称为Chargaff规则。这其实就是在高中生物中，我们学过的在DNA的结构中有一个碱基配对原则，因为DNA是双链结构，两条链上的碱基满足配对关系：A与T配对，C与G配对，于是$p_A=p_T,p_C=p_G$。

验证Chargaff规则的统计量

一个值得讨论的问题是$p_A=p_T,p_C=p_G$是否成立，用统计决策的方法建模，我们需要检验：
$$\begin{gathered} H_0:Chargaff规则不成立 \\ {H}_a:p_A=p_T,p_C=p_G \end{gathered}$$

我们可以回顾一下我们学过的假设检验工具：

根据我们需要做的假设检验，显然这是一个四总体的比例检验问题，因此我们应该用列联表。

如果不了解列联表方法，我们也可以尝试定义一个简单的统计量来验证Chargaff规则。定义$${\chi }^2=(p_A-p_T{)}^2+(p_C-p_G{)}^2$$

直观地理解一下这个统计量，在原假设下，这个统计量等于0，所以统计量的取值越小，我们越能信任原假设。

statChf = function(x){
    
    
  sum((x[, "C"] - x[, "G"])^2 + (x[, "A"] - x[, "T"])^2)
}
chfstat = statChf(ChargaffTable)
permstat = replicate(100000, {
    
    
     permuted = t(apply(ChargaffTable, 1, sample))
     colnames(permuted) = colnames(ChargaffTable)
     statChf(permuted)
})
pChf = mean(permstat <= chfstat)
pChf

## [1] 0.00014

说明
前三行定义的函数statChf作用是计算我们定义的统计量${\chi }^2$，第四行是用这个函数代入Chargaff的实验数据计算统计量${\chi }^2$的值；

第五到八行通过replicate函数对原数据做bootstrap，并用bootstrap样本计算${\chi }^2$统计量，得到${\chi }^2$的一个经验分布。第一个输入100000表示我们想得到100000组bootstrap样本，第二个输入表示我们希望用这些bootstrap样本执行{}中的语句，大概就是对每一行的比例做置换得到新的样本，然后用statChf函数计算${\chi }^2$统计量。

第九行到第十行是在根据经验分布计算检验的p-值，结果是0.00014，也就是说我们可以显著拒绝原假设，因此Chargaff规则成立。下面的柱状图表示经验分布，红线表示实验数据的${\chi }^2$统计量。

hist(permstat, breaks = 100, main = "", col = "lavender")
abline(v = chfstat, lwd = 2, col = "red")