---
title: "Part 8"
output:
  BiocStyle::html_document
---

```{r, message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")}
library("BloodCancerMultiOmics2017")
library("Biobase")
library("RColorBrewer")
library("grid")
library("ggplot2")
library("survival")
library("gtable")
library("forestplot")
library("xtable")
library("maxstat")
```

```{r echo=FALSE}
plotDir = ifelse(exists(".standalone"), "", "part08/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
```

```{r}
options(stringsAsFactors=FALSE)
```

# Survival analysis

Load the data.
```{r}
data(lpdAll, patmeta, drugs)
```

Prepare survival data.
```{r}
lpdCLL <- lpdAll[ , lpdAll$Diagnosis=="CLL"   ]

# data rearrangements
survT = patmeta[colnames(lpdCLL),]
survT[which(survT[,"IGHV"]=="U") ,"IGHV"] = 0
survT[which(survT[,"IGHV"]=="M") ,"IGHV"] = 1
survT$IGHV = as.numeric(survT$IGHV)

colnames(survT) = gsub("Age4Main", "age", colnames(survT))

survT$ibr45  <- 1-Biobase::exprs(lpdCLL)[ "D_002_4:5", rownames(survT)  ]  
survT$ide45  <- 1-Biobase::exprs(lpdCLL)[ "D_003_4:5", rownames(survT)  ] 
survT$prt45  <- 1-Biobase::exprs(lpdCLL)[ "D_166_4:5", rownames(survT)  ] 
survT$selu45 <- 1-Biobase::exprs(lpdCLL)[ "D_012_4:5", rownames(survT)  ]
survT$ever45 <- 1-Biobase::exprs(lpdCLL)[ "D_063_4:5", rownames(survT)  ]

survT$nut15 <- 1-Biobase::exprs(lpdCLL)[ "D_010_1:5", rownames(survT)  ] 
survT$dox15 <- 1-Biobase::exprs(lpdCLL)[ "D_159_1:5", rownames(survT)  ] 
survT$flu15 <- 1-Biobase::exprs(lpdCLL)[ "D_006_1:5", rownames(survT)  ] 

survT$SF3B1      <- Biobase::exprs(lpdCLL)[ "SF3B1",      rownames(survT)  ]
survT$NOTCH1     <- Biobase::exprs(lpdCLL)[ "NOTCH1",     rownames(survT)  ]
survT$BRAF       <- Biobase::exprs(lpdCLL)[ "BRAF",       rownames(survT)  ]
survT$TP53       <- Biobase::exprs(lpdCLL)[ "TP53",       rownames(survT)  ]
survT$del17p13   <- Biobase::exprs(lpdCLL)[ "del17p13",   rownames(survT)  ]
survT$del11q22.3 <- Biobase::exprs(lpdCLL)[ "del11q22.3", rownames(survT)  ]
survT$trisomy12 <-  Biobase::exprs(lpdCLL)[ "trisomy12", rownames(survT)  ]
survT$IGHV_cont <- patmeta[ rownames(survT) ,"IGHV Uppsala % SHM"]


# competinting risk endpoint fpr 
survT$compE <- ifelse(survT$treatedAfter == TRUE, 1, 0)
survT$compE <- ifelse(survT$treatedAfter == FALSE & survT$died==TRUE,
                      2, survT$compE )
survT$T7  <- ifelse(survT$compE == 1, survT$T5, survT$T6 )
```


## Univariate survival analysis

### Define forest functions

