313 lines
11 KiB
R
Executable File
313 lines
11 KiB
R
Executable File
#!/usr/bin/env Rscript
|
||
# 启动子分析增强可视化脚本 V3
|
||
# 修订内容:
|
||
# 1. 七大分类顺序修正,Other放最后,分类内细分元件
|
||
# 2. 组合图:现有进化树PDF + category_stacked结构 + 红绿渐变热图
|
||
# 3. big7柱状图细分元件类别
|
||
|
||
rm(list = ls())
|
||
|
||
cat("=== 加载R包 ===\n")
|
||
required_packages <- c("tidyverse", "ggtree", "treeio", "ape",
|
||
"patchwork", "aplot", "RColorBrewer", "ggnewscale")
|
||
|
||
for(pkg in required_packages) {
|
||
if(!require(pkg, character.only = TRUE, quietly = TRUE)) {
|
||
cat("错误: 包", pkg, "未安装\n")
|
||
quit(save = "no", status = 1)
|
||
}
|
||
}
|
||
|
||
cat("✓ 所有包加载完成\n\n")
|
||
|
||
####----命令行参数----####
|
||
args <- commandArgs(trailingOnly = TRUE)
|
||
|
||
tree_pdf_file <- args[1]
|
||
group_file <- args[2]
|
||
promoter_file <- args[3]
|
||
element_desc_file <- ifelse(length(args) >= 4, args[4], NA)
|
||
output_prefix <- ifelse(length(args) >= 5, args[5], "promoter_analysis")
|
||
|
||
tree_nwk_file <- sub("_promoter_phylogenetic_tree\\.pdf$", ".tree.nwk", tree_pdf_file)
|
||
if(!file.exists(tree_nwk_file)) {
|
||
tree_nwk_file <- sub("\\.pdf$", ".nwk", tree_pdf_file)
|
||
}
|
||
|
||
cat("=== 启动子增强可视化 V3 ===\n")
|
||
cat("进化树PDF:", tree_pdf_file, "\n")
|
||
cat("进化树NWK:", tree_nwk_file, "\n\n")
|
||
|
||
####----读取数据----####
|
||
tree <- treeio::read.newick(tree_nwk_file)
|
||
|
||
group_df <- read.delim(group_file, header = FALSE, stringsAsFactors = FALSE)
|
||
if(ncol(group_df) >= 2) colnames(group_df) <- c("label", "Subgroup")
|
||
group_df$label <- group_df$label %>% str_trim()
|
||
group_df$label_base <- group_df$label %>%
|
||
str_remove("\\.\\d+\\.T\\d+$") %>%
|
||
str_remove("\\.T\\d+$") %>%
|
||
str_remove("\\.\\d+$")
|
||
|
||
promoter <- read.delim(promoter_file, header = FALSE, stringsAsFactors = FALSE)
|
||
if(ncol(promoter) >= 8) {
|
||
promoter <- promoter[, c(1, 2, 4, 5, 6, 8)]
|
||
colnames(promoter) <- c("gene", "element", "position", "length", "strand", "description")
|
||
}
|
||
promoter$gene <- promoter$gene %>% str_trim()
|
||
|
||
# 定义7大分类和顺序(Other放最后)
|
||
big7_order <- c("Core promoter", "Light response", "Hormone signaling",
|
||
"Stress response", "TF binding", "Development", "Other")
|
||
|
||
# 元件分类函数
|
||
categorize_element <- function(x) {
|
||
if(grepl("TATA|CAAT|Inr|Y-patch", x, ignore.case = TRUE)) return("Core promoter")
|
||
if(grepl("ABRE|GA|JA|SA|Auxin|Ethylene|hormone", x, ignore.case = TRUE)) return("Hormone signaling")
|
||
if(grepl("G-box|GT1|I-box|Sp1|Box.4|GATA|light", x, ignore.case = TRUE)) return("Light response")
|
||
if(grepl("DRE|LTR|ARE|STRE|W.box|W-box|stress|drought|cold", x, ignore.case = TRUE)) return("Stress response")
|
||
if(grepl("MYB|MYC|WRKY|TF|binding", x, ignore.case = TRUE)) return("TF binding")
|
||
if(grepl("RY|AAGAA|CAT|GCN4|Development|seed|meristem", x, ignore.case = TRUE)) return("Development")
|
||
return("Other")
|
||
}
|
||
|
||
if(!is.na(element_desc_file) && file.exists(element_desc_file)) {
|
||
element_desc <- read.csv(element_desc_file, stringsAsFactors = FALSE)
|
||
element_desc$big7 <- sapply(element_desc$element, categorize_element)
|
||
element_desc$big7 <- factor(element_desc$big7, levels = big7_order)
|
||
promoter <- promoter %>%
|
||
left_join(element_desc %>% select(element, big7), by = "element")
|
||
} else {
|
||
promoter$big7 <- sapply(promoter$element, categorize_element)
|
||
promoter$big7 <- factor(promoter$big7, levels = big7_order)
|
||
}
|
||
|
||
tree_labels <- tree$tip.label
|
||
tree_labels_base <- tree_labels %>%
|
||
str_remove("\\.\\d+\\.T\\d+$") %>%
|
||
str_remove("\\.T\\d+$") %>%
|
||
str_remove("\\.\\d+$")
|
||
|
||
group_df <- group_df %>% filter(label_base %in% tree_labels_base)
|
||
promoter <- promoter %>% filter(gene %in% tree_labels_base)
|
||
|
||
cat("数据读取完成: 树", length(tree_labels), "基因, 启动子", nrow(promoter), "记录\n\n")
|
||
|
||
# 保存函数
|
||
save_plot <- function(plot, filename_base, width, height) {
|
||
pdf_file <- paste0(filename_base, ".pdf")
|
||
png_file <- paste0(filename_base, ".png")
|
||
ggsave(pdf_file, plot, width = width, height = height, device = cairo_pdf)
|
||
cat(" PDF:", pdf_file, "\n")
|
||
tryCatch({
|
||
png(filename = png_file, width = width * 100, height = height * 100, res = 100)
|
||
print(plot)
|
||
dev.off()
|
||
cat(" PNG:", png_file, "\n")
|
||
}, error = function(e) {
|
||
cat(" PNG生成失败\n")
|
||
})
|
||
}
|
||
|
||
####----可视化1: 三合一组合图----####
|
||
cat("[1/4] 绘制三合一组合图...\n")
|
||
|
||
# 1A. 尝试读取已有进化树图(如果magick可用)
|
||
p_tree_external <- NULL
|
||
if(require("magick", quietly = TRUE)) {
|
||
tryCatch({
|
||
tree_img <- image_read(tree_pdf_file)
|
||
tree_grob <- rasterGrob(tree_img, interpolate = TRUE)
|
||
p_tree_external <- ggplot() +
|
||
annotation_custom(tree_grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf) +
|
||
theme_void()
|
||
cat(" ✓ 已读取外部进化树图\n")
|
||
}, error = function(e) {
|
||
cat(" ⚠ 读取外部树图失败,将使用ggtree绘制\n")
|
||
})
|
||
} else {
|
||
cat(" ℹ magick包未安装,将使用ggtree绘制树\n")
|
||
}
|
||
|
||
# 如果外部图读取失败,使用ggtree
|
||
if(is.null(p_tree_external)) {
|
||
group_df_for_tree <- group_df %>% mutate(label = label_base)
|
||
subfam_colors <- brewer.pal(max(3, length(unique(group_df$Subgroup))), "Dark2")
|
||
subfam_colors <- subfam_colors[1:length(unique(group_df$Subgroup))]
|
||
names(subfam_colors) <- unique(group_df$Subgroup)
|
||
|
||
p_tree_external <- ggtree(tree) %<+% group_df_for_tree +
|
||
geom_tippoint(aes(color = Subgroup), size = 3) +
|
||
scale_color_manual(values = subfam_colors, name = "Subfamily") +
|
||
theme_tree2() +
|
||
theme(legend.position = "none")
|
||
}
|
||
|
||
# 1B. 结构图: category_stacked样式(显示7大分类)
|
||
structure_data <- promoter %>%
|
||
group_by(gene, big7) %>%
|
||
summarise(count = n(), .groups = "drop") %>%
|
||
complete(gene = tree_labels_base, big7, fill = list(count = 0)) %>%
|
||
group_by(gene) %>%
|
||
mutate(prop = count / sum(count)) %>%
|
||
ungroup() %>%
|
||
mutate(gene = factor(gene, levels = rev(tree_labels_base)))
|
||
|
||
big7_colors <- c(
|
||
"Core promoter" = "#e41a1c",
|
||
"Light response" = "#377eb8",
|
||
"Hormone signaling" = "#4daf4a",
|
||
"Stress response" = "#984ea3",
|
||
"TF binding" = "#ff7f00",
|
||
"Development" = "#f781bf",
|
||
"Other" = "#999999"
|
||
)
|
||
|
||
p_structure <- ggplot(structure_data, aes(x = prop, y = gene, fill = big7)) +
|
||
geom_bar(stat = "identity", position = "stack", width = 0.8) +
|
||
scale_fill_manual(values = big7_colors, name = "Category") +
|
||
labs(x = "Proportion", y = NULL) +
|
||
theme_bw() +
|
||
theme(
|
||
axis.text.y = element_blank(),
|
||
axis.ticks.y = element_blank(),
|
||
panel.grid = element_blank()
|
||
)
|
||
|
||
# 1C. 热图: 红绿渐变,小方格
|
||
top_elements <- promoter %>%
|
||
group_by(element) %>%
|
||
summarise(count = n(), .groups = "drop") %>%
|
||
arrange(desc(count)) %>%
|
||
head(25) %>%
|
||
pull(element)
|
||
|
||
heatmap_subset <- promoter %>%
|
||
filter(element %in% top_elements) %>%
|
||
group_by(gene, element) %>%
|
||
summarise(count = n(), .groups = "drop") %>%
|
||
complete(gene = tree_labels_base, element, fill = list(count = 0)) %>%
|
||
mutate(gene = factor(gene, levels = rev(tree_labels_base)))
|
||
|
||
p_heatmap_rg <- ggplot(heatmap_subset, aes(x = element, y = gene, fill = count)) +
|
||
geom_tile(color = "white", size = 0.2) +
|
||
scale_fill_gradient2(low = "green", mid = "black", high = "red",
|
||
midpoint = max(heatmap_subset$count)/2,
|
||
name = "Count") +
|
||
labs(x = "Element", y = NULL) +
|
||
theme_bw() +
|
||
theme(
|
||
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 7),
|
||
axis.text.y = element_blank(),
|
||
axis.ticks.y = element_blank(),
|
||
panel.grid = element_blank()
|
||
)
|
||
|
||
# 1D. 组合(图例放最右)
|
||
p_combined <- (p_tree_external + theme(legend.position = "none")) |
|
||
(p_structure + theme(legend.position = "none")) |
|
||
(p_heatmap_rg + theme(legend.position = "none")) |
|
||
guide_area()
|
||
|
||
p_combined <- p_combined + plot_layout(widths = c(2, 1.5, 2.5, 0.5), guides = "collect")
|
||
|
||
save_plot(p_combined, paste0(output_prefix, "_01_combined"), 24, 14)
|
||
cat(" ✓ 三合一组合图完成\n\n")
|
||
|
||
####----可视化2: 7大分类细分柱状图----####
|
||
cat("[2/4] 绘制7大分类细分柱状图...\n")
|
||
|
||
# 统计每个大分类内的元件
|
||
big7_element_stats <- promoter %>%
|
||
group_by(big7, element) %>%
|
||
summarise(count = n(), .groups = "drop") %>%
|
||
group_by(big7) %>%
|
||
arrange(desc(count)) %>%
|
||
mutate(element_rank = row_number()) %>%
|
||
ungroup()
|
||
|
||
# 每个大分类取Top 10元件
|
||
top_elements_per_big7 <- big7_element_stats %>%
|
||
group_by(big7) %>%
|
||
slice_max(order_by = count, n = 10) %>%
|
||
ungroup()
|
||
|
||
# 绘制细分柱状图(每个大分类单独排序)
|
||
p_big7_detailed <- ggplot(top_elements_per_big7,
|
||
aes(x = reorder(element, count), y = count, fill = big7)) +
|
||
geom_bar(stat = "identity", width = 0.7) +
|
||
scale_fill_manual(values = big7_colors, name = "Category") +
|
||
facet_wrap(~big7, scales = "free", ncol = 3) +
|
||
coord_flip() +
|
||
labs(x = "Element", y = "Count", title = "Top 10 Elements in Each of 7 Major Categories") +
|
||
theme_bw() +
|
||
theme(
|
||
plot.title = element_text(size = 14, face = "bold"),
|
||
axis.text.y = element_text(size = 7),
|
||
strip.background = element_rect(fill = "lightblue"),
|
||
strip.text = element_text(face = "bold")
|
||
)
|
||
|
||
save_plot(p_big7_detailed, paste0(output_prefix, "_02_big7_detailed"), 16, 12)
|
||
cat(" ✓ 7大分类细分柱状图完成\n\n")
|
||
|
||
####----可视化3: 元件热图(红绿渐变)----####
|
||
cat("[3/4] 绘制红绿渐变热图...\n")
|
||
|
||
heatmap_full <- promoter %>%
|
||
group_by(gene, element) %>%
|
||
summarise(count = n(), .groups = "drop") %>%
|
||
complete(gene = tree_labels_base, element, fill = list(count = 0)) %>%
|
||
mutate(gene = factor(gene, levels = rev(tree_labels_base)))
|
||
|
||
p_heatmap_full <- ggplot(heatmap_full, aes(x = element, y = gene, fill = count)) +
|
||
geom_tile(color = "white", size = 0.1) +
|
||
geom_text(aes(label = ifelse(count > 0, count, "")), size = 2, color = "white") +
|
||
scale_fill_gradient2(low = "darkgreen", mid = "black", high = "red",
|
||
midpoint = median(heatmap_full$count[heatmap_full$count > 0]),
|
||
name = "Count") +
|
||
labs(x = "Cis-acting Element", y = "Gene") +
|
||
theme_bw() +
|
||
theme(
|
||
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6),
|
||
axis.text.y = element_text(size = 8),
|
||
panel.grid = element_blank()
|
||
)
|
||
|
||
save_plot(p_heatmap_full, paste0(output_prefix, "_03_heatmap_rg"), 18, 12)
|
||
cat(" ✓ 红绿渐变热图完成\n\n")
|
||
|
||
####----可视化4: 类别堆叠图----####
|
||
cat("[4/4] 绘制类别堆叠图...\n")
|
||
|
||
category_props <- promoter %>%
|
||
group_by(gene, big7) %>%
|
||
summarise(count = n(), .groups = "drop") %>%
|
||
complete(gene = tree_labels_base, big7, fill = list(count = 0)) %>%
|
||
group_by(gene) %>%
|
||
mutate(prop = count / sum(count)) %>%
|
||
ungroup() %>%
|
||
mutate(gene = factor(gene, levels = rev(tree_labels_base)))
|
||
|
||
p_stacked <- ggplot(category_props, aes(x = gene, y = prop, fill = big7)) +
|
||
geom_bar(stat = "identity", position = "stack", width = 0.8) +
|
||
scale_fill_manual(values = big7_colors, name = "Category") +
|
||
coord_flip() +
|
||
labs(x = "Gene", y = "Proportion", title = "7 Major Categories Distribution by Gene") +
|
||
theme_bw() +
|
||
theme(
|
||
plot.title = element_text(size = 14, face = "bold"),
|
||
axis.text.y = element_text(size = 7),
|
||
legend.position = "right"
|
||
)
|
||
|
||
save_plot(p_stacked, paste0(output_prefix, "_04_stacked"), 12, 10)
|
||
cat(" ✓ 类别堆叠图完成\n\n")
|
||
|
||
####----完成----####
|
||
cat("===== 所有可视化完成 =====\n")
|
||
cat("输出文件:\n")
|
||
for(i in 1:4) {
|
||
cat(sprintf("%d. %s_0%d_*.pdf/png\n", i, output_prefix, i))
|
||
}
|