Expr_Boxplot_v2.r

#!/usr/bin/env Rscript

## 脚本说明
# 1) 运行环境: mamba activate r452
# 2) 应用:  查看基因表达（箱型图 + 柱状图 + 堆叠图）

# =============================================================================
# bulkRNAseq Gene Expression Boxplot / Barplot / Stacked Script
# Author: Auto-generated
# Date: 2026-07-25
# =============================================================================

suppressPackageStartupMessages({
    library(ggplot2)
    library(ggpubr)
    library(dplyr)
    library(tidyr)
    library(tibble)
    library(RColorBrewer)
    library(optparse)
})

# =============================================================================
# Option Parser (optparse)
# =============================================================================
option_list <- list(
    make_option(c("-C", "--counts"), type = "character", action = "store", 
                help = "Path to gene count matrix (genes x samples, TSV/CSV) [required]"),
    make_option(c("-M", "--meta"), type = "character", action = "store",
                help = "Path to sample metadata (sample_id, group, ...) [required]"),
    make_option(c("-g", "--genes"), type = "character", action = "store",
                help = "Comma-separated gene list or file with one gene per line [required]"),
    make_option(c("-o", "--outdir"), type = "character", action = "store", default = "./",
                help = "Output directory [default: ./]"),
    make_option(c("-p", "--outprefix"), type = "character", action = "store", default = "boxplot",
                help = "Output file prefix [default: boxplot]"),
    make_option(c("-G", "--group_col"), type = "character", action = "store", default = "group",
                help = "Column name for grouping in metadata [default: group]"),
    make_option(c("-S", "--sample_col"), type = "character", action = "store", default = "sample",
                help = "Column name for sample ID in metadata [default: sample]"),
    make_option(c("-O", "--group_order"), type = "character", action = "store", default = NULL,
                help = "Comma-separated group order (e.g., 'WT,KO,OE')"),
    make_option(c("-c", "--colors"), type = "character", action = "store", default = NULL,
                help = "Comma-separated custom colors (e.g., '#E41A1C,#377EB8')"),
    make_option(c("-T", "--test_method"), type = "character", action = "store", default = "wilcox.test",
                help = "Statistical test method: wilcox.test, t.test, anova, kruskal.test, none [default: wilcox.test]"),
    make_option(c("-m", "--comparisons"), type = "character", action = "store", default = NULL,
                help = "Custom comparisons: 'WT vs KO,KO vs OE'"),
    make_option(c("-y", "--y_label"), type = "character", action = "store", default = "Normalized Expression (log2)",
                help = "Y-axis label [default: 'Normalized Expression (log2)']"),
    make_option(c("-x", "--x_label"), type = "character", action = "store", default = NULL,
                help = "X-axis label [default: group_col]"),
    make_option(c("-t", "--title"), type = "character", action = "store", default = NULL,
                help = "Plot title [default: gene name]"),
    make_option(c("-W", "--width"), type = "double", action = "store", default = 6,
                help = "Plot width [default: 6]"),
    make_option(c("-H", "--height"), type = "double", action = "store", default = 5,
                help = "Plot height [default: 5]"),
    make_option(c("-d", "--dpi"), type = "integer", action = "store", default = 300,
                help = "Output DPI [default: 300]"),
    make_option(c("-f", "--font_size"), type = "double", action = "store", default = 12,
                help = "Base font size [default: 12]"),
    make_option(c("-j", "--no_jitter"), type = "logical", action = "store_true", default = FALSE,
                help = "Hide individual data points (jitter)"),
    make_option(c("-l", "--log_transform"), type = "logical", action = "store_true", default = FALSE,
                help = "Apply log2(count + 1) transformation"),
    make_option(c("-n", "--ncol"), type = "integer", action = "store", default = NULL,
                help = "Number of columns for faceting [default: auto]"),
    make_option(c("-F", "--format"), type = "character", action = "store", default = "pdf",
                help = "Output format: pdf, png, svg, tiff [default: pdf]"),
    make_option(c("-P", "--paired"), type = "logical", action = "store_true", default = FALSE,
                help = "Use paired test (requires pair_id column in metadata)"),
    make_option(c("-I", "--pair_col"), type = "character", action = "store", default = "pair_id",
                help = "Column name for paired samples [default: pair_id]"),
    make_option(c("-s", "--y_start"), type = "double", action = "store", default = 0,
                help = "Y-axis start value [default: 0]"),
    # ---- 新增参数 ----
    make_option(c("-N", "--genes_per_page"), type = "integer", action = "store", default = 6,
                help = "Max genes per page for faceting [default: 6]"),
    make_option(c("-R", "--row_height"), type = "double", action = "store", default = 3.5,
                help = "Row height (inch) per facet row [default: 3.5]"),
    make_option(c("-B", "--barplot"), type = "logical", action = "store_true", default = FALSE,
                help = "Also output bar plot (Mean ± SEM)"),
    make_option(c("-K", "--stacked"), type = "logical", action = "store_true", default = FALSE,
                help = "Also output stacked bar plot for multi-gene"),
    make_option(c("-E", "--stacked_percent"), type = "logical", action = "store_true", default = FALSE,
                help = "Also output percent stacked bar plot for multi-gene")
)

