用r，WT_CpG.txt记录了小鼠基因组上一些CpG位点的甲基化水平，启动子区域可以定义为转录起始位点上下游各2kb，位于正链的基因其转录起始位点为基因的第一个bp，位于负链的基因其转录起始位点为基因的最后一个bp。输出mm9.main.refGene.txt中启动子区域甲基化水平最高的前十个基因。（mean[∑freqC/(freqC+freqT)]，不区分正负链）

以下是用R语言实现的代码：

# 读取数据
WT_CpG <- read.table("WT_CpG.txt", header = TRUE)

# 读取参考基因组注释文件
refGene <- read.table("mm9.main.refGene.txt", header = TRUE)

# 定义函数，计算甲基化水平
methylation_level <- function(freqC, freqT) {
  mean(freqC / (freqC + freqT))
}

# 定义函数，获取基因的启动子区域
get_promoter_region <- function(tx_start, tx_end, strand) {
  if (strand == "+") {
    promoter_start <- tx_start - 2000
    promoter_end <- tx_start + 2000
  } else {
    promoter_start <- tx_end - 2000
    promoter_end <- tx_end + 2000
  }
  return(c(promoter_start, promoter_end))
}

# 计算每个基因的启动子区域甲基化水平
promoter_methylation <- numeric(nrow(refGene))
for (i in 1:nrow(refGene)) {
  chr <- refGene$chrom[i]
  tx_start <- refGene$txStart[i]
  tx_end <- refGene$txEnd[i]
  strand <- refGene$strand[i]
  promoter_region <- get_promoter_region(tx_start, tx_end, strand)
  
  # 提取启动子区域内的CpG位点
  if (strand == "+") {
    CpG_in_promoter <- WT_CpG[WT_CpG$chrom == chr & WT_CpG$pos >= promoter_region[1] & 
                                 WT_CpG$pos <= promoter_region[2], ]
  } else {
    CpG_in_promoter <- WT_CpG[WT_CpG$chrom == chr & WT_CpG$pos <= promoter_region[2] & 
                                 WT_CpG$pos >= promoter_region[1], ]
  }
  
  # 计算甲基化水平
  if (nrow(CpG_in_promoter) > 0) {
    promoter_methylation[i] <- methylation_level(sum(CpG_in_promoter$freqC), sum(CpG_in_promoter$freqT))
  } else {
    promoter_methylation[i] <- 0
  }
}

# 找出启动子区域甲基化水平最高的前10个基因
top10_genes <- refGene[order(promoter_methylation, decreasing = TRUE), ][1:10, ]

解释一下代码的思路：

1. 首先读取数据，包括CpG位点甲基化水平文件和参考基因组注释文件。
2. 定义一个函数，用于计算甲基化水平，输入参数为C和T的频率，输出为甲基化水平。
3. 定义一个函数，用于获取基因的启动子区域，输入参数为基因的转录起始位点、终止位点和正负链信息，输出为基因的启动子区域范围。
4. 对于每个基因，调用get_promoter_region函数获取其启动子区域范围，然后从CpG位点文件中提取该范围内的CpG位点。
5. 对于每个基因的启动子区域，调用methylation_level函数计算其甲基化水平。
6. 最后，根据启动子区域甲基化水平对所有基因进行排序，选出甲基化水平最高的前10个基因。