Expr_Heatmap.r

#!/usr/bin/env Rscript

# =============================================================================
# 脚本功能：基于数值矩阵和样本分组信息绘制ComplexHeatmap热图
# 输入：
#   - 数值矩阵：行=基因/剪接事件，列=样本，制表符分隔
#   - 分组文件：两列，制表符分隔，第一列样本名（必须与矩阵列名匹配），第二列分组标签
# 输出：
#   - PDF格式热图，上调（高值）红色，下调（低值）蓝色
# 依赖包：ComplexHeatmap, circlize, optparse, RColorBrewer
# =============================================================================

# 加载或安装所需包
# suppressPackageStartupMessages({
#   if (!require("optparse", quietly = TRUE)) install.packages("optparse", repos = "https://cloud.r-project.org/")
#   if (!require("ComplexHeatmap", quietly = TRUE)) {
#     if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
#     BiocManager::install("ComplexHeatmap")
#   }
#   if (!require("circlize", quietly = TRUE)) install.packages("circlize", repos = "https://cloud.r-project.org/")
#   if (!require("RColorBrewer", quietly = TRUE)) install.packages("RColorBrewer", repos = "https://cloud.r-project.org/")
# })

library(optparse)
library(ComplexHeatmap)
library(circlize)
library(RColorBrewer)

# -----------------------------------------------------------------------------
# 解析命令行参数
# -----------------------------------------------------------------------------
option_list <- list(
  make_option(c("-m", "--matrix"), type = "character", default = NULL,
              help = "数值矩阵文件（制表符分隔，行名=基因/事件，列名=样本）", metavar = "file"),
  make_option(c("-l", "--genelist"), type = "character", default = "all",
              help = "需要绘制基因（逗号分隔 [默认 %default]）", metavar = "string"),
  make_option(c("-s", "--showgp"), type = "character", default = "all",
              help = "热图展示的分组，逗号分隔", metavar = "string"),
  make_option(c("-g", "--group"), type = "character", default = NULL,
              help = "样本分组文件（制表符分隔，第一列样本名，第二列分组名）", metavar = "file"),
  make_option(c("-o", "--out"), type = "character", default = "./",
              help = "输出图片路径 [默认 %default]", metavar = "path"),
  make_option(c("-p", "--prefix"), type = "character", default = "no",
              help = "输出图片文件名前前缀 [默认 %default]", metavar = "string"),
  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 = "Heatmap",
              help = "热图标题 [默认 %default]", metavar = "string"),
  make_option(c("--cluster_rows"), action = "store_true", default = FALSE,
              help = "是否对行聚类 [默认 %default]"),
  make_option(c("--cluster_cols"), action = "store_true", default = FALSE,
              help = "是否对列聚类 [默认 %default]"),
  make_option(c("--show_rownames"), action = "store_true", default = FALSE,
              help = "是否显示行名 [默认 %default]"),
  make_option(c("--show_colnames"), action = "store_true", default = FALSE,
              help = "是否显示列名 [默认 %default]"),
  make_option(c("--legend_title"), type = "character", default = "Value",
              help = "图例标题 [默认 %default]", metavar = "string")
)

opt_parser <- OptionParser(option_list = option_list, description = "绘制带样本分组的ComplexHeatmap热图")
opt <- parse_args(opt_parser)

# 检查必需参数
if (is.null(opt$matrix)) {
  print_help(opt_parser)
  stop("错误：必须提供数值矩阵文件（--matrix）", call. = FALSE)
}
if (is.null(opt$genelist)) {
  print_help(opt_parser)
  stop("错误：必须提供绘制基因（--genelist）", call. = FALSE)
}
if (is.null(opt$group)) {
  print_help(opt_parser)
  stop("错误：必须提供样本分组文件（--group）", call. = FALSE)
}
prefix <- opt$prefix
print(prefix)
# -----------------------------------------------------------------------------
# 读取数据
# -----------------------------------------------------------------------------
# 读取数值矩阵（假设第一列为行名）
# 读取矩阵（不将0转为NA，自动处理分隔符）
mat_df <- tryCatch({
  read.csv(opt$matrix, header = TRUE, row.names = 1, check.names = FALSE, 
           stringsAsFactors = FALSE)
}, error = function(e) {
  stop("读取数值矩阵失败：", e$message)
})
print(head(mat_df))

