王威
目录
1. 安装加载相应的R包
library(readxl) # readxl是用来读取excel文件的 library(ggplot2) library(reshape2) # 用来调整数据格式的 library(matrixStats) library(Hmisc)
如果提示找不到相关的包，那就用以下命令安装一下，如：
install.packages("readxl")
2. 导入数据
original_data <- read_excel("Dropbox/Your_file_path/file_name1.xlsx", sheet = "Sheet2")
dt <- original_data[,2:ncol(original_data)] row.names(dt) <- original_data$Gene # 以original_data的Gene名来命名dt的行名称
因为要做线拟合图，所以还有另外一组数据
dt_line <- read_excel("Dropbox/Your_file_path/file_name2.xlsx", sheet = "Sheet3")
3. 相关作图函数
# 计算方差，并通过T检验验证显著性，control_value是对照值
# 在计算PCR基因表达量的时候，1表示没有变化，小于1则表示表达量降低
# 大于1则表示表达量升高
get_df_tb<-function(x,cotrol_value,data){
y=rowMeans(data, na.rm=TRUE)
sd=rowSds(as.matrix(data),na.rm=TRUE)
d <- data.frame(x=x,y=y)
pvalue=rep(NA,nrow(data))
for(i in 1:nrow(data)){
a <- na.omit(as.numeric(data[i,]))
pvalue[i] <- t.test(a,rep(cotrol_value,length(a)))$p.value
}
d$signif <- pvalue < 0.05
d$heigh <- d$y + sd +0.05
d$label <- rep(NA,nrow(d))
d$p <- pvalue
for(j in 1:nrow(d)){
if(d$p[j]<0.005){
d$label[j] <- "***"
}
else if(d$p[j]<0.01){
d$label[j] <- "**"
}
else if(d$p[j]<0.05){
d$label[j] <- "*"
}
else{
d$label[j] <- NA
}
}
d
}
# 为后面作图设置一个通用的格式背景，因为ggplot做出来的图; # 默认情况下是有灰底背景，白色线条网格的，接下来三个图; # 都需要把样式调整一下，所以先写一个ggplot调整样式的公共代码. theme_all=theme( axis.title.x = element_text(margin = margin(0.5,0,0,0,'cm'),size=22), axis.title.y = element_text(margin = margin(0,0.5,0,0.5,'cm'),size=22), axis.text.x = element_text(margin = margin(0.5,0,0,0,'cm'),size=18), axis.text.y = element_text(margin = margin(0,0.2,0,0,'cm'),size=18), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(size=0.3,color="#BDBCBC",linetype = 3), panel.background = element_blank(), axis.line = element_line(size=0.3,color="black") )
# 柱状图函数，df是数据表，ggplot出图后，X轴是按照每个单词首字母单词从A~Z排列的
# 而实际上我们是想要按照我们数据展示顺序来排列，所以要先指定一下排列权重顺序
# 而这个d_levels，就是顺序排列，x_title与y_title分别是坐标轴标题
# d_theme就是主题样式，后面作图的时候，会将theme_all代入其中
get_column_plot<-function(df,d_levels,d_theme, x_title,y_title){
df$x<-factor(df$x,levels=d_levels)
ggplot(df,aes(x=x,y=y,na.rm = TRUE)) +
geom_bar(stat="identity",fill="#BDBCBC") +
geom_errorbar(data=df,mapping = aes(x=x,ymin=y-sd,ymax=y+sd),width=0.1,color="black",size=0.4) +
scale_y_continuous(labels=scales::label_comma(accuracy = 0.1))+
geom_text(data=df,mapping = aes(label=label,y=heigh),size=6) +
labs(x=x_title, y=y_title)+
d_theme
}
# 箱线图 box_plot，相关参数与上述柱状图一直
get_box_plot <- function(df,d_levels,d_theme, x_title,y_title){
df$Gene <- factor(df$Gene,levels=d_levels)
ggplot(df,aes(x=Gene,y=value,na.rm = TRUE)) + geom_boxplot() +
scale_y_continuous(labels=scales::label_comma(accuracy = 0.1))+
labs(x=x_title, y=y_title) +
# geom_text(data=data.frame(), aes(x=names(df,))) +
d_theme
}
# 线性相关图
get_line_plot <- function(df,d_theme, x_title,y_title){
# 这个lm_eqn函数是为了标注出线性方程
lm_eqn <- function(df){
m <- lm(y ~ x, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(R)^2~"="~r2,
list(a = format(unname(coef(m)[1]), digits = 2),
b = format(unname(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
ggplot(data=df,aes(x=x,y=y))+
geom_smooth(method="lm",se=FALSE,color="black",formula=y~x) +
geom_point() +
labs(x=x_title, y=y_title) +
geom_text(x = 0.7, y = 2.2, label = lm_eqn(df), parse = TRUE,size=7)+
#此处x和y是线性方程标注的位置坐标，可以根据实际情况调整
d_theme
}
4. 调用函数作图
# 1.柱状图 x<-original_dt$Gene # x轴就是基因名 data1 <- get_df_tb(x,1,dt) # dt已经在前面赋值过了 get_column_plot(df,x,theme_all,"Target Genes","Gene Expression")
# 2.箱线图，此处的x已经在上面赋值过了，不需要重新赋值了，直接调用即可 # ggplot作图，需要先把数据格式转换一下 data2 <- melt(original_dt,x=Gene) get_box_plot(data2,x,theme_all,"Target Genes","Gene Expression")
转换后的数据data2，如下图
最后做出的图如下
# 3. 线性关系图
names(df_line)<-c("x","y") #因为函数里面用的数x和y，所以需要把表格的列名称重新设置一下
get_point_plot(df_line,theme_all)
5. 存在的问题
box_plot如何标注显著性，其中一个办法是用ggsignif包来标注，但是这个包需要把对照组放上去，而且画线条将对照组与实验组相比较得出显著性，而我们其实并不需要这样，我们对照的就是1，所以只需要在每个柱子上面标注一下就可以了，那其实就是一个geom_text()函数，但是这个*标注多高呢？这个问题搜了很多地方，也没找到怎么解决，下一篇文章就更新这个问题吧
