RT-pcr:定量分析自动化工具

改版说明：下面是第一版，格式是横排基因竖排样本。还有第二版，是横排样本竖排基因。

===============第一版（）====================

第一次写一个完整功能的R代码：
R版本3.5.1

改版说明：第1.1版本的算法可能有问题，下面留着该版本的内容，仅供学习R语言练习使用。正确算法，见1.2版本。

1.1 此版本算法不正确

使用说明：更改workname后的路径到对应的本地路径
更改文件名6595-6598为对应的文件名。不要有文件格式的后缀。
在路径中放入对应的excel文件
目录结构如下：

F:/guo-gene/
6595-6598.xlsx

6595-6598.xlsx格式为

原始数据的文件格式
第一列是A01-A12B01-B12……H1-H12
第二列是对应Cq值
第三列是基因的名称顺序
第四列不需要添加内容，但需要按照我的格式保存首行内容
第五列是样品名称
对应的96孔板的布局如图所示
96孔板布局

横排一行是同一个基因，纵向是样品，三个重复依次放。

install.packages("Xlsx")
install.packages("stringr")

rm(list=ls())             #清除所有环境变量
library(stringr)
library(ggplot2)

library("xlsx")
workname <- str_c("F:/guo-gene/","6595-6598")
workshop <- str_c(workname,".xlsx")
all_data <- read.xlsx(workshop,1,encoding="UTF-8") #获取所有数据

dir.create(str_c(workname,"/"))
setwd(str_c(workname,"/"))
gene_name <- all_data[,3]
gene <- gene_name[!is.na(gene_name)] #获取基因的编号列表
cq <- all_data[,2]   #获取所有的cq值
sample_name <- all_data[,4] #获取样本的列表
sample <- sample_name[!is.na(sample_name)]  #获取该次的样本编号
sim_sample <- all_data[,5]#获取样本列表，无重复
sim <- sim_sample[!is.na(sim_sample)]

#repeatnum <- 3
#reptat_sim <- paste(sim,1:3,seq="-")
g_gene <- 0
col_gene <- list()
for (gn in 1:8){
  g_gene[[gn]]  <- gene[gn]   #g-gene[g]对应基因编号
  lf <- gn*12-11
  rg <- gn*12
  col_gene[gn] <- list(cq[lf:rg])  #新列表col_gene 赋值 
  
  
}


anctin <- as.vector(unlist(col_gene[7])) 
ubq <- as.vector(unlist(col_gene[8])) 


Cq_shiyan <- list()
Cq_jizhun <- 0
Cq_chazhi <- list()
for (k in 1:6){
  ms <- 0
  ms <- as.vector(unlist(col_gene[k]))
  Cq_shiyan[k] <- list(ms - anctin) #实验样品のCq   8个
  
  Cq_jizhun[k]<- cq[k*12-11] - anctin[1] #基准样品のCq  以第一列sample为基准 6个
  
  shiyan_cq <- as.vector(unlist(Cq_shiyan[k]))
  de_cq <- 0
  for (nei in 1:12){
    de_cq[nei] <- shiyan_cq[nei]- Cq_jizhun[k]
  }
  Cq_chazhi[k] <- list(de_cq)  #ののCq
}