```{r forest}
forest <- function(Time, endpoint, title, sdrugs, split, sub) {  
  stopifnot(is.character(Time), is.character(title), is.character(split),
            is.character(endpoint), 
            all(c(Time, split, endpoint) %in% colnames(survT)),
            is.logical(survT[[endpoint]]),
            is.character(sdrugs), !is.null(names(sdrugs)))
  
  clrs <- fpColors(box="royalblue",line="darkblue", summary="royalblue")

  res <- lapply(sdrugs, function(g) { 
    drug <- survT[, g] * 10

    suse <- if (identical(sub, "none")) 
      rep(TRUE, nrow(survT)) 
    else 
      (survT[[split]] == sub)
    stopifnot(sum(suse, na.rm = TRUE) > 1)
    
    surv <- coxph(Surv(survT[,Time], survT[,endpoint]) ~ drug, subset=suse)  
    sumsu <- summary(surv) 
    c(p      = sumsu[["coefficients"]][, "Pr(>|z|)"], 
      coef   = sumsu[["coefficients"]][, "exp(coef)"], 
      lower  = sumsu[["conf.int"]][, "lower .95"], 
      higher = sumsu[["conf.int"]][, "upper .95"])
  })
  
  s <- do.call(rbind, res)
  rownames(s) <- names(sdrugs)

  tabletext <- list(c(NA, rownames(s)), append(list("p-value"),
                                               sprintf("%.4f", s[,"p"])))

  forestplot(tabletext, 
           rbind(
             rep(NA, 3), 
             s[, 2:4]), 
           page = new,
           clip = c(0.8,20), 
           xlog = TRUE, xticks = c(0.5,1, 1.5), title = title,
           col = clrs, 
           txt_gp = fpTxtGp(ticks = gpar(cex=1) ),
           new_page = TRUE)
}
```

Combine OS and TTT in one plot.
```{r forest-together}

com <- function( Time, endpoint, scaleX, sub, d, split, drug_names) {  
  
  res <- lapply(d, function(g)  { 
  
  drug <- survT[,g] * scaleX
  ## all=99, M-CLL=1, U-CLL=0
  if(sub==99) { surv <- coxph(Surv(survT[,paste0(Time)],
                                   survT[,paste0(endpoint)] == TRUE) ~ drug)} 
  if(sub<99)  { surv <- coxph(Surv(survT[,paste0(Time)],
                                   survT[,paste0(endpoint)] == TRUE) ~ drug,
                              subset=survT[,paste0(split)]==sub)}    
  
  c(summary(surv)[[7]][,5], summary(surv)[[7]][,2], 
                summary(surv)[[8]][,3], 
                summary(surv)[[8]][,4])
 })
 s <- do.call(rbind, res)
 colnames(s) <- c("p", "HR", "lower", "higher")
 rownames(s) <- drug_names
 s
}


fp <- function( sub, title, d, split, drug_names, a, b, scaleX) {  
   ttt <- com(Time="T5", endpoint="treatedAfter", sub=sub, d=d,
              split=split, drug_names=drug_names, scaleX=scaleX)
   rownames(ttt) <- paste0(rownames(ttt), "_TTT")
   
   os <-  com(Time="T6", endpoint="died", sub=sub, d=d, split=split,
              drug_names=drug_names, scaleX=scaleX)
   rownames(os) <- paste0(rownames(os), "_OS")
  
   n <- c( p=NA, HR=NA, lower=NA, higher=NA )
   nn <- t( data.frame( n ) )
   for (i in 1:(nrow(ttt)-1) ) { nn <-rbind(nn, n )  }
   rownames(nn) <- drug_names
   
   od <- order( c(seq(nrow(nn)), seq(nrow(ttt)), seq(nrow(os)) ))
      
   s <- data.frame( rbind(nn, ttt, os)[od,  ] )
   s$Name <- rownames(s)
   s$x <- 1:nrow(s)
   s$col <- rep(c("white", "black", "darkgreen"), nrow(ttt) )
   s$Endpoint <- factor( c(rep("nn", nrow(nn) ), rep("TTT", nrow(ttt) ),
                           rep("OS", nrow(os) ) )[od] )
   s$features <- "";  s[ which(s$Endpoint=="OS"),"features"] <- drug_names
   s[which(s$Endpoint=="nn"), "Endpoint"] <- "" 
   s <- rbind(s, rep(NA, 8))

   p <- ggplot(data=s ,aes(x=x, y=HR, ymin=lower, ymax=higher,
                           colour=Endpoint)) +  geom_pointrange() + 
     theme(legend.position="top", legend.text = element_text(size = 20) ) +
     scale_x_discrete(limits=s$x, labels=s$features ) +
     expand_limits(y=c(a,b)) +
     scale_y_log10(breaks=c(0.01,0.1,0.5,1,2,5,10),
                   labels=c(0.01,0.1,0.5,1,2,5,10)) +
     theme(
           panel.grid.minor = element_blank(),
           axis.title.x = element_text(size=16),
           axis.text.x = element_text(size=16, colour="black"),
           axis.title.y  = element_blank(),
           axis.text.y = element_text(size=12, colour="black"),
           legend.key = element_rect(fill = "white"),
           legend.background = element_rect(fill = "white"),
           legend.title = element_blank(),
           panel.background = element_rect(fill = "white", color="black"),
           panel.grid.major = element_blank(),
           axis.ticks.y = element_blank() 
          ) +
     coord_flip() + 
     scale_color_manual(values=c("OS"="darkgreen", "TTT"="black"),
                        labels=c("OS", "TTT", "")) +
     geom_hline(aes(yintercept=1), colour="black", size=1.5,
                linetype="dashed", alpha=0.3)  +
     annotate("text", x = 1:nrow(s)+0.5, y = s$HR+0.003,
              label = ifelse( s$p<0.001, paste0("p<","0.001"), 
                paste0("p=", round(s$p,3) ) ), colour=s$col)
   plot(p)
}
```


