---

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())

```