R语言代码快速分析PCR结果

 

 在跑完荧光定量PCR后，我们可以得到每个样本的CT值，但是如何根据CT值来比较不同组之间基因的相对表达情况呢？6年前，我就是用excel 对着公式，分别计算每个样本内参的平均ct值，然后在依次求得每个样本与自身内参之间的ct值差（ΔCt），实验组样本与对照组之间的ct值差（ΔΔCt），再用excel公式计算2^-ΔΔCt。但近几年随着R语言不断的开发及广泛应用，使得我们对数据的处理和分析变得更简单。因此，今天我将自己写的用来分析 PCR结果的R言语代码（非常实用的干货）分享给大家。

 

首先我们需要根据自己的实验目的，样本信息及基因名称按以下格式修改自己的数据表格（excel 表格）。本示例数据中一个sample 代表一个cell line, 实验共用了5个不同的cell line,所以sample 有5个。 Group 指的是实验分组（运行代码时需更换你的对照组及分组顺序）。Replicate 指的是复孔情况，第一个复孔为1，第二个复孔为2，第三个复孔为3。Gene 是你所测的基因（运行代码时注意替换你的内参基因）。运行代码前需要先按这个格式整理表格，Ct异常值或缺失值请用NA取代。另外，如果你只有1个细胞系，实验批次不一样，也可以按照我的方法计算RNA相对表达量，即同一个批次为同一个sample。 

 

Sample	Group	Replicate	Gene	Ct
NH132	Control	1	ACTB	26.10088
NH132	TNFa	1	ACTB	23.15563
NH132	TL	1	ACTB	23.56308
NH132	TD	1	ACTB	23.48223
NH172	Control	1	ACTB	23.29007
NH172	TNFa	1	ACTB	23.82369
NH172	TL	1	ACTB	24.50135
NH172	TD	1	ACTB	23.96439
NH172	L	1	ACTB	29.32579
NH172	D	1	ACTB	24.99323
NH179	Control	1	ACTB	25.30431
NH179	TNFa	1	ACTB	22.67969
NH179	TL	1	ACTB	23.30441
NH179	TD	1	ACTB	22.24917
NH179	L	1	ACTB	26.71028
NH179	D	1	ACTB	28.84222
NH198	Control	1	ACTB	32.52913
NH198	TNFa	1	ACTB	24.14458
NH198	TL	1	ACTB	26.83226
NH198	TD	1	ACTB	24.36236
NH198	L	1	ACTB	29.07011
NH198	D	1	ACTB	26.5998
NH430	Control	1	ACTB	NA
NH430	TNFa	1	ACTB	27.82982
NH132	Control	2	ACTB	25.64198
NH132	TNFa	2	ACTB	23.09304
NH132	TL	2	ACTB	23.36409

 

 分析 PCR结果的R言语代码示例如下：

 

====================================================

 

# 加载必要的库，如果你没有安装以下这些包，请自行安装

 

library(tidyverse)

 

library(readxl)

 

library(tidyverse)

 

library(readxl)

 

library(ggpubr)

 

====================================================

 

# 读取数据（这里我将我的PCR数据结果excel 文件放在同一个R project 下面，所以可以直接读取）

 

pcr_data <- read_excel("pcr.xls")

 

colnames(pcr_data)

 

 

 

# 检查数据类型并将Ct列转换为数值类型

 

pcr_data <- pcr_data %>%

 

  mutate(Ct = as.numeric(Ct))

 

 

 

# 检查转换后的数据结构

 

str(pcr_data)

 

 

 

# 检查是否有NA（缺失值）值

 

sum(is.na(pcr_data$Ct))

 

 

 

# 过滤掉NA（缺失值）值的行

 

pcr_data <- pcr_data %>%

 

  filter(!is.na(Ct))

 

 

 

# 再次检查数据结构

 

str(pcr_data)

 

====================================================

 

# 计算平均Ct值

 

pcr_data_avg <- pcr_data %>%

 

  group_by(Group, Sample, Gene) %>%

 

  summarize(Avg_Ct = mean(Ct), .groups = 'drop')

 

 

 

# 查看计算后的数据

 

head(pcr_data_avg)

 

 

 

# 检查是否所有组和样本中都有ACTB的Ct值 （这里的ACTB是内参βactin,你们如果选择GAPDH则将下面代码中的ACTB改为GAPDH）

 

actb_check <- pcr_data_avg %>%

 

  filter(Gene == "ACTB") %>%

 

  group_by(Group, Sample) %>%

 

  summarize(count = n())

 

 

 

