SherlockLung-EvolutionaryTrajectory-Analysis / Analyses / MDM2_amp_comaprison.R
MDM2_amp_comaprison.R
Raw 
#MDM2 question:
# it is not identified by timing model as not gained in enough samples
# however, samples in which it is gained, it is highly amplified
# is it...?
# a) amplified in more NSD-A samples than others
# b) amplified to a higher level in NSD-A
# c) both

#ANSWER: b.

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

NSDA_allsegs = read.table("path/to/data/SherlockLung_LUAD_2025_DN2/SherlockLung_LUAD_2025_DN2_allsegments.txt", header=T, stringsAsFactors = F, sep='\t')
NSDB_allsegs = read.table("path/to/data/Data/SherlockLung_LUAD_2025_DN3/SherlockLung_LUAD_2025_DN3_allsegments.txt", header=T, stringsAsFactors = F, sep='\t')
SD_allsegs = read.table("path/to/data/SherlockLung_LUAD_2025_DN1/SherlockLung_LUAD_2025_DN1_allsegments.txt", header=T, stringsAsFactors = F, sep='\t')

NSDA_MDM2_gains = NSDA_allsegs[NSDA_allsegs$chr == 12 & NSDA_allsegs$endpos >= start(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDA_allsegs$startpos <= end(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDA_allsegs$CNA %in% c("cGain","sGain","cAmp","sAmp"),]
NSDB_MDM2_gains = NSDB_allsegs[NSDB_allsegs$chr == 12 & NSDB_allsegs$endpos >= start(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDB_allsegs$startpos <= end(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDB_allsegs$CNA %in% c("cGain","sGain","cAmp","sAmp"),]
SD_MDM2_gains = SD_allsegs[SD_allsegs$chr == 12 & SD_allsegs$endpos >= start(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & SD_allsegs$startpos <= end(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & SD_allsegs$CNA %in% c("cGain","sGain","cAmp","sAmp"),]

length(unique(NSDA_MDM2_gains$Tumour_Name)) / length(sample_groups[[1]]) # 36%
length(unique(NSDB_MDM2_gains$Tumour_Name)) / length(sample_groups[[2]]) # 31%
length(unique(SD_MDM2_gains$Tumour_Name)) / length(sample_groups[[3]]) # 24%

p.adjust(c(fisher.test(cbind(c(length(unique(NSDA_MDM2_gains$Tumour_Name)), length(sample_groups[[1]])-length(unique(NSDA_MDM2_gains$Tumour_Name))), c(length(unique(NSDB_MDM2_gains$Tumour_Name)), length(sample_groups[[2]])-length(unique(NSDB_MDM2_gains$Tumour_Name)))))$p.value, fisher.test(cbind(c(length(unique(NSDA_MDM2_gains$Tumour_Name)), length(sample_groups[[1]])-length(unique(NSDA_MDM2_gains$Tumour_Name))), c(length(unique(SD_MDM2_gains$Tumour_Name)), length(sample_groups[[3]])-length(unique(SD_MDM2_gains$Tumour_Name)))))$p.value, fisher.test(cbind(c(length(unique(NSDB_MDM2_gains$Tumour_Name)), length(sample_groups[[2]])-length(unique(NSDB_MDM2_gains$Tumour_Name))), c(length(unique(SD_MDM2_gains$Tumour_Name)), length(sample_groups[[3]])-length(unique(SD_MDM2_gains$Tumour_Name)))))$p.value), method="fdr")
# Number of gains - FDRs:
# NSD-A vs NSD-B 0.36855123
# NSD-A vs SD    0.05426258
# NSD-B vs SD    0.24159362

MDM2_start = start(all_gene_gr[all_gene_gr$symbol == "MDM2"])
MDM2_end = end(all_gene_gr[all_gene_gr$symbol == "MDM2"])
MDM2_size = (MDM2_end - MDM2_start) + 1

NSDA_MDM2_gains$startpos = sapply(NSDA_MDM2_gains$startpos, max, MDM2_start)
NSDA_MDM2_gains$endpos = sapply(NSDA_MDM2_gains$endpos, min, MDM2_end)

NSDB_MDM2_gains$startpos = sapply(NSDB_MDM2_gains$startpos, max, MDM2_start)
NSDB_MDM2_gains$endpos = sapply(NSDB_MDM2_gains$endpos, min, MDM2_end)

SD_MDM2_gains$startpos = sapply(SD_MDM2_gains$startpos, max, MDM2_start)
SD_MDM2_gains$endpos = sapply(SD_MDM2_gains$endpos, min, MDM2_end)

all(NSDA_MDM2_gains$startpos <= MDM2_start)
all(NSDB_MDM2_gains$startpos <= MDM2_start)
all(SD_MDM2_gains$startpos <= MDM2_start)
all(NSDA_MDM2_gains$endpos >= MDM2_end)
all(NSDB_MDM2_gains$endpos >= MDM2_end)
all(SD_MDM2_gains$endpos >= MDM2_end)

get_mean_majCNs = function(gains, target_size) {
  sapply(unique(gains$Tumour_Name), function(bc) {
    gainy = gains[gains$Tumour_Name == bc,]
	mean_majCN = mean(gainy[,"nMaj1_A"])
	gain_size = (max(gainy$endpos) - min(gainy$startpos)) + 1
	size_multiplier = gain_size / target_size
	mean_majCN * size_multiplier
  })
}

NSDA_MDM2_majCN = get_mean_majCNs(NSDA_MDM2_gains, MDM2_size)
NSDB_MDM2_majCN = get_mean_majCNs(NSDB_MDM2_gains, MDM2_size)
SD_MDM2_majCN = get_mean_majCNs(SD_MDM2_gains, MDM2_size)

MDM2_gains_majCN_df = melt(list(`NSD-A`=NSDA_MDM2_majCN, `NSD-B`=NSDB_MDM2_majCN, `SD`=SD_MDM2_majCN))
colnames(MDM2_gains_majCN_df) = c("Major copy number in gains covering MDM2", "Trajectory")

p <- ggplot(MDM2_gains_majCN_df, aes(x=Trajectory, y=`Major copy number in gains covering MDM2`, fill=Trajectory), show.legend=FALSE) +
    geom_boxplot(outlier.shape = NA) +
    theme_bw(base_family = "Roboto Condensed") +
    theme(legend.position="none", panel.border=element_rect(linewidth = 1.25, fill=NA), axis.ticks.x=element_blank(), axis.ticks.y=element_blank(), text = element_text(size=20), axis.text.x = element_text(size=20), axis.text.y = element_text(size=18)) +
    geom_jitter(shape=16, position=position_jitter(0.25)) +
	scale_fill_manual(values=c("#8adbf590","#a371cb90","#5a5a6990")) +
    xlab('')
p
ggsave(file='Sherlock/New_Histology/2025/Figures/FigS3/MDM2_gains_majorCN_boxplot.png',plot = p,width = 7,height = 7,device = "png")

p.adjust(c(wilcox.test(NSDA_MDM2_majCN, NSDB_MDM2_majCN)$p.value, wilcox.test(NSDA_MDM2_majCN, SD_MDM2_majCN)$p.value, wilcox.test(NSDB_MDM2_majCN, SD_MDM2_majCN)$p.value), method = "fdr")
# Major copy number of gains - FDRs
# NSD-A vs NSD-B 0.0007495787
# NSD-A vs SD    0.0001410188
# NSD-B vs SD    0.5403193828



NSDA_MDM2_losses = NSDA_allsegs[NSDA_allsegs$chr == 12 & NSDA_allsegs$endpos >= start(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDA_allsegs$startpos <= end(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDA_allsegs$CNA %in% c("cLOH", "sLOH"),]
NSDB_MDM2_losses = NSDB_allsegs[NSDB_allsegs$chr == 12 & NSDB_allsegs$endpos >= start(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDB_allsegs$startpos <= end(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & NSDB_allsegs$CNA %in% c("cLOH","sLOH"),]
SD_MDM2_losses = SD_allsegs[SD_allsegs$chr == 12 & SD_allsegs$endpos >= start(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & SD_allsegs$startpos <= end(all_gene_gr[all_gene_gr$symbol == "MDM2"]) & SD_allsegs$CNA %in% c("cLOH","sLOH"),]

mdm2_gain_bc = c(NSDA_MDM2_gains$Tumour_Name, NSDB_MDM2_gains$Tumour_Name, SD_MDM2_gains$Tumour_Name)
mdm2_loss_bc = c(NSDA_MDM2_losses$Tumour_Name, NSDB_MDM2_losses$Tumour_Name, SD_MDM2_losses$Tumour_Name)

fisher.test(rbind(
c(length(
   setdiff(
    c(unique(NSDA_allsegs$Tumour_Name), unique(NSDB_allsegs$Tumour_Name), unique(SD_allsegs$Tumour_Name)),
	c(mdm2_gain_bc, mdm2_loss_bc)
   )
  )+1,
  length(setdiff(mdm2_gain_bc, mdm2_loss_bc))
 ),
c(length(setdiff(mdm2_loss_bc, mdm2_gain_bc)), length(intersect(mdm2_gain_bc, mdm2_loss_bc)))
))
# p = 0.09


# ecDNA analysis

ecDNA_mdm2 = read.table("path/to/data/ecDNA-MDM2.txt", header=T, stringsAsFactors = F, sep='\t')
ecDNA_mdm2_contTab = sapply(sample_groups, function(sg) { table(ecDNA_mdm2[ecDNA_mdm2$Tumor_Barcode %in% gsub("01_01", "01.01", sg),"MDM2_class"]) })
ecDNA_mdm2_amp_vs_not = rbind(ecDNA_mdm2_contTab[1,], apply(ecDNA_mdm2_contTab[2:3,], 2, sum))
any_mdm2_amp_vs_not = rbind(apply(ecDNA_mdm2_contTab[1:2,], 2, sum), ecDNA_mdm2_contTab[3,])

p.adjust(c(fisher.test(ecDNA_mdm2_amp_vs_not[,1:2])$p.value, fisher.test(ecDNA_mdm2_amp_vs_not[,c(1,3)])$p.value, fisher.test(ecDNA_mdm2_amp_vs_not[,c(2,3)])$p.value), method="fdr")
p.adjust(c(fisher.test(any_mdm2_amp_vs_not[,1:2])$p.value, fisher.test(any_mdm2_amp_vs_not[,c(1,3)])$p.value, fisher.test(any_mdm2_amp_vs_not[,c(2,3)])$p.value), method="fdr")
Powered by GitFront