# 强制转换为数值矩阵（非数值转为NA）
mat <- as.matrix(mat_df)
if (!is.numeric(mat)) {
  mat <- apply(mat_df, 2, as.numeric)
  rownames(mat) <- rownames(mat_df)
  colnames(mat) <- colnames(mat_df)
  warning("矩阵中存在非数值数据，已转换为NA，请检查原始文件")
}
if (all(is.na(mat))) stop("矩阵全为NA，无法绘图")
# 获取基因列表
genelist <- opt$genelist
if(genelist != "all"){
  genelist <- unlist(strsplit(genelist,","))
  diff <- setdiff(genelist,rownames(mat))
  if(length(diff) > 0){
    print(diff)
    genelist <- intersect(genelist,rownames(mat))
  }
  mat <- mat[genelist,]
}
print(genelist)

read_group_file <- function(file) {
  first_line <- readLines(file, n = 1)
  sep <- if (grepl("\t", first_line)) "\t" else ","
  
  # 尝试有表头
  df <- tryCatch({
    tmp <- read.table(file, header = TRUE, sep = sep, stringsAsFactors = FALSE)
    if (ncol(tmp) >= 2) tmp[, 1:2] else stop("列数不足")
  }, error = function(e) NULL)
  
  if (!is.null(df)) {
    colnames(df)[1:2] <- c("sample", "group")
    return(df)
  }
  
  # 无表头
  df <- tryCatch({
    read.table(file, header = FALSE, sep = sep, stringsAsFactors = FALSE, 
               col.names = c("sample", "group"))
  }, error = function(e) stop("分组文件解析失败"))
  df
}

# 读入样本分组文件
gp_file <- read_group_file(opt$group)
# 筛选展示的热图
if (opt$showgp != "all"){
  showgp <- unlist(strsplit(opt$showgp,","))
}else{
  showgp <- unique(gp_file$group)
}
group_df <- gp_file[which(gp_file$group %in% showgp),]
# 设定分组顺序
group_df$group <- factor(group_df$group, levels = showgp)
group_df <- group_df[order(group_df$group),]
print(group_df)

# 检查矩阵列名是否都在分组文件中
samples_in_mat <- colnames(mat)
missing <- setdiff(samples_in_mat, group_df$sample)
if (length(missing) > 0) {
    message("⚠️ 矩阵中的样本在分组文件中缺失: ", paste(missing, collapse = ", "))
} else {
    message("✅ 所有矩阵样本都在分组文件中")
}

mat <- mat[,intersect(group_df$sample, samples_in_mat), drop = FALSE]
samples_in_mat <- colnames(mat)

print(head(mat))
print(samples_in_mat)

# 定义行归一化函数：将每一行缩放到 [-1, 1]
row_scale_to_minus1_1 <- function(mat) {
      row_names <- rownames(mat)
      col_names <- colnames(mat)
      scaled <- t(apply(mat, 1, function(x) {
        min_x <- min(x, na.rm = TRUE)
        max_x <- max(x, na.rm = TRUE)
        if (min_x == max_x) return(rep(0, length(x)))
        else return(2 * (x - min_x) / (max_x - min_x) - 1)
      }))
      rownames(scaled) <- row_names
      colnames(scaled) <- col_names
      return(scaled)
}

# -----------------------------------------------------------------------------
# 准备热图注释
# -----------------------------------------------------------------------------
# 定义分组颜色（自动为每个分组分配颜色）
groups <- levels(group_df$group)
n_groups <- length(groups)
if (n_groups == 2){
    group_colors <- setNames(c("#2166AC", "#B2182B"), groups)
}else if(n_groups <= 8) {
    group_colors <- setNames(brewer.pal(n_groups, "Set1"), groups)
} else {
    group_colors <- setNames(colorRampPalette(brewer.pal(8, "Set1"))(n_groups), groups)
}
print(group_colors)

# 创建列注释
col_ha <- HeatmapAnnotation(
    Group = group_df$group,
    col = list(Group = group_colors),
    annotation_name_side = "right",  # 注释名称位置
    annotation_legend_param = list(
        Group = list(title = "Group", title_gp = gpar(fontsize = 10))
    )
)

# -----------------------------------------------------------------------------
# 绘制热图
# -----------------------------------------------------------------------------
# 设置行聚类、列聚类参数
cluster_rows <- opt$cluster_rows
cluster_cols <- opt$cluster_cols

# 对于行聚类，若行数过少（<2）则自动禁用
if (nrow(mat) < 2) cluster_rows <- FALSE
if (ncol(mat) < 2) cluster_cols <- FALSE

# 为使正负值区分明显，若数据有正有负，以0为中心；否则以数据范围中心为界
if (min(mat, na.rm = TRUE) < 0 & max(mat, na.rm = TRUE) > 0) {
  col_fun <- colorRamp2(c(min(mat, na.rm = TRUE), 0, max(mat, na.rm = TRUE)), 
                        c("#2166AC", "white", "#B2182B"))  # 蓝色 -> 白色 -> 红色
} else {
  # 数据全为非负或全为非正时，使用连续渐变
  if (min(mat, na.rm = TRUE) >= 0) {
    col_fun <- colorRamp2(c(min(mat, na.rm = TRUE), max(mat, na.rm = TRUE)), 
                          c("#2166AC", "#B2182B"))  # 低值蓝，高值红
  } else {
    col_fun <- colorRamp2(c(min(mat, na.rm = TRUE), max(mat, na.rm = TRUE)), 
                          c("#B2182B", "#2166AC"))  # 低值红，高值蓝（负值越低越红）
  }
}

