我在一个拷贝区域中有4种基因型的观察计数样本。我想做的是计算这些基因型的等位基因频率,然后使用R中的卡方检验这些频率与预期值25%:25%:25%:25%的显著偏差。
到目前为止,我得到了:
> a <- c(do.call(rbind, strsplit(as.character(gdr18[1,9]), ",")), as.character(gdr18[1,8]))
> a
[1] "27" "30" "19" "52"
接下来我得到的总数:
> sum <- as.numeric(a[1]) + as.numeric(a[2]) + as.numeric(a[3]) + as.numeric(a[4])
> sum
[1] 128
现在的频率:
> af1 <- as.numeric(a[1])/sum
> af2 <- as.numeric(a[2])/sum
> af3 <- as.numeric(a[3])/sum
> af4 <- as.numeric(a[4])/sum
> af1
[1] 0.2109375
> af2
[1] 0.234375
> af3
[1] 0.1484375
> af4
[1] 0.40625
我现在迷路了。我想知道af1、af2、af3和af4是否显著偏离0.25、0.25、0.25和0.25
如何在R中执行此操作?
谢谢,Adrian
编辑:
好吧,我正在按照建议尝试chisq.test((:
> p <- c(0.25,0.25,0.25,0.25)
> chisq.test(af, p=p)
Chi-squared test for given probabilities
data: af
X-squared = 0.146, df = 3, p-value = 0.9858
Warning message:
In chisq.test(af, p = p) : Chi-squared approximation may be incorrect
警告信息试图告诉我什么?为什么近似值不正确?
为了测试这种方法,我选择了远低于预期的0.25:
> af=c(0.001,0.200,1.0,0.5)
> chisq.test(af, p=p)
Chi-squared test for given probabilities
data: af
X-squared = 1.3325, df = 3, p-value = 0.7214
Warning message:
In chisq.test(af, p = p) : Chi-squared approximation may be incorrect
在这种情况下,H0仍然没有被拒绝,即使这些值与预期的0.25值相差甚远。
observed <- c(27,30,19,52)
chisq.test(observed)
这表明这样的频率或比这更极端的频率将仅在大约0.03%的时间内偶然出现(p=0.0003172(
如果你的零假设不是四个类别之间的25:25:25:25分布,而是说问题是这些数据是否与3:3:1:9的预期有很大差异,你需要明确计算预期频率:
expected <- sum(observed)*c(3,3,1,9)/16
chisq.test(observed,p=c(3,3,1,9),rescale.p=TRUE)