1. 准备文件

/home/gengxin/ZLbud_1.methratio.txt

2. 运行命令

Rscript meth_depth_analysis.R /home/gengxin/ZLbud_1.methratio.txt /home/gengxin/result

#!/usr/bin/env Rscript

# 功能：筛选8条染色体（LG01-LG08），浅色系配色绘制甲基化水平柱状图并输出汇总表
# 运行方式：Rscript meth_depth_analysis.R <输入methratio.txt文件路径> <输出文件前缀>

# 1. 获取并校验命令行参数
args <- commandArgs(trailingOnly = TRUE)
if (length(args) != 2) {
  stop("使用格式：Rscript meth_depth_analysis.R <输入文件路径> <输出文件前缀>\n示例：Rscript meth_depth_analysis.R /home/gengxin/ZLbud_1.methratio.txt /home/gengxin/result")
}

# 提取输入输出路径
input_file <- args[1]
output_prefix <- args[2]

# 2. 静默加载依赖包
suppressPackageStartupMessages({
  library(tidyverse)
})

# 3. 读取并预处理甲基化数据
tryCatch({
  methratio <- read_delim(
    file = input_file,
    delim = "\t",
    escape_double = FALSE,
    trim_ws = TRUE
  ) %>%
    dplyr::select(
      chr = `#chromsome`,
      loci,
      strand,
      context,
      C_count,
      CT_count,
      methRatio,
      MethContext,
      FiveContext = `5context`
    ) %>%
    mutate(context = if_else(context == "CG", "CpG", context)) %>%
    mutate(context = factor(context, levels = c("CpG", "CHG", "CHH")))
}, error = function(e) {
  stop(paste("读取输入文件失败：", e$message))
})

# ============= 核心修改1：筛选LG01-LG08这8条染色体的数据 =============
message(paste("原始数据总行数：", nrow(methratio)))
# 定义目标染色体列表
target_chrs <- c("LG01", "LG02", "LG03", "LG04", "LG05", "LG06", "LG07", "LG08")
# 筛选数据（仅保留8条染色体）
methratio_filtered <- methratio %>%
  filter(chr %in% target_chrs)
# 校验筛选后数据
if (nrow(methratio_filtered) == 0) {
  stop("筛选后无有效数据，请检查染色体名称是否与数据一致！")
}
message(paste("8条染色体筛选后数据总行数：", nrow(methratio_filtered)))

# 4. 执行数据汇总（基于筛选后的数据）
# 4.1 全局汇总（无分组，仅8条染色体）
methratio_global <- summarise(
  methratio_filtered,
  count = n(),
  methRatio = round(mean(methRatio, na.rm = TRUE), digits = 3)
)
message("8条染色体全局甲基化汇总结果：")
print(methratio_global)

# 4.2 按chr+context分组汇总（仅8条染色体）
methratio_sum <- methratio_filtered %>%
  group_by(chr, context) %>%
  summarise(
    count = n(),
    methRatio = round(mean(methRatio, na.rm = TRUE), digits = 3),
    .groups = "drop"
  )
# 预览前6行分组汇总结果
message("\n8条染色体分组甲基化汇总前6行：")
print(head(methratio_sum))

# 5. 保存汇总表（无表头）
# 保存分组汇总表
sum_file <- paste0(output_prefix, "_methratio_group_sum.txt")
write_delim(methratio_sum, file = sum_file, delim = "\t", col_names = FALSE)
message(paste("\n8条染色体分组汇总表已保存至：", sum_file))

# 保存全局汇总表
global_sum_file <- paste0(output_prefix, "_methratio_global_sum.txt")
write_delim(methratio_global, file = global_sum_file, delim = "\t", col_names = FALSE)
message(paste("8条染色体全局汇总表已保存至：", global_sum_file))

# ============= 核心修改2：定义8种浅色系配色（适配8条染色体） =============
light_colors <- c(
  "#F8C8DC",  # 浅粉
  "#C1E1C1",  # 浅绿
  "#B0E0E6",  # 浅蓝
  "#eeeaeb",  # 浅黄
  "#E6E6FA",  # 浅紫
  "#FFE4B5",  # 浅橙
  "#cee9dc",  # 浅薄荷绿
  "#FFE4E1"   # 浅珊瑚粉
)

# 6. 绘制柱状图（应用浅色系配色，基于8条染色体数据）
plot_file <- paste0(output_prefix, "_methRatio_barplot.pdf")
p <- ggplot(data = methratio_sum, aes(x = chr, y = methRatio)) +
  geom_col(
    aes(fill = chr),
    position = "dodge",
    color = "white",  # 白色边框，增强浅色柱状图辨识度
    width = 0.7,
    show.legend = F
  ) +
  # 应用自定义浅色系配色
  scale_fill_manual(values = light_colors) +
  geom_text(
    aes(y = methRatio * 1.05, label = methRatio),
    angle = 90,
    vjust = 0.5,
    size = 3
  ) +
  labs(x = NULL, y = 'Methylation Levels', fill = NULL) +
  facet_wrap(~context, scales = "free") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

# 保存图片（适配8条染色体分面展示）
ggsave(
  filename = plot_file,
  plot = p,
  width = 14,
  height = 8,
  dpi = 300
)
message(paste("8条染色体浅色系甲基化柱状图已保存至：", plot_file))

# 7. 运行完成提示
message("\n所有分析完成！")

3. 结果文件