Expr_VolcanoPlot.r

#!/usr/bin/env Rscript

# =============================================================================
# 脚本功能：绘制火山图，支持差异分析结果表格
# 输入：差异分析结果文件（制表符分隔，含基因名、log2FC、p值/padj列）
# 过滤条件：|log2FC| >= 阈值 且 显著性 <= 阈值
# 配色：上调红色，下调蓝色，不显著灰色
# 标注：用户指定基因列表，与显著上调和显著下调基因取交集后标注
# =============================================================================

# 加载所需包
suppressPackageStartupMessages({
  if (!require("optparse", quietly = TRUE)) install.packages("optparse", repos = "https://cloud.r-project.org/")
  if (!require("ggplot2", quietly = TRUE)) install.packages("ggplot2", repos = "https://cloud.r-project.org/")
  if (!require("ggrepel", quietly = TRUE)) install.packages("ggrepel", repos = "https://cloud.r-project.org/")
})

library(optparse)
library(ggplot2)
library(ggrepel)

# -----------------------------------------------------------------------------
# 命令行参数定义
# -----------------------------------------------------------------------------
option_list <- list(
  make_option(c("-i", "--input"), type = "character", default = NULL,
              help = "差异分析结果文件（制表符分隔，必须包含基因名/特征名列、logFC列、p值或padj列）", metavar = "file"),
  make_option(c("--gene_col"), type = "character", default = "gene",
              help = "基因名列名 [默认 %default]", metavar = "colname"),
  make_option(c("--sig_only"), type = "logical", default = FALSE,
              help = "只标注显著分子 [默认 %default]", action = "store_true"),
  make_option(c("--logFC_col"), type = "character", default = "log2FoldChange",
              help = "log2FoldChange列名 [默认 %default]", metavar = "colname"),
  make_option(c("--pvalue_col"), type = "character", default = "pvalue",
              help = "p值列名（当不启用 --use_padj 时使用） [默认 %default]", metavar = "colname"),
  make_option(c("--padj_col"), type = "character", default = "padj",
              help = "校正后p值列名（当启用 --use_padj 时使用） [默认 %default]", metavar = "colname"),
  make_option(c("--use_padj"), action = "store_true", default = FALSE,
              help = "是否使用校正后p值（padj/FDR）进行筛选 [默认 %default]"),
  make_option(c("--fc_threshold"), type = "numeric", default = 1,
              help = "log2FC绝对值阈值 [默认 %default]", metavar = "numeric"),
  make_option(c("--p_threshold"), type = "numeric", default = 0.05,
              help = "p值或padj阈值 [默认 %default]", metavar = "numeric"),
  make_option(c("--up_color"), type = "character", default = "red",
              help = "上调基因颜色 [默认 %default]", metavar = "color"),
  make_option(c("--down_color"), type = "character", default = "blue",
              help = "下调基因颜色 [默认 %default]", metavar = "color"),
  make_option(c("--ns_color"), type = "character", default = "gray",
              help = "不显著基因颜色 [默认 %default]", metavar = "color"),
  make_option(c("--label_file"), type = "character", default = NULL,
              help = "需要标注的基因列表文件（每行一个基因名）", metavar = "file"),
  make_option(c("--label_list"), type = "character", default = NULL,
              help = "需要标注的基因列表（逗号分隔，如 GeneA,GeneB,GeneC）", metavar = "string"),
  make_option(c("--comparison"), type = "character", default = "",
              help = "对比组信息，将显示在X轴下方", metavar = "string"),
  make_option(c("-o", "--output"), type = "character", default = "./",
              help = "输出图片路径 [默认 %default]", metavar = "path"),
  make_option(c("-w", "--width"), type = "numeric", default = 8,
              help = "图片宽度（英寸） [默认 %default]", metavar = "numeric"),
  make_option(c("-H", "--height"), type = "numeric", default = 7,
              help = "图片高度（英寸） [默认 %default]", metavar = "numeric"),
  make_option(c("-t", "--title"), type = "character", default = "Volcano Plot",
              help = "图片标题 [默认 %default]", metavar = "string"),
  make_option(c("--xlab"), type = "character", default = "Log2 Fold Change",
              help = "X轴标签 [默认 %default]", metavar = "string"),
  make_option(c("--ylab"), type = "character", default = "-Log10(P-value)",
              help = "Y轴标签 [默认 %default]", metavar = "string"),
  make_option(c("--point_size"), type = "numeric", default = 1.5,
              help = "点大小 [默认 %default]", metavar = "numeric"),
  make_option(c("--label_size"), type = "numeric", default = 3,
              help = "标注字体大小 [默认 %default]", metavar = "numeric"),
  make_option(c("--max_overlaps"), type = "integer", default = 10,
              help = "ggrepel最大重叠数 [默认 %default]", metavar = "integer")
)

opt_parser <- OptionParser(option_list = option_list, description = "绘制火山图（支持差异分析结果）")
opt <- parse_args(opt_parser)

# 检查必需参数
if (is.null(opt$input)) {
  print_help(opt_parser)
  stop("错误：必须提供差异分析结果文件（--input）", call. = FALSE)
}

# -----------------------------------------------------------------------------
# 赋值
# -----------------------------------------------------------------------------
input <- opt$input
gene_col <- opt$gene_col
sig_only <- opt$sig_only
logFC_col <- opt$logFC_col
pvalue_col <- opt$pvalue_col
padj_col <- opt$padj_col
use_padj <- opt$use_padj
fc_threshold <- opt$fc_threshold
p_threshold <- opt$p_threshold
up_color <- opt$up_color
down_color <- opt$down_color
ns_color <- opt$ns_color
label_file <- opt$label_file
label_list <- opt$label_list
comparison <- opt$comparison
output <- opt$output
width <- opt$width
height <- opt$height
title <- opt$title
xlab <- opt$xlab
ylab <- opt$ylab
point_size <- opt$point_size
label_size <- opt$label_size
max_overlaps <- opt$max_overlaps