### Forest plot for genetic factors

```{r Fig5A, fig.path=plotDir, fig.width=5.55, fig.height=( (1+8*1.2) ), dev = c("png", "pdf"), warning=FALSE}
#FIG# S27
d <- c("SF3B1", "NOTCH1", "BRAF", "TP53", "del17p13", "del11q22.3",
       "trisomy12", "IGHV")
drug_names <- c("SF3B1", "NOTCH1", "BRAF", "TP53", "del17p13", "del11q22.3",
                "Trisomy12" ,"IGHV")

fp(sub=99, d=d, drug_names=drug_names, split="IGHV", title="", a=0, b=10,
   scaleX=1)
```

### Forest plot for drug responses

```{r Fig5B, fig.path=plotDir, fig.width=6.0, fig.height=( (1+8*1.2) ), dev = c("png", "pdf"), warning=FALSE}
#FIG# 6A
d <- c("flu15", "nut15", "dox15", "ibr45", "ide45", "prt45", "selu45",
       "ever45")
drug_names <- c("Fludarabine",  "Nutlin-3", "Doxorubicine", "Ibrutinib",
                "Idelalisib", "PRT062607 HCl", "Selumetinib" ,"Everolimus")

fp(sub=99, d=d, drug_names=drug_names, split="TP53", title="", a=0, b=5,
   scaleX=10)
```

## Kaplan-Meier curves 

### Genetics factors

```{r SFig_genetics1, fig.path=plotDir, fig.width = 8.46, fig.height = 4.5, dev = c("png", "pdf")}
#FIG# S27 left (top+bottom)
par(mfcol=c(1,2))

for (fac in paste(c("IGHV", "TP53"))) {
survplot( Surv(survT$T5, survT$treatedAfter == TRUE)  ~ as.factor(survT[,fac]), 
   snames=c("wt", "mut"),
   lwd=1.5, cex.axis = 1, cex.lab=1, col= c("darkmagenta", "dodgerblue4"),
   show.nrisk = FALSE,
   legend.pos = FALSE, stitle = "", hr.pos= "topright",
   main = paste(fac), 
   xlab = 'Time (Years)', ylab = 'Time to treatment')
}
```

```{r SFig_genetics2, fig.path=plotDir, fig.width = 8.46, fig.height = 4.5, dev = c("png", "pdf")}
#FIG# 6B left
#FIG# S27 right (top+bottom)
par(mfcol=c(1,2))

for (fac in paste(c("IGHV", "TP53"))) {
survplot( Surv(survT$T6, survT$died == TRUE)  ~ as.factor(survT[,fac]), 
   snames=c("wt", "mut"),
   lwd=1.5, cex.axis = 1.0, cex.lab=1.0, col= c("darkmagenta", "dodgerblue4"),
   show.nrisk = FALSE,
   legend.pos = FALSE, stitle = "", hr.pos= "bottomleft",
   main = paste(fac), 
   xlab = 'Time (Years)', ylab = 'Overall survival')
}  
```


