R语言之于定量PCR分析

这里以几位小朋友的练手实验数据为例，对20个样本的Gapdh、Cyp2e1、Atf4三个基因在两块96孔板上做了qPCR。更多信息在后续分析过程中注解。使用的是LightCycler® 96定量PCR仪，使用Roche配套软件进行分析获取Ct值并将两块板子的数据批量导入到了同一个文件中（Result Batch Export）。由于后续设置样本名等更容易，导出数据之前并未对每个孔设置样本名分组等信息。

首先加载需要用到的R包，并且读入数据

options(stringsAsFactors=FALSE)
library(reshape)               # 用于数据变形
library(referenceIntervals)    # 使用其dixon.outliers函数检验异常值
library(plyr)
library(ggplot2)
 
dat <- read.delim("ftp://ftp.biolab.wang/data/LC96qPCRdata.txt")

看一下数据的情况

这么多列数据，有用的并不多，File.Name、Position、Cq三列是需要保留的。 这里需要根据自己的情况选择读入的文件和数据。

根据实验记录，将P1的1-60孔为Gapdh基因，P2的1-60孔为Cyp2e1基因，P1的61-90及P2的61-90为Atf4基因，对应的小鼠编号都是1-20号，除了采样时的时候1号鼠和4号鼠装管颠倒了。
添加样本名和基因名比较关键，需要根据自己的实际情况添加对应好。

# 取出所需的列，并且将行的顺序调好
dat <- dat[c(1:60, 1:60 + 96, 61:90, 61:90 + 96), c("File.Name", "Position", "Cq")]
 
# 增加一列记录对应鼠的编号，M01-M20，每个鼠同一个基因占3个孔
dat$Mouse.ID <- sprintf("M%02d", rep(1:20, each=3))  
 
# 由于M01和M04鼠颠倒，这里还需要加一步。如果不理解如下一行，可以先用dat[dat$Mouse.ID %in% c("M01", "M04"), ]看看
dat[dat$Mouse.ID %in% c("M01", "M04"), ]$Mouse.ID <- rep(c("M04", "M01"), each=3)
 
# 增加基因名
dat$Gene.Name <- rep(c("Gapdh", "Cyp2e1", "Atf4"), each=60)
 
# 有小朋友们初学，加样的时候Atf4基因14号鼠的和15号鼠的加反了，因此这里需要加一步
dat[dat$Gene.Name =="Atf4" & dat$Mouse.ID %in% c("M14", "M15"), ]$Mouse.ID <- rep(c("M15", "M14"), each=3)

需要对3个重复的孔进修过滤和求平均，每连续的3个孔属于同一只鼠的同一个基因，给这些孔对应ID设为D1-D3

# 添加一列，记录重复的ID
dat$Dup.ID <- c("D1", "D2", "D3")
 
# 如果Cq一列不是数值类型，要当心，如果这一列是factor而直接转为数值的话，会出错
if(class(dat$Cq) != "numeric"){
    dat$Cq <- as.numeric(as.character(dat$Cq))
}
 
# 对于没有数值的孔，将它设置为40
if(sum(is.na(dat$Cq)) > 0){
    dat$Cq[is.na(dat$Cq)] <- 40
}
 
# 对数据进行变形，使得同一个鼠的同一个基因的三个孔的数据在同一行，方便处理
Dat <- data.frame(cast(dat, Gene.Name + Mouse.ID ~Dup.ID, value="Cq"))
 
# 这里可以用head(Dat)看看变形之后的效果
 
# 先看看数据质量如何，是否有的3个重复孔之中存在异常值，使用dixon.outliers函数，自己看看相关说明
Dat$OL <- apply(Dat[ , c("D1", "D2", "D3")], 1, function(x){
    if(sd(x) == 0) {
        return("")
    } else {
        return(paste(dixon.outliers(x)$outliers, collapse=","))
    }
})
 
Dat[Dat$OL != "", ]  #这样查看具有异常值的列
 
# 考虑到异常值的情况下来求Ct的平均值
# 使用dixon.outliers函数来查看三个重复中的异常值
# 如果存在异常值，并且异常值与其它两个值的平均值差距大于0.5，那么就把异常值去掉之后再求平均值
Dat$Ct <- apply(Dat[ , c("D1", "D2", "D3")], 1, mean)
 
# 如果存在异常值，那么需要考虑有异常值的情况
if(sum(Dat$OL != "") > 0){
    Dat[Dat$OL != "", ]$Ct <- apply(Dat[Dat$OL != "", c("D1", "D2", "D3")], 1, function(x){
        Ou <- dixon.outliers(x)$outliers      # 异常值
        Ou <- as.numeric(as.character(Ou))    # 有时候返回的异常值是字符型的，如果还报错，看看下一行是不是也要转换
        Su <- dixon.outliers(x)$subset        # 除去异常值之后剩下的数值
        if(abs(Ou - mean(Su)) > 0.5){
            return(mean(Su))                  # 如果异常值与非异常值的平均数之差大于0.5，则去掉异常值
        } else {
            return(mean(x))                   # 否则还是使用全部3个数的平均值
        }
    })
}

