PlotScripts/Expr_VolcanoPlot.r at main · wchenqi/PlotScripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
#!/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("--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)
}

# -----------------------------------------------------------------------------
# 读取数据
# -----------------------------------------------------------------------------
data <- tryCatch({
  read.csv(opt$input, header = TRUE, stringsAsFactors = FALSE)
}, error = function(e) {
  stop("读取输入文件失败：", e$message)
})
print(head(data))
# 检查必需列是否存在
required_cols <- c(opt$gene_col, opt$logFC_col)
if (opt$use_padj) {
  sig_col <- opt$padj_col
} else {
  sig_col <- opt$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[[opt$gene_col]]
logFC <- data[[opt$logFC_col]]
p_sig <- data[[sig_col]]

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

# -----------------------------------------------------------------------------
# 定义显著性分组
# -----------------------------------------------------------------------------
is_sig <- abs(logFC) >= opt$fc_threshold & p_sig <= opt$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(opt$label_file)) {
  label_genes <- tryCatch({
    as.character(read.csv(opt$label_file, header = FALSE, stringsAsFactors = FALSE)[[1]])
  }, error = function(e) {
    warning("读取标注文件失败：", e$message)
    return(c())
  })
}
if (!is.null(opt$label_list)) {
  label_genes <- c(label_genes, trimws(unlist(strsplit(opt$label_list, ","))))
}
label_genes <- unique(label_genes)

# 取交集：只标注既在label_genes中的基因
label_keep <- intersect(label_genes, plot_df$Gene)
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" = opt$up_color, "Down" = opt$down_color, "NS" = opt$ns_color)

# 创建ggplot对象
p <- ggplot(plot_df, aes(x = logFC, y = negLogP, color = Group)) +
            geom_point(alpha = 0.5, size = opt$point_size) +
            scale_color_manual(values = color_map) +
            geom_hline(yintercept = -log10(opt$p_threshold), linetype = "dashed", color = "grey40") +
            geom_vline(xintercept = c(-opt$fc_threshold, opt$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 = 15,
                            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 <- opt$output
if (!dir.exists(out_dir)) {
  dir.create(out_dir, recursive = TRUE)
}
ggsave(file.path(out_dir,"VolcanoPlot.pdf"), plot = p, width = opt$width, height = opt$height)

write.csv(plot_df,file.path(out_dir,"VolcanoPlot.csv"),row.names=TRUE,quote=FALSE)

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