library(ggplot2)
library(hrbrthemes)
library(viridis)
library(ggrepel)

mergedseg = read.table("path/to/data/SherlockLung_AllHQ_LUAD_2025_mergedseg_G3.txt", header=T, stringsAsFactors=F, sep='\t')

all_s_s4 = landscape[landscape$Smoking_SBS4 == "Smoker_SBS4","Tumor_Barcode"]
all_ns_ns4 = landscape[landscape$Smoking_SBS4 == "NonSmoker_NonSBS4","Tumor_Barcode"]

s_s4 = intersect(all_s_s4, unique(mergedseg$Tumour_Name))
ns_ns4 = intersect(all_ns_ns4, unique(mergedseg$Tumour_Name))

all_events = unique(mergedseg$ID)
ns_events = unique(mergedseg[mergedseg$Tumour_Name %in% ns_ns4,"ID"])
s_events = unique(mergedseg[mergedseg$Tumour_Name %in% s_s4,"ID"])

# MAKE SURE sample_groups is ordered such that NSD-A is 1st, NSD-B is 2nd, and SD is 3rd item

sl_ns = intersect(ns_ns4, sample_groups[[3]])
nsl_ns = intersect(ns_ns4, c(sample_groups[[1]],sample_groups[[2]]))

sl_s = intersect(s_s4, sample_groups[[3]])
nsl_s = intersect(s_s4, c(sample_groups[[1]],sample_groups[[2]]))

# get data for volcano plots
get_vp_data = function(events, sample_set_1, sample_set_2, mergedseg) {
  contTabs = lapply(events, function(ev) {
    samples_with_event = mergedseg[mergedseg$ID == ev,"Tumour_Name"]
    num_ss1_with_event = length(intersect(sample_set_1, samples_with_event))
    num_ss1_without_event = length(sample_set_1) - num_ss1_with_event
    num_ss2_with_event = length(intersect(sample_set_2, samples_with_event))
    num_ss2_without_event = length(sample_set_2) - num_ss2_with_event
    rbind(
      c(num_ss1_without_event, num_ss1_with_event),
	  c(num_ss2_without_event, num_ss2_with_event)
    )
  })
  names(contTabs) = events
  
  ORs = sapply(contTabs, function(x) { fisher.test(x)$estimate[['odds ratio']] })
  ps = sapply(contTabs, function(x) { fisher.test(x)$p.value })
  fdrs = p.adjust(ps, method="fdr")

  l2or = log2(ORs)

  hiRange = round(max(l2or[-which(is.infinite(l2or))]))
  loRange = round(min(l2or[-which(is.infinite(l2or))]))
  range_limit = max(abs(hiRange), abs(loRange))
  exl_ss2_val = range_limit + 1
  exl_ss1_val = -exl_ss2_val

  l2or[which(l2or == Inf)] = exl_ss2_val
  l2or[which(l2or == -Inf)] = exl_ss1_val

  event_names = apply(unique(mergedseg[mergedseg$Tumour_Name %in% c(sample_set_1, sample_set_2),c("CNA","ID")]), 1, function(x) { x[1] = substr(x[1], 2, nchar(x[1])); x[2] = gsub("dMut_", "", x[2]); paste(x[1], x[2], sep='_')})

  vp_data = as.data.frame(cbind(`event`=event_names, `log2OR`=l2or, `FDR`=fdrs, `highlight`=rep("none", length(event_names))))
  vp_data[,'log2OR'] = as.numeric(vp_data[,'log2OR'])
  vp_data[,'FDR'] = as.numeric(vp_data[,'FDR'])

  vp_data$highlight[vp_data$log2OR > log2(3/2) & vp_data$FDR < 0.05] <- "UP"
  vp_data$highlight[vp_data$log2OR < log2(2/3) & vp_data$FDR < 0.05] <- "DOWN"

  vp_data$vp_label <- NA
  vp_data$vp_label[vp_data$highlight != "none"] <- vp_data$event[vp_data$highlight != "none"]
  vp_data$event_type <- NA
  vp_data$event_type[grep("^Mut_", vp_data$event)] <- "Mut"
  vp_data$event_type[grep("^LOH_", vp_data$event)] <- "Loss"
  vp_data$event_type[grep("^Gain_", vp_data$event)] <- "Gain"
  vp_data$event_type[grep("^HD_", vp_data$event)] <- "HD"

  vp_data$vpd$vp_label = event_ids[gsub("dGain_|dLOH_|dHD_", "", paste0("d", vp_data$vp_label))]

  vp_data
}