### Drug responses

Drug responses were dichotomized using maximally selected rank statistics. The analysis is also perforemd within subgroups: TP53 wt/ mut. and IGHV wt/ mut. 

```{r KM, echo=FALSE}
km <- function(drug, split, title, t, hr, c) { 
  stopifnot(is.character(drug), length(drug)==1, is.character(title),
            length(title)==3)

  surv <- survT[!(is.na(survT[,split])), ]
  k <- Biobase::exprs(lpdCLL)[ drug, rownames(surv) ]
  
  ms5 <- maxstat.test(Surv(T5, treatedAfter)  ~ k, 
                             data = surv,
                             smethod = "LogRank",
                             minprop = 0.2, 
                             maxprop = 0.8, 
                             alpha = NULL)
  ms6 <- maxstat.test(Surv(T6, died) ~ k, 
                             data = surv,
                             smethod = "LogRank",
                             minprop = 0.2, 
                             maxprop = 0.8, 
                             alpha = NULL)

  
  # median & TTT
  if (c=="med" & t=="TTT") {    
   surv$cutA <- ifelse(
     k >= median( k[which(!(is.na(surv$T5)) ) ] ), "weak", "good")
   surv$cutM <- ifelse(
     k >= median( k[ which( surv[,paste0(split)]==1 & !(is.na(surv$T5)) ) ],
                  na.rm=TRUE ), "weak", "good") 
   surv$cutU <- ifelse(
     k >= median( k[ which( surv[,paste0(split)]==0 & !(is.na(surv$T5)) ) ],
                  na.rm=TRUE ), "weak", "good")
   
  }
  
  # median & OS
  if (c=="med" & t=="OS") {    
   surv$cutA <- ifelse(k >= median(k), "weak", "good")
   surv$cutM <- ifelse(k >= median( k[ which( surv[,paste0(split)]==1 ) ] ),
                       "weak", "good") 
   surv$cutU <- ifelse(k >= median( k[ which( surv[,paste0(split)]==0 ) ] ),
                       "weak", "good")
  }
  
  #TTT & maxstat
  if (c=="maxstat" & t=="TTT") {    
   surv$cutA <- surv$cut5 <- ifelse(k >= ms5$estimate, "weak", "good")
   surv$cutM <- surv$cut5 <- ifelse(k >= ms5$estimate, "weak", "good") 
   surv$cutU <- surv$cut5 <- ifelse(k >= ms5$estimate, "weak", "good")
  }
  
  #OS & maxstat
  if (c=="maxstat" & t=="OS") {    
   surv$cutA <- surv$cut5 <- ifelse(k >= ms6$estimate, "weak", "good")
   surv$cutM <- surv$cut5 <- ifelse(k >= ms6$estimate, "weak", "good") 
   surv$cutU <- surv$cut5 <- ifelse(k >= ms6$estimate, "weak", "good")
  }

  drName <- toCaps(drugs[stripConc(drug), "name"])
  
   sp <- function(...)
    survplot(..., 
      lwd = 3, cex.axis = 1.2, cex.lab = 1.5,
      col= c("royalblue", "darkred"), show.nrisk = FALSE,
      legend.pos = FALSE, stitle = "",
      hr.pos=ifelse(hr=="bl", "bottomleft", "topright" ),
      xlab = 'Time (Years)') 
    
  if (t=="TTT") {
         yl <- "Fraction w/o treatment"
         if (c=="med"){
           cat(sprintf("%s median-cutpoint for TTT: %5.2g\n",
                       drName, median(k) ) ) } else 
                       { cat(sprintf("%s cutpoint for TTT: %5.2g\n", drName,
                         ms5$estimate )) }
         
         sp(Surv(surv$T5, surv$treatedAfter) ~ surv$cutA,
            subset = rep(TRUE, nrow(surv)), ylab = yl, main = drName)
         sp(Surv(surv$T5, surv$treatedAfter) ~ surv$cutM,
            subset = surv[, split]==1, ylab = yl,
            main = paste(drName, title[1], title[3])) 
         sp(Surv(surv$T5, surv$treatedAfter) ~ surv$cutU,
            subset = surv[ ,split]==0, ylab = yl,
            main = paste(drName, title[1], title[2])) }
   # OS  
   else {
         yl <- "Fraction overall survival"
         if (c=="med"){
           cat(sprintf("%s median-cutpoint for OS: %5.2g\n",
                       drName, median(k) ) ) } else {
                         cat(sprintf("%s cutpoint for OS: %5.2g\n",
                                     drName, ms6$estimate ))}
         sp(Surv(surv$T6, surv$died) ~ surv$cutA,
            subset = rep(TRUE, nrow(surv)), ylab = yl, main = drName)
         sp(Surv(surv$T6, surv$died) ~ surv$cutM,
            subset = surv[, split]==1, ylab = yl,
            main = paste(drName, title[1], title[3])) 
         sp(Surv(surv$T6, surv$died) ~ surv$cutU,
            subset = surv[ ,split]==0, ylab = yl,
            main = paste(drName, title[1], title[2]))
         }
}
```


