生物学的现代统计方法与应用第一讲列联表1：验证Chargaff规则（碱基配对规则）

问题描述：Chargaff规则
验证Chargaff规则的统计量

问题描述：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 p_A=p_T,p_C=p_G pA=pT,pC=pG。

验证Chargaff规则的统计量

一个值得讨论的问题是 p A = p T , p C = p G p_A=p_T,p_C=p_G pA=pT,pC=pG是否成立，用统计决策的方法建模，我们需要检验：
H 0 : C h a r g a f f 规则不成立 H a : p A = p T , p C = p G H_0:Chargaff规则不成立\\ H_a:p_A = p_T, p_C = p_G H0:Chargaff规则不成立Ha:pA=pT,pC=pG

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

总体	检验均值	检验比例
单总体	Z检验、T检验	proportional z检验
两总体	Z检验、T检验	proportional z检验
多总体	ANOVA F检验	列联表卡方检验

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

如果不了解列联表方法，我们也可以尝试定义一个简单的统计量来验证Chargaff规则。定义 χ 2 = ( p A − p T ) 2 + ( p C − p G ) 2 \chi^2=(p_A-p_T)^2+(p_C-p_G)^2 χ2=(pA−pT)2+(pC−pG)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作用是计算我们定义的统计量 χ 2 \chi^2 χ2，第四行是用这个函数代入Chargaff的实验数据计算统计量 χ 2 \chi^2 χ2的值；

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

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

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