parser <- OptionParser(option_list = option_list,
                       usage = "Usage: %prog [options]
Description: Draw boxplot/barplot/stacked for bulkRNAseq gene expression")
args <- parse_args(parser)

# =============================================================================
# Helper Functions
# =============================================================================

# Error bar: mean +- SEM
mean_se <- function(x) {
    n <- length(x)
    m <- mean(x, na.rm = TRUE)
    se <- sd(x, na.rm = TRUE) / sqrt(n)
    data.frame(y = m, ymin = m - se, ymax = m + se)
}

# =============================================================================
# Load Data
# =============================================================================

# Read count matrix
counts_file <- args$counts
if (grepl("\\.csv$", counts_file, ignore.case = TRUE)) {
    counts <- read.csv(counts_file, row.names = 1, check.names = FALSE, stringsAsFactors = FALSE)
} else {
    counts <- read.delim(counts_file, row.names = 1, check.names = FALSE, stringsAsFactors = FALSE)
}

# Read metadata
meta_file <- args$meta
if (grepl("\\.csv$", meta_file, ignore.case = TRUE)) {
    meta <- read.csv(meta_file, stringsAsFactors = FALSE)
} else {
    meta <- read.delim(meta_file, stringsAsFactors = FALSE)
}

# Parse gene list
if (file.exists(args$genes)) {
    gene_list <- readLines(args$genes)
    gene_list <- trimws(gene_list)
    gene_list <- gene_list[gene_list != ""]
} else {
    gene_list <- strsplit(args$genes, ",")[[1]]
    gene_list <- trimws(gene_list)
}

# Check genes exist
missing_genes <- setdiff(gene_list, rownames(counts))
if (length(missing_genes) > 0) {
    warning("Genes not found in count matrix: ", paste(missing_genes, collapse = ", "))
    gene_list <- intersect(gene_list, rownames(counts))
}

if (length(gene_list) == 0) {
    stop("No valid genes found in count matrix!")
}

# Subset counts
counts_sub <- counts[gene_list, , drop = FALSE]

# =============================================================================
# Prepare Data
# =============================================================================

# Match samples
common_samples <- intersect(colnames(counts_sub), meta[[args$sample_col]])
if (length(common_samples) == 0) {
    stop("No matching samples between count matrix and metadata!")
}

counts_sub <- counts_sub[, common_samples, drop = FALSE]
meta <- meta[meta[[args$sample_col]] %in% common_samples, ]

# Ensure sample order matches
meta <- meta[match(colnames(counts_sub), meta[[args$sample_col]]), ]

# Long format
df_long <- counts_sub %>%
    as.data.frame() %>%
    rownames_to_column("Gene") %>%
    pivot_longer(cols = -Gene, names_to = "Sample", values_to = "Expression")

# Add group info
df_long[[args$group_col]] <- meta[[args$group_col]][match(df_long$Sample, meta[[args$sample_col]])]

# Log transform if requested
if (args$log_transform) {
    df_long$Expression <- log2(df_long$Expression + 1)
    if (args$y_label == "Normalized Expression (log2)") {
        args$y_label <- "log2(Expression + 1)"
    }
}

# Set group order
if (!is.null(args$group_order)) {
    group_order <- strsplit(args$group_order, ",")[[1]]
    group_order <- trimws(group_order)
    df_long[[args$group_col]] <- factor(df_long[[args$group_col]], levels = group_order)
    meta[[args$group_col]] <- factor(meta[[args$group_col]], levels = group_order)
} else {
    df_long[[args$group_col]] <- factor(df_long[[args$group_col]])
}

# Remove NA groups
df_long <- df_long %>% filter(!is.na(!!sym(args$group_col)))

# =============================================================================
# Color Palette (Group)
# =============================================================================

n_groups <- length(unique(df_long[[args$group_col]]))

if (!is.null(args$colors)) {
    custom_colors <- strsplit(args$colors, ",")[[1]]
    custom_colors <- trimws(custom_colors)
    if (length(custom_colors) >= n_groups) {
        fill_colors <- custom_colors[1:n_groups]
    } else {
        warning("Not enough custom colors provided, using default palette")
        fill_colors <- brewer.pal(max(n_groups, 3), "Set1")[1:n_groups]
    }
} else {
    if (n_groups <= 9) {
        fill_colors <- brewer.pal(9, "Set1")[1:n_groups]
    } else {
        fill_colors <- colorRampPalette(brewer.pal(9, "Set1"))(n_groups)
    }
}

# =============================================================================
# Color Palette (Gene) for stacked plots
# =============================================================================

n_genes_total <- length(gene_list)
if (n_genes_total <= 12) {
    gene_colors <- brewer.pal(12, "Set3")[1:n_genes_total]
} else {
    gene_colors <- colorRampPalette(brewer.pal(12, "Set3"))(n_genes_total)
}
names(gene_colors) <- gene_list

# =============================================================================
# Statistical Comparisons
# =============================================================================

if (!is.null(args$comparisons)) {
    comp_pairs <- strsplit(args$comparisons, ",")[[1]]
    comp_pairs <- lapply(comp_pairs, function(x) {
        pair <- strsplit(trimws(x), " vs |VS|vs\\.")[[1]]
        trimws(pair)
    })
    comp_pairs <- comp_pairs[sapply(comp_pairs, length) == 2]
} else {
    groups <- levels(df_long[[args$group_col]])
    if (length(groups) >= 2) {
        comp_pairs <- combn(groups, 2, simplify = FALSE)
    } else {
        comp_pairs <- list()
    }
}

# =============================================================================
# Common Theme
# =============================================================================

common_theme <- function() {
    theme_bw(base_size = args$font_size) +
    theme(
        plot.title = element_text(hjust = 0.5, face = "bold", size = args$font_size + 2),
        axis.text.x = element_text(angle = 90, hjust = 1, size = args$font_size),
        axis.text.y = element_text(size = args$font_size - 1),
        axis.title = element_text(face = "bold", size = args$font_size + 1),
        legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(color = "black", fill = NA, size = 0.8)
    )
}

# =============================================================================
# Plotting Functions
# =============================================================================

# ---- Single Gene Boxplot ----
plot_box_single <- function(gene_name, df_gene) {
    y_min_data <- min(df_gene$Expression, na.rm = TRUE)
    y_max_data <- max(df_gene$Expression, na.rm = TRUE)
    y_start <- args$y_start
    if (y_min_data < y_start) {
        y_start <- min(0, y_min_data * 0.95)
    }

    p <- ggplot(df_gene, aes(x = .data[[args$group_col]], y = Expression, fill = .data[[args$group_col]])) +
        geom_boxplot(width = 0.5, outlier.shape = NA, alpha = 0.7, lwd = 0.6) +
        stat_summary(fun = mean, geom = "point", shape = 23, size = 2.5,
                     fill = "white", color = "black", stroke = 0.8) +
        stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2, color = "black", lwd = 0.6)

    if (!args$no_jitter) {
        p <- p + geom_jitter(width = 0.15, size = 1.5, alpha = 0.6,
                             shape = 21, color = "black", stroke = 0.3)
    }

    if (args$paired && args$pair_col %in% colnames(meta)) {
        df_gene$pair_id <- meta[[args$pair_col]][match(df_gene$Sample, meta[[args$sample_col]])]
        p <- p + geom_line(aes(group = pair_id), color = "gray50", alpha = 0.4, linetype = "dashed")
    }

    # ---- 统计检验部分 ----
    if (args$test_method != "none" && length(comp_pairs) > 0) {
        if (args$test_method %in% c("wilcox.test", "t.test")) {
            p <- p + stat_compare_means(
                comparisons = comp_pairs,
                method = args$test_method,
                paired = args$paired,
                tip.length = 0.02,
                bracket.size = 0.5,
                label = "p.format",
                size = 3
            )
        } else if (args$test_method %in% c("anova", "kruskal.test")) {
            p <- p + stat_compare_means(
                method = args$test_method,
                label.y = max(df_page$Expression, na.rm = TRUE) * 1.1,
                size = 3
            )
        }
    }

    plot_title <- ifelse(is.null(args$title), gene_name, args$title)
    x_lab <- ifelse(is.null(args$x_label), args$group_col, args$x_label)

    p <- p + scale_fill_manual(values = fill_colors) +
        labs(title = plot_title, x = x_lab, y = args$y_label) +
        coord_cartesian(ylim = c(y_start, y_max_data * 1.15)) +
        common_theme()

    return(p)
}