# get grid for creating volcano plot
get_vp_grid = function(vp_data, ss1_name, ss2_name, label=FALSE, down_colour="#8BD6FFBB", up_colour="#A944C3BB", nsig_colour="#FFFFFF88") {
    extent = round(max(abs(vp_data$log2OR))+0.5)
	xbreaks = -extent:extent
	xticklabs = as.character(xbreaks)
    xticklabs[1] = paste0("exclusive in\n", ss1_name)
	xticklabs[length(xticklabs)] = paste0("exclusive in\n", ss2_name)
    p <- ggplot(data=vp_data, aes(x=log2OR, y=-log10(FDR), fill=highlight, label=vp_label, shape=event_type)) + 
        scale_x_continuous(breaks=xbreaks, minor_breaks=NULL, labels=xticklabs, limits=c(-extent, extent), expand=c(0, 0)) +
        scale_y_continuous(expand=c(0,0,0,1)) +
        theme_bw(base_family="Roboto Condensed") +
        geom_vline(xintercept=c(log2(2/3), log2(3/2)), col="green", linetype="dashed") +
        geom_hline(yintercept=-log10(0.05), col="green", linetype="dashed") +
        coord_cartesian(clip = 'off') +
        scale_fill_manual(values=c(down_colour, nsig_colour, up_colour)) +
        scale_shape_manual(values=c(24,4,25,21)) +
        labs(x="log2(odds ratio)") +
        geom_point(size=2.5)
    if (label) {
        p = p + geom_text_repel(max.overlaps=Inf, size=3)
    }
    g <- grid::grid.force(ggplotGrob(p))
    pos <- g$layout[g$layout$name == "panel", 1:5]
    g <- gtable::gtable_add_grob(x=ggplotGrob(p),  grobs=grid::getGrob(g, "point", grep=TRUE), t=pos$t, l=pos$l, z=max(g$layout$z) + 1L, clip=FALSE)
    g
}

ms = mergedseg
msl = read.table("path/to/data/SherlockLung_AllHQ_LUAD_2025_LOH_mergedsegs.txt", header=T, stringsAsFactors = F, sep='\t')
msl$Tumour_Name = gsub("01_01", "01.01", msl$Tumour_Name)
msl_compare = apply(msl, 2, function(x) { gsub("^ *", "", x) })
ms_compare = apply(ms[,1:10], 2, function(x) { gsub("^ *", "", x) })
msl_in_ms = apply(msl_compare, 1, function(x) { any(apply(ms_compare, 1, function(y) { identical(x, y) })) })
sLOH_to_add = msl[!msl_in_ms,]
sLOH_to_add = cbind(sLOH_to_add, `ID`=apply(sLOH_to_add, 1, function(x) { gsub(" *", "", paste(paste0("chr", x["chr"]), x["startpos"], x["endpos"], sep="_")) }), `no.chr.bearing.mut`=as.character(NA))
mergedseg_extended = rbind(ms, sLOH_to_add)


# All LUAD: NSD vs SD

vpd = get_vp_data(all_events, c(sample_groups[[1]], sample_groups[[2]]), sample_groups[[3]], mergedseg_extended)
fn = "path/to/figures/Fig2/LUAD_all_samples_event_trajectory_preferences_redblack_mse.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", label=FALSE, down_colour="#c25050BB", up_colour="#5a5a69BB"))
dev.off()

# commands to aid with labelling
vpd = vpd[order(vpd$log2OR, decreasing=T),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl <- vpd$vp_label[i]; vph <- vpd$highlight[i]; c(vpl, vph) }))
t(t(conversions[conversions[,3] == "UP",2]))
vpd = vpd[order(vpd$log2OR),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "DOWN",2]))

# output table
tab_s1 = vpd[,1:3]
tab_s1[,1] = event_ids[gsub("LOH_|Gain_|HD_", "", gsub("Mut_", "dMut_", tab_s1[,1]))]
colnames(tab_s1) = c("Event", "log2(OR)", "FDR")
write.table(tab_s1, file="path/to/figures/Fig2/All_LUAD_NSD_vs_SD_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)


# SBS4-negative NS: NSD vs SD

vpd = get_vp_data(ns_events, nsl_ns, sl_ns, mergedseg_extended)
fn = "path/to/figures/Fig2/LUAD_NS_event_trajectory_preferences_redblack_mse.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", label=FALSE, down_colour="#c25050BB", up_colour="#5a5a69BB"))
dev.off()

