Untitled
unknown
plain_text
17 days ago
3.2 kB
3
Indexable
Never
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) # 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, max(res$log_pvalue[res$padj <= sigthresh]) + 1) + # Set upper y limit based on max significant points geom_vline(xintercept = c(-lfcthresh, lfcthresh), linetype="dashed", color="grey") + geom_hline(yintercept = -log10(sigthresh), 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
Leave a Comment