PlotScripts/Venn_plot.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
#!/usr/bin/env Rscript

# =============================================================================
# 脚本功能：基于输入文件绘制Venn图
# 输入文件格式：两列，tab分隔，无表头
#   第一列：集合名称（label）
#   第二列：逗号分隔的集合元素（例如：gene1,gene2,gene3）
#   Control   geneA,geneB,geneC,geneD
#   KO        geneB,geneC,geneE,geneF
#   KD        geneA,geneC,geneG

# 输出：Venn图保存到输出目录（默认PDF格式）
# 依赖包：ggvenn, ggplot2, optparse
# 使用示例: Rscript venn_plot.R -i sets.txt -o ./results -f pdf -w 8 -H 8
# =============================================================================

# 加载或安装必要的包
source("/data/med-wangcq/01CondaEnv/02Git_repo/00MyGit_wchenqi/Multi_omics/RNA-seq/00Functions/loadPackages_v4.r")
# library(optparse)
# library(ggvenn)
# library(ggplot2)

# 解析命令行参数
option_list <- list(
  make_option(c("-i", "--infile"), type = "character", default = NULL,
              help = "输入文件路径（两列，tab分隔，第一列label，第二列逗号分隔的元素; 无列名）", metavar = "file"),
  make_option(c("-o", "--outdir"), type = "character", default = ".",
              help = "输出目录 [默认 %default]", metavar = "path"),
  make_option(c("-f", "--format"), type = "character", default = "pdf",
              help = "输出格式：pdf, png, svg 等 [默认 %default]", metavar = "format"),
  make_option(c("-p", "--prefix"), type = "character", default = "no",
              help = "指定输出文件前缀 [默认 %default]"),
  make_option(c("-w", "--width"), type = "numeric", default = 8,
              help = "图片宽度（英寸） [默认 %default]", metavar = "numeric"),
  make_option(c("-H", "--height"), type = "numeric", default = 8,
              help = "图片高度（英寸） [默认 %default]", metavar = "numeric")
)

opt_parser <- OptionParser(option_list = option_list,
                           description = "绘制基于集合列表的Venn图")
opt <- parse_args(opt_parser)

prefix <- opt$prefix
infile <- opt$infile
outdir <- opt$outdir
format <- opt$format
width <- opt$width
height <- opt$height

cat("prefix: ",prefix,"\n")
cat("infile: ",infile,"\n")
cat("outdir: ",outdir,"\n")
cat("format: ",format,"\n")
cat("width: ",width,"\n")
cat("height: ",height,"\n")

# 检查必要参数
if (is.null(infile)) {
  print_help(opt_parser)
  stop("错误：必须提供输入文件（--infile）", call. = FALSE)
}

# 创建输出目录（如果不存在）
if (!dir.exists(outdir)) {
  dir.create(outdir, recursive = TRUE)
}

# 读取输入文件
message("正在读取输入文件：", infile)
input_data <- tryCatch({
  read.table(infile, sep = "\t", header = FALSE, stringsAsFactors = FALSE)
}, error = function(e) {
  stop("无法读取输入文件，请检查文件路径和格式：", e$message)
})

# 验证列数
if (ncol(input_data) < 2) {
  stop("输入文件必须至少包含两列（tab分隔）")
}

# 构建集合列表
sets <- list()
for (i in 1:nrow(input_data)) {
  label <- trimws(input_data[i, 1])
  elements_raw <- trimws(input_data[i, 2])
  # 处理空行或空元素
  if (elements_raw == "" || is.na(elements_raw)) {
    elements <- c()
  } else {
    # 按逗号分割，并去除每个元素前后的空格
    elements <- trimws(unlist(strsplit(elements_raw, ",")))
    # 去除空字符串
    elements <- elements[elements != ""]
  }
  label <- paste0(label,"_",length(elements))
  sets[[label]] <- elements
}

# 检查集合数量
n_sets <- length(sets)
if (n_sets < 2) {
  stop("至少需要2个集合才能绘制Venn图，当前集合数：", n_sets)
}
if (n_sets > 5) {
  warning("Venn图通常适用于2-5个集合，当前集合数为 ", n_sets,
          "，建议改用UpSet图。但仍将尝试绘制，可能视觉效果不佳。")
}

message("检测到 ", n_sets, " 个集合：", paste(names(sets), collapse = ", "))

# 计算每个集合的大小
set_sizes <- sapply(sets, length)
message("集合大小：")
for (name in names(set_sizes)) {
  message("  ", name, ": ", set_sizes[name])
}

# 绘制Venn图
message("正在绘制Venn图...")
venn_plot <- ggvenn(sets,
                    fill_color = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00"),
                    stroke_size = 0.5,
                    set_name_size = 5,
                    text_size = 4)

# 保存图片
outfile <- file.path(outdir, paste0(prefix,"_venn_plot.", format))
ggsave(filename = outfile, plot = venn_plot,
       width = width, height = height, dpi = 300)
message("Venn图已保存至：", outfile)

# 可选：同时输出集合交集的详细信息（文本文件）
intersection_info <- file.path(outdir, paste0(prefix,"_venn_intersections.txt"))
cat("Venn图集合交集统计\n", file = intersection_info)
cat("生成时间：", as.character(Sys.time()), "\n\n", file = intersection_info, append = TRUE)
cat("集合列表：\n", file = intersection_info, append = TRUE)
for (name in names(sets)) {
  cat(sprintf("  %s: %d 个元素\n", name, length(sets[[name]])), file = intersection_info, append = TRUE)
}

# 如果有2-3个集合，输出两两/三三交集的具体元素
if (n_sets <= 3) {
  cat("\n交集详细信息：\n", file = intersection_info, append = TRUE)
  # 生成所有非空子集的组合（排除空集）
  all_combinations <- list()
  for (k in 1:n_sets) {
    combos <- combn(names(sets), k, simplify = FALSE)
    for (combo in combos) {
      intersect_elements <- Reduce(intersect, sets[combo])
      if (length(intersect_elements) > 0) {
        combo_str <- paste(combo, collapse = " ∩ ")
        cat(sprintf("\n%s (%d个元素):\n", combo_str, length(intersect_elements)),
            file = intersection_info, append = TRUE)
        cat(paste(intersect_elements, collapse = ", "), "\n", file = intersection_info, append = TRUE)
      }
    }
  }
}
message("交集统计信息已保存至：", intersection_info)

message("完成！")