# 查看缺失的ACTB值

 

actb_check %>%

 

  filter(count == 0)

 

 

 

# 过滤掉缺失ACTB值的样本（这里的ACTB是内参βactin,你们如果选择GAPDH则将下面代码中的ACTB改为GAPDH）

 

pcr_data_avg <- pcr_data_avg %>%

 

  group_by(Group, Sample) %>%

 

  filter(any(Gene == "ACTB")) %>%

 

  ungroup()

 

====================================================

 

# 计算ΔCt值（这里的ACTB是内参βactin,你们如果选择GAPDH则将下面代码中的ACTB改为GAPDH）

 

pcr_data_avg <- pcr_data_avg %>%

 

  group_by(Group, Sample) %>%

 

  mutate(Delta_Ct = Avg_Ct - Avg_Ct[Gene == "ACTB"]) %>%

 

  ungroup() %>%

 

  filter(Gene != "ACTB") # 过滤掉ACTB行，因为我们只关注目标基因

 

 

 

# 计算对照组每个样本/每个基因的平均ΔCt值

 

control_mean_Delta_Ct <- pcr_data_avg %>%

 

  filter(Group == "Control") %>%

 

  group_by(Sample,Gene) %>%

 

  summarize(Mean_Delta_Ct_Control = mean(Delta_Ct), .groups = 'drop')

 

 

 

# 查看计算对照组的平均ΔCt值结果

 

head(control_mean_Delta_Ct)

 

head(pcr_data_avg)

====================================================

 

# 合并数据，确保列后缀正确

 

pcr_data_avg <- pcr_data_avg %>%

 

  inner_join(control_mean_Delta_Ct, by = c("Sample", "Gene"))

 

 

 

# 检查合并后的数据结构

 

head(pcr_data_avg)

 

 

 

# 确认列名

 

colnames(pcr_data_avg)

 

 

 

# 计算Delta_Delta_Ct值

 

pcr_data_avg <- pcr_data_avg %>%group_by(Sample,Gene) %>%

 

  mutate(Delta_Delta_Ct = Delta_Ct - Mean_Delta_Ct_Control)

 

 

 

# 计算相对表达量

 

pcr_data_avg <- pcr_data_avg %>%

 

  mutate(Relative_Expression = 2^(-Delta_Delta_Ct))

 

 

 

#确定 Group 列的顺序 （根据你自己的组进行排序，一般对照组放在第一个）

 

pcr_data_avg$Group <- factor(pcr_data_avg$Group, levels = c("Control", "TNFa", "TL", "TD", "L", "D"))

 

====================================================

 

# 绘制调整后的相对表达量图

 

p_adjusted <- ggplot(pcr_data_avg, aes(x = Gene, y = Relative_Expression, fill = Group)) +

 

  geom_bar(stat = "identity", position = "dodge") +

 

  labs(title = "Relative Expression of Target Genes",

 

       x = "Gene",

 

       y = "Relative Expression") +

 

  theme_minimal()

 

 

 

# 保存图表为 JPEG 图片

 

file_path_adjusted <- "my_plot_adjusted.jpg"

 

ggsave(file = file_path_adjusted, plot = p_adjusted, device = "jpeg", width = 8, height = 6, units = "in")

 

 

 

 

 

# 绘制箱型图

 

p_adjusted <- ggplot(pcr_data_avg, aes(x = Gene, y = Relative_Expression, fill = Group)) +

 

  geom_boxplot(outlier.shape = NA, alpha = 0.6) + # 绘制箱型图，隐藏异常值

 

  geom_jitter(position = position_jitter(width = 0.2), alpha = 0.6) + # 添加散点图

 

  labs(title = "Relative Expression of Target Genes",

 

       x = "Gene",

 

       y = "Relative Expression") +

 

  theme_minimal()

 

 

 

# 打印图形

 

print(p_adjusted)

====================================================

 

图1. 目的基因的表达量统计

 

以上是R代码分析PCR数据的全部内容。希望可以帮助到大家，如果感兴趣我们可以相互交流学习。

 

当然在数据分析前你需要检查自己的数据，并去除异常值。

 

请关注微信公众号，更多精彩内容实时更新中

                                                               

---

本网站发布所有原创内容，版权归属尖端生物及相关版权方，内容仅供学术交流，如有侵权请联系删除。未经授权的转载是侵权行为，版权方保留追究法律责任的权利。投稿，转载或版权问题，请联系:submit@advanced-biotech.cn;商务合作请联系:cc@advanced-biotech.cn