Time to next treatment (maxstat).
```{r KM-TTT-maxstat, fig.path=plotDir, fig.width = 10, fig.height = 3.3, dev = c("png", "pdf")}
par(mfrow=c(1,3), mar=c(5,5,2,0.9))
km(drug = "D_006_1:5", split = "TP53", t="TTT",
   title=c("(TP53", "wt)", "mut)"),  hr="tr", c="maxstat")
km(drug = "D_159_1:5", split = "TP53", t="TTT",
   title=c("(TP53", "wt)", "mut)"), hr="tr",  c="maxstat")
km(drug = "D_010_1:5", split = "TP53", t="TTT",
   title=c("(TP53", "wt)", "mut)"), hr="tr",  c="maxstat")  

km(drug = "D_002_4:5", split = "IGHV", t="TTT",
   title=c("(IGHV",  "wt)" , "mut)"), hr="tr", c="maxstat" )
km(drug = "D_003_4:5", split = "IGHV", t="TTT",
   title=c("(IGHV",  "wt)" , "mut)"), hr="tr", c="maxstat" )
km(drug = "D_166_4:5", split = "IGHV", t="TTT",
   title=c("(IGHV",  "wt)" , "mut)"), hr="tr", c="maxstat" ) 
```

Overall survival (maxstat).
```{r KM-OS-maxstat, fig.path=plotDir, fig.width = 10, fig.height = 3.3, dev = c("png", "pdf")}
par(mfrow=c(1,3), mar=c(5,5,2,0.9))
km(drug = "D_006_1:5", split = "TP53", t="OS",
   title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat")

#FIG# 6B right
#FIG# 6C
km(drug = "D_159_1:5", split = "TP53", t="OS", # doxorubicine
   title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat" )

#FIG# 6B middle
km(drug = "D_010_1:5", split = "TP53", t="OS", # nutlin-3
   title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat" )

km(drug = "D_002_4:5", split = "IGHV", t="OS",
   title=c("(IGHV",  "wt)" , "mut)"), hr="bl", c="maxstat" )
km(drug = "D_003_4:5", split = "IGHV", t="OS",
   title=c("(IGHV",  "wt)" , "mut)"), hr="bl", c="maxstat" )
km(drug = "D_166_4:5", split = "IGHV", t="OS",
   title=c("(IGHV",  "wt)" , "mut)"), hr="bl", c="maxstat" ) 
```


## Multivariate Cox-model

