数据分析：RT-qPCR分析及R语言绘图

介绍

转录组分析是一种用于研究细胞或组织中所有RNA分子的表达水平的高通量技术。完成转录组分析后，科学家们通常需要通过定量实时聚合酶链式反应（qRT-PCR）来验证二代测序（Next-Generation Sequencing, NGS）结果的可靠性。这是因为qRT-PCR是一种精确的定量方法，可以用来验证特定基因的表达水平。

荧光定量PCR（Quantitative Real-Time PCR）是一种利用荧光信号来实时监测PCR扩增过程的技术。它允许研究者在PCR反应进行时实时检测DNA的累积量，从而实现对基因表达水平的定量分析。

在进行相对定量分析时，常用的方法之一是双标曲线法（也称为标准曲线法或绝对定量法）。这种方法的基本步骤如下：

1. 标准曲线的构建：首先，需要通过一系列已知浓度的标准品（通常是目标基因的克隆DNA）进行PCR扩增，以获得一系列的Ct值（阈值循环数，即PCR扩增过程中荧光信号首次超过阈值的循环次数）。然后，将这些Ct值对数转换后与相应的DNA浓度绘制成图，形成标准曲线。
2. 样本的Ct值测定：接下来，对实验样本进行qRT-PCR，记录目标基因的Ct值。
3. 相对定量计算：利用标准曲线，根据样本的Ct值计算出样本中目标基因的相对浓度。这通常涉及到将样本的Ct值转换为DNA浓度，然后与标准品的浓度进行比较。
4. 数据归一化：由于qRT-PCR可能会受到实验操作和样本制备的影响，因此需要使用一个或多个内参基因（通常是表达水平相对稳定的基因）来归一化数据，以消除这些潜在的变异。
5. 结果表达：最终，研究者会以目标基因相对于内参基因的表达水平来表达结果，通常是以2的幂次方来表示倍数变化。

通过这种方法，研究者可以验证NGS结果的准确性，并进一步探索基因表达的调控机制。ΔCt法，2^{-ΔΔCt法(Livak法)，用参照基因的2}-ΔΔCt法(Livak法):

qRT-PCR介绍及计算公式

该部分引用自下方参考链接1

qRT-PCR原理

以基因的cDNA为模板进行PCR扩增，在PCR扩增过程中，通过收集荧光信号，对PCR进程进行实时检测。由于在PCR扩增的指数时期，模板的Ct值和该模板的起始拷贝数存在线性关系，所以可以定量。

Ct值

Ct值的含义是：每个反应管内的荧光信号达到设定的域值时所经历的循环数 (cycle)。 qRT-PCR在扩增的时候都会有平台期，在平台期之前，PCR 扩增就是简单的指数增长，也就是 1 变 2，2 变 4，4 变 8 …扩增。数学形式就是 2 的 ct 次方，到了平台期所有基因扩增的数目是一致的，而唯一有区别的则是 ct 值的不同。所以不难推断出 ct 值越小，反应扩增到达平台期所需循环数越少，目的基因起始含量越高。这里可以得到公式：

计算 -ΔΔCt：内参基因分为对照组和处理组内参基因

1. 先计算对照组和处理组的内参基因Ct的均值： $$Mea{n}_{内参基因}=mean(对照组或处理组内参基因)$$

2. 计算对照组待检测目的基因减去对照组内参基因的平均Ct值：$$\Delta C{t}_{对照组目的基因i}=C{t}_{对照组目的基因i}-C{t}_{对照组内参基因的平均值}$$

3. 计算处理组待检测目的基因减去处理组内参基因的平均Ct值：$$\Delta C{t}_{处理组目的基因i}=C{t}_{处理组目的基因i}-C{t}_{处理组内参基因的平均值}$$

4. 计算基于对照组的-ΔΔCt，处理组待检测目的基因的ΔCt减去对照组待检测基因的ΔCt的平均值：$$-\Delta \Delta C{t}_{处理组目的基因i}=\Delta C{t}_{处理组目的基因i}-\Delta C{t}_{对照组目的基因i的平均值}$$

5. 相对表达量计算，也就是相对于对照组: 2^-ΔΔct: $$2^{-(-\Delta \Delta Ct)}$$

6. 条形图或相关性点图可视化结果

R代码

加载R包

knitr::opts_chunk$set(warning = F, message = F)
library(dplyr)
library(tibble)
library(ggplot2)
library(xlsx)
library(Rmisc)

R函数

