#!/usr/bin/env Rscript
chooseCRANmirror(ind=1)  # 选择第一个镜像，您可以根据实际情况选择
# 定义一个函数来检查包是否存在，如果不存在则安装，存在则加载
check_and_install <- function(package_name) {
  if (!requireNamespace(package_name, quietly = TRUE)) {
    message(paste(package_name, "is not installed. Installing now..."))
    if (package_name == "GenomicRanges" || package_name == "IRanges") {
      # 这两个包属于Bioconductor，需要使用BiocManager安装
      if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
      BiocManager::install(package_name, character.only = TRUE)
    } 
  } else {
    message(paste(package_name, "is already installed."))
  }
  library(package_name, character.only = TRUE)
}

# 检查并安装或加载包
check_and_install("IRanges")
check_and_install("GenomicRanges")

# 自定义函数，接受文件路径作为参数
# 自定义函数，接受文件路径作为参数
run_analysis <- function(data1_file, raw1_file, out_put_file) {
  # 读取数据
  data1 <- read.csv(file = data1_file, header = TRUE, sep = ",", stringsAsFactors = FALSE, comment.char = "#", check.names = FALSE)
  data1$chr <- paste0("chr", data1$chr)
  GR1 <- makeGRangesFromDataFrame(data1)
  
  raw1 <- read.table(file = raw1_file, header = TRUE, stringsAsFactors = FALSE, comment.char = "#", check.names = FALSE)
  names(raw1)[1:3] <- c("chr", "start", "end")
  sGR1 <- makeGRangesFromDataFrame(raw1)
  
  # findOverlaps
  hit <- findOverlaps(GR1, sGR1, type = "any", minoverlap = 2)
  qry <- data1[hit@from, ]
  sbj <- raw1[hit@to, ]
  res1 <- cbind(qry, sbj)
  colnames(res1) <- make.unique(colnames(res1))
  res1$chr.1 <- NULL
  res1$start.1 <- NULL
  res1$end.1 <- NULL
  res1$cM <- NULL
  
  custom_aggregate <- function(x) {
    list(
      mean = mean(x),
      max = max(x)
    )
  }
  
  res2 <- aggregate(cMperMb ~ ., data = res1, FUN = custom_aggregate)
  temp <- as.data.frame(lapply(as.data.frame(res2$cMperMb), unlist))
  res2 <- cbind(res2, temp)
  res2$cMperMb <- NULL
  
  # 写入文件时保留原始列名
  write.table(res2, out_put_file, sep = "\t", quote = FALSE, row.names = FALSE, col.names = TRUE)
}

# 获取命令行参数
args <- commandArgs(trailingOnly = TRUE)

# 检查参数数量
if (length(args) < 3) {
  stop("Usage: Rscript assigin_RecRate_for_hapSetSelected.R <data1.csv> <raw1.txt> <out_put_file>", call. = FALSE)
}

# 调用函数
run_analysis(args[1], args[2], args[3])