#构建计算出图函数calc_expression,传参数样本名称
calc_expression <- function(sam_name,id){
  sum_single <- list()
  result <- list()
  for (num in 1:6){
    Cq_cha <- as.vector(unlist(Cq_chazhi[num])) 
    result[num] <- list(2^-(Cq_cha))
    
    single <- NULL
    for (i in 1:4){
      lfi <- 3*i-2
      rgi <- 3*i
      single[i] <- mean(result[[num]][lfi:rgi])      #3个技术重复合并计算平均值
      
    }
    
   
    sum_single[[num]] <- list(single)
   
  }
  
  #if (false){
   # cq12 <- 1:12
    #for (i in cq12) {
     # result[i] <- 2^-(Cq_chazhi[i])  #最终的相对表达量,共12个样品的Cq(三个技术重复分开计算)
    #}
  #}
    
    sam_name <- as.character(sam_name)
    img_sam_name <- str_c("g",sam_name)  #拼接字符串，制定输出的文件名为:g15798
    img_title <- str_c(sam_name," gene expression") #拼接字符串，图的标题是15798 gene expression
    y_lim <- as.vector(unlist(sum_single[id]))
    sim_resultframe <- data.frame(sim,y_lim)   #合并三个技术重复后的数据框
    print(mode(sim))
    img_sam <- ggplot(sim_resultframe,aes(x=sim,y=y_lim,label=sam_name))+labs(x="sample",y="expression")+geom_bar(stat = "identity") #合并三个技术重复的数据，生成表达量分布图
    img_sam <-  img_sam + ggtitle(img_title)+theme(plot.title = element_text(hjust=0.5)) #添加标题同时居中
    
    #end_result <- data.frame(result,sample)  #构建最终的数据框，用于画图。
    
    #ggplot(end_result,aes(x=sample,y=result))+geom_bar(stat = "identity") #画出该基因的直方图，横轴是样本编号，纵轴是相对表达量。
    
    
    img_png <- 0
    img_png <- str_c(img_sam_name,".png")
    
    #return(img_sam)
    ggsave(img_png,img_sam)
    
    
  }
  
 for (gene_id in 1:6){
  p <- calc_expression(gene[gene_id],gene_id)

}

save.image()

1.2 改版后的算法

96孔板的布局是 1.1版本中的布局。
下机数据格式整理为文件first-all.xlsx
文件first-all.xlsx的内容格式如下：

first-all格式

前几列是直接把下机的Cq粘贴过来就行。次截图未完整显示前6列内容

#此版本是用于qRT-PCR的结果的计算
#适用于横向是基因，最后一行是内参。因为我有2个内参，运算时使用最后一个内参ubq。纵向是4个样本，3个技术重复。
#first-all.xlsx是我的运算的最初的数据，参考这个格式修改你自己的数据。
#更改first-all为你的文件名，最终会生成一个总的all.csv 该csv需要自己略作调整
#其他参数，不需要调整


install.packages("Xlsx")
install.packages("stringr")

rm(list=ls())             #清除所有环境变量
library(stringr)
library(ggplot2)

library("xlsx")
workname <- str_c("F:/guo-gene/","first-all")
workshop <- str_c(workname,".xlsx")
all_data <- read.xlsx(workshop,1,encoding="UTF-8") #获取所有数据
#dir.create(str_c(workname,"/"))
setwd(str_c(workname,"/"))

gene_name <- all_data[,7]
gene <- gene_name[!is.na(gene_name)] #获取基因的编号列表
sample_name <- all_data[,8] #获取样本的列表
sample <- sample_name[!is.na(sample_name)]  #获取该次的样本编号


# 基因编号 12899  15798  39040  39301  45261  47761  anctin UBQ
#样品编号 6575 6576 6577 6578 6579 6580 6581 6582 6583 6584 6585 6586 6587 6588 6589 6590 6591 6592 6593 6594 6595 6596 6597 6598
Cq <- list()
for (i in 1:6){
  Cq[i] <- list(all_data[,i])  #Cq[i] 是第i列Cq值
}



