分組的圖
# ==============================
# 分組 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. 循環(huán)處理每個(gè) 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"
? # 過(guò)濾非標(biāo)準(zhǔn)染色體
? 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]
? # 導(dǎo)出處理后的 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
)
# 美化：白色背景 + 坐標(biāo)軸清晰 + 字體調(diào)整
p <- p +
? ? theme_bw(base_size = 14) +? # 白色背景，基礎(chǔ)字體 14
? ? theme(
? ? ? ? panel.grid.major = element_line(color = "grey90"), # 主網(wǎng)格線淺灰
? ? ? ? panel.grid.minor = element_blank(),? ? ? ? ? ? ? ? # 去掉次網(wǎng)格線
? ? ? ? axis.text = element_text(color = "black"),? ? ? ? # 坐標(biāo)軸文字黑色
? ? ? ? axis.title = element_text(face = "bold"),? ? ? ? ? # 坐標(biāo)軸標(biāo)題加粗
? ? ? ? legend.position = "right",? ? ? ? ? ? ? ? ? ? ? ? # 圖例放右邊
? ? ? ? legend.title = element_text(face = "bold")? ? ? ? # 圖例標(biāo)題加粗
? ? )
print(p)
# 4. 輸出 PDF
pdf(file.path(folder_path, "Sham_vs_SNI_metagene_plot.pdf"), width=8, height=6)
print(p)
dev.off()
六個(gè)樣本不分組的圖
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) 循環(huán)處理每個(gè) 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"
? # 過(guò)濾非標(biāo)準(zhǔn)染色體、長(zhǎng)度>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
? # 導(dǎo)出處理后的 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) 讀取轉(zhuǎn)錄組數(shù)據(jù)庫(kù)
# ==============================
txdb_file <- system.file("extdata", "mm10_toy.sqlite", package="Guitar")
txdb <- loadDb(txdb_file)
# ==============================
# 4) 構(gòu)建 GuitarPlot 對(duì)象
# ==============================
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")
最后美化了一波
# ========== 三種風(fēng)格對(duì)比 ==========
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)
用的圖是
# 極簡(jiǎn)主題（去掉淡灰網(wǎng)格線）
p_minimal <- p +
? theme_minimal(base_size = 14) +
? theme(
? ? panel.grid.major = element_blank(),? # 去掉主網(wǎng)格線
? ? panel.grid.minor = element_blank()? # 去掉次網(wǎng)格線
? )
ggsave("guitar_minimal.pdf", p_minimal, width = 7, height = 5)