# commands to aid with labelling
vpd = vpd[order(vpd$log2OR, decreasing=T),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "UP",2]))
vpd = vpd[order(vpd$log2OR),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "DOWN",2]))

# output table
tab_s2 = vpd[,1:3]
tab_s2[,1] = event_ids[gsub("LOH_|Gain_|HD_", "", gsub("Mut_", "dMut_", tab_s2[,1]))]
colnames(tab_s2) = c("Event", "log2(OR)", "FDR")
write.table(tab_s2, file="path/to/figures/Fig2/SBS4negNS_NSD_vs_SD_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)


# SBS4-positive smokers: NSD vs SD

vpd <- get_vp_data(s_events, nsl_s, sl_s, mergedseg_extended)
fn = "path/to/figures/Fig2/S-LUAD_event_trajectory_preferences_redblack_mse.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_vp_grid(get_vp_data(s_events, nsl_s, sl_s, mergedseg_extended), "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", label=FALSE, down_colour="#c25050BB", up_colour="#5a5a69BB"))
dev.off()

# commands to aid with labelling
vpd = vpd[order(vpd$log2OR, decreasing=T),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "UP",2]))
vpd = vpd[order(vpd$log2OR),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "DOWN",2]))

# output table
tab_s3 = vpd[,1:3]
tab_s3[,1] = event_ids[gsub("LOH_|Gain_|HD_", "", gsub("Mut_", "dMut_", tab_s3[,1]))]
colnames(tab_s3) = c("Event", "log2(OR)", "FDR")
write.table(tab_s3, file="path/to/figures/Fig2/SBS4posS_NSD_vs_SD_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)



# All LUAD: NSD-A vs NSD-B

nsd_events = unique(mergedseg_extended[mergedseg_extended$Tumour_Name %in% c(sample_groups[[1]], sample_groups[[2]]),"ID"])

vpd <- get_vp_data(nsd_events, sample_groups[[1]], sample_groups[[2]], mergedseg_extended)
fn = "path/to/figures/Fig2/LUAD_NS_AvB_event_trajectory_preferences_rothcols_mse.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_vp_grid(vpd, "NS-LUAD\ndominant A", "NS-LUAD\ndominant B", label=FALSE, down_colour="#8adbf5BB", up_colour="#9c66c7BB"))
dev.off()

# commands to aid with labelling
vpd = vpd[order(vpd$log2OR, decreasing=T),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "UP",2]))
vpd = vpd[order(vpd$log2OR),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "DOWN",2]))

# output table
tab_s4 = vpd[,1:3]
tab_s4[,1] = event_ids[gsub("LOH_|Gain_|HD_", "", gsub("Mut_", "dMut_", tab_s4[,1]))]
colnames(tab_s4) = c("Event", "log2(OR)", "FDR")
write.table(tab_s4, file="path/to/figures/Fig2/NSDA_vs_NSDB_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)


# All LUAD: Smokers vs NS

vpd = get_vp_data(all_events, c(nsl_ns, sl_ns), c(sl_s, nsl_s), mergedseg)
fn = "path/to/figures/Fig2/LUAD_NS_S_event_trajectory_preferences.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_vp_grid(vpd, "NS-LUAD", "S-LUAD", label=FALSE, down_colour="#3DB738BB", up_colour="#F3D871BB"))
dev.off()

# commands to aid with labelling
vpd = vpd[order(vpd$log2OR, decreasing=T),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "UP",2]))
vpd = vpd[order(vpd$log2OR),]
to_include = which(!is.na(vpd$vp_label)); conversions = t(sapply(to_include, function(i) { vpl = vpd$vp_label[i]; vpe <- gsub("^[^_]+_", "", vpl); eid <- event_ids[vpe]; vph <- vpd$highlight[i]; c(vpl, eid, vph) }))
t(t(conversions[conversions[,3] == "DOWN",2]))