#大循环遍历6个数据
out_put <- c(1:4)
for (big in 1:6){
  
  #实验材料的Cq-对应的内参
  
  d_cq <- list()
  for (i in 1:8){
    d_cq[i] <- list(Cq[[big]][(i*12-11):(i*12)])
  }
  anctin <- as.vector(unlist(d_cq[7]))
  ubq <- as.vector(unlist(d_cq[8]))
  
  #shiyan_cq是d_cq的每行减去对应的内参
  shiyan_cq <- list()
  for (i in 1:6){
    shiyan_cq[[i]] <- as.vector(unlist(d_cq[[i]]))-ubq
  }
  
  
  #average d_cq 合并三个技术重复  avg_dcq是合并后的平均值的list
  vol <- 0
  avg_dcq <- list()
  for (i in 1:6){
    for (r in 1:4){
      vol[r] <- mean(shiyan_cq[[i]][(3*r-2):(3*r)])
    }
    avg_dcq[i] <- list(vol)
  }
  
  
  #dd_cq是所有的样品的d_cq-avg_cq[[1]][1]  设置左侧第一个样品第一个基因为基准，其他的都减去这个样本
  dd_cq <- list()
  ddcq <- 0
  for (i in 1:6){
    for (k in 1:12){
      ddcq[k] <- shiyan_cq[[i]][k]- avg_dcq[[1]][1]
    }
    dd_cq[i] <- list(ddcq)
  }
  
  
  #2^-dd_cq 
  dd2_cq <- list()
  dd2cq <- 0
  for (i in 1:6) {
    for (k in 1:12){
      dd2cq[k] <- 2^(-dd_cq[[i]][k])
    }
    dd2_cq[i] <- list(dd2cq)
  }
  
  
  #avg_dd2_cq  2^-dd_cq 三个重复的平均值
  avg_dd2_cq <- list()
  avgdd2cq <- 0
  for (i in 1:6){
    for (k in 1:4) {
      avgdd2cq[k] <- mean(dd2_cq[[i]][(k*3-2):(k*3)])
    }
    avg_dd2_cq[i] <- list(avgdd2cq)
  }
  avg_dd2_cq[[7]] <- list(sample[(big*4-3):(big*4)])
  #输出数据框
  out_put <- data.frame(out_put,avg_dd2_cq)
  #filename <- str_c(big,".csv")
  #write.csv(avg_dd2_cq,file=filename)   #此处是分批次输出对应的结果

}
write.csv(out_put,file = 'all.csv')   #输出所有结果到一个csv

===================第2版==========================

此版本是横排是样本名称，竖排是基因编号。内参必须在最后。
最终会新建你的文件的同名的文件夹，并在该文件夹输出同名的csv文件。
此版本包括两种：

2.1 第一种是有2个基因+1个内参

96孔板布局如下图：

版本2.1

我的2-4的文件格式是：

数据格式

使用前修改

修改2处 2-4改为你对应的文件名，最后一处是你输出的文件名为2-4.csv

修改替换39589,14341为对应的按照顺序的基因编号，ubq改为对应的内参名

修改你的下机数据格式为我的2-4的参考格式

实现计算代码：

#使用说明文件
#本R程序是自动化分析qRT-pcr的结果的
#针对的是2个基因，3个技术重复，8个样品。
#如果你的样品数目比8大，请仔细修改程序中的8出现的地方。小于8则无所谓
#数据格式

#使用前修改
#修改2处  2-4改为你对应的文件名，最后一处是你输出的文件名为2-4.csv
#修改替换39589,14341为对应的按照顺序的基因编号，ubq改为对应的内参名





install.packages("xlsx")
install.packages("stringr")

rm(list=ls())             #清除所有环境变量
library(stringr)
library(ggplot2)

library("xlsx")
workname <- str_c("F:/guo-gene/","2-2")
workshop <- str_c(workname,".xlsx")
all_data <- read.xlsx(workshop,1,encoding="UTF-8") #获取所有数据

dir.create(str_c(workname,"/"))
setwd(str_c(workname,"/"))
gene_name <- all_data[,2]
gene <- gene_name[!is.na(gene_name)] #获取基因的编号列表
cq <- all_data[,1]   #获取所有的cq值
cq <- cq[!is.na(cq)]
sample_name <- all_data[,3] #获取样本的列表
sample <- sample_name[!is.na(sample_name)]  #获取该次的样本编号
sim_sample <- all_data[,3]#获取样本列表，无重复
sim <- sim_sample[!is.na(sim_sample)]
gene_list_len <- length(gene)  #获取基因的列表长度，用以确定运算位置。最后一个样品是内参
#repeatnum <- 3
#reptat_sim <- paste(sim,1:3,seq="-")


sample_gene <- list()
col_gene <- list()

ubq <- list()
g_14341 <- list()
g_39589 <- list()
#初始化list的每一个向量
for (seq in 1:3){
  g_39589[[seq]] <- 1
  g_14341[[seq]] <- 1
  ubq[[seq]] <- 1
}

#获取每一个基因的对应材料的Cq值
for (gn in 1:8){
  pos <- c(0,gene_list_len*3,gene_list_len*6,gene_list_len*9,gene_list_len*12,gene_list_len*15,gene_list_len*18,gene_list_len*21)
  for (sa1 in 1:3){
    print(cq[pos[gn]+sa1])
    
    g_39589[[sa1]][gn] <- cq[pos[gn]+sa1]  #g_39589[1]代表39589基因的第一个技术重复的8个Cq值
    g_14341[[sa1]][gn] <- cq[pos[gn]+sa1+3]
    ubq[[sa1]][gn] <- cq[pos[gn]+sa1+6]
  }
}