print(input)
print(gene_col)
print(sig_only)
print(logFC_col)
print(pvalue_col)
print(padj_col)
print(use_padj)
print(fc_threshold)
print(p_threshold)
print(label_file)
print(label_list)
print(comparison)
print(output)
print(width)
print(height)
print(title)
print(max_overlaps)
print(point_size)     # 1.5
print(label_size)     # 3
print(up_color)       # "red"
print(down_color)     # "blue"
print(ns_color)       # "gray"
print(xlab)           # "Log2 Fold Change"
print(ylab)           # "-Log10(P-value)"
# -----------------------------------------------------------------------------
# 读取数据
# -----------------------------------------------------------------------------
data <- tryCatch({
  read.csv(input, header = TRUE, stringsAsFactors = FALSE)
}, error = function(e) {
  stop("读取输入文件失败：", e$message)
})
print(head(data))
# 检查必需列是否存在
required_cols <- c(gene_col, logFC_col)
if (use_padj) {
  sig_col <- padj_col
} else {
  sig_col <- pvalue_col
}
required_cols <- c(required_cols, sig_col)
missing_cols <- setdiff(required_cols, colnames(data))
if (length(missing_cols) > 0) {
  stop("输入文件中缺少以下必需列：", paste(missing_cols, collapse = ", "))
}

# 提取数据
genes <- data[[gene_col]]
logFC <- data[[logFC_col]]
p_sig <- data[[sig_col]]

# 计算 -log10(p值)
negLogP <- -log10(p_sig)

# -----------------------------------------------------------------------------
# 定义显著性分组
# -----------------------------------------------------------------------------
is_sig <- abs(logFC) >= fc_threshold & p_sig <= p_threshold
is_up <- is_sig & logFC > 0
is_down <- is_sig & logFC < 0
group <- ifelse(is_up, "Up", ifelse(is_down, "Down", "NS"))

# 创建绘图数据框
plot_df <- data.frame(Gene = genes,
                      logFC = logFC,
                      negLogP = negLogP,
                      Group = group,
                      stringsAsFactors = FALSE
                    )

# -----------------------------------------------------------------------------
# 处理标注基因：与显著差异基因取交集
# -----------------------------------------------------------------------------
label_genes <- c()
if (!is.null(label_file)) {
  label_genes <- tryCatch({
    as.character(read.csv(label_file, header = FALSE, stringsAsFactors = FALSE)[[1]])
  }, error = function(e) {
    warning("读取标注文件失败：", e$message)
    return(c())
  })
}
if (!is.null(label_list)) {
  label_genes <- c(label_genes, trimws(unlist(strsplit(label_list, ","))))
}
label_genes <- unique(label_genes)

# 取交集：只标注既在label_genes中的基因
if(sig_only){
  label_keep <- intersect(label_genes, plot_df$Gene[which(plot_df$Group != "NS")])
}else{
  label_keep <- label_genes
}

plot_df$Label <- ifelse(plot_df$Gene %in% label_keep, plot_df$Gene, "")
top_label <- plot_df[plot_df$Label != "",]
# 输出标注的基因数量
message("标注基因总数（用户提供）：", length(label_genes))
message("显著差异基因中可标注的基因数：", length(label_keep))
if (length(label_keep) > 0) {
  message("标注的基因：", paste(label_keep, collapse = ", "))
}

# -----------------------------------------------------------------------------
# 绘制火山图
# -----------------------------------------------------------------------------
# 设置颜色映射
color_map <- c("Up" = up_color, "Down" = down_color, "NS" = ns_color)

# 创建ggplot对象
p <- ggplot(plot_df, aes(x = logFC, y = negLogP, color = Group)) +
            geom_point(alpha = 0.5, size = point_size) +
            scale_color_manual(values = color_map) +
            geom_hline(yintercept = -log10(p_threshold), linetype = "dashed", color = "grey40") +
            geom_vline(xintercept = c(-fc_threshold, fc_threshold), linetype = "dashed", color = "grey40")
# 添加文本标注（仅当有可标注基因时）
if (nrow(plot_df[plot_df$Label != "", ]) > 0) {
    p <- p + geom_text_repel(data = top_label,
                            aes(x = logFC, y = negLogP, label = Label),
                            color = "black", size = 2,
                            inherit.aes = FALSE, max.overlaps = as.numeric(max_overlaps),
                            box.padding = 0.5, min.segment.length = 0.1) +
             labs(x = "logFC", y = "-log10(FDR)") +
             theme_bw() + theme(legend.position = "none")
} else {
  message("警告：没有可标注的显著差异基因，将跳过文本标注")
}

# -----------------------------------------------------------------------------
# 保存图片
# -----------------------------------------------------------------------------
# 创建输出目录
out_dir <- output
if (!dir.exists(out_dir)) {
  dir.create(out_dir, recursive = TRUE)
}
filename <- gsub(" ","_",title)
ggsave(file.path(out_dir,paste0(filename,".pdf")), plot = p, width = width, height = height)
write.csv(plot_df,file.path(out_dir,paste0(filename,".csv")),row.names=TRUE,quote=FALSE)

message("火山图已保存至：", output)