pkgs <- c("fs", "futile.logger", "configr", "stringr", "ggpubr", "ggthemes",
"jhtools", "glue", "ggsci", "cowplot", "patchwork", "tidyverse",
"circlize", "grid", "gplots", "ggrepel", "ComplexHeatmap",
"GenomicRanges", "jhuanglabRNAseq", "ggh4x", "gghalves"
)
for (pkg in pkgs){
suppressPackageStartupMessages(library(pkg, character.only = T))
}
setwd("/cluster/home/jhuang/projects/mm/analysis/meta/human/manuscript")
# ======================== circos plot of overall genomic alterations =======================
band_data <- read.delim("rds/cytoBand_hg38.txt", header = F)
colnames(band_data) <- c("chr", "start", "end", "band", "stain_value")
band_data <- band_data[nchar(band_data$chr) <= 5, ]
m_cnv_band <- read_rds("rds/merged_cnv_band_circos.rds")
m_cnv_band_fil <- m_cnv_band[!grepl("C$", m_cnv_band$sample_id), ]
sample_count <- length(unique(m_cnv_band_fil$sample_id))
stat_cnv_band <- m_cnv_band_fil %>% group_by(chr, band, bivalue) %>% summarise(prop = n() / sample_count)
stat_cnv_band <- left_join(stat_cnv_band, band_data[, 1:4])
gain_df <- stat_cnv_band[stat_cnv_band$bivalue %in% 1, ]
loss_df <- stat_cnv_band[stat_cnv_band$bivalue %in% -1, ]
gain_df$value <- gain_df$prop
loss_df$value <- loss_df$prop
gain_df <- gain_df[, c("chr", "start", "end", "value")]
loss_df <- loss_df[, c("chr", "start", "end", "value")]
# load mut data (may alter later)
merged_mutations <- read_rds("rds/merged_mutations_circos.rds")
merged_mutations <- merged_mutations[!grepl("C$", merged_mutations$sample_id), ]
sample_count2 <- length(unique(merged_mutations$sample_id))
stat_mut <- merged_mutations %>% group_by(chr, start, end) %>% summarise(value = n() / sample_count2)
mutation_df <- as.data.frame(stat_mut)
mutation_df <- mutation_df[mutation_df$value %notin% min(mutation_df$value), ]
mutation_df <- mutation_df[mutation_df$value %notin% min(mutation_df$value), ]
col_fun <- colorRamp2(
c(min(mutation_df$value, na.rm = TRUE),
median(mutation_df$value, na.rm = TRUE),
max(mutation_df$value, na.rm = TRUE)),
c("steelblue", "yellow2", "tomato")
)
# load label data
label_df <- read_rds("rds/white_list_label_circos.rds") %>%
dplyr::filter(!chr %in% c("chrX"))
legend_obj <- ComplexHeatmap::Legend(
col_fun = col_fun,
title = "Mutation Frequency",
direction = "horizontal",
legend_width = unit(4, "cm"),
title_position = "topcenter"
)
# plot
pdf("figure2_circos_plot_hg38.pdf", width = 8, height = 8)
circos.clear()
circos.par(start.degree = 90, "canvas.ylim" = c(-1.25, 1.25))
circos.initializeWithIdeogram(
species = "hg38",
chromosome.index = paste0("chr", c(1:22)),
plotType = "none"
)
circos.genomicLabels(
label_df,
labels.column = 4,
side = "outside",
connection_height = mm_h(5),
line_col = "grey30",
cex = 0.5,
labels_height = mm_h(3)
)
circos.clear()
par(new = TRUE)
circos.par(
start.degree = 90,
gap.degree = 1,
track.margin = c(0.01, 0.01),
canvas.ylim = c(-1.5, 1.5)
)
# ======= initialization =======
circos.initializeWithIdeogram(
species = "hg38",
chromosome.index = paste0("chr", c(1:22)),
plotType = c("ideogram", "labels"),
ideogram.height = 0.04,
track.height = 0.02
)
# ======= Track 3: mutation =======
circos.genomicTrack(mutation_df,
track.height = 0.1,
bg.border = NA,
ylim = c(0, 1),
panel.fun = function(region, value, ...) {
circos.segments(region[[1]], 0, region[[2]], 1,
col = col_fun(value[[1]]), lwd = 1)
})
# ======= Track 1: cn gain =======
circos.genomicTrack(gain_df,
track.height = 0.12,
bg.border = NA,
ylim = c(0, max(gain_df$value, na.rm = TRUE)),
panel.fun = function(region, value, ...) {
circos.rect(xleft = region[[1]],
ybottom = 0,
xright = region[[2]],
ytop = value[[1]],
col = "#E04477",
border = NA)
})
# ======= Track 2: cn loss =======
circos.genomicTrack(loss_df,
track.height = 0.12,
bg.border = NA,
ylim = c(0, max(loss_df$value, na.rm = TRUE)),
panel.fun = function(region, value, ...) {
circos.rect(xleft = region[[1]],
ybottom = value[[1]],
xright = region[[2]],
ytop = 0,
col = "#55CC33",
border = NA)
})
# ======= Track 4: label =======
# ======= center legend =======
pushViewport(viewport(x = 0.5, y = 0.65, width = 0.1, height = 0.1, just = c("center", "center")))
grid.text("Legend", y = 0.65, gp = gpar(fontsize = 10, fontface = "bold"))
grid.text("Red: Amplification", y = 0.55, gp = gpar(col = "#E04477", fontsize = 8))
grid.text("Green: Deletion", y = 0.45, gp = gpar(col = "#55CC33", fontsize = 8))
grid.text("Labeled: Selected genes", y = 0.35, gp = gpar(fontsize = 8))
popViewport()
pushViewport(viewport(x = 0.5, y = 0.4))
grid.draw(legend_obj)
popViewport()
dev.off()