迈出脚步
体验世界

R语言(R language)分析qPCR数据

目录

  1. 安装并加载相应的R包
  2. 导入数据
  3. 相关作图函数
  4. 调用函数作图
  5. 存在的问题

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")

original_data

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")

column_img

 

# 2.箱线图,此处的x已经在上面赋值过了,不需要重新赋值了,直接调用即可
# ggplot作图,需要先把数据格式转换一下
data2 <- melt(original_dt,x=Gene)
get_box_plot(data2,x,theme_all,"Target Genes","Gene Expression")

转换后的数据data2,如下图

data2

最后做出的图如下

box_img

# 3. 线性关系图
names(df_line)<-c("x","y") #因为函数里面用的数x和y,所以需要把表格的列名称重新设置一下
get_point_plot(df_line,theme_all)

liner_img

 

5. 存在的问题

box_plot如何标注显著性,其中一个办法是用ggsignif包来标注,但是这个包需要把对照组放上去,而且画线条将对照组与实验组相比较得出显著性,而我们其实并不需要这样,我们对照的就是1,所以只需要在每个柱子上面标注一下就可以了,那其实就是一个geom_text()函数,但是这个*标注多高呢?这个问题搜了很多地方,也没找到怎么解决,下一篇文章就更新这个问题吧

http://xzh.i3geek.com
赞(0) 喜欢我
未经允许不得转载:王威 » R语言(R language)分析qPCR数据
分享到: 更多 (0)

评论 抢沙发

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

微信扫一扫打赏