--- a
+++ b/vignettes/src/part08.Rmd
@@ -0,0 +1,595 @@
+---
+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())
+```