# output table
tab_s5 = vpd[,1:3]
tab_s5[,1] = event_ids[gsub("LOH_|Gain_|HD_", "", gsub("Mut_", "dMut_", tab_s5[,1]))]
colnames(tab_s5) = c("Event", "log2(OR)", "FDR")
write.table(tab_s5, file="path/to/figures/Fig2/NS-LUAD_S-LUAD_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)



# Signatures for supplementary figure

sample_set_1 = c(sample_groups[[1]], sample_groups[[2]])
sample_set_2 = sample_groups[[3]]

id_sigs  = read.table("Sherlock/New_Histology/2025/Data/Signatures/Activities_Sherlock_ID83_10k_normalization_COSMIC_decomposition_v3.4.txt", header=T, stringsAsFactors =F, sep='\t')
cn_sigs  = read.table("Sherlock/New_Histology/2025/Data/Signatures/Activities_Sherlock_CN_COSMIC_decomposition_v3.4.txt", header=T, stringsAsFactors =F, sep='\t')
sbs_sigs = read.table("Sherlock/New_Histology/2025/Data/Signatures/Activities_Sherlock_SBS288_10k_normalization_COSMIC_decomposition_v3.4.txt", header=T, stringsAsFactors =F, sep='\t')
dbs_sigs = read.table("Sherlock/New_Histology/2025/Data/Signatures/Activities_Sherlock_DBS78_10k_normalization_COSMIC_decomposition_v3.4.txt", header=T, stringsAsFactors =F, sep='\t')
sv_sigs  = read.table("Sherlock/New_Histology/2025/Data/Signatures/Activities_Sherlock_SV_COSMIC_decomposition_v3.4.txt", header=T, stringsAsFactors =F, sep='\t')

sigs_list = list(sbs_sigs, dbs_sigs, id_sigs, cn_sigs, sv_sigs)

get_signature_vp_data = function(sample_set_1, sample_set_2) {

  contTabs = lapply(sigs_list, function(sigs) {
  
    rownames(sigs) = sigs[,1]
    sigs = sigs[,-1]
    sig_presence = sigs
    sig_presence[sig_presence > 0] = 1
	
    sig_type_contTabs = lapply(1:ncol(sig_presence), function(sig_i) {
      num_ss1_with_event = sum(sig_presence[sample_set_1,sig_i])
      num_ss1_without_event = length(sample_set_1) - num_ss1_with_event
      num_ss2_with_event = sum(sig_presence[sample_set_2,sig_i])
      num_ss2_without_event = length(sample_set_2) - num_ss2_with_event
      contTab = rbind(
        c(num_ss1_without_event, num_ss1_with_event),
        c(num_ss2_without_event, num_ss2_with_event)
      )
	  colnames(contTab) = c("Absent","Present")
	  rownames(contTab) = c("Set1","Set2")
	  contTab
	})
	names(sig_type_contTabs) = colnames(sig_presence)
	sig_type_contTabs
  })
  names(contTabs) = c("SBS","DBS","ID","CN","SV")
  
  all_contTabs = Reduce(c, contTabs)
  
  ORs = sapply(all_contTabs, function(x) { fisher.test(x)$estimate[['odds ratio']] })
  ps = sapply(all_contTabs, function(x) { fisher.test(x)$p.value })
  fdrs = p.adjust(ps, method="fdr")

  l2or = log2(ORs)

  hiRange = round(max(l2or[-which(is.infinite(l2or))]))
  loRange = round(min(l2or[-which(is.infinite(l2or))]))
  range_limit = max(abs(hiRange), abs(loRange))
  exl_ss2_val = range_limit + 1
  exl_ss1_val = -exl_ss2_val

  l2or[which(l2or == Inf)] = exl_ss2_val
  l2or[which(l2or == -Inf)] = exl_ss1_val

  event_names = unlist(lapply(contTabs, names))

  vp_data = as.data.frame(cbind(`event`=event_names, `log2OR`=l2or, `FDR`=fdrs, `highlight`=rep("none", length(event_names))))
  vp_data[,'log2OR'] = as.numeric(vp_data[,'log2OR'])
  vp_data[,'FDR'] = as.numeric(vp_data[,'FDR'])

  vp_data$highlight[vp_data$log2OR > log2(3/2) & vp_data$FDR < 0.05] <- "UP"
  vp_data$highlight[vp_data$log2OR < log2(2/3) & vp_data$FDR < 0.05] <- "DOWN"

  vp_data$vp_label <- NA
  vp_data$vp_label[vp_data$highlight != "none"] <- vp_data$event[vp_data$highlight != "none"]
  vp_data$event_type <- gsub('[0-9]+.*', '', vp_data$event)

  vp_data
}

get_signature_vp_grid = function(vp_data, ss1_name, ss2_name, label=FALSE, down_colour="#8BD6FFBB", up_colour="#A944C3BB", nsig_colour="#FFFFFF88") {
    extent = round(max(abs(vp_data$log2OR))+0.5)
	xbreaks = -extent:extent
	xticklabs = as.character(xbreaks)
    xticklabs[1] = paste0("exclusive in\n", ss1_name)
	xticklabs[length(xticklabs)] = paste0("exclusive in\n", ss2_name)
    p <- ggplot(data=vp_data, aes(x=log2OR, y=-log10(FDR), fill=highlight, label=vp_label, shape=event_type)) + 
	    scale_x_continuous(breaks=xbreaks, minor_breaks=NULL, labels=xticklabs, limits=c(-extent, extent), expand=c(0, 0)) +
        scale_y_continuous(expand=c(0,0,0,1)) +
        theme_bw(base_family="Roboto Condensed") +
        geom_vline(xintercept=c(log2(2/3), log2(3/2)), col="green", linetype="dashed") +
        geom_hline(yintercept=-log10(0.05), col="green", linetype="dashed") +
        coord_cartesian(clip = 'off') +
        scale_fill_manual(values=c(down_colour, nsig_colour, up_colour)) +
        scale_shape_manual(values=c(24,22,23,21,25)) +
        labs(x="log2(odds ratio)") +
        geom_point(size=2.5)
    if (label) {
        p = p + geom_text_repel(max.overlaps=Inf, size=3)
    }
    g <- grid::grid.force(ggplotGrob(p))
    pos <- g$layout[g$layout$name == "panel", 1:5]
    g <- gtable::gtable_add_grob(x=ggplotGrob(p),  grobs=grid::getGrob(g, "point", grep=TRUE), t=pos$t, l=pos$l, z=max(g$layout$z) + 1L, clip=FALSE)
    g
}

vpd = get_signature_vp_data(c(sample_groups[[1]], sample_groups[[2]]), sample_groups[[3]])
grid::grid.draw(get_signature_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", label=TRUE, down_colour="#c25050BB", up_colour="#5a5a69BB"))
fn = "path/to/supplementary/figures/LUAD_all_samples_signature_trajectory_preferences.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_signature_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", down_colour="#c25050BB", up_colour="#5a5a69BB"))
dev.off()