Cq_shiyan_39589 <- list()
Cq_shiyan_14341 <- list()
#の实验Cq  Cq_shiyan
for (i in 1:3){
 Cq_shiyan_39589[i] <- list(g_39589[[i]] - ubq[[i]])
 Cq_shiyan_14341[i] <- list(g_14341[[i]]-ubq[[i]])
}


#average Cq
avg_39589 <- (Cq_shiyan_39589[[1]]+Cq_shiyan_39589[[2]]+Cq_shiyan_39589[[3]])/3
avg_14341 <- (Cq_shiyan_14341[[1]]+Cq_shiyan_14341[[2]]+Cq_shiyan_14341[[3]])/3

#ののCq
dd_39589 <- list()
dd_14341 <- list()
for (i in 1:3){
  dd_39589[[i]] <- 1
  dd_14341[[i]] <- 1
}
for (k in 1:8){
  for (i in 1:3) {
    dd_39589[[i]][k] <- Cq_shiyan_39589[[i]][k] - avg_39589[k]
    dd_14341[[i]][k] <- Cq_shiyan_14341[[i]][k] - avg_14341[k]
  }
}

#2^-ののCq
dd2_39589 <- list()
dd2_14341 <- list()
for (i in 1:3){
  dd2_39589[[i]] <- 1
  dd2_14341[[i]] <- 1
}
for (k in 1:8){
  for (i in 1:3) {
    dd2_39589[[i]][k] <- 2^(avg_39589[1]-Cq_shiyan_39589[[i]][k])
    dd2_14341[[i]][k] <- 2^(avg_14341[1]-Cq_shiyan_14341[[i]][k])
  }
}


#2^-ののCq 的平均值
avg2_39589 <- (dd2_39589[[1]]+dd2_39589[[2]]+dd2_39589[[3]])/3
avg2_14341 <- (dd2_14341[[1]]+dd2_14341[[2]]+dd2_14341[[3]])/3



frame_39589 <- data.frame(sim,avg2_39589)
frame_14341 <- data.frame(sim,avg2_14341)


t <- merge(frame_39589,frame_14341,all = TRUE,sort = TRUE)
write.csv(t,file='2-2.csv')

2.2 第二种是3个基因+1个内参

96孔板布局如下图：

3gene+1内参

我的3-1文件的下机数据的格式：

数据格式

使用前修改3-1为你对应的文件名，不包括后缀。

修改替换42886，18342，12391，ubq为你实际位置对应的基因编号。内参必须在最后。

文件格式为普通xlsx格式。具体数据位置见下图。

实现代码：

#使用说明文件
#本R程序是自动化分析qRT-pcr的结果的
#针对的是3个基因，3个技术重复，8个样品。
#如果你的样品数目比8大，请仔细修改程序中的8出现的地方。小于8则无所谓


#使用前修改3-1为你对应的文件名，不包括后缀。
#修改替换42886，18342，12391，ubq为你实际位置对应的基因编号。内参必须在最后。
#文件格式为普通xlsx格式。具体数据位置见下图。


#install.packages("Xlsx")
#install.packages("stringr")

rm(list=ls())             #清除所有环境变量
library(stringr)
library(ggplot2)

library("xlsx")
workname <- str_c("F:/guo-gene/","3-1")
workshop <- str_c(workname,".xlsx")
all_data <- read.xlsx(workshop,1,encoding="UTF-8") #获取所有数据

dir.create(str_c(workname,"/"))
setwd(str_c(workname,"/"))
gene_name <- all_data[,2]
gene <- gene_name[!is.na(gene_name)] #获取基因的编号列表
cq <- all_data[,1]   #获取所有的cq值
cq <- cq[!is.na(cq)]
sim_sample <- all_data[,3]#获取样本列表，无重复
sim <- sim_sample[!is.na(sim_sample)]
gene_list_len <- length(gene)  #获取基因的列表长度，用以确定运算位置。最后一个样品是内参
#repeatnum <- 3
#reptat_sim <- paste(sim,1:3,seq="-")


