qPCR结果绘图-详细注释

# 清除先前数据并设置工作目录
rm(list=ls())
setwd("D:/great_work/R_wd/qPCR")

# 读取数据
# txt文件
# 第一列为“group”
# 第二列为内参Ct（列名为内参基因名）
# 第三列为基因Ct（靶基因名）
# “Tab”分隔不同列
data <- read.table("qPCR-LMF.txt",head=T,sep = "\t") #经典读取txt，以第一行为列名，tab为分隔符取出数据

#基础计算
library(magrittr)
library(dplyr)
#计算ΔCt（dCt），并添加一列
#transform() 函数的作用是为数据框增加新的列或修改已有列的值
data2 <- transform(data, dCt = data[,3] - data[,2])
#计算对照组dCt均值，N取决于第一组的数据量
#如果“对照组的”ΔCt间变化较大，建议用几何平均值，更好地处理异常值
#用于计算ddCt
M <- mean(data2$dCt[1:6])
#计算
data3 <- transform(data2, quantity = 2^(-(data2$dCt-M)))#计算绘图数据ddCt

#设置绘图顺序
data4 <- data3[,c(1,5)] %>% # 只要group和ddCt两列
  mutate(treat = factor(rep(c("Leaf", "M", "F"), c(6, 3, 6)))) %>%# 自定义绘图时横坐标顺序
  arrange(treat) %>%
  mutate(group = factor(group, levels = c("Leaf", "M", "F")))

#拟合模型，计算P值
#进行单因素方差分析（ANOVA）并使用 Tukey's HSD 方法进行多重比较
a.aov <- aov(quantity~group, data = data4)#aov函数用于计算 ddCt 在 group 变量的不同水平之间的方差，生成一个 aov 对象。
library(multcomp)
multcomp_res <- glht(a.aov, linfct = mcp(group = "Tukey"))#glht函数用于对 aov 对象进行多重比较，并计算出每个组之间的显著性差异。
multcomp_summary <- summary(multcomp_res)#summary函数用于汇总多重比较的结果，包括每个组之间的比较结果、置信区间、p 值等
multcomp_summary

#从多重比较结果中提取P值和比较组并调整为输入格式
#提取比较组
comparisons <- names(multcomp_summary$test$coefficients)
comparisons <- lapply(comparisons, function(x) strsplit(x, " - ")[[1]])
comparisons
#提取P值
pvalues <- as.numeric(multcomp_summary$test$pvalues)
p.format <- function(x) {
  ifelse(x < 0.001, "***", 
         ifelse(x < 0.01, "**", 
                ifelse(x < 0.05, "*", "ns")))
}
star_labels <- c(p.format(pvalues))
star_labels

#计算sd等统计学特征
library(Rmisc)
count_ana <- summarySE(data4, measurevar = "quantity", #可以考虑psych包的describeBy() 函数
                       groupvars = "group")
count_ana <- transform(count_ana, sem = sd/N)#增加sem列，适用小样本绘图
count_ana

#设置颜色方案
library(RColorBrewer)
colors <- brewer.pal(3, "Set2")#group的数量对应选取方案中颜色的数量

#绘图函数
library(ggplot2)
library(ggsignif)
p <- ggplot(count_ana, aes(x = group, y = quantity, fill = group)) +
     geom_bar(position = position_dodge(0.1), 
           width = 0.5, stat = "identity") +
  
     scale_fill_manual(values = colors) +
  
     geom_errorbar(aes(ymin = quantity - sem, ymax = quantity + sem),#sem适用小样本，sd适用大样本
                position = position_dodge(0.6), width = 0.25) +

     labs(x = NULL,#x轴可以设置标题
       y = paste(colnames(data)[3], "mRNA expression (Related to GAPC1)"), # 利用paste函数自动添加目标基因
       fill = "Group") + 
  
     scale_y_continuous(expand = c(0, 0),# 纵坐标的刻度设置需要数据范围修改
                     limits = c(0, 20),
                     breaks = seq(0, 20, by = 2)) + 
  
     theme_classic() +
  
     theme(axis.title = element_text(size = 13),#x轴标签的大小
        axis.text.x = element_text(vjust = 0.55,#x轴标签距离坐标轴的距离
                                   angle = 0))#x轴标签旋转的角度

for (i in 1:length(star_labels)) {
  p <- p + 
    geom_signif(comparisons = comparisons[i],
                annotations = star_labels[i],
                y_position = max(count_ana$quantity) + i*1.5 + 1.5)
}

p