分组的图
# ==============================
# 分组 metagene plot 脚本（Sham vs SNI）
# ==============================
library(Guitar)
library(rtracklayer)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(GenomicRanges)
# 1. 定义文件夹和文件
folder_path <- "F:/data/ac4crip"
sham_files <- c("Sham2_summits.bed", "Sham3_summits.bed", "Sham4_summits.bed")
sni_files <- c("SNI1_summits.bed", "SNI2_summits.bed", "SNI4_summits.bed")
all_files <- c(sham_files, sni_files)
all_groups <- c(rep("Sham", length(sham_files)), rep("SNI", length(sni_files)))
# 2. 循环处理每个 BED 文件
bed_paths <- c()
for (i in seq_along(all_files)) {
file_name <- all_files[i]
group_name <- all_groups[i]
bed_path <- file.path(folder_path, file_name)
bed <- import(bed_path, format="bed")
# 添加 chr 前缀
bed_ucsc <- GRanges(
seqnames = paste0("chr", seqnames(bed)),
ranges = ranges(bed),
strand = strand(bed),
mcols = mcols(bed)
)
# chrMT -> chrM
seqlevels(bed_ucsc)[seqlevels(bed_ucsc) == "chrMT"] <- "chrM"
# 过滤非标准染色体
valid_chr <- seqlevels(TxDb.Mmusculus.UCSC.mm10.knownGene)
valid_chr <- valid_chr[grep("^chr[0-9XYM]+$", valid_chr)]
bed_filtered <- bed_ucsc[seqnames(bed_ucsc) %in% valid_chr]
# 导出处理后的 BED 文件
output_file <- file.path(folder_path, sub("summits.bed", "summits_UCSC_final.bed", file_name))
export(bed_filtered, con = output_file)
bed_paths <- c(bed_paths, output_file)
}
# 3. 绘制分组 metagene plot
p <- GuitarPlot(
txTxdb = TxDb.Mmusculus.UCSC.mm10.knownGene,
stBedFiles = bed_paths,
headOrtail = FALSE,
enableCI = FALSE,
mapFilterTranscript = TRUE,
pltTxType = "mrna",
stGroupName = all_groups
)
# 美化：白色背景 + 坐标轴清晰 + 字体调整
p <- p +
theme_bw(base_size = 14) + # 白色背景，基础字体 14
theme(
panel.grid.major = element_line(color = "grey90"), # 主网格线浅灰
panel.grid.minor = element_blank(), # 去掉次网格线
axis.text = element_text(color = "black"), # 坐标轴文字黑色
axis.title = element_text(face = "bold"), # 坐标轴标题加粗
legend.position = "right", # 图例放右边
legend.title = element_text(face = "bold") # 图例标题加粗
)
print(p)
# 4. 输出 PDF
pdf(file.path(folder_path, "Sham_vs_SNI_metagene_plot.pdf"), width=8, height=6)
print(p)
dev.off()
六个样本不分组的图
library(Guitar)
library(rtracklayer)
library(GenomicRanges)
# ==============================
# 1) 文件路径和样本
# ==============================
folder_path <- "F:/data/"
sample_files <- c(
"Sham2_summits.bed", "Sham3_summits.bed", "Sham4_summits.bed",
"SNI1_summits.bed", "SNI2_summits.bed", "SNI4_summits.bed"
)
sample_names <- c("Sham2","Sham3","Sham4","SNI1","SNI2","SNI4")
bed_paths <- c()
# ==============================
# 2) 循环处理每个 BED 文件
# ==============================
for (i in seq_along(sample_files)) {
bed_path <- file.path(folder_path, sample_files[i])
bed <- import(bed_path, format="bed")
# 添加 chr 前缀
bed_ucsc <- GRanges(
seqnames = paste0("chr", seqnames(bed)),
ranges = ranges(bed),
strand = strand(bed),
mcols = mcols(bed)
)
seqlevels(bed_ucsc)[seqlevels(bed_ucsc)=="chrMT"] <- "chrM"
# 过滤非标准染色体、长度>0、非 NA
valid_chr <- seqlevels(TxDb.Mmusculus.UCSC.mm10.knownGene)
valid_chr <- valid_chr[grep("^chr[0-9XYM]+$", valid_chr)]
bed_filtered <- bed_ucsc[
seqnames(bed_ucsc) %in% valid_chr &
width(bed_ucsc) > 0 &
!is.na(start(bed_ucsc)) &
!is.na(end(bed_ucsc))
]
# 去掉 mcols NA
if (length(bed_filtered) > 0) {
mcols(bed_filtered) <- mcols(bed_filtered)[
rowSums(is.na(as.data.frame(mcols(bed_filtered))))==0, , drop=FALSE
]
}
if (length(bed_filtered) == 0) next
# 导出处理后的 BED 文件
output_file <- file.path(folder_path, sub("summits.bed","summits_UCSC_final.bed",sample_files[i]))
export(bed_filtered, con=output_file)
bed_paths <- c(bed_paths, output_file)
}
# ==============================
# 3) 读取转录组数据库
# ==============================
txdb_file <- system.file("extdata", "mm10_toy.sqlite", package="Guitar")
txdb <- loadDb(txdb_file)
# ==============================
# 4) 构建 GuitarPlot 对象
# ==============================
guitar_obj <- GuitarPlot(
txTxdb = txdb,
stBedFiles = bed_paths, # 必须是文件路径
stGroupName = sample_names,
pltTxType = "mrna",
headOrtail = FALSE,
enableCI = FALSE,
mapFilterTranscript = TRUE
)
# ==============================
# 5) 输出 PDF
# ==============================
library(ggplot2)
ggsave("6samples_metagene_plot.pdf", plot = guitar_obj, width = 10, height = 6, device = "pdf")
dev.off()
cat("✅ 6 samples metagene plot saved to PDF successfully!\n")
最后美化了一波
# ========== 三种风格对比 ==========
p_bw <- p + theme_bw() +
theme(panel.grid = element_line(color = "grey90"),
plot.title = element_text(hjust = 0.5, size = 14))
p_classic <- p + theme_classic() +
theme(axis.text = element_text(size = 12),
axis.title = element_text(size = 14))
p_minimal <- p + theme_minimal() +
theme(panel.grid = element_line(color = "grey85", linetype = "dashed"),
axis.text = element_text(size = 12))
# ========== 输出 ==========
# 分别保存三张图
ggsave("guitar_bw.pdf", p_bw, width = 7, height = 5)
ggsave("guitar_classic.pdf", p_classic, width = 7, height = 5)
ggsave("guitar_minimal.pdf", p_minimal, width = 7, height = 5)
用的图是
# 极简主题（去掉淡灰网格线）
p_minimal <- p +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_blank(), # 去掉主网格线
panel.grid.minor = element_blank() # 去掉次网格线
)
ggsave("guitar_minimal.pdf", p_minimal, width = 7, height = 5)