vpd = vpd[order(vpd$log2OR),]
tab_s6 = vpd[,1:3]
colnames(tab_s6) = c("Event", "log2(OR)", "FDR")
write.table(tab_s6, file="path/to/supplementary/figures/AllLUAD_signature_trajectory_preference_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)

vpd = get_signature_vp_data(nsl_ns, sl_ns)
grid::grid.draw(get_signature_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", label=TRUE, down_colour="#c25050BB", up_colour="#5a5a69BB"))
fn = "path/to/supplementary/figures/LUAD_NS_signature_trajectory_preferences.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_signature_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", down_colour="#c25050BB", up_colour="#5a5a69BB"))
dev.off()

vpd = vpd[order(vpd$log2OR),]
tab_s7 = vpd[,1:3]
colnames(tab_s7) = c("Event", "log2(OR)", "FDR")
write.table(tab_s7, file="path/to/supplementary/figures/SBS4-NS_LUAD_signature_trajectory_preference_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)

vpd = get_signature_vp_data(nsl_s, sl_s)
grid::grid.draw(get_signature_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", label=TRUE, down_colour="#c25050BB", up_colour="#5a5a69BB"))
fn = "path/to/supplementary/figures/LUAD_S_signature_trajectory_preferences.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_signature_vp_grid(vpd, "NS-LUAD\ndominant\ntrajectories", "S-LUAD\ndominant\ntrajectory", down_colour="#c25050BB", up_colour="#5a5a69BB"))
dev.off()

vpd = vpd[order(vpd$log2OR),]
tab_s8 = vpd[,1:3]
colnames(tab_s8) = c("Event", "log2(OR)", "FDR")
write.table(tab_s8, file="path/to/supplementary/figures/SBS4+S_LUAD_signature_trajectory_preference_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)

vpd = get_signature_vp_data(sample_groups[[1]], sample_groups[[2]])
grid::grid.draw(get_signature_vp_grid(vpd, "NSD-A", "NSD-B", label=TRUE))
fn = "path/to/supplementary/figures/LUAD_NSD_A_B_signature_trajectory_preferences.png"
png(fn, width=2200, height=1500, res=300)
grid::grid.draw(get_signature_vp_grid(vpd, "NS-LUAD\ndominant A", "NS-LUAD\ndominant B"))
dev.off()

vpd = vpd[order(vpd$log2OR),]
tab_s9 = vpd[,1:3]
colnames(tab_s9) = c("Event", "log2(OR)", "FDR")
write.table(tab_s9, file="path/to/supplementary/figures/NSDA_vs_NSDB_signature_trajectory_preference_log2OR_FDR.txt", quote=FALSE, sep='\t', col.names=TRUE, row.names=FALSE)