1. 准备文件

/home/gengxin/03.BatMeth2/genes.bed5

TEs.bed5

文件示例图：

2. 运行代码

需要修改的地方：

文件路径、选择运行的染色体号、图片标题

# =================== 加载包 ===================
library(dplyr)
library(ggplot2)
library(readr)

# =================== 1. 读取文件 ===================
genes <- read.delim("/home/gengxin/03.BatMeth2/genes.bed5", header = F, sep = "\t")
colnames(genes) <- c("chr","start","end","id","strand")

tes <- read.delim("TEs.bed5", header = F, sep = "\t")
colnames(tes) <- c("chr","start","end","type","strand")

# 只保留 1号染色体
genes <- filter(genes, chr == "LG08")
tes   <- filter(tes,   chr == "LG08")

window_size <- 100000  # 100kb窗口

# =================== 2. 计算密度 ===================
calc_density <- function(bed, window_size) {
  bed %>%
    group_by(chr) %>%
    do({
      chr_max <- max(.$end)
      windows <- seq(0, chr_max + window_size, by = window_size)
      count <- hist(.$start, breaks = windows, plot = FALSE)$count
      data.frame(
        chr = unique(.$chr),
        pos = windows[-1],
        density = count
      )
    })
}

gene_density <- calc_density(genes, window_size)
te_density <- calc_density(tes, window_size)
comb <- inner_join(gene_density, te_density, by=c("chr","pos"), suffix=c("_gene","_te"))

# =================== 3. 绘图：加上标题 ===================
ggplot(comb) +
  geom_line(aes(x=pos/1e6, y=density_gene, color="Gene"), linewidth=1) +
  geom_line(aes(x=pos/1e6, y=density_te/10, color="TE"), linewidth=1) +
  
  # 双轴刻度
  scale_y_continuous(
    name = "Gene Density",
    sec.axis = sec_axis(~ . *10, name = "TE Density")
  ) +
  
  # 颜色与图例
  scale_color_manual(
    name = NULL,
    values = c("Gene"="#377EB8", "TE"="#E41A1C")
  ) +
  
  labs(
    x = "Position (Mb)",
    title = "Chromosome 1"  # 标题
  ) +
  
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = "top",
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold")  # 标题居中、加粗
  )

# 保存
ggsave("chr1_gene_TE_density_with_title.pdf", width=12, height=5)

3. 结果