get_qPCR <- function(dataset=dat,
                     ref_gene="GAPDH",
                     control_group="6H NC",
                     grp=c("6H M1")){
  
  # dataset=dat                   # 初始数据
  # ref_gene="GAPDH"              # 内参基因名字
  # control_group="6H NC"         # 对照组
  # grp=c("6H M1")                # 实验组排序
  
  
  if(!any(is.element(colnames(dataset), c("Sample_Name", "Target_Name", "CT")))){
    stop("Check the sheet's colnames")
  }
  sampleid <- c("Sample_Name", "Target_Name", "CT")
  dat <- dataset %>% select(sampleid)
  
  # step1: 计算对照组和处理组的内参基因平均值
  dat_ref_gene <- dat %>% filter(Target_Name == ref_gene) 
  ref_gene_mean <- dat_ref_gene %>% group_by(Sample_Name) %>%
    dplyr::summarise(CT_ref_mean = mean(CT))
  
  # step2: 计算对照组和处理组待检测目的基因减去对应分组的内参基因的平均Ct值
  dat_gene <- dat %>% filter(Target_Name != ref_gene) 
  dat_gene_merge <- dat_gene %>% inner_join(ref_gene_mean, by = "Sample_Name")
  dat_gene_merge$CT_delta <- with(dat_gene_merge, CT - CT_ref_mean) 
  
  dat_control <- dat_gene_merge %>% filter(Sample_Name == control_group) %>%
    group_by(Sample_Name, Target_Name) %>%
    dplyr::summarise(Delta_CT_control_mean=mean(CT_delta)) %>% 
    dplyr::rename(Sample_Name_control=Sample_Name)
  dat_treat <- dat_gene_merge %>% filter(Sample_Name != control_group) %>%
    # group_by(Sample_Name, Target_Name) %>%
    # dplyr::summarise(Delta_CT_treat_mean=mean(CT_delta)) %>% 
    dplyr::rename(Sample_Name_treat=Sample_Name)
  
  # step3: 计算对照组检测基因的平均Δ值
  dat_double_delta <- inner_join(dat_treat, dat_control,
                                 by = "Target_Name")
  dat_double_delta$CT_delta_delta <- with(dat_double_delta, CT_delta - Delta_CT_control_mean)
  
  # step4: 基于对照组检测基因的平均Δ值，计算实验组的2-ΔΔCt值
  dat_double_delta$qPCR <- 2^-(dat_double_delta$CT_delta_delta) 
  
  # step5: 条形图或相关性散点图可视化
  dat_plot <- dat_double_delta %>% 
    dplyr::rename(Sample_Name=Sample_Name_treat) %>%
    dplyr::select(Sample_Name, Target_Name, qPCR) 
  dat_plot_bar <- Rmisc::summarySE(dat_plot, measurevar = "qPCR", 
                                   groupvars = c("Sample_Name", "Target_Name")) %>%
    mutate(Sample_Name=factor(Sample_Name, levels = grp),
           Target_Name=factor(Target_Name)) %>% 
    group_by(Sample_Name, Target_Name) %>%
    mutate(ylimit=(qPCR+sd)) %>%
    ungroup()
  
  dat_plot_bar_ymax <- dat_plot_bar %>% 
    group_by(Target_Name) %>% 
    summarise_at(vars(ylimit), max)
  
  # dat_plot_range <- dat_plot %>% group_by(Sample_Name, Target_Name) %>%
  #   summarise(ymin=min(qPCR), ymax=max(qPCR))
  # setting y axis scale
  y_group <- c()
  y_scale <- c()
  for(i in 1:nrow(dat_plot_bar_ymax)){
    y_group <- c(y_group, rep(as.character(dat_plot_bar_ymax$Target_Name[i]), 2))
    y_scale <- c(y_scale, c(0, ceiling(dat_plot_bar_ymax$ylimit[i])))
  }
  blank_data <- data.frame(Target_Name = y_group, 
                           Sample_Name = 1, 
                           qPCR = y_scale)
  
  # step6: visualization
  pl <- ggplot(dat_plot_bar, aes(x=Sample_Name, weight=qPCR))+
    geom_hline(aes(yintercept = qPCR), color = "gray")+
    geom_bar(color = "black", width = .4, position = "dodge")+
    geom_errorbar(aes(ymin = qPCR, ymax = qPCR + se), 
                  width = 0.25, size = 0.5, position = position_dodge(0.7))+
    labs(x="", y=expression(paste(log[2], " fold change in expression")))+ 
    geom_blank(data = blank_data, aes(x = Sample_Name, y = qPCR))+
    expand_limits(y = 0)+
    scale_y_continuous(expand = c(0, 0))+
    facet_wrap(. ~ Target_Name, scales = "free")+
    theme_bw()+
    theme(axis.title = element_text(face = "bold", color = "black", size = 14),
          axis.text = element_text(color = "black", size = 10),
          axis.text.x = element_text(angle = 60, hjust = 1, face = "bold"),
          text = element_text(size = 10, color = "black", family="serif"),
          panel.grid = element_blank(),
          legend.position = "right",
          legend.key.height = unit(0.6, "cm"),
          legend.text = element_text(face = "bold", color = "black", size = 10),
          strip.text = element_text(face = "bold", size = 14))
  res <- list(dat=dat_double_delta, plot=pl)
  return(res)  
}

读取数据

单个样本三个技术重复，检验不同的目的基因扩增效率。数据下载链接：

百度云盘链接: https://pan.baidu.com/s/1W4Uvoy3Z7i9s8dWGOOkKWw
提取码: vh4s

dat <- read.xlsx("qPCR.xlsx", sheetIndex = 1)
head(dat)

计算结果

qPCR_res <- get_qPCR(dataset=dat,
                     ref_gene="GAPDH",
                     control_group="6H NC",
                     grp=c("6H M1"))
DT::datatable(qPCR_res$dat)

可视化结果

qPCR_res$plot

结果： IL-1B 和INOS基因相比NC组而言，其含量越多

参考

1. qRT-PCR相对定量计算详解

2. geom_lines in different facet