话说当年我在暨大工作的时候,学生都很害怕汇报RT-PCR的结果,因为每次一讲到用了t test做了检验,底下某小老板,就一顿狂喷,说学生瞎搞,这不对那不对,但从没告诉正确的姿势,学生们都怕了他。于是我写了这篇文章,普及了一下相关的知识,RT-PCR那点东西,全在这里。
RT-PCR相对定量有两个模型,一是 模型,一个是扩增效率校正模型。这里我们先讨论简单的模型: 模型,在这一模型中,假定扩增效率为2,即每个PCR cycle,产物倍增,由以下公式给出:
其中ΔΔCt = ΔΔCtreated - ΔΔCcontrol, ΔΔCtreated和ΔΔCcontrol分别是treatment组和control组中目标基因和参照基因的Ct差,即ΔΔCt = ΔΔCtarget - ΔΔCreference.
扩增效率肯定是在(1,2)之间,理想状态下才是2,而理想状态通常是不存在的。这个模型的假设点其实是treatment组和control组的扩增效率一样,至于是不是2,比2小多少,并不影响后续的统计分析。
我们来看以下一份数据,4组重复实验,用RT-PCR测了gene01和gene02的表达量,HK代表house keeping gene,即参照。以下数值为原始的Ct值。
ct <->->
sample = rep(rep(1:4, each=3), 2),
treatment = rep(c('Control','Treated'), each=12),
gene01 = c(23.22,23.34,23.12,24.06,24.15,24.15,23.18,23.13,
23.10,24.78,24.45,24.67,23.11,22.99,23.10,22.77,
22.99,23.06,23.73,24.01,23.80,23.73,23.83,23.73),
gene02 = c(29.08,29.04,29.39,28.23,28.01,28.12,28.79,28.43,
28.49,31.37,30.74,31.09,27.11,27.24,27.37,25.52,
25.72,25.52,27.43,26.73,26.65,27.96,27.84,27.98),
HK = c(19.68,19.69,19.80,19.95,19.93,19.97,19.93,20.02,20.27,
19.93,19.88,19.90,20.61,19.98,20.57,19.68,19.95,
19.85,20.27,20.08,20.07,20.10,20.07,20.10)
)
print(ct)
## sample treatment gene01 gene02 HK
## 1 1 Control 23.22 29.08 19.68
## 2 1 Control 23.34 29.04 19.69
## 3 1 Control 23.12 29.39 19.80
## 4 2 Control 24.06 28.23 19.95
## 5 2 Control 24.15 28.01 19.93
## 6 2 Control 24.15 28.12 19.97
## 7 3 Control 23.18 28.79 19.93
## 8 3 Control 23.13 28.43 20.02
## 9 3 Control 23.10 28.49 20.27
## 10 4 Control 24.78 31.37 19.93
## 11 4 Control 24.45 30.74 19.88
## 12 4 Control 24.67 31.09 19.90
## 13 1 Treated 23.11 27.11 20.61
## 14 1 Treated 22.99 27.24 19.98
## 15 1 Treated 23.10 27.37 20.57
## 16 2 Treated 22.77 25.52 19.68
## 17 2 Treated 22.99 25.72 19.95
## 18 2 Treated 23.06 25.52 19.85
## 19 3 Treated 23.73 27.43 20.27
## 20 3 Treated 24.01 26.73 20.08
## 21 3 Treated 23.80 26.65 20.07
## 22 4 Treated 23.73 27.96 20.10
## 23 4 Treated 23.83 27.84 20.07
## 24 4 Treated 23.73 27.98 20.10
每个sample,测了三次技术重复,我们使用平均值来提高精度。
require(plyr)
## Loading required package: plyr
ct <- ddply(ct,="" .(sample,="" treatment),="" function(x)="">->
data.frame(gene01=mean(x$gene01),
gene02=mean(x$gene02), HK=mean(x$HK)))
print(ct)
## sample treatment gene01 gene02 HK
## 1 1 Control 23.22667 29.17000 19.72333
## 2 1 Treated 23.06667 27.24000 20.38667
## 3 2 Control 24.12000 28.12000 19.95000
## 4 2 Treated 22.94000 25.58667 19.82667
## 5 3 Control 23.13667 28.57000 20.07333
## 6 3 Treated 23.84667 26.93667 20.14000
## 7 4 Control 24.63333 31.06667 19.90333
## 8 4 Treated 23.76333 27.92667 20.09000
按照ΔΔCt 模型,我们先计算ΔCt,然后计算ΔΔCt,最终计算出Ratio。
delta.ct <- ddply(ct,="" .(sample,="" treatment),="" function(x)="">->
data.frame(gene01=x$gene01-x$HK,
gene02=x$gene02-x$HK))
print(delta.ct)
## sample treatment gene01 gene02
## 1 1 Control 3.503333 9.446667
## 2 1 Treated 2.680000 6.853333
## 3 2 Control 4.170000 8.170000
## 4 2 Treated 3.113333 5.760000
## 5 3 Control 3.063333 8.496667
## 6 3 Treated 3.706667 6.796667
## 7 4 Control 4.730000 11.163333
## 8 4 Treated 3.673333 7.836667
## Delta Delta Ct
dd.ct <- subset(delta.ct,="" treatment="='Treated',">->
select=c(gene01, gene02)) -
subset(delta.ct,
treatment=='Control',
select=c(gene01, gene02))
print(dd.ct)
## gene01 gene02
## 2 -0.8233333 -2.593333
## 4 -1.0566667 -2.410000
## 6 0.6433333 -1.700000
## 8 -1.0566667 -3.326667
ratio <->->
print(ratio)
## gene01 gene02
## 2 1.769490 6.034915
## 4 2.080120 5.314743
## 6 0.640232 3.249010
## 8 2.080120 10.032899
这个ratio值计算出,经过扩增后,目的基因在处理组中的产量是对照组的多少倍。
计算出ratio后,很多人就拿它来做t检验。这种做法是不对的,因为ratio的分布明显是右偏的。从理论上看,PCR的扩增可以分成以下三个阶段,指数扩增、线性扩增、平台期。
当试剂量足够的时候,就处于指数扩增阶段,只有当PCR试剂不足时,才进入线性扩增期,最后当试剂被耗尽时,产物不增长,位于平台期。
在我们的实验中,通常PCR反应是处于指数增长期。我们将其进行对数处理,指数增长就会变成线性增长,这样数据的右偏现象就会有明显减弱,当然我不敢说对数转换后,它就能呈正态分布,但起码对数处理后数据的分布比原始ratio数据更接近正态了。
只有偏态不明显的情况下,才可以用t检验进行分析。
## T test for gene01
t.test(log(ratio[,1], base=2), mu = 0)
##
## One Sample t-test
##
## data: log(ratio[, 1], base = 2)
## t = 1.4009, df = 3, p-value = 0.2558
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.729139 1.875806
## sample estimates:
## mean of x
## 0.5733333
## T test for gene02
t.test(log(ratio[,2], base=2), mu = 0)
##
## One Sample t-test
##
## data: log(ratio[, 2], base = 2)
## t = 7.5039, df = 3, p-value = 0.004904
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 1.44405 3.57095
## sample estimates:
## mean of x
## 2.5075
如果我们要使用T检验,就得做对数处理,细心的人应该会发现ratio的计算就是指数运算,ratio取对数之后,就是−ΔΔCt,所以我们可以直接使用它来做检验。
## T test for gene02
t.test(-dd.ct[,2])
##
## One Sample t-test
##
## data: -dd.ct[, 2]
## t = 7.5039, df = 3, p-value = 0.004904
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 1.44405 3.57095
## sample estimates:
## mean of x
## 2.5075
另一种方法,使用ANOVA进行统计,这也是参数统计方法,在我们这种小样本的情况下,同样对数据分布有严格的要求,我们不能基于ratio做统计,要基于ΔCt值。
print(delta.ct)
## sample treatment gene01 gene02
## 1 1 Control 3.503333 9.446667
## 2 1 Treated 2.680000 6.853333
## 3 2 Control 4.170000 8.170000
## 4 2 Treated 3.113333 5.760000
## 5 3 Control 3.063333 8.496667
## 6 3 Treated 3.706667 6.796667
## 7 4 Control 4.730000 11.163333
## 8 4 Treated 3.673333 7.836667
summary(aov(gene01 ~ treatment, data=delta.ct))
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 0.6574 0.6574 1.687 0.242
## Residuals 6 2.3385 0.3898
summary(aov(gene02 ~ treatment, data=delta.ct))
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 12.575 12.575 9.963 0.0197 *
## Residuals 6 7.573 1.262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA也是比较常用的,而且可以应用于多组实验和多个因素中。但是如果按照上面的代码来做,还不如用T test,因为实验是成对(tretated vs control)进行的,而这个信息在上面的ANOVA分析中,被扔了。
实际上 −ΔΔCt 的单样本T检验,相当于ΔΔCt 的双样本T检验。
t.test(gene02 ~ treatment, delta.ct, paired=T)
##
## Paired t-test
##
## data: gene02 by treatment
## t = 7.5039, df = 3, p-value = 0.004904
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.44405 3.57095
## sample estimates:
## mean of the differences
## 2.5075
而这里用ANOVA,没有用到成对的信息,计算出来的p值和unpaired T test就基本相等。
t.test(gene01 ~ treatment, delta.ct, paired=F)
##
## Welch Two Sample t-test
##
## data: gene01 by treatment
## t = 1.2988, df = 5.2396, p-value = 0.2482
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.5460188 1.6926855
## sample estimates:
## mean in group Control mean in group Treated
## 3.866667 3.293333
t.test(gene02 ~ treatment, delta.ct, paired=F)
##
## Welch Two Sample t-test
##
## data: gene02 by treatment
## t = 3.1565, df = 5.064, p-value = 0.02476
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.473167 4.541833
## sample estimates:
## mean in group Control mean in group Treated
## 9.319167 6.811667
做为成对T检验的扩展,使用方差分析进行多组实验的比较,需要用到的是重复测量方差分析(repeated-measures ANOVA)。
delta.ct[,1] = factor(delta.ct[,1])
delta.ct[,2] = factor(delta.ct[,2])
summary(aov(gene01 ~ treatment+Error(sample/treatment), data=delta.ct))
##
## Error: sample
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 3 1.333 0.4445
##
## Error: sample:treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 0.6574 0.6574 1.962 0.256
## Residuals 3 1.0050 0.3350
summary(aov(gene02 ~ treatment+Error(sample/treatment), data=delta.ct))
##
## Error: sample
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 3 6.903 2.301
##
## Error: sample:treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 12.57 12.575 56.31 0.0049 **
## Residuals 3 0.67 0.223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
这样算出来的p值就和成对T检验几乎一样。
我们计算出来的ratio,是经过PCR扩增之后产物的比值,这是经过指数扩增的,而如果我们取对数,相当于消除了指数增长的效果,−ΔΔCt 其实就是原始量的差值。
所以ratio的比较,相当于是扩增后的比较,而log2ratio 的比较,则相当于是扩增之前的比较。这两种数据当然都是可以拿来做统计分析的,只不过T test等参数检验方法对数据分布有要求,不适合于ratio值的统计分析而已。如果我们使用非参数检验,则不管是ratio值也好,log2ratio值也好,结果都是一样的。非参数检验并不使用数据的原始值,而是使用其排序值来做检验,我们并不能确定对数处理后的数据就呈正态分布,还有就是通常我们的实验重复次数是非常少的,在这种情况下,使用非参数检验,其实也是一种保守的、不容易犯错的方法。
## Wilcoxon Rank Sum test for gene02, original data
wilcox.test( ratio[,2], mu=1)
##
## Wilcoxon signed rank test
##
## data: ratio[, 2]
## V = 10, p-value = 0.125
## alternative hypothesis: true location is not equal to 1
## Wilcoxon Rank Sum test for gene02, log2 transformed data
wilcox.test( log(ratio[,2], base=2), mu=0)
##
## Wilcoxon signed rank test
##
## data: log(ratio[, 2], base = 2)
## V = 10, p-value = 0.125
## alternative hypothesis: true location is not equal to 0
说它不容易犯错,是因为它对数据分布没要求;说它保守,是因为在计算p-value上,灵敏度没有T test这么高,p-value会大一些。
生物学的数据通常比较messy,non-normal,重复次数又少。在我们不能确定数据分布是否正态的情况下,使用参数检验,很可能结果bias太大,使用非参数检验,又觉得power不够大,那么我们可以使用permutation test。这不失为一种好方法,通过随机变换两组数据的标签,观察它们的均值差,我们将得到均值差的分布,将我们实际的测量值对比于随机得到的差值分布,算出随机得到的差值大于实际的测量值的概率。
perm.test <- function(d1,="" d2,="" niter="">->10000) {
if (length(d2) == 1)
d2 <- rep(d2,="">->length(d1))
m <->->
pooledData <- c(d1,="">->
n <->->length(d1)
meanDiff <->->
for (i in 1:nIter) {
idx <->->1:length(pooledData), n, replace=FALSE)
d1 <->->
d2 <->->
meanDiff[i] <- mean(d2)="" -="">->
}
p <->->abs(meanDiff) >= abs(m))
return(p)
}
perm.test(ratio[,1],1)
## [1] 0.1388
perm.test(-dd.ct[,1], 0)
## [1] 0.1454
perm.test(ratio[,2],1)
## [1] 0.0288
perm.test(-dd.ct[,2], 0)
## [1] 0.0266
从结果上看,用ratio还是用−ΔΔCt 差不多,这表明原始数据的分布完全不关事。
没有吃没有穿
自有那敌人送上前
没有枪没有炮
敌人给我们造
----游击队之歌
这个方法的奇妙之处,就在于通过bootstrap,产生针对我们统计量的一个背景分布,再由这个背景分布计算出p-value。
至于画图,我推荐使用−ΔΔCt 来画,这个CONTROL是0,不用画。如果要用ratio来画,你觉得CONTROL每次都是画一个mean=1,sd=0的柱子,有意义吗?anyway,我下面的画图代码,同样也适用于拿来画ratio值。
BTW:这里图中标的p值,使用的是T检验算出来的p值。
require(reshape2)
## Loading required package: reshape2
dd.ctm <->->
## No id variables; using all as measure variables
print(dd.ctm)
## variable value
## 1 gene01 0.8233333
## 2 gene01 1.0566667
## 3 gene01 -0.6433333
## 4 gene01 1.0566667
## 5 gene02 2.5933333
## 6 gene02 2.4100000
## 7 gene02 1.7000000
## 8 gene02 3.3266667
dd.cts <- ddply(dd.ctm,="" .(variable),="">->function(x)
data.frame(mean=mean(x$value), sd=sd(x$value)))
print(dd.cts)
## variable mean sd
## 1 gene01 0.5733333 0.8185353
## 2 gene02 2.5075000 0.6683222
require(ggplot2)
## Loading required package: ggplot2
p <->->ggplot(dd.cts, aes(variable, mean, fill=variable, width=.5))+
geom_bar(stat='identity', position=position_dodge(width=.8), colour='black')+
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2, position=position_dodge(width=.8))
p <->->p+scale_fill_manual(values=c('white', 'grey80'))
p <->->p+theme_bw()+theme(legend.position='none') + xlab('')+ylab(expression(paste('-', Delta, Delta, 'Ct')))
p <->->p+theme(axis.text.x=element_text(face='bold', size=14), axis.text.y=element_text(face='bold', size=14), axis.title.y=element_text(size=18, face='bold'))
p <->->p+annotate('text', x=2,y=3.5, label='p <>, size=4)
print(p)
再来一个箱式图,我其实觉得箱式图比柱状图要好一些,这里我除了用箱式图画数据分布之后,把原始的数据点也给画了。
p2 <- ggplot(dd.ctm,="" aes(variable,="" value,="" fill="variable," colour="">->
p2 <- p2+geom_boxplot(alpha=".3," width="">->'y', stackdir='center', dotsize=.5)
p2 <- p2+theme_bw()+theme(legend.position="">->'none') + xlab('')+ylab(expression(paste('-', Delta, Delta, 'Ct')))
p2 <- p2+theme(axis.text.x="">->'bold', size=14), axis.text.y=element_text(face='bold', size=14), axis.title.y=element_text(size=18, face='bold'))
p2 <->->'text', x=2,y=3.5, label='p <>, size=4)
print(p2)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
据我的观察,生物坑里大家对T检验特别有爱,我甚至看到用T检验去分析免疫组化数据的,而且还发在Nature子刊上,我都吐槽无力了。如果你喜欢T检验,你要注意数据分布,那么在分析RT-PCR数据上,请用−ΔΔCt作为统计量,或者将你算完的ratio值取对数也是一样的。用ANOVA时,也必须同样注意。
在这个重复次数少,数据又非常messy的生物世界里,保守点的做法是使用非参数统计,即使你使用ratio这种严重右偏的数据,你的统计分析照样没问题,多学点总是好的,凡事不能太依赖于T检验,因为在很多其它的场景里,正如这里的ratio值,你不能简单套用T检验。
Permutation test其实才是工具箱里不可或缺的一件,应用场景和非参一样广,同时灵敏度又比非参高,你值得拥有。
联系客服