介紹

做完轉(zhuǎn)錄組分析之后货矮，一般都要求做qRT-PCR來驗證二代測序得到的轉(zhuǎn)錄本表達(dá)是否可靠佃乘。熒光定量PCR是一種相對表達(dá)定量的方法薪贫，他的計算方法有很多艘狭，常用的相對定量數(shù)據(jù)分析方法有雙標(biāo)曲線法，ΔCt法囚霸，2^{-ΔΔCt法(Livak法)腰根，用參照基因的2}-ΔΔCt法(Livak法):

qRT-PCR介紹及計算公式

該部分引用自下方參考鏈接1

qRT-PCR原理

以基因的cDNA為模板進(jìn)行PCR擴(kuò)增，在PCR擴(kuò)增過程中拓型，通過收集熒光信號额嘿，對PCR進(jìn)程進(jìn)行實時檢測。由于在PCR擴(kuò)增的指數(shù)時期劣挫，模板的Ct值和該模板的起始拷貝數(shù)存在線性關(guān)系册养，所以可以定量。

Ct值

計算 -ΔΔCt：內(nèi)參基因分為對照組和處理組內(nèi)參基因

1. 先計算對照組和處理組的內(nèi)參基因Ct的均值：

2. 計算對照組待檢測目的基因減去對照組內(nèi)參基因的平均Ct值：

3. 計算處理組待檢測目的基因減去處理組內(nèi)參基因的平均Ct值：

4. 計算基于對照組的-ΔΔCt捡絮，處理組待檢測目的基因的ΔCt減去對照組待檢測基因的ΔCt的平均值：

5. 相對表達(dá)量計算，也就是相對于對照組: 2^-ΔΔct:

6. 條形圖或相關(guān)性點圖可視化結(jié)果

R代碼

加載R包

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

R函數(shù)

get_qPCR <- function(dataset=dat,
                     ref_gene="GAPDH",
                     control_group="6H NC",
                     grp=c("6H M1")){
  
  # dataset=dat                   # 初始數(shù)據(jù)
  # ref_gene="GAPDH"              # 內(nèi)參基因名字
  # 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: 計算對照組和處理組的內(nèi)參基因平均值
  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: 計算對照組和處理組待檢測目的基因減去對應(yīng)分組的內(nèi)參基因的平均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: 條形圖或相關(guān)性散點圖可視化
  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)  
}

讀取數(shù)據(jù)

單個樣本三個技術(shù)重復(fù)福稳，檢驗不同的目的基因擴(kuò)增效率

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

計算結(jié)果

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

可視化結(jié)果

qPCR_res$plot

結(jié)果： IL-1B 和INOS基因相比NC組而言，其含量越多

同一基因多分組結(jié)果圖

get_qPCR <- function(dataset=dat,
                     ref_gene="Gadph",
                     control_group="2d_Control",
                     grp=c("2d_100uM")){
  
  # dataset=dat                   # 初始數(shù)據(jù)
  # ref_gene="Gadph"              # 內(nèi)參基因名字
  # control_group="2d_Control"    # 對照組
  # grp=c("2d_100uM")             # 實驗組排序
  
  
  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: 計算對照組和處理組的內(nèi)參基因平均值
  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: 計算對照組和處理組待檢測目的基因減去對應(yīng)分組的內(nèi)參基因的平均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: 條形圖或相關(guān)性散點圖可視化
  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))
  
  # step6: visualization
  pl <- ggplot(dat_plot_bar, aes(x=Sample_Name, weight=qPCR))+
    geom_hline(yintercept = seq(0, round(max(dat_plot_bar$qPCR), 1), 0.2), color = "gray")+
    geom_hline(yintercept = 1, color = "black", linetype = 2, size = 1)+ 
    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")))+    
    scale_y_continuous(breaks = seq(0, round(max(dat_plot_bar$qPCR), 1), 0.2),
        expand = c(0, 0),
        limits = c(0, round(max(dat_plot_bar$qPCR), 1)+round(max(dat_plot_bar$sd), 1)))+
    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)  
}

dat <- read.xlsx("qPCR.xlsx", sheetName = "6d")
# 6 days 
qPCR_res <- get_qPCR(dataset=dat,
                     ref_gene="Gapdh",
                     control_group="6d_control",
                     grp=c("6d_100", "6d_300", "6d_1000"))

# 計算結(jié)果
DT::datatable(qPCR_res$dat)

# 可視化結(jié)果
qPCR_res$plot

參考

1. qRT-PCR相對定量計算詳解

2. geom_lines in different facet