PlotScripts/PathView.r at main · wchenqi/PlotScripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#!/usr/bin/env Rscript

#### 脚本说明:
# 1) mamba activate r452
# 2) 基于基因列表画pathview图

## 下载全部通路KGML文件
# 安装并加载KEGGREST包
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("KEGGREST")
## mamba install bioconda::bioconductor-keggrest

library(KEGGREST)
library(pathview)
library(SBGNview)

# 定义物种和对应的文件夹
species_list <- list(
  hsa = list(name = "Human", dir = "./kegg_local_data_hsa"),
  mmu = list(name = "Mouse", dir = "./kegg_local_data_mmu")
)

# 循环处理每个物种
for (sp_code in names(species_list)) {
  cat("\n===== 开始处理物种:", sp_code, "=====\n")

  # 创建该物种的专用文件夹
  sp_dir <- species_list[[sp_code]]$dir
  if (!dir.exists(sp_dir)) {
    dir.create(sp_dir, recursive = TRUE)
  }

  # 切换到该文件夹（或直接使用绝对路径，但keggGet会自动下载到当前工作目录，所以需要setwd）
  original_wd <- getwd()
  setwd(sp_dir)

  # 获取该物种所有通路ID
  pathways <- keggList("pathway", sp_code)
  pathway_ids <- sub("path:", "", names(pathways))

  # 下载每个通路的KGML文件
  for (pid in pathway_ids) {
    message("下载 ", sp_code, " 通路: ", pid)
    kgml <- keggGet(pid, option = "kgml")
    writeLines(kgml, con = paste0(pid, ".kgml"))
    Sys.sleep(1)  # 避免请求过快
  }

  setwd(original_wd)
  cat("物种", sp_code, "下载完成，文件保存在:", sp_dir, "\n")
}

## 使用本地文件绘图
# 人类通路绘图
pathview(gene.data = human_data, pathway.id = "hsa04910",
         species = "hsa", kegg.dir = "./kegg_local_data_hsa/")

# 小鼠通路绘图
pathview(gene.data = mouse_data, pathway.id = "mmu04910",
         species = "mmu", kegg.dir = "./kegg_local_data_mmu/")


### 或者使用SBGNview作为替代
## 支持KEGG、Reactome、PANTHER、MetaCyc、SMPDB、MetaCrop等 7大主流数据库及用户自定义通路
# 安装: BiocManager::install(c("SBGNview", "SBGNview.data"))

library(SBGNview)
library("org.Mm.eg.db")
library(clusterProfiler)

## 设置参数
changeID <- "no"        # 默认no, 或者设置转换组 SYMBOL_ENTREZID
gg <- ""                #
species <- "mouse"      # 默认"mouse", 或者"human", 对应替换基因名要用的数据库


if (changeID != "no"){
    # 替换基因名
    ENTREZID <- bitr(geneID = gg,
                    fromType = "SYMBOL",
                    toType  = "ENTREZID",
                    OrgDb   = org.Mm.eg.db,
                    drop = FALSE)

    if(species == "mouse"){
        translist <- mapIds(org.Mm.eg.db, keys = features, keytype = ENSEMBL, column= SYMBOL)
    }else if(species == "human"){
        translist <- mapIds(org.Hs.eg.db, keys = features, keytype = ENSEMBL, column= SYMBOL)
    }

    # 模拟你的数据（gene是Entrez ID，value是logFC）
    my_gene_data <- c(1236, 1000, 1234, 10000)   # Entrez IDs
    names(my_gene_data) <- c(1236, 1000, 1234, 10000)
    my_gene_data <- log2(my_gene_data)

    # 使用SBGNview直接绘图
    SBGNview(gene.data = my_gene_data,
            pathway.id = "hsa04910",  # 这里是示例ID：人源的胰岛素通路
            species = "hsa",
            out.suffix = "my_analysis_hsa")
}