所有bug暫且不表许布，只記錄有用的
環(huán)境&包，因?yàn)槭切码娔X，全是新裝的
if (!require("BiocManager", quietly = TRUE))
? ? install.packages("BiocManager")
BiocManager::install(c(
? ? "GenomicRanges",? ? ? # 用于處理基因組范圍數(shù)據(jù)
? ? "rtracklayer",? ? ? ? # 用于導(dǎo)入/導(dǎo)出BED、GTF等文件
? ? "TxDb.Mmusculus.UCSC.mm10.knownGene",? # mm10基因組注釋
? ? "Guitar"? ? ? ? ? ? ? # 用于繪制metagene plot
))
install.packages(c(
? ? "tidyverse",? ? ? ? ? # 數(shù)據(jù)整理和可視化
? ? "ggplot2",? ? ? ? ? ? # 高級(jí)繪圖
? ? "data.table",? ? ? ? ? # 高效數(shù)據(jù)處理
? ? "devtools"? ? ? ? ? ? # 用于從GitHub安裝包
))
library(GenomicRanges)
library(rtracklayer)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(Guitar)
library(tidyverse)
SNI的，主要bug是mm10.fa的染色質(zhì)序號(hào)沒有chr前綴先舷，
BED文件中的線粒體染色體名稱為?MT，而?TxDb?中使用?M滓侍，需要手動(dòng)修正
以及存在非標(biāo)準(zhǔn)染色體
報(bào)錯(cuò)信息大概是 1: In .merge_two_Seqinfo_objects(x, y) : Each of the 2 combined objects has sequence levels not in the other: - in 'x': chrGL456233.2, chrJH584304.1, chrMT, chrMU069435.1
# 1. 讀取BED文件
bed_path_sni <- "F:/XZHMU/huangyue/motif/RIP_SNI_peaks_summits.bed"?
bed_sni <- import(bed_path_sni, format = "bed")? # 明確指定文件格式為bed
# 2. 修正染色體名稱（添加chr前綴）
bed_ucsc_sni <- GRanges(
? seqnames = paste0("chr", seqnames(bed_sni)),
? ranges = ranges(bed_sni),
? strand = strand(bed_sni),
? mcols = mcols(bed_sni)
)
# 3. 統(tǒng)一線粒體染色體名稱（chrMT -> chrM）
seqlevels(bed_ucsc_sni)[seqlevels(bed_ucsc_sni) == "chrMT"] <- "chrM"
# 4. 過濾非標(biāo)準(zhǔn)染色體（僅保留chr1-chr19, chrX, chrY, chrM）
valid_chr <- seqlevels(TxDb.Mmusculus.UCSC.mm10.knownGene)
valid_chr <- valid_chr[grep("^chr[0-9XYM]+$", valid_chr)]
bed_filtered_sni <- bed_ucsc_sni[seqnames(bed_ucsc_sni) %in% valid_chr]
# 5. 導(dǎo)出為新的BED文件
export(bed_filtered_sni, "F:/XZHMU/huangyue/motif/RIP_SNI_peaks_summits_UCSC_final.bed")
# 6. 運(yùn)行GuitarPlot
p_sni <- GuitarPlot(
? txTxdb = TxDb.Mmusculus.UCSC.mm10.knownGene,
? stBedFiles = "F:/XZHMU/huangyue/motif/RIP_SNI_peaks_summits_UCSC_final.bed",?
? headOrtail = FALSE,
? enableCI = FALSE,
? mapFilterTranscript = TRUE,
? pltTxType = "mrna",
? stGroupName = "SNI"
)
# 7. 顯示結(jié)果
print(p_sni)
SH
# 1. 讀取BED文件bed_path_sh <- "F:/XZHMU/huangyue/motif/RIP_sh_peaks_summits.bed"?
bed_sh <- import(bed_path_sh, format = "bed")? # 明確指定文件格式為bed
# 2. 修正染色體名稱（添加chr前綴）
bed_ucsc_sh <- GRanges(
? seqnames = paste0("chr", seqnames(bed_sh)),
? ranges = ranges(bed_sh),
? strand = strand(bed_sh),
? mcols = mcols(bed_sh)
)
# 3. 統(tǒng)一線粒體染色體名稱（chrMT -> chrM）
seqlevels(bed_ucsc_sh)[seqlevels(bed_ucsc_sh) == "chrMT"] <- "chrM"
# 4. 過濾非標(biāo)準(zhǔn)染色體（僅保留chr1-chr19, chrX, chrY, chrM）
valid_chr <- seqlevels(TxDb.Mmusculus.UCSC.mm10.knownGene)
valid_chr <- valid_chr[grep("^chr[0-9XYM]+$", valid_chr)]
bed_filtered_sh <- bed_ucsc_sh[seqnames(bed_ucsc_sh) %in% valid_chr]
# 5. 導(dǎo)出為新的BED文件
export(bed_filtered_sh, "F:/XZHMU/huangyue/motif/RIP_sh_peaks_summits_UCSC_final.bed")
# 6. 運(yùn)行GuitarPlot
p_sh <- GuitarPlot(
? txTxdb = TxDb.Mmusculus.UCSC.mm10.knownGene,
? stBedFiles = "F:/XZHMU/huangyue/motif/RIP_sh_peaks_summits_UCSC_final.bed",? # 去掉多余的空格
? headOrtail = FALSE,
? enableCI = FALSE,
? mapFilterTranscript = TRUE,
? pltTxType = "mrna",
? stGroupName = "sh"
)
# 7. 顯示結(jié)果
print(p_sh)
SHAM
# 1. 讀取BED文件
bed_path_sham <- "F:/XZHMU/huangyue/motif/RIP_Sham_peaks_summits.bed"?
bed_sham <- import(bed_path_sham)
# 2. 修正染色體名稱（添加chr前綴）
bed_ucsc_sham <- GRanges(
? seqnames = paste0("chr", seqnames(bed_sham)),
? ranges = ranges(bed_sham),
? strand = strand(bed_sham),
? mcols = mcols(bed_sham)
)
# 3. 統(tǒng)一線粒體染色體名稱（chrMT -> chrM）
seqlevels(bed_ucsc_sham)[seqlevels(bed_ucsc_sham) == "chrMT"] <- "chrM"
# 4. 過濾非標(biāo)準(zhǔn)染色體（僅保留chr1-chr19, chrX, chrY, chrM）
valid_chr <- seqlevels(TxDb.Mmusculus.UCSC.mm10.knownGene)
valid_chr <- valid_chr[grep("^chr[0-9XYM]+$", valid_chr)]
bed_filtered_sham <- bed_ucsc_sham[seqnames(bed_ucsc_sham) %in% valid_chr]
# 5. 導(dǎo)出為新的BED文件
export(bed_filtered_sham, "F:/XZHMU/huangyue/motif/RIP_sham_peaks_summits_UCSC_final.bed")
# 6. 運(yùn)行GuitarPlot
p_sham <- GuitarPlot(
? txTxdb = TxDb.Mmusculus.UCSC.mm10.knownGene,
? stBedFiles = "F:/XZHMU/huangyue/motif/RIP_sham_peaks_summits_UCSC_final.bed",? # 傳入文件路徑
? headOrtail = FALSE,
? enableCI = FALSE,
? mapFilterTranscript = TRUE,
? pltTxType = "mrna",
? stGroupName = "sham"
)
# 7. 顯示結(jié)果
print(p_sham)
合并兩個(gè)圖
# 讀取兩組數(shù)據(jù)
bed_path_sham <- "F:/XZHMU/huangyue/motif/RIP_sham_peaks_summits_UCSC_final.bed"
bed_path_sni <- "F:/XZHMU/huangyue/motif/RIP_SNI_peaks_summits_UCSC_final.bed"
bed_sham <- import(bed_path_sham, format = "bed")
bed_sni <- import(bed_path_sni, format = "bed")
# 添加分組信息
mcols(bed_sham)$group <- "sham"
mcols(bed_sni)$group <- "SNI"
# 合并兩組數(shù)據(jù)
combined_bed <- c(bed_sham, bed_sni)
# 讀取轉(zhuǎn)錄組數(shù)據(jù)庫
txdb_file <- system.file("extdata", "mm10_toy.sqlite", package = "Guitar")
txdb <- loadDb(txdb_file)
# 運(yùn)行 GuitarPlot
p_combined <- GuitarPlot(
? txTxdb = txdb,
? stBedFiles = list(bed_sham, bed_sni),? # 使用列表傳遞兩組數(shù)據(jù)
? stGroupName = c("sham", "SNI"),? # 分組名稱
? colors = c("blue", "red"),? # 為不同分組指定顏色
? pltTxType = "mrna",
? headOrtail = FALSE,
? enableCI = FALSE,
? mapFilterTranscript = TRUE
)
# 顯示結(jié)果
print(p_combined)
