應(yīng)實驗室同學(xué)的需求，為了減輕他們計算結(jié)果的工作量，來編寫了以下的R語言代碼，利用單因素方差分析（多重比較分析）分析qPCR的Ct值，并繪制表達(dá)量柱狀圖。
- 方差分析用于兩個及兩個以上樣本均數(shù)差別的顯著性檢驗。
-
輸入的文件格式輸入的文件
rm(list=ls())
options(scipen = 200) # 取消科學(xué)記數(shù)法
library(ggplot2)
library(ggsci)
library(tidyverse)
library(ggpubr)
#讀取數(shù)據(jù)
rt_ct <- data.table::fread(file = "./Desktop/Ct-values.csv", data.table = F)
# 整理數(shù)據(jù),
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))) %>% # 為了自定義在繪圖是橫坐標(biāo)按照自己想要的分組順序來繪制，因為默認(rèn)按照字母順序排列
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)
# 統(tǒng)計分析, 計算標(biāo)準(zhǔn)差、標(biāo)準(zhǔn)誤以及95%的置信區(qū)間
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的值要根據(jù)需要標(biāo)記顯著性的兩組之間差值確定,比較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函數(shù)自動添加目標(biāo)基因，避免下次分析另外基因所需的手動修改
fill = "Group") +
scale_y_continuous(expand = c(0, 0),
limits = c(0, 1.2),
breaks = seq(0, 1.2, by = .2)) + # 縱坐標(biāo)的刻度設(shè)置需要根據(jù)每次的數(shù)據(jù)來更改
theme_classic() +
theme(axis.title = element_text(size = 13),
axis.text.x = element_text(vjust = 0.55,
angle = 45))
-
分析結(jié)果柱狀圖帶顯著性的柱狀圖
