一个处理qpcr数据的程序(第二版,未debug)

2020-06-14  本文已影响0人  一只烟酒僧
######################################################## 
#-------------------------------------------------------
# Topic:进行qRTpcr运算
# Author:Wang Haiquan
# Date:Sat Jun 13 21:28:16 2020
# Mail:mg1835020@smail.nju.edu.cn
#-------------------------------------------------------
########################################################
#列是基因重复,第一个基因为对照
#行是样本名
set.seed(1111)
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
neican_name<-"gene1"
control="sample1"
##############################################################################
#将复孔取平均值,同时顺带把差值大于0.5的复孔删除,同时支持不同复孔的实验
mydata<-read.table("clipboard",sep = "\t",row.names = 1,header = T)
index_gene<-factor(colnames(mydata))
#按照基因把数据框拆成列表
mydata<-split(as.data.frame(t(mydata)),index_gene)
mydata<-lapply(mydata,t)
mydata<-lapply(mydata,function(x){a=ncol(x);if(a==1){return(x)}else{
  ct=as.data.frame(matrix(0,ncol = 1,nrow = nrow(x)))
  colnames(ct)<-colnames(x)[1]
  rownames(ct)<-rownames(x)
  for (i in 1:nrow(x)) {
    temp<-x[i,]
    temp<-temp[order(temp)]
    temp1<-c()
    for (j in 2:length(temp)) {
      temp2<-temp[j]-temp[j-1]<0.5
      temp1=c(temp1,temp2)
    }
    if(Reduce(sum,temp1)==length(temp)-1){temp_ct=mean(temp)}else if(Reduce(sum,temp1)==0){
      print(paste("我认为",rownames(x)[i],":",colnames(x)[1],"有点问题",sep=""));temp_ct=0}else{
        index_ct<-which(temp1==T); temp_ct=mean(c(temp[index_ct],temp[index_ct+1]))
      }
    ct[i,1]<-temp_ct
  }
  return(ct)
}})
ct<-Reduce(cbind,mydata)
ct<-t(ct)
ct<-as.data.frame(ct)
#样本靶基因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-x[index]}))
}
#指数运算
power<-2^-double_deltact
power

上一篇下一篇

猜你喜欢

热点阅读