一个处理qPCR数据的程序(未debug版本,慎用!)

2020-06-13  本文已影响0人  一只烟酒僧
######################################################## 
#-------------------------------------------------------
# Topic:进行qRTPCR运算
# Author:Wang Haiquan
# Date:Sat Jun 13 21:28:16 2020
# Mail:mg1835020@smail.nju.edu.cn
#-------------------------------------------------------
########################################################
#计算过程
#1、将复孔取平均值,同时顺带把差值大于0.5的复孔删除
#2、样本靶基因ct-样本内参ct,求出Δct,只支持一个内参
#3、样本Δct-对照Δct,可设置多个对照,也可不设置对照
#4、指数运算
########################################################
#模拟数据
#列是基因重复,第一个基因为对照
#行是样本名
mydata<-as.data.frame(matrix(abs(rnorm(384,0,1)),nrow = 32))
mydata
colnames(mydata)<-rep(paste("gene",1:4,sep = ""),each=3)
rownames(mydata)<-paste("sample",1:32,sep = "")
##############################################################################
#输入
mydata=mydata #按照列为gene+复孔重复,行为样本名,必须使用三个复孔重复
neican_name<-"gene1" #作为内参的基因名,仅支持设置一个内参基因
control="sample1"  #作为对照的样本名,可设置多个对照,也可不设置对照,只计算deltact
##############################################################################
#从粘贴板上读取数据
mydata<-read.table("clipboard",sep = "\t",header = T,row.names = 1)
ct<-as.data.frame(matrix(NA,ncol=nrow(mydata),nrow=ncol(mydata)/3))
#将复孔取平均值,同时顺带把差值大于0.5的复孔删除
for (i in 1:nrow(mydata)) {
  for (j in 1:(ncol(mydata)/3)) {
    #取1,4,7.。。等差数列的列
    index<-(3*j-2):(3*j)
    temp<-mydata[i,index]
    temp<-temp[order(temp)]
    temp<-as.numeric(temp)
    temp2<-temp[3]-temp[2]
    temp3<-temp[2]-temp[1]
    #如果差值都两两大于0.5,那直接返回0,如果只有一个大于0.5,则返回另两个值的均值,否则返回三个值的均值
    if(temp2>0.5&&temp3>0.5){print(paste("你要注意!我觉得",rownames(mydata)[i],":",colnames(mydata)[3*j-2],"不能要!",sep=""));temp_ct=NA}else{
      if(temp2>0.5){temp_ct=mean(c(temp[2],temp[1]))}
      if(temp3>0.5){temp_ct=mean(c(temp[3],temp[2]))}
      if(temp3<=0.5&&temp2<=0.5){temp_ct=mean(c(temp[1],temp[2],temp[3]))}
      
    }
    ct[j,i]<-temp_ct
  }
}
colnames(ct)<-rownames(mydata)
rownames(ct)<-unique(colnames(mydata))
#样本靶基因ct-样本内参ct,求出Δct,只支持一个内参

index<-which(rownames(ct)==neican_name)
rownames_deltact<-rownames(ct)
deltact<-sapply(ct,function(x){x-x[index]},simplify = T)
rownames(deltact)<-rownames_deltact

#样本Δct-对照Δct,这里只支持一个对照
if(is.null(control)){double_deltact=deltact}else{
index<-which(colnames(deltact)==control)
double_deltact<-t(apply(deltact,1,function(x){x-mean(x[index])}))
}
#指数运算
power<-2^-double_deltact


上一篇下一篇

猜你喜欢

热点阅读