# ---- Multi-gene Boxplot (one page) ----
plot_box_page <- function(page_genes, df_page, page_num, total_pages) {
    n_genes_page <- length(page_genes)
    ncol_facet <- ifelse(is.null(args$ncol), ceiling(sqrt(n_genes_page)), args$ncol)
    n_row <- ceiling(n_genes_page / ncol_facet)

    p <- ggplot(df_page, aes(x = .data[[args$group_col]], y = Expression,
                             fill = .data[[args$group_col]])) +
        geom_boxplot(width = 0.5, outlier.shape = NA, alpha = 0.7, lwd = 0.6) +
        stat_summary(fun = mean, geom = "point", shape = 23, size = 2,
                     fill = "white", color = "black") +
        stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2,
                     color = "black", lwd = 0.6)

    if (!args$no_jitter) {
        p <- p + geom_jitter(width = 0.15, size = 1.2, alpha = 0.5,
                             shape = 21, color = "black", stroke = 0.3)
    }

    if (args$test_method != "none" && n_groups == 2) {
        p <- p + stat_compare_means(method = args$test_method, paired = args$paired,
                                    label = "p.format", size = 3)
    }

    p <- p + facet_wrap(~ Gene, scales = "free_y", ncol = ncol_facet) +
        scale_fill_manual(values = fill_colors) +
        labs(x = ifelse(is.null(args$x_label), args$group_col, args$x_label),
             y = args$y_label) +
        theme_bw(base_size = args$font_size) +
        theme(
            strip.background = element_rect(fill = "lightgray", color = "black"),
            strip.text = element_text(face = "bold", size = args$font_size),
            axis.text.x = element_text(angle = 90, hjust = 1, size = args$font_size - 1),
            legend.position = "none",
            panel.grid.major.x = element_blank(),
            panel.grid.minor = element_blank()
        )

    plot_width <- max(args$width, ncol_facet * 3)
    plot_height <- max(args$height, n_row * args$row_height)
    plot_width <- min(plot_width, 49)
    plot_height <- min(plot_height, 49)

    suffix <- ifelse(total_pages > 1, paste0("_page", page_num), "")
    out_file <- file.path(args$outdir, paste0(args$outprefix, "_multi_gene", suffix, ".", args$format))
    ggsave(out_file, p, width = plot_width, height = plot_height, dpi = args$dpi, limitsize = FALSE)
    cat("Boxplot page saved to:", out_file, "\n")
}

