R语言分析qPCR结果
2020-09-05 本文已影响0人
焱黎
应实验室同学的需求，为了减轻他们计算结果的工作量，来编写了以下的R语言代码，利用单因素方差分析（多重比较分析）分析qPCR的Ct值，并绘制表达量柱状图。
- 方差分析用于两个及两个以上样本均数差别的显著性检验。
-
输入的文件格式
输入的文件
rm(list=ls())
options(scipen = 200) # 取消科学记数法
library(ggplot2)
library(ggsci)
library(tidyverse)
library(ggpubr)
#读取数据
rt_ct <- data.table::fread(file = "./Desktop/Ct-values.csv", data.table = F)
# 整理数据,
rt_ct_q <- rt_ct %>%
mutate(aver_ct = rt_ct[,3]-rt_ct[,2]) # 计算ΔCt ，并添加一列
rt_ct_q <- rt_ct_q %>%
mutate(quantity = c(2^(-(rt_ct_q$aver_ct[1]-rt_ct_q$aver_ct[1])), # 计算2^(-ΔΔCt)值
2^(-(rt_ct_q$aver_ct[2]-rt_ct_q$aver_ct[2])),
2^(-(rt_ct_q$aver_ct[3]-rt_ct_q$aver_ct[3])),
2^(-(rt_ct_q$aver_ct[4]-rt_ct_q$aver_ct[1])),
2^(-(rt_ct_q$aver_ct[5]-rt_ct_q$aver_ct[2])),
2^(-(rt_ct_q$aver_ct[6]-rt_ct_q$aver_ct[3])),
2^(-(rt_ct_q$aver_ct[7]-rt_ct_q$aver_ct[1])),
2^(-(rt_ct_q$aver_ct[8]-rt_ct_q$aver_ct[2])),
2^(-(rt_ct_q$aver_ct[9]-rt_ct_q$aver_ct[3
]))
))
rt_ct_q <- rt_ct_q[,c(1,5)] %>% # 只要group和quantity两列
mutate(treat = factor(rep(seq(1,3), each = 3))) %>% # 为了自定义在绘图是横坐标按照自己想要的分组顺序来绘制，因为默认按照字母顺序排列
arrange(treat) %>%
mutate(group = factor(group, levels = c("control", "treat1", "treat2")))
# 单因素方差分析-Multiple comparisons
a.aov <- aov(quantity~group, data = rt_ct_q)
library(multcomp)
multcomp_res <- glht(a.aov, linfct = mcp(group = "Tukey"))
summary(multcomp_res)
# 统计分析, 计算标准差、标准误以及95%的置信区间
library(Rmisc)
count_ana <- summarySE(rt_ct_q, measurevar = "quantity",
groupvars= "group")
count_ana
# 绘制图
ggplot(count_ana, aes(x = group, y = quantity, fill = group)) +
geom_bar(position = position_dodge(0.6),
width = 0.5, stat = "identity") +
scale_fill_manual(values = c("#000000" ,"#696969", "#A9A9A9")) +
geom_errorbar(aes(ymin = quantity - sd, ymax = quantity + sd),
position = position_dodge(0.6), width = 0.25) +
geom_signif(annotations = c("***"), y_position = c(1.1), # 显著性星号由summary(multcomp_res)来得知
xmin = c(1), xmax = c(3), # xmax的值要根据需要标记显著性的两组之间差值确定,比较control与treat1就改成2。
tip_length = c(c(0.05, 0.85)), vjust = -1) +
labs(x = NULL,
y = paste(colnames(rt_ct)[3], "mRNA expression (Related to GAPDH)"), # 利用paste函数自动添加目标基因，避免下次分析另外基因所需的手动修改
fill = "Group") +
scale_y_continuous(expand = c(0, 0),
limits = c(0, 1.2),
breaks = seq(0, 1.2, by = .2)) + # 纵坐标的刻度设置需要根据每次的数据来更改
theme_classic() +
theme(axis.title = element_text(size = 13),
axis.text.x = element_text(vjust = 0.55,
angle = 45))
-
分析结果柱状图
带显著性的柱状图