我们看看具有异常值的行的平均值计算情况

很显然，一些过于偏离的异常值被排除了，没有用于计算平均的Ct值。

这里Gapdh是参照基因，把Gapdh拿出来作为一列

# Gapdh的数据单独拿出来
Dat1 <- Dat[Dat$Gene.Name == "Gapdh", c("Mouse.ID", "Ct")]
names(Dat1) <- c("Mouse.ID", "Cr")   # r 代表 reference gene
 
# 和其它两个目标基因的数据merge在一起
DAT <- merge(Dat[Dat$Gene.Name != "Gapdh", c("Gene.Name", "Mouse.ID", "Ct")], Dat1, by="Mouse.ID", all.x=TRUE)
 
# 求ΔCt，用基因的Ct减去参照基因的Ct
DAT$dCt <- DAT$Ct -DAT$Cr
 
# 接下来秋ΔΔCt， 这个针对每个基因单独求，因此先求出各个基因最大的ΔCt
DAT1 <- ddply(DAT, .(Gene.Name), summarize, Cm = max(dCt))   # Cm表示最大的dCt
 
# 再把各个基因对应的最大的ΔCt与各个样本的ΔCt对应merge在一起
DAT <- merge(DAT, DAT1, by="Gene.Name", all.x=TRUE)
 
# 求ΔΔCt
DAT$ddCt <- DAT$dCt - DAT$Cm   #
 
# 通过ΔΔCt来求每个样本对应的相对表达量
DAT$Fc <- 2^(-DAT$ddCt)   # Fc表示Fold Change，这里是相对表达量

这样，相对的表达量就计算出来了。

到此，特别针对qPCR的分析就结束了，接下来对这些数据的分析就和别的数据一样了。可以根据自己的实验设计进行各种检验。

求出了每个小鼠各个基因的相对表达量，接下来就要看看各个组之间是不是具有差异，那么再把小鼠对应到相应的处理组采样时间等等上边去。

# 首先根据实际情况，构建一个数据库记录每只小鼠对应的各种信息
Mouse_Info <- data.frame(
    Mouse.ID = sprintf("M%02d", 1:20),
    Treatment = rep(c("Nothing", "SAL_ZT00", "SAL_ZT12", "ALC_ZT00", "ALC_ZT12"), time=4),
    Sample.Time = rep(c("ZT00", "ZT04", "ZT08", "ZT12"), each=5))
 
# 然后把小鼠的信息与基因表达信息merge到一起
Data <- merge(Mouse_Info, DAT, by="Mouse.ID")
 
# 我们先来使用boxplot来看看数据的情况
ggplot(Data) +
    geom_boxplot(aes(x=Sample.Time, y=Fc)) +
    geom_point(aes(x=Sample.Time, y=Fc, colour=Treatment)) +
    scale_y_continuous(limits=c(0, 25), expand=c(0,0), breaks=seq(0, 25, 5), label=seq(0, 25, 5))+
    ylab("Relative Gene Expression")+
    facet_grid(.~Gene.Name)+
    theme_bw()+
    theme(legend.position="top",
          axis.text=element_text(size=12),
          axis.title=element_text(size=15))

可以看到SAL_ZT12 ZT08采样的鼠有些与众不同。 由于整个实验并没有生物学重复，因此没法知道这个小鼠是否真的异常。

把这两个值去掉之后看看：

可以看出来，还是有一定的趋势的，只是由于样本量，实验过程等等各种因素，或者整个趋势就不是很明显等等的原因，并不能下什么结论。

比如，这里我们可以使用ANOVA分析看看采用时间和处理类型是否对表达有影响。

fit <- aov(Fc ~ Treatment, Data[Data$Gene.Name == "Cyp2e1", ])
summary(fit)
TukeyHSD(fit)
 
fit <- aov(Fc ~ Sample.Time, Data[Data$Gene.Name == "Cyp2e1", ])
summary(fit)
TukeyHSD(fit)
 
# Two way ANOVA
fit <- aov(Fc ~ Treatment + Sample.Time, Data[Data$Gene.Name == "Cyp2e1", ])
summary(fit)
TukeyHSD(fit)

结果看起来似乎不能说有显著的影响。

高级分析自然是多种多样的，需要根据自己具体的实验去使用相应的统计方法。 这里不展开了。

本网页已被访问396次。