# =============================================================================
# Generate Plots
# =============================================================================

if (!dir.exists(args$outdir)) {
    dir.create(args$outdir, recursive = TRUE)
}

# ---- Single Gene ----
if (length(gene_list) == 1) {
    df_gene <- df_long %>% filter(Gene == gene_list[1])

    # Boxplot
    p_box <- plot_box_single(gene_list[1], df_gene)
    out_box <- file.path(args$outdir, paste0(args$outprefix, "_", gene_list[1], ".", args$format))
    ggsave(out_box, p_box, width = args$width, height = args$height, dpi = args$dpi, limitsize = FALSE)
    cat("Boxplot saved to:", out_box, "\n")

# ---- Multi Gene (paginated) ----
} else {
    n_genes <- length(gene_list)
    genes_per_page <- args$genes_per_page
    n_pages <- ceiling(n_genes / genes_per_page)

    cat("Total genes:", n_genes, "| Genes per page:", genes_per_page, "| Total pages:", n_pages, "\n")

    for (page in 1:n_pages) {
        start_idx <- (page - 1) * genes_per_page + 1
        end_idx <- min(page * genes_per_page, n_genes)
        page_genes <- gene_list[start_idx:end_idx]
        df_page <- df_long %>% filter(Gene %in% page_genes)

        # Boxplot page
        plot_box_page(page_genes, df_page, page, n_pages)

    }

    # Stacked plots (use all genes, single figure)
    if (args$stacked) {
        plot_stacked(df_long)
    }
    if (args$stacked_percent) {
        plot_stacked_percent(df_long)
    }
}

# =============================================================================
# Save Data Table
# =============================================================================

out_table <- file.path(args$outdir, paste0(args$outprefix, "_data.tsv"))
write.table(df_long, out_table, sep = "\t", row.names = FALSE, quote = FALSE)
cat("Data table saved to:", out_table, "\n")

cat("\nDone!\n")