IC50, drc

ic50=function(
  x=data.frame(),log=TRUE,base=10,plicated=2,sep="..",model=LL.4(),AUCmode='mean',
  file=NA,apch="T"){
  if(nrow(x)==0&ncol(x)==0){
    cat("Drug data is not available",sep='\n')
    cat("ic50(x, log, base, plicated, sep,model, AUCmode, file=apch)",sep='\n')
    cat("1. x : Drug numeric matrix 1st col should be dose concentration",sep='\n')
    cat("       from 2nd column, remaining columns should be drug response data(numeric)",sep='\n')
    cat("For example of data input format")
    print(data.frame(con=c(0.5,2,8,32),drugA1=c(95,80,75,40),drugA2=c(96,88,73,35),drugB1=c(87,60,50,43),drugB2=c(99,55,43,37)))
    cat("2. log : Whether dose concentration is log scaled or not; Defalut 'TRUE'",sep='\n')
    cat("3. base : if log scaled is true, it works. default is 10",sep='\n')
    cat("4. plicated : Determine how many does this experiment conducted Default 2","          Same experiment should be located serially. ex. A1, A2, A3",sep='\n')
    cat("5. sep : For function strsplit, it determines spliting character",sep='\n')
    cat("         Drug will be extracted from header part",sep='\n')
    cat("6. model : Log Logistic model, default LL.4() [Four parameter]",sep='\n')
    cat("7. AUCmode : There are 'min','max','mean',and 'median' mode to calculrate trapezoid AUC value",sep='\n')
    #cat("8. norm(default=T) : Whether response should be fitted and start from 100%",sep='\n')
    cat("8. file : PDF file name. The file will be stored at current Workdir",sep='\n')
    cat("9. apch(default=TRUE) : It is logical and character vector to decide whether last summary figure points should be labelled by ascii code",sep='\n')
    cat("                        If you set true, the points will be ascii code according to their orders(i.e first drug will have !).",sep='\n')
    cat("                        I combined asccii letter with hangul(Korean) so the number of labels is 308",sep='\n')
    cat("                        If you set give numeric value or character value",sep='\n')
    cat("                        The last figure points are converted according to your pch style",sep='\n')
  } else {pdf(paste0(file,'.pdf'))
    print(paste0('Image.pdf will be stored at ', getwd()))
    print("Drug data is imported");print("Attach drc library")
    library("drc");print("drc attached")
    if(log){
      print(paste("Drug data is converted from log to exponential scale.", "Base =",base,sep=" "))
      con=base^x[,1];print('exponentilization finished')
    }else{
      x=as.data.frame(x,stringsAsFactors = F)
      print("Drug concentration is not log, No conversion")
      con=x[,1]};tmp=x[,-1];tmp1=tmp[,1:plicated]
     
      #drug list
      print('Extracting unique drug name')
      cols=1:length(colnames(tmp))
      for(i in 1:length(colnames(tmp))){
        cols[i]=unlist(strsplit(colnames(tmp)[i],fixed="T",sep))[1]
      }
      tmp=tmp[,-c(1:plicated)]
      cols=unique(cols)
      for(i in 1:length(cols)){
        if(length(grep(pattern=cols[i],cols,ignore.case=TRUE))>1){
          cols[max(grep(pattern=cols[i],cols,ignore.case=TRUE))]=NA
        }}
      cols=na.omit(cols)
      print('unique drug name list finished')
      print('Make table for three columns; concentration, response, and label')
     
      tmp1[,plicated+1]=1;rownames(tmp1)=paste0('V',1:nrow(tmp1))
      colnames(tmp1)=c(paste('TRY',1:plicated,sep=""),'label')
      for(i in 1:(ncol(tmp)/plicated)){
        tmp2=tmp[,1:plicated]
        tmp=tmp[,-c(1:plicated)]
        tmp2[,plicated+1]=i+1
        colnames(tmp2)=c(paste('TRY',1:plicated,sep=""),'label')
        tmp1=rbind(tmp1,tmp2)
      } # tmp1 is m x ? matrix
      tmp2=tmp1[,c(1,plicated+1)]
      colnames(tmp2)=c('response','label')
      for(i in 2:(ncol(tmp1)-1)){
        tmp3=tmp1[,c(i,plicated+1)]
        colnames(tmp3)=c('response','label')
        tmp2=rbind(tmp2, tmp3)
      }
      colnames(tmp2)=c('response','label')
      tmp2=tmp2[order(tmp2$label, decreasing=FALSE),]
      tmp2$conc=con
      IC50=data.frame(
        Drug = cols,
        estiIC50 = NA, trapezoid_AUC=NA, drc_AUC=NA,
        absIC50=NA, LoBound=NA, UpBound=NA, hill=NA
      );print("table read to be used")
      for(i in 1:max(tmp2$label)){ #73
        tryCatch({ #124i=59
          tmp3=tmp2[which(tmp2$label==i),]
          q=drm(response~conc, label, data=tmp3,fct=model)
           # 2016Mar18 AUC calculation
            hillslope=q$coefficients[1]
            bottom=q$coefficients[2];top=q$coefficients[3];mid=mean(c(bottom,top));EC50=q$coefficients[4]
            fxdrc=function(r){bottom + ((top-bottom)/(1+(r/EC50)^hillslope))}
            IC50[i,2]=EC50 # relative IC50
            IC50[i,5]=(((top-50)/(50-bottom))^(1/hillslope))*EC50 # Absolute IC50
            IC50[i,6]=bottom;IC50[i,7]=top;IC50[i,8]=hillslope
          #if(IC50[i,2]>max(tmp3$conc)|IC50[i,2]<min(tmp3$conc)){
           # IC50[i,2]=NA} # Whether calculated IC50 is out of from drug region
          #rn=fx(c(runif(100,min=min(con),max=max(con)),min(con)*0.9,max(con)*1.1))
          #if(any(bottom>rn)&any(top<rn)){}else{# Whether bottom and top is within drc curve
                                               # I decided to give 20% points more range to original drug range
            #IC50[i,2]=NA} # should satisfy top<rn and bottom>rn or NA will be reported
              if(log(min(tmp2$conc),base)<0){
                abs(log(min(tmp2$conc),base))->K
                IC50[i,4]=integrate(fxdrc, lower=log(min(tmp2$conc),base)+K,upper=log(max(tmp2$conc),base)+K)$value
              } else {
                IC50[i,4]=integrate(fxdrc, lower=log(min(tmp2$conc),base),upper=log(max(tmp2$conc),base))$value
              };rm(q)}, error=function(e){})
           # trapezoid AUC
            tmp4=tmp3
            if(AUCmode=='mean'){
             tmp4=aggregate(tmp4,by=list(tmp4$conc),FUN = mean)}else{
            if(AUCmode=='median'){
              tmp4=aggregate(tmp4, by=list(tmp4$conc),FUN=median)}else{
            if(AUCmode=='max'){
              tmp4=aggregate(tmp4, by=list(tmp4$conc),FUN=max)}else{
            if(AUCmode=='min'){
              tmp4=aggregate(tmp4,by=list(tmp4$conc),FUN=min)}}}}
            size=1:(nrow(tmp4)-1)
            tmp4$conc=log(tmp4$conc,base);tmp4=tmp4[order(tmp4$conc),]
            for(j in 1:(nrow(tmp4)-1)){
               size[j]=(tmp4$conc[j+1]-tmp4$conc[j])*(tmp4$response[j+1]+tmp4$response[j])/2
             };IC50[i,3]=sum(size)
            # plotting
             plot(tmp4$conc,tmp4$response,main=IC50[i,1],xlim=c(min(log(tmp3$conc,base)),max(log(tmp3$conc,base))),"b",
                 ylim=c(0,max(tmp3$response)*1.05),xaxt="n", ylab="Response",xlab="Dose",las=1,pch=20)
            abline(v=log(IC50[i,2],base),col=2)
            axis(side=1, at=tmp4$conc, labels=format(tmp4$Group.1,digits =2,nsmall= 3),cex.axis=0.7)
              for(p in 1:nrow(tmp4)){
                arrows(tmp4$conc[p],min(tmp3$response[which(tmp3$conc==tmp4$Group.1[p])]),
                       tmp4$conc[p],max(tmp3$response[which(tmp3$conc==tmp4$Group.1[p])]), code=3, length=0.05, angle=90)}
            if(exists('fxdrc')){
              tmp5=tmp4
              tmp5=tmp5[,c(1,2)];colnames(tmp5)=c('x','y');tmp5[,2]=fxdrc(tmp5[,1])
              xa=log(tmp5$x,base);lxa=length(xa)-1
              for(iq in 1:lxa){me=c(xa[iq],xa[iq+1]);xa=append(xa,c((min(me)/3)+max(me)*2/3,mean(me),((min(me)*2/3)+max(me)/3)))};xa=sort(xa)
              lines(xa,fxdrc(base^xa), col=2, lty=2);rm(fxdrc);rm(bottom);rm(top);rm(mid);rm(EC50)}
            if(is.na(IC50[i,2])){}else{
              mtext(text=format(IC50[i,2],digits=3,nsmall=3),at=log(IC50[i,2],base),side=1,cex = 0.5) 
            }
            mtext(text=paste('trapezoid=',format(IC50[i,3],digits=3,nsmall=3),sep=""),side=3,at=max(tmp4$conc)*0.9,cex=0.7)
            mtext(text=paste("absIC50=",format(IC50[i,5],digits=3,nsmall=3),sep=""),side=3,at=max(tmp4$conc)*0.95,line=0.7,cex=0.7)
            mtext(text=paste('AUC=',format(IC50[i,4],digits=3,nsmall=3),sep=""),side=3,at=max(tmp4$conc)*0.96,cex=0.7,line=1.4)
            mtext(text=paste("IC50=",format(IC50[i,2],digits=3,nsmall=3),sep=""),
                  side=3,at=max(tmp4$conc)*0.9,cex=0.7, line=2.1)
            abline(h=50,col=2,lty=3);abline(v=log(IC50[i,5],base),col=2,lty=3)
            if(is.na(IC50[i,5])){}else{
              mtext(text=format(IC50[i,5],digits=3,nsmall=3),side=1,at=log(IC50[i,5],base),cex=0.5)
            }
            # making trapezoidal drug response curve from all integrated drugs
            if(i==max(tmp2$label)){
              plot(NULL,ylim=c(min(tmp2$response),max(tmp2$response)),xlim=c(min(log(tmp2$conc,base)),max(log(tmp2$conc,base))),
                    xlab="Drug",ylab="Response",main="All drug response", xaxt="n")
              axis(side=1,at=log(tmp2$conc,base),labels = format(tmp2$conc,digits = 2, nsmall=3))
            wq=c();for(y in 1:max(tmp2$label)){
              k=tmp2[which(tmp2$label==y),];if(AUCmode=='mean'){k2=aggregate(k,by=list(k$conc),FUN = mean)}else{
                                            if(AUCmode=='median'){k2=aggregate(k, by=list(k$conc),FUN=median)}else{
                                            if(AUCmode=='max'){k2=aggregate(k, by=list(k$conc),FUN=max)}else{
                                            if(AUCmode=='min'){k2=aggregate(k,by=list(k$conc),FUN=min)}}}}
              if(apch){sty=asha[y]}
              if(is.numeric(apch)){sty=apch}
              if(is.character(apch)){sty=apch}
              points(log(k2$conc,base),k2$response,"b",
                     col=colorRampPalette(c("blue","cyan","skyblue","yellow","red"))(max(tmp2$label))[y],pch=sty)
              wq=c(wq,colorRampPalette(c("blue","cyan","skyblue","yellow","red"))(max(tmp2$label))[y])}
                legend('topright',legend=cols,pch = sty,col = wq);dev.off()}
                #for(yo in 1:nrow(tmp2)){
                 # arrows(log(k2$conc[yo],base),min(k$response[which(k$conc==k2$Group.1[yo])]),
                  #       log(k2$conc[yo],base),max(k$response[which(k$conc==k2$Group.1[yo])]), code=3, length=0.05, angle=90)}
             
            }
              }
      return(IC50)}