Expr_Boxplot_v2.r

#!/usr/bin/env Rscript

## 脚本说明
# 1) 运行环境: mamba activate r452
# 2) 应用:  查看基因表达

# =============================================================================
# bulkRNAseq Gene Expression Boxplot Script
# Author: Auto-generated
# Date: 2026-06-03
# =============================================================================

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]")
)

parser <- OptionParser(option_list = option_list,
                       usage = "Usage: %prog [options]
Description: Draw boxplot 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
# =============================================================================

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)
    }
}

# =============================================================================
# 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()
    }
}

# =============================================================================
# Plotting Function
# =============================================================================

plot_gene <- function(gene_name, df_gene) {

    # Y轴起始值设置
    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

    # 如果数据最小值小于y_start，自动调整
    if (y_min_data < y_start) {
        y_start <- min(0, y_min_data * 0.95)
    }

    # FIX: 使用 .data[[var]] 替代已弃用的 aes_string()
    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
        )

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

    # Add paired lines if requested
    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"
        )
    }

    # Statistical tests
    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.5
            )
        } else if (args$test_method %in% c("anova", "kruskal.test")) {
            p <- p + stat_compare_means(
                method = args$test_method,
                label.y = max(df_gene$Expression, na.rm = TRUE) * 1.1,
                size = 3.5
            )
        }
    }

    # Labels and theme
    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)) +
        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)
        )

    return(p)
}

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

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

# Single gene or multiple genes
if (length(gene_list) == 1) {
    df_gene <- df_long %>% filter(Gene == gene_list[1])
    p <- plot_gene(gene_list[1], df_gene)

    out_file <- file.path(args$outdir, paste0(args$outprefix, "_", gene_list[1], ".", args$format))

    ggsave(out_file, p, width = args$width, height = args$height, dpi = args$dpi)
    cat("Plot saved to:", out_file, "\n")

} else {
    # Multiple genes - facet
    # FIX: 使用 .data[[var]] 替代已弃用的 aes_string()
    p <- ggplot(df_long, 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)
    }

    # Add stats per facet if only 2 groups
    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)
    }

    # Facet
    ncol_facet <- ifelse(is.null(args$ncol), ceiling(sqrt(length(gene_list))), args$ncol)

    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()
        )

    out_file <- file.path(args$outdir, paste0(args$outprefix, "_multi_gene.", args$format))

    # Adjust dimensions for multi-panel
    n_row <- ceiling(length(gene_list) / ncol_facet)
    plot_width <- max(args$width, ncol_facet * 3)
    plot_height <- max(args$height, n_row * 3.5)
    
    # FIX: 限制最大尺寸不超过 49 英寸，并添加 limitsize = FALSE 防止报错
    plot_width <- min(plot_width, 49)
    plot_height <- min(plot_height, 49)

    ggsave(out_file, p, width = plot_width, height = plot_height, dpi = args$dpi, limitsize = FALSE)
    cat("Multi-gene plot saved to:", out_file, "\n")
}

# =============================================================================
# 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")