# 构建热图对象
ht <- Heatmap(mat,
              name = opt$legend_title,
              col = col_fun,
              row_title = "Features",
              column_title = opt$title,
              column_title_gp = gpar(fontsize = 10, fontface = "bold"),
              cluster_rows = cluster_rows,
              cluster_columns = cluster_cols,
              show_row_names = opt$show_rownames,
              show_column_names = opt$show_colnames,
              top_annotation = col_ha,
              heatmap_legend_param = list(
                title_gp = gpar(fontsize = 10, fontface = "bold"),
                labels_gp = gpar(fontsize = 9)
              ),
              # 可选：调整行名、列名大小
              row_names_gp = gpar(fontsize = 9),
              column_names_gp = gpar(fontsize = 9, angle = 45, hjust = 1))

# 在 mat 筛选完基因后，绘制热图之前，调用此函数
mat_scaled <- row_scale_to_minus1_1(mat)
print(head(mat_scaled))
# 为缩放热图单独设置颜色映射（范围固定为 -1 到 1）
if (min(mat_scaled, na.rm = TRUE) < 0 & max(mat_scaled, na.rm = TRUE) > 0) {
  col_fun <- colorRamp2(c(min(mat_scaled, na.rm = TRUE), 0, max(mat_scaled, na.rm = TRUE)), 
                        c("#2166AC", "white", "#B2182B"))  # 蓝色 -> 白色 -> 红色
} else {
  # 数据全为非负或全为非正时，使用连续渐变
  if (min(mat_scaled, na.rm = TRUE) >= 0) {
    col_fun <- colorRamp2(c(min(mat_scaled, na.rm = TRUE), max(mat_scaled, na.rm = TRUE)), 
                          c("#2166AC", "#B2182B"))  # 低值蓝，高值红
  } else {
    col_fun <- colorRamp2(c(min(mat_scaled, na.rm = TRUE), max(mat_scaled, na.rm = TRUE)), 
                          c("#B2182B", "#2166AC"))  # 低值红，高值蓝（负值越低越红）
  }
}

ht_scale <- Heatmap(mat_scaled,
                    name = opt$legend_title,
                    col = col_fun,
                    row_title = "Features",
                    column_title = opt$title,
                    column_title_gp = gpar(fontsize = 10, fontface = "bold"),
                    cluster_rows = cluster_rows,
                    cluster_columns = cluster_cols,
                    show_row_names = opt$show_rownames,
                    show_column_names = opt$show_colnames,
                    top_annotation = col_ha,
                    heatmap_legend_param = list(
                      at = c(-1, 0, 1),
                      labels = c("-1", "0", "1"),
                      title = opt$legend_title,
                      title_gp = gpar(fontsize = 10, fontface = "bold"),
                      labels_gp = gpar(fontsize = 9)
                    ),
                    row_names_gp = gpar(fontsize = 9),
                    column_names_gp = gpar(fontsize = 9, angle = 45, hjust = 1))

# -----------------------------------------------------------------------------
# 保存图片
# -----------------------------------------------------------------------------
# 创建输出目录（如果不存在）
out_dir <- opt$out
if (!dir.exists(out_dir)) {
  dir.create(out_dir, recursive = TRUE)
}

if(prefix != "no"){
  pdf(file.path(out_dir,paste0(prefix,"_Heatmap.pdf")), width = opt$width, height = opt$height)
  draw(ht, annotation_legend_side = "right", heatmap_legend_side = "right")
  invisible(dev.off())

  pdf(file.path(out_dir,paste0(prefix,"_Heatmap_scale.pdf")), width = opt$width, height = opt$height)
  draw(ht_scale, annotation_legend_side = "right", heatmap_legend_side = "right")
  invisible(dev.off())
}else{
  pdf(file.path(out_dir,"Heatmap.pdf"), width = opt$width, height = opt$height)
  draw(ht, annotation_legend_side = "right", heatmap_legend_side = "right")
  invisible(dev.off())

  pdf(file.path(out_dir,"Heatmap_scale.pdf"), width = opt$width, height = opt$height)
  draw(ht_scale, annotation_legend_side = "right", heatmap_legend_side = "right")
  invisible(dev.off())
}
message("热图已成功保存至：", out_dir)