######################################################################### ###### Vulvar ONT HPV integrant methylation ############################# ###### Adapted from Vanessa's github #################################### ###### AUTHOR: MICHELLE NG ############################################## ###### DATE: OCT 2026 ################################################### ######################################################################### ### https://github.com/vanessa-porter/Porter-ONT-HPV-FigureScripts2023/blob/master/scripts/Fig5abcde-hpvMethylation.R ### ---------------------------------------------------------- ### LOAD LIBRARIES ### ---------------------------------------------------------- ## Set Library Path .libPaths(c("/projects/marralab/vporter_prj/R/x86_64-centos7-linux-gnu-library/4.3/", "/projects/marralab/mng_prj/R/x86_64-centos7-linux-gnu-library/4.3/", "/gsc/software/linux-x86_64-centos7/R-4.3.2/")) ## Load Libraries library(reshape2) library(ggplot2) library(tidyr) library(stringr) library(dplyr) library(ggsci) library(cowplot) library(valr) # caused error: no package called valr library(pheatmap) library(cutr) # caused error: no package called cutr --> need this one install.packages("fs", version="2.1", lib="/projects/marralab/mng_prj/R/x86_64-centos7-linux-gnu-library/4.3/") install.packages("devtools", lib="/projects/marralab/mng_prj/R/x86_64-centos7-linux-gnu-library/4.3/") devtools::install_github("moodymudskipper/cutr", lib="/projects/marralab/mng_prj/R/x86_64-centos7-linux-gnu-library/4.3/") ### ---------------------------------------------------------- ### READ IN DATAFRAMES ### ---------------------------------------------------------- # Event locations events <- read.delim("intSitesSummaryAll.txt", header = TRUE) events$id <- paste0(events$sample, "/", events$event) events <- events[events$event != "none",] # Event categories summary <- read.delim("integrationTypes.txt", header = T) summary$id <- paste0(summary$sample, "/", summary$event) ### ---------------------------------------------------------- ### HPV GENE POSITION INFO ### ---------------------------------------------------------- # HPV Gene Info hpvAll <- read.delim("/projects/marralab/vporter_prj/ref/hpv_gff/HPV_all.gff", header = F) hpv <- hpvAll[hpvAll$V3 == "gene",] hpvLengths <- hpvAll[hpvAll$V4 == 1 & hpvAll$V3 == "region",] hpv_info <- strsplit(hpv$V9, ";") hpv_info <- lapply(hpv_info, function(x){ x <- x[grep("Name=", x)] x <- gsub("Name=", "", x) return(x) }) hpv$gene <- unlist(hpv_info) hpv$gene <- gsub("ORF putative ", "", hpv$gene) hpv <- hpv[,-c(2,3,6,7,8,9)] colnames(hpv) <- c("genome", "start", "end","gene") # Add Reading Frame hpv$reading.frame <- ifelse(hpv$gene %in% c("E7","E1","E5","L2"), 1, ifelse(hpv$gene %in% c("E6","E2","L1"), 2, 3)) hpv$position <- hpv$start + (hpv$end - hpv$start)/2 ## Create HPV Position Plots plot1 <- lapply(c("HPV16", "HPV18"), function(i) { ggplot(hpv[hpv$genome == i,], aes(xmin = start, xmax = end, ymin = reading.frame - 0.4, ymax = reading.frame + 0.4)) + geom_rect(aes(fill = gene), colour = "black") + scale_fill_lancet() + geom_text(mapping = aes(x = position, y = reading.frame, label = gene), size = 3, fontface = "bold")+ scale_x_continuous(breaks = c(seq(0,8000,by=1000)), limits = c(0,8000)) + theme(axis.line.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(), axis.title=element_blank(), panel.background=element_blank(), panel.border=element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.background=element_blank(), axis.line.x = element_line(size = 0.5), axis.text.x = element_text(size = 12, colour = "black"), legend.position = "none") }) ### ---------------------------------------------------------- ### INTEGRATION EVENT METHYLATION DATA ### ---------------------------------------------------------- # Set file paths htmcp_path <- "callONTIntegration/output/" # Find event methylation files eventMethFiles <- grep("event_methyl_freq.tsv", list.files(htmcp_path, recursive = TRUE, full.names = TRUE), value = TRUE) # Extract library names from file paths ontLib <- gsub(paste0(htmcp_path, "/|/methylation/event_methyl_freq.tsv"), "", eventMethFiles) # Import ont sites eventMeth <- lapply(eventMethFiles, read.delim, header = TRUE, sep = "\t") names(eventMeth) <- ontLib eventMeth <- dplyr::bind_rows(eventMeth, .id = "library") # Subset HPV data eventMethHPV <- eventMeth[grep("HPV", eventMeth$chromosome),] eventMethHPV <- eventMethHPV[eventMethHPV$chromosome %in% c("HPV16", "HPV18"),] # Find and remap HPV genes eventMethHPV$hpv.gene <- NA eventMethHPV$E6.start <- NA for (i in 1:nrow(eventMethHPV)) { gene <- hpv$gene[hpv$start < eventMethHPV$start[i] & hpv$end > eventMethHPV$start[i] & hpv$genome == eventMethHPV$chromosome[i]] eventMethHPV$hpv.gene[i] <- ifelse(length(gene) > 0, paste(gene, collapse = ","), ifelse(eventMethHPV$start[i] < hpv$start[hpv$genome == eventMethHPV$chromosome[i] & hpv$gene == "E6"] | eventMethHPV$start[i] > hpv$end[hpv$genome == eventMethHPV$chromosome[i] & hpv$gene == "L1"], "LCR", NA)) eventMethHPV$E6.start[i] <- ifelse(eventMethHPV$start[i] < hpv$end[hpv$genome == eventMethHPV$chromosome[i] & hpv$gene == "L1"], eventMethHPV$start[i] - hpv$start[hpv$genome == eventMethHPV$chromosome[i] & hpv$gene == "E6"], -(hpv$start[hpv$genome == eventMethHPV$chromosome[i] & hpv$gene == "E6"] + (hpvLengths$V5[hpvLengths$V1 == eventMethHPV$chromosome[i]] - eventMethHPV$start[i]))) } rm(i, gene) # Reorder factor levels eventMethHPV$chromosome <- factor(eventMethHPV$chromosome, levels = c("HPV16", "HPV18", "HPV45", "HPV82", "HPV52", "HPV58", "HPV33", "HPV31","HPV68", "HPV73","HPV35")) ### ---------------------------------------------------------- ### UNINTEGRATED HPV METHYLATION ### ---------------------------------------------------------- # unintegrated samples unint_samples <- c("GTB_9887","GTB_9931") # Find event methylation files hpvFiles <- grep("hpv_methylation_frequency.bed", list.files(htmcp_path, recursive = TRUE, full.names = TRUE), value = TRUE) # Extract library names from file paths hpvMethNames <- gsub(paste0(htmcp_path, "/|/methylation/hpv_methylation_frequency.bed"), "", hpvFiles) # Import hpvMeth <- lapply(hpvFiles, read.delim, header = FALSE, sep = "\t") names(hpvMeth) <- hpvMethNames hpvMeth <- dplyr::bind_rows(hpvMeth, .id = "library") hpvMeth <- hpvMeth[hpvMeth$V4 == "m",] hpvMeth <- hpvMeth[grep("HPV", hpvMeth$chromosome),] hpvMeth <- hpvMeth[hpvMeth$library %in% unint_samples, c(1:3,13:14,11:12)] colnames(hpvMeth) <- c("library","chromosome","start","methylated","unmethylated","total.calls","perc.methylated") hpvMeth$integration.type <- "no_integration_detected" hpvMeth$id <- hpvMeth$library hpvMeth$event <- "event0" # Find and remap HPV genes hpvMeth$hpv.gene <- NA hpvMeth$E6.start <- NA for (i in 1:nrow(hpvMeth)) { gene <- hpv$gene[hpv$start < hpvMeth$start[i] & hpv$end > hpvMeth$start[i] & hpv$genome == hpvMeth$chromosome[i]] hpvMeth$hpv.gene[i] <- ifelse(length(gene) > 0, paste(gene, collapse = ","), ifelse(hpvMeth$start[i] < hpv$start[hpv$genome == hpvMeth$chromosome[i] & hpv$gene == "E6"] | hpvMeth$start[i] > hpv$end[hpv$genome == hpvMeth$chromosome[i] & hpv$gene == "L1"], "LCR", NA)) hpvMeth$E6.start[i] <- ifelse(hpvMeth$start[i] < hpv$end[hpv$genome == hpvMeth$chromosome[i] & hpv$gene == "L1"], hpvMeth$start[i] - hpv$start[hpv$genome == hpvMeth$chromosome[i] & hpv$gene == "E6"], -(hpv$start[hpv$genome == hpvMeth$chromosome[i] & hpv$gene == "E6"] + (hpvLengths$V5[hpvLengths$V1 == hpvMeth$chromosome[i]] - hpvMeth$start[i]))) } rm(i, gene) ### ---------------------------------------------------------- ### GET THE METHYLATION DATA ADJACENT TO HPV INTEGRATION ### ---------------------------------------------------------- # Get a list of files containing "methylregions.bed" files <- grep("methylregions.bed", list.files("callONTIntegration/output/", recursive = TRUE, full.names = TRUE), value = TRUE) # Extract sample names from file paths names1 <- gsub("callONTIntegration/output//|/methylation/regions|/methylregions.bed", "", files) # Read data from files into a list mrList <- lapply(files, read.delim, header = FALSE) names(mrList) <- names1 # Combine all data into one data frame and add an "id" column mr <- dplyr::bind_rows(mrList, .id = "id") # Add event categories and other information to the methyl regions data frame mr$integration.type <- summary$int.type[match(mr$id, summary$id)] mr$hpv.type <- summary$hpv.type[match(mr$id, summary$id)] # Subset event locations data frame to only the events in the methyl regions data frame eventsSub <- events[events$id %in% mr$id,] # Add 5' and 3' positions to the event locations data frame pos <- events %>% dplyr::filter(event != "none") %>% dplyr::group_by(id) %>% dplyr::summarise(up_pos = min(pos), down_pos = max(pos)) pos$HPVdir <- NA for (i in 1:nrow(pos)) { pos$HPVdir[i] <- ifelse(events$hpv.pos[events$pos == pos$up_pos[i]] > events$hpv.pos[events$pos == pos$down_pos[i]], "forward", "reverse") } rm(i) # Add 5' and 3' positions and HPV direction to the methyl regions data frame mr$up_hpv_pos <- pos$up_pos[match(mr$id, pos$id)] mr$down_hpv_pos <- pos$down_pos[match(mr$id, pos$id)] mr$HPVdir <- pos$HPVdir[match(mr$id, pos$id)] # Add HPV distance to the methyl regions data frame mr$hpv.dist <- ifelse(grepl("5", mr$V12), mr$V2 - mr$up_hpv_pos, mr$V2 - mr$down_hpv_pos) ################# # switch the direction to match the translation direction mrTln <- mr mrTln$hpv.dist <- dplyr::case_when( mr$HPVdir == "forward" & mr$V12 == "5_upstream" ~ -(abs(mr$hpv.dist)), mr$HPVdir == "reverse" & mr$V12 == "5_upstream" ~ abs(mr$hpv.dist), mr$HPVdir == "forward" & mr$V12 == "3_downstream" ~ abs(mr$hpv.dist), mr$HPVdir == "reverse" & mr$V12 == "3_downstream" ~ -(abs(mr$hpv.dist)), mr$HPVdir == "forward" & mr$V12 == "5_downstream" ~ -(abs(mr$hpv.dist)), mr$HPVdir == "reverse" & mr$V12 == "5_downstream" ~ abs(mr$hpv.dist), mr$HPVdir == "forward" & mr$V12 == "3_upstream" ~ abs(mr$hpv.dist), mr$HPVdir == "reverse" & mr$V12 == "3_upstream" ~ -(abs(mr$hpv.dist)), TRUE ~ mr$hpv.dist ) ################# ## BREAK ######## ## can't install cutr package for smart_cut() ## try installing locally and exporting workspace ## bin the methlation into 500 bp bins mrTln$bins <- smart_cut(mrTln$hpv.dist, seq(from = -5000, to = 5000, by = 500), labels = ~paste0(.y[1],':',.y[2]-1), simplify = FALSE) allbins <- smart_cut(seq(from = -5000, to = 5000, by = 500), seq(from = -5000, to = 5000, by = 500), labels = ~paste0(.y[1],':',.y[2]-1), simplify = FALSE) # make an average value for each bin mrAve <- mrTln %>% dplyr::group_by(id,bins,V12, integration.type, hpv.type) %>% dplyr::summarise(average_methyl = mean(V8)) mrAve$bins <- factor(mrAve$bins, levels = levels(allbins))