```{r extract}
extractSome <- function(x) {
  sumsu <- summary(x)
  data.frame(
    `p-value`      = 
      sprintf("%6.3g", sumsu[["coefficients"]][, "Pr(>|z|)"]),
    `HR`           = 
      sprintf("%6.3g", signif( sumsu[["coefficients"]][, "exp(coef)"], 2) ), 
    `lower 95% CI` = 
      sprintf("%6.3g", signif( sumsu[["conf.int"]][, "lower .95"], 2) ),
    `upper 95% CI` = 
      sprintf("%6.3g", signif( sumsu[["conf.int"]][, "upper .95"], 2),
              check.names = FALSE) )
}
```

Define covariates and effects.
```{r covariates, echo=FALSE}
survT$age <- survT$age/10
survT$IC50beforeTreatment <- ifelse(survT$IC50beforeTreatment==TRUE, 1, 0)
survT$IGHVwt <- ifelse(survT$IGHV==1, 0, 1)

survT$flu15 <- survT$flu15*10
survT$dox15 <- survT$dox15*10

survT$ibr45 <- survT$ibr45*10
survT$ide45 <- survT$ide45*10
survT$prt45 <- survT$prt45*10
```


### Chemotherapies

#### Fludarabine
```{r}
surv1 <- coxph(
  Surv(T6, died) ~  
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    as.factor(TP53) +
    IGHVwt +
    flu15,       # continuous
    #dox15 +     # continuous
    #flu15:TP53,
    #TP53:dox15,
  data = survT )
extractSome(surv1)

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv1)$n, summary(surv1)[6] ) )
```

```{r echo=FALSE, results='hide', eval=FALSE}
write(print(xtable(extractSome(surv1))), file=paste0(plotDir,"flu_MD.tex"))
```

#### Doxorubicine

```{r}
surv2 <- coxph(
  Surv(T6, died) ~   #as.factor(survT$TP53) , data=survT )
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    as.factor(TP53) +
    IGHVwt +
    #flu15 +    # continuous
    dox15 ,     # continuous
    #flu15:TP53 ,
    #TP53:dox15,
  data = survT )
extractSome(surv2)

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv2)$n, summary(surv2)[6] ) )
```

```{r echo=FALSE, results='hide', eval=FALSE}
write(print(xtable(extractSome(surv2))), file=paste0(plotDir,"dox_MD.tex"))
```

### Targeted therapies

#### Ibrutinib TTT

```{r}
surv4 <- coxph(
  Surv(T5, treatedAfter) ~ 
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    IGHVwt +
    ibr45 +
    IGHVwt:ibr45,
  data = survT )

extractSome(surv4)

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv4)$n, summary(surv4)[6] ) )
```

```{r echo=FALSE, results='hide', eval=FALSE}
write(print(xtable(extractSome(surv4))), file=paste0(plotDir,"ibr_TTT.tex"))
```

#### Idelalisib TTT

```{r}
surv6 <- coxph(
  Surv(T5, treatedAfter) ~ 
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    IGHVwt +
    ide45 +
    IGHVwt:ide45,
  data = survT )

extractSome(surv6)

cat(sprintf("%s patients considerd in the model; number of events %1g\n",
            summary(surv6)$n, summary(surv6)[6] ) )
```

```{r echo=FALSE, results='hide', eval=FALSE}
write(print(xtable(extractSome(surv6))), file=paste0(plotDir,"ide_TTT.tex"))
```

#### PRT062607 HCl TTT
```{r}
surv8 <- coxph(
  Surv(T5, treatedAfter) ~ 
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    IGHVwt +
    prt45 +
    IGHVwt:prt45,
  data = survT )

extractSome(surv8)

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv8)$n, summary(surv8)[6] ) )
```

```{r echo=FALSE, results='hide', eval=FALSE}
write(print(xtable(extractSome(surv8))), file=paste0(plotDir,"prt_TTT.tex"))
```


```{r, include=!exists(".standalone"), eval=!exists(".standalone")}
sessionInfo()
```

```{r, message=FALSE, warning=FALSE, include=FALSE}
rm(list=ls())
```
