library(survival)
library(ggplot)


dat <- read.csv("") #file name


names(dat)<-c("x","tmb.level","Gender","tmb","smoke","ethnicity.y","status.y","time.met.y","age.60","tissue.y","log.tmb","log.tmb.squared")

#recode factor variables as appropriate
tmb.levs <- c("Low", "Intermediate", "High", "Very High")
dat$tmb.level <- factor(tmb.levs[dat$tmb.level], levels=tmb.levs)
dat$smoke <- factor(c("No","Yes")[dat$smoke+1])
dat$tissue.y <- relevel(dat$tissue.y, "Other")
dat$age.60 <- factor(c(" <60", ">=60")[dat$age.60+1])
dat$X <- factor(dat$X)
dat$ethnicity.y <- sub("other", "Other", dat$ethnicity.y)
dat$ethnicity.y <- factor(sub("African american", "African-American", dat$ethnicity.y))
dat$ethnicity.y <- relevel(dat$ethnicity.y, "White")
dat$status.y <- dat$status.y - 1

facts <- c("age.60", "Gender", "ethnicity.y", "smoke", "tissue.y", "tmb.level")

unis <- list()
logr <- list()
chisqs <- matrix(NA, nrow=length(facts), ncol=length(facts))
rownames(chisqs) <- facts
colnames(chisqs) <- facts

for(f in facts){
  logr[[f]] <- survfit(formula(paste0("Surv(time.met.y, status.y) ~ ", f)), dat)
  unis[[f]] <- coxph(formula(paste0("Surv(time.met.y, status.y) ~ ", f)), dat)
  for(g in facts){
    if(match(f, facts) < match(g, facts)){
      chisqs[f,g] <- chisq.test(dat[,f], dat[,g])$p.value
    }
  }
}
unis[["poly.tmb"]] <- coxph(formula(Surv(time.met.y, status.y) ~ poly(log(tmb),2)), dat)

CIs <- lapply(unis, function(x)round(exp(cbind(x$coef, confint(x))),2))
pvals <- lapply(unis, function(x)x$pval)


#Repeating the above for a multivariable adjusted curve averaging the effect of these other covariables
full.mod <- coxph(formula(paste0("Surv(time.met.y, status.y) ~ ", paste(c("age.60", "ethnicity.y", "smoke", "tissue.y", "tmb.level"), collapse=" + "))), dat)

round(summary(full.mod)$conf.int,2)[,-2]
round(summary(update(full.mod, . ~ . - age.60 - smoke))$conf.int,2)[,-2]
full.mod.cont <- coxph(formula(paste0("Surv(time.met.y, status.y) ~ ", paste(c("age.60", "ethnicity.y", "smoke", "tissue.y", "poly(log(tmb),2)"), collapse=" + "))), dat)

xs <- seq(min(dat$tmb), max(dat$tmb), len=1000)
ys.uni <- exp(predict(unis[["poly.tmb"]], data.frame(tmb=xs)))
ys.multi <- predict(full.mod.cont, data.frame(age.60=" <60", ethnicity.y = "White", smoke = "No", tissue.y="Other", tmb=xs))
ys.multi <- exp(ys.multi - mean(predict(full.mod.cont, data.frame(age.60=" <60", ethnicity.y = "White", smoke = "No", tissue.y="Other", tmb=dat$tmb))))


par(cex=1.5)
plot(log(xs), ys.uni, type="l", main="Hazard ratio vs. log(TMB)", xlab="log(TMB)", ylab="Hazard ratio",
     col="red", lwd=2)
lines(log(xs), ys.multi, col="blue", lwd=2)
tmb.cuts <- c(1, 5.5, 19.5, 49.5, 860)
abline(v=log(tmb.cuts[2:4]), lty=2)
text(0, 0.275, "TMB level:", adj=c(0.05, -1))
text(log(tmb.cuts[-1]) - diff(log(tmb.cuts)/2), 0.2, c("Low", "Inter.", "High", "Very high"), pos=3)
legend("topright", c("Univariate", "Multivariate"), col=c("red", "blue"), lwd=2, bg="white",cex=0.8)


#############

#Coding the patient-specific point curves stratified by tumor type
library("ggsci") 
Color<-brewer.pal(n=9, name = "Set1")[c(1,2,3,4,5,7,9)]
tumor.type<-c("Brain","Breast","Colorectal","Heme","Lung","Skin")
d.figure1b<-univariate.x
d.figure1b$Tissue<-factor(d.figure1b$tissue,levels = c("Brain","Breast","Colorectal","Heme","Lung","Skin","Other"))
cox5<-coxph(Surv( time.met ,status) ~ smoke + Tissue + ethnicity + log.tmb  + log.tmb.squared , 
            data = d.figure1b)

p.fig1C<-ggplot(data = d.figure1b, 
                aes(x = tmb,y=exp(predict(cox5,d.figure1b)),colour = Tissue) )  + 
  geom_point(aes(x = tmb,y=exp(predict(cox5,d.figure1b)),colour = Tissue),size = 4,shape = 18) + 
  scale_x_continuous(trans='log10',breaks = c(0,1,5,10,50,100,500,1000),minor_breaks = FALSE,name = "Tumor Mutational Burden") +
  scale_y_continuous(breaks = c(0,0.5,1,1.5,2,2.5,3,3.5),name = "Hazard Ratio") + expand_limits(x=1,y=0) +
  scale_color_manual(values = Color) + 
  theme_minimal() +
  theme(axis.text.x = element_text(size = 20,colour = "black"),
        axis.title.x = element_text(size = 25,colour = "black"),
        axis.text.y = element_text(size = 20,colour = "black"),
        axis.title.y = element_text(size = 25,colour = "black"),
        legend.title = element_text(size = 25,colour = "black"),
        legend.text = element_text(size = 20,colour = "black")
  )

# ggsave(filename = "",plot = p.fig1C,device = "tiff",dpi = 300,
       # width = 10,height = 8)