survival analysis function

#' survival analysis and get Log-rank p.value
#' 
#' This function is to get Log-rank p-value and Kaplan-meier plot for survival analysis.
#' @param class (default=NULL,factor) : Input vector, 1st row drug response, 2nd mutation
#' @param survival_time (default=NULL,numeric) : Input vector, put them as numeric.
#' @param status (default=NULL,numeric) : Input vector, 0(alive) or 1(dead)
#' @param plot (default=T,logical) : True or False to draw the KM plot
#' @param colors (default=NULL,character) : Colors for each groups
#' @param xlab (default=NULL,character) : Put time scale (days, months, years....)
#' @param title (default=NULL,character) : Put title of plot
#' @param export (default=T,character) : True or False to export the results of survival analysis.
#' @param scientific (default=F) : True or False whether p-value should be scientific notation because of too significant p-value.
#' @param nlarge (default=5) : If scientific is true, this parameter is used to decide whether p-value is too significant or not.\cr
#' 5 means (10^-5). If a p-value is more significant than the nlarge, the p-value will computed as scientific notation (ex. 5E-6)
#' @keywords Log-rank test
#' @return
#' surv_formula : formula used in survfit and survdiff
#' survfit : survfit(surv_formula)
#' p.value : Log-rank p value survdiff(surv_formula)
#' @export
#' @examples
#' class=as.factor(paste0('GPC',rep(c(1,2,3),each=50)))
#' survival_time=c(rnorm(50,10),rnorm(50,100,sd=5),rnorm(50,50,sd=10))
#' status=sample(c(0,1),150,replace=T)
#' #Drawing only
#' survival_analysis(class=cls,survival_time=surv_time,status=status,plot=T,
#' colors=c('skyblue','tomato','darkgreen'),
#' xlab='Months',title='example',export=F,scientific=T)
#' #Result and plotting
#' survival_analysis(class=cls,survival_time=surv_time,status=status,plot=T,
#' colors=c('skyblue','tomato','darkgreen'),
#' xlab='Months',title='Result and plotting',export=T)
  #-----------------------------------------
  # draw plot
  #-----------------------------------------
    survival_analysis=function(class=NULL,survival_time=NULL,status=NULL,plot=T,colors=NULL,
                               xlab=NULL,title=NULL,export=T,scientific=F,nlarge=5){
      frms=survival::Surv(survival_time,status)~class
      os=survival::survfit(frms)
      # Calculate Log-rank p.value
        pv=survival::survdiff(frms)
        pv=pchisq(pv$chisq,length(pv$n)-1,lower.tail = F)
      # Plotting?
        if(plot){
          plot(os,ylab='Survival probability',xlab=xlab,main=title,mark.time=T,cex=5,lwd=5,col=colors)
          # p-value formating
            if(scientific){
              if(pv<10^(-nlarge)){
                pv2=toupper(format(pv,scientific=T,nsmall=3,digits=3))
                pv2=gsub(pv2,pattern = '-0',replacement = '-')
              }else{
                pv2=format(pv,scientific=F,nsmall=3,digits=3)
              }
            }else{
              pv2=format(pv,scientific=F,nsmall=3,digits=3)
            }
          mtext(paste0('Log-rank p=',pv2),adj=1,side=3)
          # Legend
            legend.txt=gsub(paste0(names(os$strata),'(N=',os$n,')'),pattern='class=',replacement = "",fixed=T)
            legend('topright',legend = legend.txt,lwd=5,col=colors)
        }
      # Calculate hazard ratio
        # Get a subtype showing the best prognosis (longest median survival)
          median.os=quantile(os)$quantile[,2]
          median.os=sort(median.os,decreasing = T) # Sort the classes according to median survival period
          sorted.su=gsub(names(median.os),pattern = 'class=',replacement = '',fixed = T)
        # Change class name
          class2=as.factor(multi_sub(as.character(class),pattern = sorted.su,
                                     replacement = 1:length(sorted.su),exact=T))
          frms2=survival::Surv(survival_time,status)~class2
          hazard.ratio=coxph(frms2)
          hazard.ratio2=exp(hazard.ratio$coefficients)
        # Return the hazard ratio
          hazard.ratio3=c(1,hazard.ratio2)
          names(hazard.ratio3)=sorted.su
       
      # Export result?
        if(export){return(list('surv_formula'=frms,'survfit'=os,'p.value'=pv,'hazard.ratio'=hazard.ratio3))}
    }