library(riskRegression)
library(pammtools)
library(adjustedCurves)
library(survival)
library(survminer)

# NOTE: this script requires dn_landscape to have been generated
# Code for this is at the top of metric_boxplots.R

dss = read.csv("path/to/data/sherlock_20240909_derived_stage_and_survival.csv", header=T, stringsAsFactors = F)
dss = cbind(dss, `Broad_Stage`=gsub("[ABC123]", "", dss$STAGE_DERIVED))

# add survival information to landscape matrix with de novo group information
surv_dn_landscape = cbind(dn_landscape, `Stage`=dss[rownames(dn_landscape),"Broad_Stage"])
surv_dn_landscape = cbind(surv_dn_landscape, `Derived_Survival`=dss[rownames(dn_landscape),"SURVIVAL_TIME_WEEKS_DERIVED"])
surv_dn_landscape = cbind(surv_dn_landscape, `event`=1)
surv_dn_landscape[,"DN_group"] = factor(surv_dn_landscape[,"DN_group"], levels=c("NSD-A","NSD-B","SD"))

# fit Cox proportional hazards model on all LUAD, adjusting for stage, sex (stored under the column name "Gender"), and age
fit.coxph <- coxph(Surv(Derived_Survival, event)~ Stage + Gender + Age + DN_group, data = surv_dn_landscape, x=TRUE)
# fit adjusted survival curves for all LUAD, comparing de novo groups
adjsurv <- adjustedsurv(data=surv_dn_landscape,
                        variable="DN_group",
                        ev_time="Derived_Survival",
						event="event",
                        method="direct",
                        outcome_model=fit.coxph,
                        conf_int=TRUE,
						bootstrap=TRUE)
png(file='path/to/supplementary/figures/DN_Adjusted_Survival.png',width = 7,height = 7,units="in",res=300)
plot(adjsurv, custom_colors = c("#8adbf590","#a371cb90","#5a5a6990"), conf_int=TRUE, xlab="Time (weeks)", legend.position="none") + theme(axis.text = element_text(size = 15))
dev.off()

# select only SBS4- subjects who have not smoked
surv_dn_landscape_ns = surv_dn_landscape[surv_dn_landscape$Smoking_SBS4 == "NonSmoker_NonSBS4",]

# fit Cox proportional hazards model on NS-LUAD, adjusting for stage, sex (stored under the column name "Gender"), and age
fit.coxph <- coxph(Surv(Derived_Survival, event)~ Stage + Gender + Age + DN_group, data = surv_dn_landscape_ns, x=TRUE)
# fit adjusted survival curves for NS-LUAD, comparing de novo groups
adjsurv_ns <- adjustedsurv(data=surv_dn_landscape_ns,
                        variable="DN_group",
                        ev_time="Derived_Survival",
						event="event",
                        method="direct",
                        outcome_model=fit.coxph,
                        conf_int=TRUE,
						bootstrap=TRUE)
png(file='path/to/supplementary/figures/DN_NS_Adjusted_Survival.png',width = 7,height = 7,units="in",res=300)
plot(adjsurv_ns, custom_colors = c("#8adbf590","#a371cb90","#5a5a6990"), conf_int=TRUE, xlab="Time (weeks)", legend.position="none") + theme(axis.text = element_text(size = 15))
dev.off()

# convert self-reported smoking status column to factor
surv_dn_landscape[,"Smoking"] = factor(surv_dn_landscape[,"Smoking"], levels=unique(surv_dn_landscape[,"Smoking"]))

# fit Cox proportional hazards model on all LUAD, adjusting for stage, sex (stored under the column name "Gender"), and age
fit.coxph <- coxph(Surv(Derived_Survival, event)~ Stage + Gender + Age + Smoking, data = surv_dn_landscape, x=TRUE)
# fit adjusted survival curves for all LUAD, comparing self-reported smoking status (smoked vs never smoked)
adjsurv_by_smoking <- adjustedsurv(data=surv_dn_landscape,
                        variable="Smoking",
                        ev_time="Derived_Survival",
						event="event",
                        method="direct",
                        outcome_model=fit.coxph,
                        conf_int=TRUE,
						bootstrap=TRUE)
png(file='path/to/supplementary/figures/Smoking_Adjusted_Survival.png',width = 7,height = 7,units="in",res=300)
plot(adjsurv_by_smoking, custom_colors = c("#F3D871BB", "#3DB738BB"), conf_int=TRUE, xlab="Time (weeks)", legend.position="none") + theme(axis.text = element_text(size = 15))
dev.off()

# compare 5-year adjusted survival of all LUAD by de novo group
act <- adjusted_curve_test(adjsurv, to=260)
# multiple testing correction
p.adjust(unlist(sapply(act[1:3], function(x) { x["p_value"] })), method="fdr")

# compare 5-year adjusted survival of NS-LUAD by de novo group
act <- adjusted_curve_test(adjsurv_ns, to=260)
# multiple testing correction
p.adjust(unlist(sapply(act[1:3], function(x) { x["p_value"] })), method="fdr")

# compare 5-year adjusted survival of all LUAD by smoking status
adjusted_curve_test(adjsurv_by_smoking, to=260)