g_42886 <- list()
g_18342 <- list()
g_12391 <- list()
ubq <- list()

#初始化list的每一个向量
for (seq in 1:3){
  g_42886[[seq]] <- 1
  g_18342[[seq]] <- 1
  g_12391[[seq]] <- 1
  ubq[[seq]] <- 1
}

#获取每一个基因的对应材料的Cq值
for (gn in 1:8){
  pos <- c(0,gene_list_len*3,gene_list_len*6,gene_list_len*9,gene_list_len*12,gene_list_len*15,gene_list_len*18,gene_list_len*21)
  for (sa1 in 1:3){
    print(cq[pos[gn]+sa1])
    
    g_42886[[sa1]][gn] <- cq[pos[gn]+sa1]  #g_42886[1]代表42886基因的第一个技术重复的8个Cq值
    g_18342[[sa1]][gn] <- cq[pos[gn]+sa1+3]
    g_12391[[sa1]][gn] <- cq[pos[gn]+sa1+6]
    ubq[[sa1]][gn] <- cq[pos[gn]+sa1+9]
  }
}


Cq_shiyan_42886 <- list()
Cq_shiyan_18342 <- list()
Cq_shiyan_12391 <- list()
#の实验Cq  Cq_shiyan
for (i in 1:3){
  Cq_shiyan_42886[i] <- list(g_42886[[i]] - ubq[[i]])
  Cq_shiyan_18342[i] <- list(g_18342[[i]]-ubq[[i]])
  Cq_shiyan_12391[i] <- list(g_12391[[i]]-ubq[[i]])
}


#average Cq
avg_42886 <- (Cq_shiyan_42886[[1]]+Cq_shiyan_42886[[2]]+Cq_shiyan_42886[[3]])/3
avg_18342 <- (Cq_shiyan_18342[[1]]+Cq_shiyan_18342[[2]]+Cq_shiyan_18342[[3]])/3
avg_12391 <- (Cq_shiyan_12391[[1]]+Cq_shiyan_12391[[2]]+Cq_shiyan_12391[[3]])/3
#ののCq
dd_42886 <- list()
dd_18342 <- list()
dd_12391 <- list()
for (i in 1:3){
  dd_42886[[i]] <- 1
  dd_18342[[i]] <- 1
  dd_12391[[i]] <- 1
}
for (k in 1:8){
  for (i in 1:3) {
    dd_42886[[i]][k] <- Cq_shiyan_42886[[i]][k] - avg_42886[k]
    dd_18342[[i]][k] <- Cq_shiyan_18342[[i]][k] - avg_18342[k]
    dd_12391[[i]][k] <- Cq_shiyan_12391[[i]][k] - avg_12391[k]
  }
}

#2^-ののCq
dd2_42886 <- list()
dd2_18342 <- list()
dd2_12391 <- list()
for (i in 1:3){
  dd2_42886[[i]] <- 1
  dd2_18342[[i]] <- 1
  dd2_12391[[i]] <- 1
}
for (k in 1:8){
  for (i in 1:3) {
    dd2_42886[[i]][k] <- 2^(avg_42886[1]-Cq_shiyan_42886[[i]][k])
    dd2_18342[[i]][k] <- 2^(avg_18342[1]-Cq_shiyan_18342[[i]][k])
    dd2_12391[[i]][k] <- 2^(avg_12391[1]-Cq_shiyan_12391[[i]][k])
  }
}


#2^-ののCq 的平均值
avg2_42886 <- (dd2_42886[[1]]+dd2_42886[[2]]+dd2_42886[[3]])/3
avg2_18342 <- (dd2_18342[[1]]+dd2_18342[[2]]+dd2_18342[[3]])/3
avg2_12391 <- (dd2_12391[[1]]+dd2_12391[[2]]+dd2_12391[[3]])/3


frame_42886 <- data.frame(sim,avg2_42886)
frame_18342 <- data.frame(sim,avg2_18342)
frame_12391 <- data.frame(sim,avg2_12391)

t <- merge(frame_42886,frame_18342)
y <- merge(t,frame_12391)
write.csv(y,file='3-1.csv')