Untitled

# Custom volcano plot function with horizontal lines for labels
volcanoplot_gg <- function(res, lfcthresh=0.5, sigthresh=0.05, main="DAPs_G_vs_P", label_genes=NULL, label_genes_2=NULL) {
    res$color <- ifelse(res$padj > sigthresh | abs(res$log2FoldChange) < lfcthresh, "#D6D6D6",
                        ifelse(res$log2FoldChange > 0,  "#C02628", "#2E417A"))
    res$log_pvalue <- -log10(res$pvalue)
    filter=res[which(res$padj < 0.05 & abs(res$log2FoldChange) > 0.5), ]

    # Function to get the most significant points for labeling
    get_most_significant <- function(df, genes) {
        df <- subset(df, real_name %in% genes & padj < sigthresh & abs(log2FoldChange) > lfcthresh)
        df <- df[order(df$padj), ]  # Order by significance
        df <- df[!duplicated(df$real_name), ]  # Keep only the most significant entry per gene
        df$nudge_x <- ifelse(df$log2FoldChange > 0, 10, -10)  # Nudge to the right for positive, left for negative
        return(df)
    }

    # Base plot
    p <- ggplot(res, aes(x=log2FoldChange, y=log_pvalue, color=color)) +
        geom_point(size=1.5) +
        scale_color_identity() +
        theme_minimal() +
        theme(
            plot.title = element_text(hjust=0.5, face="bold", size=14),  # Center and bold title
            panel.border = element_rect(color = "black", fill = NA, size = 0.5),  # Thinner frame
            panel.grid = element_blank(),  # Remove grid
            axis.line = element_line(color = "black"),
            axis.text = element_text(size=12),  # Larger margin numbers
            axis.title = element_text(size=14),  # Larger axis titles
            legend.position = "none"  # Remove legend
        ) +
        labs(title=main, x="log2(Fold Change)", y="-log10(p-value)") +
        xlim(-5, 5) +
        ylim(0, 40) +
        geom_vline(xintercept = c(-lfcthresh, lfcthresh), linetype="dashed", color="grey") +
        geom_hline(yintercept = -log10(max(filter$pvalue)), linetype="dashed", color="grey")


    # Add gene labels for the first list
    if (!is.null(label_genes)) {
        label_data <- get_most_significant(res, label_genes)
        p <- p + geom_text_repel(
                data=label_data, aes(label=real_name),
                size=3, color="black", segment.color="black", max.overlaps=Inf, force=10,
                box.padding = unit(1, "lines"), point.padding = unit(0.1, "lines"), 
                nudge_x = label_data$nudge_x, nudge_y = 0
            )
    }

    # Add gene labels for the second list with grey color
    if (!is.null(label_genes_2)) {
        label_data_2 <- get_most_significant(res, label_genes_2)
        p <- p + geom_text_repel(
                data=label_data_2, aes(label=real_name),
                size=3, color="grey", segment.color="grey", max.overlaps=Inf, force=10,
                box.padding = unit(1, "lines"), point.padding = unit(0.1, "lines"), 
                nudge_x = label_data_2$nudge_x, nudge_y = 0
            )
    }

    return(p)
}

# Define genes to label
label_genes <- c("SMAD3","FOS","ITGAM","CSF1","TANK")
label_genes_2 <- c()
# Create the plot
p <- volcanoplot_gg(daps_P_vs_S_ann, main="DAPs PT vs P", lfcthresh=0.5, sigthresh=0.05, label_genes=label_genes, label_genes_2=label_genes_2)
p
# Save the plot as PNG
# ggsave("figures/PT_vs_P_ATAC_VP_genes.png", plot = p, width = 8, height = 6, units = "in", dpi = 150)