From ac55c1f01b9138bb2b7649da35b8aa1c9e9903be Mon Sep 17 00:00:00 2001 From: Presisitence Date: Wed, 19 Aug 2026 22:06:38 +0800 Subject: [PATCH 1/2] Add plugin rgraph RNA-seq downstream analysis and publication figures via local Rscript (DESeq2/edgeR/limma, volcano, heatmap, GSEA, WGCNA). --- plugins/Presisitence/rgraph/LICENSE | 21 + plugins/Presisitence/rgraph/README.md | 37 ++ plugins/Presisitence/rgraph/mcp.json | 11 + plugins/Presisitence/rgraph/plugin.json | 22 + plugins/Presisitence/rgraph/pyproject.toml | 23 + .../rgraph/rgraph_toolkit/__init__.py | 6 + .../Presisitence/rgraph/rgraph_toolkit/cli.py | 34 ++ .../rgraph/rgraph_toolkit/runner.py | 291 ++++++++++++ .../Presisitence/rgraph/rscripts/_common.R | 167 +++++++ .../Presisitence/rgraph/rscripts/build_gmt.R | 57 +++ .../rgraph/rscripts/correlation_heatmap.R | 34 ++ plugins/Presisitence/rgraph/rscripts/diff.R | 102 +++++ .../rgraph/rscripts/distribution.R | 53 +++ plugins/Presisitence/rgraph/rscripts/enrich.R | 82 ++++ .../rgraph/rscripts/enrich_circos.R | 94 ++++ .../Presisitence/rgraph/rscripts/gene_bar.R | 47 ++ .../rgraph/rscripts/gene_correlation.R | 44 ++ .../Presisitence/rgraph/rscripts/go_plot.R | 63 +++ plugins/Presisitence/rgraph/rscripts/gsea.R | 49 ++ .../rgraph/rscripts/gsea_mountain.R | 46 ++ .../Presisitence/rgraph/rscripts/heatmap.R | 73 +++ .../rgraph/rscripts/interaction_diff.R | 55 +++ .../rgraph/rscripts/interactive_volcano.R | 42 ++ .../Presisitence/rgraph/rscripts/kegg_plot.R | 44 ++ plugins/Presisitence/rgraph/rscripts/kmeans.R | 86 ++++ .../rgraph/rscripts/multilevel_scatter.R | 43 ++ .../Presisitence/rgraph/rscripts/network.R | 87 ++++ .../rgraph/rscripts/norm_matrix.R | 36 ++ .../Presisitence/rgraph/rscripts/normalize.R | 37 ++ .../rgraph/rscripts/pathway_gene_heatmap.R | 70 +++ .../rgraph/rscripts/pathway_trend.R | 56 +++ plugins/Presisitence/rgraph/rscripts/pca.R | 56 +++ plugins/Presisitence/rgraph/rscripts/ppi.R | 50 +++ .../Presisitence/rgraph/rscripts/quadrant.R | 66 +++ plugins/Presisitence/rgraph/rscripts/sankey.R | 43 ++ .../rgraph/rscripts/segmented_volcano.R | 46 ++ plugins/Presisitence/rgraph/rscripts/ssgsea.R | 40 ++ .../rgraph/rscripts/stacked_volcano.R | 62 +++ plugins/Presisitence/rgraph/rscripts/venn.R | 49 ++ .../rgraph/rscripts/virtual_knockout.R | 69 +++ .../Presisitence/rgraph/rscripts/volcano.R | 61 +++ plugins/Presisitence/rgraph/rscripts/wgcna.R | 125 ++++++ plugins/Presisitence/rgraph/server.py | 421 ++++++++++++++++++ .../rgraph/skills/rgraph/SKILL.md | 19 + .../rgraph/tests/data/gene_count.csv | 11 + .../rgraph/tests/data/sample_group.csv | 7 + 46 files changed, 3037 insertions(+) create mode 100644 plugins/Presisitence/rgraph/LICENSE create mode 100644 plugins/Presisitence/rgraph/README.md create mode 100644 plugins/Presisitence/rgraph/mcp.json create mode 100644 plugins/Presisitence/rgraph/plugin.json create mode 100644 plugins/Presisitence/rgraph/pyproject.toml create mode 100644 plugins/Presisitence/rgraph/rgraph_toolkit/__init__.py create mode 100644 plugins/Presisitence/rgraph/rgraph_toolkit/cli.py create mode 100644 plugins/Presisitence/rgraph/rgraph_toolkit/runner.py create mode 100644 plugins/Presisitence/rgraph/rscripts/_common.R create mode 100644 plugins/Presisitence/rgraph/rscripts/build_gmt.R create mode 100644 plugins/Presisitence/rgraph/rscripts/correlation_heatmap.R create mode 100644 plugins/Presisitence/rgraph/rscripts/diff.R create mode 100644 plugins/Presisitence/rgraph/rscripts/distribution.R create mode 100644 plugins/Presisitence/rgraph/rscripts/enrich.R create mode 100644 plugins/Presisitence/rgraph/rscripts/enrich_circos.R create mode 100644 plugins/Presisitence/rgraph/rscripts/gene_bar.R create mode 100644 plugins/Presisitence/rgraph/rscripts/gene_correlation.R create mode 100644 plugins/Presisitence/rgraph/rscripts/go_plot.R create mode 100644 plugins/Presisitence/rgraph/rscripts/gsea.R create mode 100644 plugins/Presisitence/rgraph/rscripts/gsea_mountain.R create mode 100644 plugins/Presisitence/rgraph/rscripts/heatmap.R create mode 100644 plugins/Presisitence/rgraph/rscripts/interaction_diff.R create mode 100644 plugins/Presisitence/rgraph/rscripts/interactive_volcano.R create mode 100644 plugins/Presisitence/rgraph/rscripts/kegg_plot.R create mode 100644 plugins/Presisitence/rgraph/rscripts/kmeans.R create mode 100644 plugins/Presisitence/rgraph/rscripts/multilevel_scatter.R create mode 100644 plugins/Presisitence/rgraph/rscripts/network.R create mode 100644 plugins/Presisitence/rgraph/rscripts/norm_matrix.R create mode 100644 plugins/Presisitence/rgraph/rscripts/normalize.R create mode 100644 plugins/Presisitence/rgraph/rscripts/pathway_gene_heatmap.R create mode 100644 plugins/Presisitence/rgraph/rscripts/pathway_trend.R create mode 100644 plugins/Presisitence/rgraph/rscripts/pca.R create mode 100644 plugins/Presisitence/rgraph/rscripts/ppi.R create mode 100644 plugins/Presisitence/rgraph/rscripts/quadrant.R create mode 100644 plugins/Presisitence/rgraph/rscripts/sankey.R create mode 100644 plugins/Presisitence/rgraph/rscripts/segmented_volcano.R create mode 100644 plugins/Presisitence/rgraph/rscripts/ssgsea.R create mode 100644 plugins/Presisitence/rgraph/rscripts/stacked_volcano.R create mode 100644 plugins/Presisitence/rgraph/rscripts/venn.R create mode 100644 plugins/Presisitence/rgraph/rscripts/virtual_knockout.R create mode 100644 plugins/Presisitence/rgraph/rscripts/volcano.R create mode 100644 plugins/Presisitence/rgraph/rscripts/wgcna.R create mode 100644 plugins/Presisitence/rgraph/server.py create mode 100644 plugins/Presisitence/rgraph/skills/rgraph/SKILL.md create mode 100644 plugins/Presisitence/rgraph/tests/data/gene_count.csv create mode 100644 plugins/Presisitence/rgraph/tests/data/sample_group.csv diff --git a/plugins/Presisitence/rgraph/LICENSE b/plugins/Presisitence/rgraph/LICENSE new file mode 100644 index 0000000..9637979 --- /dev/null +++ b/plugins/Presisitence/rgraph/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 Presisitence + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/plugins/Presisitence/rgraph/README.md b/plugins/Presisitence/rgraph/README.md new file mode 100644 index 0000000..d89561d --- /dev/null +++ b/plugins/Presisitence/rgraph/README.md @@ -0,0 +1,37 @@ +# rgraph + +RNA-seq **downstream** analysis and figures for MiniMax Code. Parameterized R scripts +(ggplot2, pheatmap, clusterProfiler, edgeR, limma, WGCNA, …) are exposed as MCP tools and +rendered by the user's local **Rscript** to png + pdf. + +The Plugin consumes the user's own count/FPKM tables. It does not ship genomes or experimental matrices. +`tests/data/` is a tiny synthetic `g1`–`g10` table for smoke tests. + +## Try it + +```text +I have gene_count.csv and sample_group.csv. Run DESeq2 for treatment vs control, then draw a +volcano plot (label top 8 genes, padj < 0.05, |log2FC| > 1) and a clustered heatmap of DEGs. +``` + +Expected result: the agent calls `rgraph_env`, then `rgraph_diff(method="deseq2")` and +`rgraph_volcano` / `rgraph_heatmap`. Output png+pdf paths are returned. DEG calls use **padj** +by default. If R or a package is missing, the tool returns install commands instead of crashing. + +## Requirements + +- Python 3.10+ and [uv](https://docs.astral.sh/uv/) on PATH. +- R with `Rscript` on PATH, or set `RGRAPH_RSCRIPT` to the Rscript executable. +- Common R packages: ggplot2, pheatmap, edgeR or DESeq2 or limma, clusterProfiler as needed. + Missing packages are reported with CRAN/Bioconductor install lines. +- Windows, macOS, and Linux. + +## Data and network + +- Default analyses are local: user CSVs in, png/pdf out. No telemetry. +- `rgraph_ppi` may contact STRING (`string-db.org`) when the user asks for a PPI edge table. +- No credentials in the package. + +## License + +MIT. See [LICENSE](LICENSE). diff --git a/plugins/Presisitence/rgraph/mcp.json b/plugins/Presisitence/rgraph/mcp.json new file mode 100644 index 0000000..1cdbe22 --- /dev/null +++ b/plugins/Presisitence/rgraph/mcp.json @@ -0,0 +1,11 @@ +{ + "$schema": "https://agent-plugins.org/schemas/1.0.0/mcp.schema.json", + "mcpServers": { + "rgraph": { + "type": "stdio", + "command": "uv", + "args": ["run", "server.py"], + "cwd": "${PLUGIN_ROOT}" + } + } +} diff --git a/plugins/Presisitence/rgraph/plugin.json b/plugins/Presisitence/rgraph/plugin.json new file mode 100644 index 0000000..2db2c7d --- /dev/null +++ b/plugins/Presisitence/rgraph/plugin.json @@ -0,0 +1,22 @@ +{ + "$schema": "https://agent-plugins.org/schemas/1.0.0/plugin.schema.json", + "name": "rgraph", + "version": "0.1.0", + "description": "RNA-seq downstream analysis and publication figures via local Rscript (DESeq2/edgeR/limma, volcano, heatmap, GSEA, WGCNA).", + "author": { + "name": "Presisitence", + "url": "https://github.com/Presisitence" + }, + "homepage": "https://github.com/Presisitence/rgraph-mcp", + "repository": "https://github.com/Presisitence/rgraph-mcp", + "license": "MIT", + "keywords": [ + "minimax-code", + "plugin", + "mcp", + "rnaseq", + "deseq2", + "ggplot2", + "wgcna" + ] +} diff --git a/plugins/Presisitence/rgraph/pyproject.toml b/plugins/Presisitence/rgraph/pyproject.toml new file mode 100644 index 0000000..3808ce4 --- /dev/null +++ b/plugins/Presisitence/rgraph/pyproject.toml @@ -0,0 +1,23 @@ +[project] +name = "rgraph-mcp" +version = "0.1.0" +description = "RNAseq 下游可视化/分析 MCP 服务:把一套 R (ggplot2/pheatmap/clusterProfiler/WGCNA...) 出图脚本参数化封装为 MCP 工具,由 Rscript 引擎渲染" +readme = "README.md" +requires-python = ">=3.10" +license = "MIT" +dependencies = [ + "mcp[cli]>=1.2.0", +] + +[project.scripts] +rgraph-cli = "rgraph_toolkit.cli:main" + +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[tool.hatch.build.targets.wheel] +packages = ["rgraph_toolkit"] + +[tool.hatch.build.targets.wheel.force-include] +"rscripts" = "rgraph_toolkit/rscripts" diff --git a/plugins/Presisitence/rgraph/rgraph_toolkit/__init__.py b/plugins/Presisitence/rgraph/rgraph_toolkit/__init__.py new file mode 100644 index 0000000..3677a2f --- /dev/null +++ b/plugins/Presisitence/rgraph/rgraph_toolkit/__init__.py @@ -0,0 +1,6 @@ +"""rgraph_toolkit: 把一套 RNAseq 下游 R 出图脚本参数化封装为 MCP 工具。""" +from __future__ import annotations + +from . import runner + +__all__ = ["runner"] diff --git a/plugins/Presisitence/rgraph/rgraph_toolkit/cli.py b/plugins/Presisitence/rgraph/rgraph_toolkit/cli.py new file mode 100644 index 0000000..c87694e --- /dev/null +++ b/plugins/Presisitence/rgraph/rgraph_toolkit/cli.py @@ -0,0 +1,34 @@ +"""命令行自检入口:`rgraph-cli` 检查 R 引擎与关键包是否就绪。""" +from __future__ import annotations + +import sys + +from . import runner + +_KEY_PKGS = [ + "ggplot2", "dplyr", "tidyr", "stringr", "pheatmap", "ggrepel", + "DESeq2", "edgeR", "limma", "clusterProfiler", "enrichplot", + "pathview", "circlize", "ComplexHeatmap", "WGCNA", "VennDiagram", + "GSVA", "fgsea", "ggpubr", "patchwork", "reshape2", +] + + +def main() -> int: + rscript = runner.find_rscript() + print(f"Rscript: {rscript or '未找到 (设置 RGRAPH_RSCRIPT 或安装 R)'}") + print(f"rscripts 目录: {runner.rscripts_dir()}") + if not rscript: + return 1 + status = runner.check_packages(_KEY_PKGS, rscript=rscript) + print("\nR 包状态:") + for pkg, ok in status.items(): + print(f" [{'x' if ok else ' '}] {pkg}") + missing = [p for p, ok in status.items() if not ok] + if missing: + print("\n缺失包安装建议:") + print(" " + runner._install_hint(missing)) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/plugins/Presisitence/rgraph/rgraph_toolkit/runner.py b/plugins/Presisitence/rgraph/rgraph_toolkit/runner.py new file mode 100644 index 0000000..911353d --- /dev/null +++ b/plugins/Presisitence/rgraph/rgraph_toolkit/runner.py @@ -0,0 +1,291 @@ +"""Rscript 执行引擎:定位 R、序列化参数、运行 rscripts/ 下的 R 模板并解析结果。 + +设计要点 +-------- +* **忠实 R 出图**:不在 Python 里重画,而是驱动一套清洗过的参数化 R 脚本(ggplot2/pheatmap/ + clusterProfiler/WGCNA 等),产出与课程一致风格的 png/pdf。 +* **零额外 R 依赖传参**:参数以「R 源文件」形式落盘(`params <- list(...)`),R 端 `source()` + 即可,不依赖 jsonlite 等包。 +* **优雅降级**:定位不到 Rscript → 返回可直接手动运行的脚本与命令;R 端缺包 → 返回缺失包名与 + 安装建议,而不是抛出晦涩的 R 报错。 +* **可复现**:以 `--no-init-file --no-site-file` 运行,屏蔽用户 .Rprofile/Rprofile.site + 的自动加载横幅,避免污染输出。 +""" +from __future__ import annotations + +import os +import shutil +import subprocess +import tempfile +from pathlib import Path +from typing import Any + +# --------------------------------------------------------------------------- +# Rscript 定位 +# --------------------------------------------------------------------------- + +_COMMON_DIRS = [ + r"C:\Program Files\R", + r"C:\Program Files (x86)\R", +] + + +def find_rscript() -> str | None: + """按优先级定位 Rscript:环境变量 RGRAPH_RSCRIPT → PATH → 常见安装目录/注册表。""" + env = os.environ.get("RGRAPH_RSCRIPT") + if env and Path(env).exists(): + return env + + on_path = shutil.which("Rscript") or shutil.which("Rscript.exe") + if on_path: + return on_path + + # 常见安装根目录下搜 R-*/bin/Rscript.exe + for root in _COMMON_DIRS: + p = Path(root) + if p.name.lower() == "bin" and (p / "Rscript.exe").exists(): + return str(p / "Rscript.exe") + if p.exists(): + for sub in sorted(p.glob("R-*"), reverse=True): + cand = sub / "bin" / "Rscript.exe" + if cand.exists(): + return str(cand) + + # Windows 注册表 + try: + import winreg # type: ignore + + for hive in (winreg.HKEY_LOCAL_MACHINE, winreg.HKEY_CURRENT_USER): + try: + with winreg.OpenKey(hive, r"SOFTWARE\R-core\R") as k: + install_path, _ = winreg.QueryValueEx(k, "InstallPath") + cand = Path(install_path) / "bin" / "Rscript.exe" + if cand.exists(): + return str(cand) + except OSError: + continue + except Exception: # noqa: BLE001 + pass + return None + + +# --------------------------------------------------------------------------- +# 路径 +# --------------------------------------------------------------------------- + +def rscripts_dir() -> Path: + """定位 rscripts 目录(源码布局为包的同级;wheel 布局为包内)。""" + here = Path(__file__).resolve().parent + for cand in (here.parent / "rscripts", here / "rscripts"): + if cand.is_dir(): + return cand + return here.parent / "rscripts" + + +# --------------------------------------------------------------------------- +# Python -> R 字面量序列化 +# --------------------------------------------------------------------------- + +def _r_str(s: str) -> str: + s = s.replace("\\", "\\\\").replace('"', '\\"').replace("\n", "\\n").replace("\t", "\\t") + return f'"{s}"' + + +def to_r_literal(value: Any) -> str: + """把 Python 值转成 R 字面量表达式。""" + if value is None: + return "NULL" + if isinstance(value, bool): + return "TRUE" if value else "FALSE" + if isinstance(value, (int, float)): + if isinstance(value, float) and (value != value): # NaN + return "NA" + return repr(value) + if isinstance(value, str): + return _r_str(value) + if isinstance(value, (list, tuple)): + if not value: + return "c()" + return "c(" + ", ".join(to_r_literal(v) for v in value) + ")" + if isinstance(value, dict): + parts = [] + for k, v in value.items(): + name = str(k) + key = name if name.isidentifier() else f"`{name}`" + parts.append(f"{key} = {to_r_literal(v)}") + return "list(" + ", ".join(parts) + ")" + # 兜底:字符串化 + return _r_str(str(value)) + + +def write_params_file(params: dict[str, Any], path: Path) -> None: + """把参数写成 `params <- list(...)` 的 R 源文件。""" + body = to_r_literal(params) + path.write_text(f"params <- {body}\n", encoding="utf-8") + + +# --------------------------------------------------------------------------- +# 执行 +# --------------------------------------------------------------------------- + +_START_FLAGS = ["--no-init-file", "--no-site-file", "--no-restore"] + + +def build_command(script_name: str, params_file: Path, rscript: str | None = None) -> list[str] | None: + rscript = rscript or find_rscript() + if not rscript: + return None + script = rscripts_dir() / script_name + common = rscripts_dir() / "_common.R" + # trailing args: [params_file, _common.R 路径],R 端 source(args[[2]]) 再读 args[[1]] + return [rscript, *_START_FLAGS, str(script), str(params_file), str(common)] + + +def run_script( + script_name: str, + params: dict[str, Any], + *, + outdir: str | None = None, + timeout: float = 600.0, +) -> dict[str, Any]: + """运行一个 R 模板脚本,返回结构化结果。 + + 返回字段: + status: ok | missing_packages | r_not_found | error + outputs: 产物文件绝对路径列表(R 端通过 `RGRAPH_OUTPUT:` 标记回传) + packages: 缺失的 R 包(status=missing_packages 时) + log: R 端 stdout+stderr(截断) + command / script / params_file: 便于手动复现 + """ + if outdir: + Path(outdir).mkdir(parents=True, exist_ok=True) + params = {**params, "outdir": outdir} + + rscript = find_rscript() + # 参数文件落在 outdir(便于复现),否则临时目录 + param_dir = Path(outdir) if outdir else Path(tempfile.mkdtemp(prefix="rgraph_")) + params_file = param_dir / f".rgraph_params_{script_name.replace('.R', '')}.R" + write_params_file(params, params_file) + + script_path = rscripts_dir() / script_name + if not script_path.exists(): + return {"status": "error", "error": f"R 模板不存在: {script_path}", + "script": str(script_path)} + + cmd = build_command(script_name, params_file, rscript) + if cmd is None: + common = rscripts_dir() / "_common.R" + manual = f'"" {" ".join(_START_FLAGS)} "{script_path}" "{params_file}" "{common}"' + return { + "status": "r_not_found", + "message": ("未找到 Rscript。请安装 R 或设置环境变量 RGRAPH_RSCRIPT 指向 Rscript(.exe)。" + "已生成可手动运行的脚本与参数文件。"), + "script": str(script_path), + "params_file": str(params_file), + "manual_command": manual, + } + + try: + proc = subprocess.run( + cmd, capture_output=True, text=True, timeout=timeout, + encoding="utf-8", errors="replace", stdin=subprocess.DEVNULL, + ) + except subprocess.TimeoutExpired: + return {"status": "error", "error": f"R 运行超时(>{timeout}s)", "command": cmd} + + stdout = proc.stdout or "" + stderr = proc.stderr or "" + combined = stdout + ("\n" + stderr if stderr else "") + + outputs: list[str] = [] + missing: list[str] = [] + for line in combined.splitlines(): + line = line.strip() + if line.startswith("RGRAPH_OUTPUT:"): + outputs.append(line[len("RGRAPH_OUTPUT:"):].strip()) + elif line.startswith("RGRAPH_MISSING_PACKAGES:"): + pk = line[len("RGRAPH_MISSING_PACKAGES:"):].strip() + missing = [x.strip() for x in pk.replace(";", ",").split(",") if x.strip()] + + log_tail = combined if len(combined) < 6000 else combined[-6000:] + + if missing: + return { + "status": "missing_packages", + "packages": missing, + "install_hint": _install_hint(missing), + "script": str(script_path), + "params_file": str(params_file), + "log": log_tail, + } + + if proc.returncode != 0: + return { + "status": "error", + "returncode": proc.returncode, + "error": "R 脚本执行失败,见 log。", + "script": str(script_path), + "params_file": str(params_file), + "command": cmd, + "log": log_tail, + } + + return { + "status": "ok", + "outputs": outputs, + "n_outputs": len(outputs), + "script": str(script_path), + "params_file": str(params_file), + "log": log_tail, + } + + +# --------------------------------------------------------------------------- +# 环境自检 / 安装建议 +# --------------------------------------------------------------------------- + +# CRAN vs Bioconductor 归属,用于给出正确的安装命令 +_BIOC = { + "DESeq2", "edgeR", "limma", "clusterProfiler", "enrichplot", "pathview", + "ComplexHeatmap", "GSVA", "GSEABase", "fgsea", "WGCNA", "AnnotationDbi", + "GO.db", "org.Hs.eg.db", "org.Mm.eg.db", "impute", "preprocessCore", "GENIE3", +} + + +def _install_hint(pkgs: list[str]) -> str: + cran = [p for p in pkgs if p not in _BIOC] + bioc = [p for p in pkgs if p in _BIOC] + lines = [] + if cran: + quoted = ", ".join(f'"{p}"' for p in cran) + lines.append(f'install.packages(c({quoted}))') + if bioc: + quoted = ", ".join(f'"{p}"' for p in bioc) + lines.append('if (!require("BiocManager")) install.packages("BiocManager")') + lines.append(f'BiocManager::install(c({quoted}))') + return " ; ".join(lines) + + +def check_packages(pkgs: list[str], rscript: str | None = None) -> dict[str, bool]: + """查询若干 R 包是否已安装。""" + rscript = rscript or find_rscript() + if not rscript: + return {p: False for p in pkgs} + vec = ", ".join(f'"{p}"' for p in pkgs) + code = f'i<-rownames(installed.packages()); for(x in c({vec})) cat(x,"=",x %in% i,"\\n")' + try: + proc = subprocess.run( + [rscript, *_START_FLAGS, "-e", code], + capture_output=True, text=True, timeout=120, + encoding="utf-8", errors="replace", stdin=subprocess.DEVNULL, + ) + except Exception: # noqa: BLE001 + return {p: False for p in pkgs} + result: dict[str, bool] = {} + for line in (proc.stdout or "").splitlines(): + if "=" in line: + name, _, val = line.partition("=") + result[name.strip()] = val.strip().upper().startswith("TRUE") + for p in pkgs: + result.setdefault(p, False) + return result diff --git a/plugins/Presisitence/rgraph/rscripts/_common.R b/plugins/Presisitence/rgraph/rscripts/_common.R new file mode 100644 index 0000000..afd6b49 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/_common.R @@ -0,0 +1,167 @@ +# ============================================================================= +# _common.R — rgraph-mcp 共享 R 库 +# 所有出图/分析模板都 source() 本文件,统一:参数读取、包预检、样本表规整、 +# 出版级主题、配色、png+pdf 导出、产物回传标记。 +# +# 相比原课程脚本的改进: +# * 不再 setwd() 到硬编码路径 / 不 rm(list=ls());一切走参数与相对/绝对入参。 +# * 统一 group / group_name 列名混用问题(二者缺一自动补齐)。 +# * 缺包时输出 RGRAPH_MISSING_PACKAGES 标记并以状态码 3 退出,交由上层友好提示。 +# * 产物统一通过 RGRAPH_OUTPUT 标记回传,Python 侧据此收集图片路径。 +# ============================================================================= + +options(stringsAsFactors = FALSE, warn = 1) + +# ---- 参数读取 --------------------------------------------------------------- +rgraph_load_params <- function() { + args <- commandArgs(trailingOnly = TRUE) + if (length(args) < 1) stop("缺少参数文件路径 (params file)") + e <- new.env() + sys.source(args[[1]], envir = e) + if (!exists("params", envir = e)) stop("参数文件未定义 params 列表") + get("params", envir = e) +} + +rgraph_opt <- function(p, key, default = NULL) { + v <- p[[key]] + if (is.null(v)) default else v +} + +`%||%` <- function(a, b) if (is.null(a)) b else a + +# ---- 包预检 ----------------------------------------------------------------- +rgraph_need <- function(pkgs) { + inst <- rownames(installed.packages()) + miss <- setdiff(pkgs, inst) + if (length(miss) > 0) { + cat("RGRAPH_MISSING_PACKAGES:", paste(miss, collapse = ","), "\n") + quit(save = "no", status = 3) + } + invisible(TRUE) +} + +rgraph_library <- function(pkgs) { + rgraph_need(pkgs) + for (p in pkgs) suppressMessages(suppressWarnings(library(p, character.only = TRUE))) + invisible(TRUE) +} + +# ---- 产物回传 / 目录 -------------------------------------------------------- +rgraph_emit <- function(path) { + cat("RGRAPH_OUTPUT:", normalizePath(path, winslash = "/", mustWork = FALSE), "\n") +} + +rgraph_ensure_dir <- function(path) { + if (!dir.exists(path)) dir.create(path, recursive = TRUE, showWarnings = FALSE) + invisible(path) +} + +# ---- 读入:样本分组表(规整 group / group_name / TvsC) -------------------- +rgraph_read_sample_group <- function(path) { + sg <- read.csv(path, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + if (!"sample_name" %in% names(sg)) stop("sample_group 缺少 sample_name 列") + sg <- sg[!is.na(sg$sample_name) & sg$sample_name != "", , drop = FALSE] + has_g <- "group" %in% names(sg) + has_gn <- "group_name" %in% names(sg) + if (has_g && !has_gn) sg$group_name <- sg$group + if (has_gn && !has_g) sg$group <- sg$group_name + if (!has_g && !has_gn) { sg$group <- sg$sample_name; sg$group_name <- sg$sample_name } + if (!"TvsC" %in% names(sg)) sg$TvsC <- NA_character_ + sg +} + +# ---- 读入:表达矩阵 (gene_id + 各样本列) ----------------------------------- +rgraph_read_matrix <- function(path) { + df <- read.csv(path, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + if (!"gene_id" %in% names(df)) { + names(df)[1] <- "gene_id" # 容错:首列视作 gene_id + } + df +} + +# 校验 sample_group 中样本是否都在矩阵列里 +rgraph_check_samples <- function(mat, sample_list) { + miss <- setdiff(sample_list, colnames(mat)) + if (length(miss) > 0) { + stop(paste0("以下样本不在表达矩阵中: ", paste(miss, collapse = ", "))) + } + invisible(TRUE) +} + +# ---- 出版级主题 ------------------------------------------------------------- +theme_rgraph <- function(base_size = 12) { + ggplot2::theme_minimal(base_size = base_size) + + ggplot2::theme( + panel.grid = ggplot2::element_blank(), + panel.border = ggplot2::element_rect(colour = "black", fill = NA), + axis.text = ggplot2::element_text(color = "black"), + axis.title = ggplot2::element_text(color = "black"), + axis.ticks = ggplot2::element_line(color = "black", linewidth = 0.4), + axis.ticks.length = ggplot2::unit(0.1, "cm"), + legend.position = "right", + plot.title = ggplot2::element_text(hjust = 0.5, colour = "black") + ) +} + +# ---- 配色 ------------------------------------------------------------------- +.RGRAPH_PAL_COURSE <- c("#f03e3e", "#FF9933", "#CCCCFF", "#00CCCC", "#00CC66", + "#009999", "#d4d445", "#0066CC", "#FFCCFF") +# Okabe-Ito 色盲友好 +.RGRAPH_PAL_CB <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", + "#D55E00", "#CC79A7", "#999999", "#000000") + +rgraph_palette <- function(n, name = "course") { + base <- switch(name, + colorblind = .RGRAPH_PAL_CB, + course = .RGRAPH_PAL_COURSE, + .RGRAPH_PAL_COURSE) + if (n <= length(base)) base[seq_len(n)] else grDevices::colorRampPalette(base)(n) +} + +# 火山图/上下调三色(默认色盲友好红-蓝,可切换课程红-绿) +rgraph_updown_colors <- function(scheme = "rb") { + if (identical(scheme, "rg")) { + c(Up = "red", Down = "green", No = "grey80") + } else { + c(Up = "#d73027", Down = "#4575b4", No = "grey80") + } +} + +# ---- 保存:ggplot 对象 → png+pdf 并回传 ------------------------------------ +rgraph_save <- function(plot, outdir, base, width = 8, height = 6, dpi = 300, + formats = c("png", "pdf")) { + rgraph_ensure_dir(outdir) + for (fmt in formats) { + f <- file.path(outdir, paste0(base, ".", fmt)) + ggplot2::ggsave(f, plot = plot, width = width, height = height, dpi = dpi, bg = "white") + rgraph_emit(f) + } + invisible(TRUE) +} + +# ---- 保存:base graphics (pheatmap/plotGOgraph/pathview) → png+pdf -------- +rgraph_save_base <- function(outdir, base, draw_fn, width = 8, height = 6, dpi = 300, + formats = c("png", "pdf")) { + rgraph_ensure_dir(outdir) + for (fmt in formats) { + f <- file.path(outdir, paste0(base, ".", fmt)) + if (fmt == "png") { + grDevices::png(f, width = width, height = height, units = "in", res = dpi) + } else { + grDevices::pdf(f, width = width, height = height) + } + draw_fn() + grDevices::dev.off() + rgraph_emit(f) + } + invisible(TRUE) +} + +# ---- 写出 CSV 并回传 -------------------------------------------------------- +rgraph_write_csv <- function(df, outdir, base, row.names = FALSE) { + rgraph_ensure_dir(outdir) + f <- file.path(outdir, paste0(base, ".csv")) + write.csv(df, f, row.names = row.names) + rgraph_emit(f) + invisible(f) +} diff --git a/plugins/Presisitence/rgraph/rscripts/build_gmt.R b/plugins/Presisitence/rgraph/rscripts/build_gmt.R new file mode 100644 index 0000000..cec383b --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/build_gmt.R @@ -0,0 +1,57 @@ +# build_gmt.R —— 为 GSEA 构建 gmt 基因集 +# params: outdir, type("go"|"kegg"), orgdb(如 org.Hs.eg.db / 自建 org.Xxx.eg.db), +# source("orgdb"|"online"), keytype("GID"|"ENSEMBL"), kegg_species("hsa"), +# min_count(5), name(默认 GSEA_GO / GSEA_KEGG) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "stringr", "AnnotationDbi", p$orgdb)) +odb <- get(p$orgdb) + +type <- rgraph_opt(p, "type", "go") +min_count <- rgraph_opt(p, "min_count", 5) +keytype <- rgraph_opt(p, "keytype", "GID") + +to_gmt <- function(df, term_col, gene_col) { + df <- na.omit(df[, c(term_col, gene_col)]) + names(df) <- c("term", "gene") + df$term <- gsub(":", "_", df$term) + df <- df %>% add_count(term, name = "n") %>% filter(n > min_count) %>% dplyr::select(term, gene) + wide <- df %>% group_by(term) %>% mutate(row = row_number()) %>% + tidyr::pivot_wider(names_from = row, values_from = gene) %>% ungroup() + wide <- wide %>% mutate(desc = "-", .after = 1) + wide +} + +if (type == "go") { + g <- AnnotationDbi::select(odb, keys = AnnotationDbi::keys(odb, keytype = keytype), + keytype = keytype, columns = c("GO")) + gcol <- intersect(c(keytype, "GID", "ENSEMBL"), names(g))[1] + wide <- to_gmt(g, "GO", gcol) + nm <- rgraph_opt(p, "name", "GSEA_GO") +} else { + if (rgraph_opt(p, "source", "orgdb") == "online") { + sp <- rgraph_opt(p, "kegg_species", "hsa") + options(timeout = 600) + pe <- read.delim(sprintf("https://rest.kegg.jp/link/pathway/%s", sp), sep = "\t", + col.names = c("ENTREZID", "KEGGID")) + pe$ENTREZID <- str_remove(pe$ENTREZID, paste0(sp, ":")) + pe$KEGGID <- str_remove(pe$KEGGID, "path:") + map <- AnnotationDbi::select(odb, keys = AnnotationDbi::keys(odb, keytype = "ENTREZID"), + keytype = "ENTREZID", columns = c("ENSEMBL")) + pe <- merge(pe, na.omit(map), by = "ENTREZID") + wide <- to_gmt(pe, "KEGGID", "ENSEMBL") + } else { + g <- AnnotationDbi::select(odb, keys = AnnotationDbi::keys(odb, keytype = keytype), + keytype = keytype, columns = c("Pathway")) + gcol <- intersect(c(keytype, "GID"), names(g))[1] + wide <- to_gmt(g, "Pathway", gcol) + } + nm <- rgraph_opt(p, "name", "GSEA_KEGG") +} + +rgraph_ensure_dir(p$outdir) +f <- file.path(p$outdir, paste0(nm, ".gmt")) +write.table(wide, f, col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE, na = "") +rgraph_emit(f) +cat("RGRAPH_DONE: build_gmt", type, "genesets=", nrow(wide), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/correlation_heatmap.R b/plugins/Presisitence/rgraph/rscripts/correlation_heatmap.R new file mode 100644 index 0000000..4c11960 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/correlation_heatmap.R @@ -0,0 +1,34 @@ +# correlation_heatmap.R —— 样本间 Pearson 相关性热图 +# params: fpkm(路径), sample_group(路径), outdir, dpi(默认300), width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "ggplot2", "scales", "reshape2")) + +fpkm <- rgraph_read_matrix(p$fpkm) +sg <- rgraph_read_sample_group(p$sample_group) +rgraph_check_samples(fpkm, sg$sample_name) + +mat <- fpkm[, sg$sample_name, drop = FALSE] +mat[is.na(mat)] <- 0 +logm <- log2(mat + 1) + +cor_mat <- cor(logm, use = "pairwise.complete.obs") +cor_df <- reshape2::melt(cor_mat) + +low <- rgraph_opt(p, "low_color", "white") +high <- rgraph_opt(p, "high_color", "blue") + +pl <- ggplot2::ggplot(cor_df, ggplot2::aes(x = Var2, y = Var1, fill = value)) + + ggplot2::geom_tile(color = "white", linewidth = 0.5) + + ggplot2::scale_fill_gradient(low = low, high = high, name = expression(R^2)) + + ggplot2::geom_text(ggplot2::aes(label = round(value, 3)), color = "black", size = 3) + + ggplot2::labs(x = NULL, y = NULL, title = "Pearson correlation between samples") + + theme_rgraph() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + + ggplot2::coord_fixed() + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "correlation"), + width = rgraph_opt(p, "width", 6), height = rgraph_opt(p, "height", 6), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: correlation_heatmap samples=", ncol(mat), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/diff.R b/plugins/Presisitence/rgraph/rscripts/diff.R new file mode 100644 index 0000000..4432f38 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/diff.R @@ -0,0 +1,102 @@ +# diff.R —— 差异表达分析(DESeq2 / edgeR / edgeR_norep / limma) +# params: count, sample_group, outdir, method("deseq2"|"edger"|"edger_norep"|"limma"), +# sig_metric("padj"|"pvalue",默认padj), pcut(0.05), log2fc(1), +# min_mean_count(1), dispersion(edger_norep用,默认0.1) +# 输出: 01.Dse2_result.csv, 02.Deg_all.csv, 03.Deg_down.csv, 04.Deg_up.csv +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() + +method <- rgraph_opt(p, "method", "deseq2") +sig_metric <- rgraph_opt(p, "sig_metric", "padj") +pcut <- rgraph_opt(p, "pcut", 0.05) +log2fc <- rgraph_opt(p, "log2fc", 1) +min_mean <- rgraph_opt(p, "min_mean_count", 1) + +rgraph_library(c("dplyr")) +count <- rgraph_read_matrix(p$count) +sg <- rgraph_read_sample_group(p$sample_group) +rgraph_check_samples(count, sg$sample_name) + +# 实验/对照组:优先 TvsC,其次取前两组 +if (any(!is.na(sg$TvsC)) && all(c("treatment", "control") %in% sg$TvsC)) { + treatment_g <- unique(sg$group[sg$TvsC == "treatment"])[1] + control_g <- unique(sg$group[sg$TvsC == "control"])[1] +} else { + gs <- unique(sg$group) + if (length(gs) < 2) stop("分组不足两组,无法差异分析") + treatment_g <- gs[1]; control_g <- gs[2] +} + +mat <- count[, c("gene_id", sg$sample_name)] +rownames(mat) <- mat$gene_id; mat$gene_id <- NULL +mat[is.na(mat)] <- 0 +grp <- factor(sg$group[match(colnames(mat), sg$sample_name)]) + +# ---- 各方法产出 stats(gene_id,log2FoldChange,pvalue,padj) 与 norm(矩阵) ----- +if (method == "deseq2") { + rgraph_library(c("DESeq2")) + coldata <- data.frame(row.names = colnames(mat), group = grp) + dds <- DESeq2::DESeqDataSetFromMatrix(round(as.matrix(mat)), coldata, design = ~group) + dds <- DESeq2::DESeq(dds) + norm <- as.data.frame(DESeq2::counts(dds, normalized = TRUE)) + res <- as.data.frame(DESeq2::results(dds, contrast = c("group", treatment_g, control_g))) + stats <- data.frame(gene_id = rownames(res), log2FoldChange = res$log2FoldChange, + pvalue = res$pvalue, padj = res$padj) +} else if (method %in% c("edger", "edger_norep")) { + rgraph_library(c("edgeR")) + dge <- edgeR::DGEList(counts = as.matrix(mat), group = grp) + dge <- edgeR::calcNormFactors(dge, method = "TMM") + if (method == "edger_norep") { + disp <- rgraph_opt(p, "dispersion", 0.1) + dge <- edgeR::estimateDisp(dge, trend = "none", tagwise = FALSE) + dge$common.dispersion <- disp + dge$trended.dispersion <- rep(disp, nrow(dge)) + dge$tagwise.dispersion <- rep(disp, nrow(dge)) + } else { + dge <- edgeR::estimateDisp(dge) + } + et <- edgeR::exactTest(dge, pair = c(control_g, treatment_g)) # logFC = treatment/control + norm <- as.data.frame(edgeR::cpm(dge, normalized.lib.sizes = TRUE)) + tt <- edgeR::topTags(et, n = Inf, adjust.method = "BH", sort.by = "none")$table + stats <- data.frame(gene_id = rownames(tt), log2FoldChange = tt$logFC, + pvalue = tt$PValue, padj = tt$FDR) +} else if (method == "limma") { + rgraph_library(c("limma")) + design <- model.matrix(~0 + grp) + colnames(design) <- make.names(levels(grp)) + v <- limma::voom(as.matrix(mat), design) + norm <- as.data.frame(2^v$E) + fit <- limma::lmFit(v, design) + ct <- limma::makeContrasts(contrasts = paste0(make.names(treatment_g), "-", make.names(control_g)), + levels = design) # 修正为 treatment - control(原课程写反) + fit2 <- limma::eBayes(limma::contrasts.fit(fit, ct)) + tt <- limma::topTable(fit2, number = Inf) + stats <- data.frame(gene_id = rownames(tt), log2FoldChange = tt$logFC, + pvalue = tt$P.Value, padj = tt$adj.P.Val) +} else { + stop(paste0("未知 method: ", method)) +} + +norm <- data.frame(gene_id = rownames(norm), norm, row.names = NULL, check.names = FALSE) +result <- merge(norm, stats, by = "gene_id", all = TRUE) + +# ---- 过滤低表达 + 去 NA ---------------------------------------------------- +ns <- length(sg$sample_name); sc <- 2:(ns + 1) +result <- result[rowSums(result[, sc] == 0) < ns * 0.5, ] +result <- result[rowMeans(result[, sc]) >= min_mean, ] +result <- result[!is.na(result$pvalue) & !is.na(result$padj), ] + +# ---- DEG 判定(默认 padj,可切 pvalue)------------------------------------- +metric <- result[[sig_metric]] +deg_all <- result[metric < pcut & abs(result$log2FoldChange) > log2fc, ] +deg_up <- deg_all[deg_all$log2FoldChange > log2fc, ] +deg_down <- deg_all[deg_all$log2FoldChange < -log2fc, ] + +rgraph_write_csv(result, p$outdir, "01.Dse2_result") +rgraph_write_csv(deg_all, p$outdir, "02.Deg_all") +rgraph_write_csv(deg_down, p$outdir, "03.Deg_down") +rgraph_write_csv(deg_up, p$outdir, "04.Deg_up") +cat("RGRAPH_DONE: diff", method, "| treatment=", treatment_g, "control=", control_g, + "| all=", nrow(result), "DEG=", nrow(deg_all), "up=", nrow(deg_up), + "down=", nrow(deg_down), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/distribution.R b/plugins/Presisitence/rgraph/rscripts/distribution.R new file mode 100644 index 0000000..335a362 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/distribution.R @@ -0,0 +1,53 @@ +# distribution.R —— 表达分布:箱线图 / 小提琴图 / 密度曲线 +# params: fpkm, sample_group, outdir, kind("box"|"violin"|"density"), +# palette, drop_zero_frac(默认0.5), dpi(300), width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2")) + +kind <- rgraph_opt(p, "kind", "box") +fpkm <- rgraph_read_matrix(p$fpkm) +sg <- rgraph_read_sample_group(p$sample_group) +rgraph_check_samples(fpkm, sg$sample_name) + +mat <- fpkm[, sg$sample_name, drop = FALSE] +# 去除 0 占比过高的行 +zf <- rgraph_opt(p, "drop_zero_frac", 0.5) +keep <- rowSums(mat == 0 | is.na(mat)) < ncol(mat) * zf +mat <- mat[keep, , drop = FALSE] +mat[is.na(mat)] <- 0 + +long <- as.data.frame(mat) |> + tidyr::pivot_longer(cols = dplyr::everything(), names_to = "Sample", values_to = "FPKM") |> + dplyr::mutate(logFPKM = log2(FPKM + 1)) +long$Group <- sg$group_name[match(long$Sample, sg$sample_name)] +long <- long |> + dplyr::mutate(Sample = factor(Sample, levels = unique(sg$sample_name)), + Group = factor(Group, levels = unique(sg$group_name))) + +cols <- rgraph_palette(length(unique(sg$group_name)), rgraph_opt(p, "palette", "course")) + +if (kind == "density") { + pl <- ggplot2::ggplot(long, ggplot2::aes(x = logFPKM, color = Group, group = Sample)) + + ggplot2::geom_density(linewidth = 0.6, alpha = 0.6) + + ggplot2::scale_color_manual(values = cols) + + ggplot2::labs(title = "Expression density", x = "log2(fpkm+1)", y = "Density") + + theme_rgraph() + w <- rgraph_opt(p, "width", 10); h <- rgraph_opt(p, "height", 6) +} else { + base <- ggplot2::ggplot(long, ggplot2::aes(x = Sample, y = logFPKM, fill = Group, colour = Group)) + + ggplot2::scale_color_manual(values = cols) + + ggplot2::scale_fill_manual(values = cols) + + ggplot2::labs(title = "Gene expression distribution", x = "Sample", y = "log2(fpkm+1)") + + theme_rgraph() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + geom <- if (kind == "violin") ggplot2::geom_violin(alpha = 0.5, linewidth = 0.6) + else ggplot2::geom_boxplot(alpha = 0.5, linewidth = 0.6, outlier.size = 0.4) + pl <- base + geom + w <- rgraph_opt(p, "width", 8); h <- rgraph_opt(p, "height", 8) +} + +nm <- rgraph_opt(p, "name", switch(kind, violin = "violinplot", density = "density", "boxplot")) +rgraph_save(pl, p$outdir, nm, width = w, height = h, dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: distribution", kind, "genes=", nrow(mat), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/enrich.R b/plugins/Presisitence/rgraph/rscripts/enrich.R new file mode 100644 index 0000000..aed52a2 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/enrich.R @@ -0,0 +1,82 @@ +# enrich.R —— GO / KEGG 富集分析(支持模式物种 OrgDb 与非模式物种自建基因集) +# params: +# gene_list(路径,列 gene_id[,log2FoldChange]), outdir, type("go"|"kegg") +# orgdb(如 "org.Hs.eg.db"), id_type("ENSEMBL"|"SYMBOL"|"ENTREZID"|"GID") +# ont("all"|"BP"|"CC"|"MF") # GO +# go_source("orgdb"|"custom") # custom: 从 OrgDb 的 GO 列自建 TERM2GENE(非模式) +# kegg_source("online"|"gmt"), kegg_species(如 "hsa"), gmt(路径), kegg_info(路径,可选) +# pAdjustMethod("BH"), pvalueCutoff(1), qvalueCutoff(1) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("clusterProfiler", "dplyr", "stringr")) + +type <- rgraph_opt(p, "type", "go") +padj_m <- rgraph_opt(p, "pAdjustMethod", "BH") +pcut <- rgraph_opt(p, "pvalueCutoff", 1) +qcut <- rgraph_opt(p, "qvalueCutoff", 1) +id_type <- rgraph_opt(p, "id_type", "ENSEMBL") + +gl <- read.csv(p$gene_list, stringsAsFactors = FALSE, check.names = FALSE) +genes_in <- gl$gene_id + +to_entrez <- function(ids) { + rgraph_library(c(p$orgdb)) + suppressWarnings(clusterProfiler::bitr(ids, fromType = id_type, toType = "ENTREZID", + OrgDb = p$orgdb))$ENTREZID +} + +if (type == "go") { + go_source <- rgraph_opt(p, "go_source", "orgdb") + ont <- rgraph_opt(p, "ont", "all") + if (go_source == "custom") { + rgraph_library(c(p$orgdb, "AnnotationDbi")) + odb <- get(p$orgdb) + g2go <- AnnotationDbi::select(odb, keys = AnnotationDbi::keys(odb), columns = c("GO")) + g2go <- na.omit(g2go[, c("GO", intersect(c("GID", "ENSEMBL"), names(g2go))[1])]) + names(g2go) <- c("term", "gene") + er <- clusterProfiler::enricher(genes_in, TERM2GENE = g2go[, c("term", "gene")], + pAdjustMethod = padj_m, pvalueCutoff = pcut, qvalueCutoff = qcut) + res <- er@result + } else { + rgraph_library(c(p$orgdb)) + genes <- if (id_type == "ENTREZID") genes_in else to_entrez(genes_in) + er <- clusterProfiler::enrichGO(gene = genes, OrgDb = p$orgdb, keyType = "ENTREZID", + ont = ont, pAdjustMethod = padj_m, + pvalueCutoff = pcut, qvalueCutoff = qcut) + res <- er@result + names(res)[names(res) == "ID"] <- "GOID" + } + names(res)[names(res) == "p.adjust"] <- "padj" + rgraph_write_csv(res, p$outdir, rgraph_opt(p, "name", "GO_enrich")) + cat("RGRAPH_DONE: enrich go terms=", nrow(res), "\n") +} else { + kegg_source <- rgraph_opt(p, "kegg_source", "online") + if (kegg_source == "gmt") { + t2g <- clusterProfiler::read.gmt(p$gmt) + er <- clusterProfiler::enricher(genes_in, TERM2GENE = t2g, pAdjustMethod = padj_m, + pvalueCutoff = pcut, qvalueCutoff = qcut) + res <- er@result + names(res)[names(res) == "p.adjust"] <- "padj" + } else { + sp <- rgraph_opt(p, "kegg_species", "hsa") + options(timeout = 600) + pe <- read.delim(sprintf("https://rest.kegg.jp/link/pathway/%s", sp), sep = "\t", + col.names = c("gene", "KEGGID")) + pe$gene <- str_remove(pe$gene, paste0(sp, ":")) + pe$KEGGID <- str_remove(pe$KEGGID, "path:") + t2g <- pe[, c("KEGGID", "gene")] + genes <- if (id_type == "ENTREZID") genes_in else to_entrez(genes_in) + er <- clusterProfiler::enricher(genes, TERM2GENE = t2g, pAdjustMethod = padj_m, + pvalueCutoff = pcut, qvalueCutoff = qcut) + res <- er@result + names(res)[names(res) == "ID"] <- "KEGGID" + names(res)[names(res) == "p.adjust"] <- "padj" + if (!is.null(p$kegg_info)) { + ki <- read.csv(p$kegg_info, stringsAsFactors = FALSE) + res <- merge(ki, res, by = "KEGGID", all.y = TRUE) + } + } + rgraph_write_csv(res, p$outdir, rgraph_opt(p, "name", "KEGG_enrich")) + cat("RGRAPH_DONE: enrich kegg pathways=", nrow(res), "\n") +} diff --git a/plugins/Presisitence/rgraph/rscripts/enrich_circos.R b/plugins/Presisitence/rgraph/rscripts/enrich_circos.R new file mode 100644 index 0000000..607faee --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/enrich_circos.R @@ -0,0 +1,94 @@ +# enrich_circos.R —— 富集分析圈图(circlize 多轨道:分类/背景基因/前景基因(上下调)/RichFactor) +# 输入: enrich(GO_enrich.csv 或 KEGG_enrich.csv), 需列: pvalue,RichFactor,Bg_gene,Count[,Up,Down] +# params: enrich, outdir, type("go"|"kegg"), topN(每类10), split_count(有Up/Down时按上下调拆分), +# show_axis(T), ifLog(F), circ1/2/3.color, dpi(300), width(8),height(8), name +# 说明: 图例需 ComplexHeatmap(缺失则跳过图例,圈图主体仍由 circlize 正常绘制)。 +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "circlize")) +has_ch <- requireNamespace("ComplexHeatmap", quietly = TRUE) + +type <- rgraph_opt(p, "type", "go") +group_col <- if (type == "kegg") "Category" else "ONTOLOGY" +id_col <- if (type == "kegg") "KEGGID" else "GOID" +topN <- rgraph_opt(p, "topN", 10) +show_axis <- rgraph_opt(p, "show_axis", TRUE) +ifLog <- rgraph_opt(p, "ifLog", FALSE) +c1 <- rgraph_opt(p, "circ1_color", c("#f7cb16", "#65c3fc", "#bfe046")) +c2 <- rgraph_opt(p, "circ2_color", c("#fee5d9", "#fb6a4a")) +c3 <- rgraph_opt(p, "circ3_color", c("#68a6e3", "#99CCFF")) + +df <- read.csv(p$enrich, stringsAsFactors = FALSE, check.names = FALSE) +for (col in c(group_col, id_col, "pvalue", "RichFactor", "Bg_gene", "Count")) + if (!col %in% names(df)) stop(paste0("富集表缺少列: ", col)) +split_count <- rgraph_opt(p, "split_count", all(c("Up", "Down") %in% names(df))) +df$.grp <- df[[group_col]]; df$.id <- df[[id_col]] +df$neg_log10_pvalue <- -log10(df$pvalue) + +dt <- df %>% dplyr::group_by(.grp) %>% dplyr::slice_min(pvalue, n = topN, with_ties = FALSE) %>% + dplyr::arrange(dplyr::desc(RichFactor), .by_group = TRUE) %>% dplyr::ungroup() +main.col <- c1[as.numeric(as.factor(dt$.grp))] +p_max <- ceiling(max(dt$neg_log10_pvalue)) +color_assign <- circlize::colorRamp2(0:p_max, grDevices::colorRampPalette(c2)(p_max + 1)) +bgmax <- max(dt$Bg_gene, na.rm = TRUE) + +df1 <- data.frame(id = dt$.id, start = 0, end = bgmax) +df2 <- data.frame(id = dt$.id, start = 0, end = if (show_axis) dt$Bg_gene else bgmax, + Bg_gene = dt$Bg_gene, col = color_assign(dt$neg_log10_pvalue)) +if (split_count) { + tempLong <- if (ifLog) log10(bgmax + 1) else bgmax + frac <- dt$Up / pmax(dt$Up + dt$Down, 1) + df3 <- rbind( + data.frame(id = dt$.id, start = 0, end = frac * tempLong, count = dt$Up, col = c3[1]), + data.frame(id = dt$.id, start = frac * tempLong, end = tempLong, count = dt$Down, col = c3[2])) + df3$count <- ifelse(df3$count == 0, "", df3$count) +} else { + df3 <- data.frame(id = dt$.id, start = 0, end = if (show_axis) dt$Count else bgmax, + count = dt$Count, col = c3[1]) +} +df4 <- data.frame(id = dt$.id, start = 0, end = bgmax, ratio = dt$RichFactor, col = main.col) +if (ifLog) { df1$end <- log10(df1$end + 1); df2$end <- log10(df2$end + 1) + if (!split_count) df3$end <- log10(df3$end + 1); df4$end <- log10(df4$end + 1) } + +draw_fn <- function() { + par(omi = c(0.1, 0.1, 0.1, 1.5)) + circlize::circos.par(track.margin = c(0.01, 0.01)) + circlize::circos.genomicInitialize(df1, plotType = "none") + circlize::circos.trackPlotRegion(ylim = c(0, 1), track.height = 0.08, bg.border = NA, bg.col = main.col, + panel.fun = function(x, y) { + circlize::circos.text(mean(circlize::get.cell.meta.data("xlim")), 0.5, + circlize::get.cell.meta.data("sector.index"), cex = rgraph_opt(p, "circ1_size", 0.45), + facing = "bending.inside", niceFacing = TRUE) }) + if (show_axis) for (si in circlize::get.all.sector.index()) + circlize::circos.axis(h = "top", labels.cex = rgraph_opt(p, "axis_size", 0.3), sector.index = si, + track.index = 1, major.at = pretty(c(0, max(df1$end, na.rm = TRUE)), n = 5), labels.facing = "clockwise") + circlize::circos.genomicTrack(df2, ylim = c(0, 1), track.height = 0.1, bg.border = "white", + panel.fun = function(region, value, ...) { + circlize::circos.genomicRect(region, value, ytop = 0, ybottom = 1, col = value[, 2], border = NA, ...) + circlize::circos.genomicText(region, value, y = 0.5, labels = value[, 1], adj = c(0.5, 0.5), + cex = rgraph_opt(p, "circ2_size", 0.5), ...) }) + circlize::circos.genomicTrack(df3, ylim = c(0, 1), track.height = 0.1, bg.border = "white", + panel.fun = function(region, value, ...) { + circlize::circos.genomicRect(region, value, ytop = 0, ybottom = 1, col = value[, 2], border = NA, ...) + circlize::circos.genomicText(region, value, y = 0.5, labels = value[, 1], adj = c(0.5, 0.5), + cex = rgraph_opt(p, "circ3_size", 0.3), ...) }) + circlize::circos.genomicTrack(df4, ylim = c(0, max(dt$RichFactor, na.rm = TRUE)), track.height = 0.35, + bg.border = "white", bg.col = "grey97", + panel.fun = function(region, value, ...) + circlize::circos.genomicRect(region, value, ytop = 0, ybottom = value[, 1], col = value[, 2], border = NA, ...)) + circlize::circos.clear() + if (has_ch) { + ml <- ComplexHeatmap::Legend(labels = unique(dt$.grp), type = "points", pch = 15, + legend_gp = grid::gpar(col = c1), title = group_col, size = grid::unit(3, "mm")) + lp <- ComplexHeatmap::Legend(col_fun = circlize::colorRamp2(round(seq(0, p_max, length.out = 6), 0), + grDevices::colorRampPalette(c2)(6)), title = "-log10(pvalue)") + ComplexHeatmap::draw(ComplexHeatmap::packLegend(ml, lp), + x = grid::unit(1, "snpc") * 0.85, y = grid::unit(1, "snpc") * 0.5, just = "left") + } +} + +rgraph_save_base(p$outdir, rgraph_opt(p, "name", paste0(toupper(type), "_loop")), draw_fn, + width = rgraph_opt(p, "width", 8), height = rgraph_opt(p, "height", 8), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: enrich_circos", type, "terms=", nrow(dt), "legend=", has_ch, "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/gene_bar.R b/plugins/Presisitence/rgraph/rscripts/gene_bar.R new file mode 100644 index 0000000..776daf0 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/gene_bar.R @@ -0,0 +1,47 @@ +# gene_bar.R —— 单基因表达柱形图(对照 vs 处理,均值±SEM,显著性星号) +# params: deg(Deg_all.csv, 含 pvalue 与各样本归一化列), sample_group(需 TvsC 与 group/group_name), +# gene_list(CSV, 列 gene_id) 或 genes(基因向量), outdir, dpi(300), width(4),height(5) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2")) + +deg <- read.csv(p$deg, check.names = FALSE) +sg <- rgraph_read_sample_group(p$sample_group) +if (!is.null(p$genes)) genes <- p$genes else genes <- read.csv(p$gene_list, check.names = FALSE)$gene_id + +deg$difference <- dplyr::case_when(deg$pvalue < 0.001 ~ "***", deg$pvalue < 0.01 ~ "**", + deg$pvalue < 0.05 ~ "*", TRUE ~ "") +ctrl <- sg$sample_name[sg$TvsC == "control"] +trt <- sg$sample_name[sg$TvsC == "treatment"] +if (length(ctrl) == 0 || length(trt) == 0) { # 无 TvsC 时按前两组 + gs <- unique(sg$group); trt <- sg$sample_name[sg$group == gs[1]]; ctrl <- sg$sample_name[sg$group == gs[2]] +} +ctrl_name <- sg$group_name[match(ctrl[1], sg$sample_name)] +trt_name <- sg$group_name[match(trt[1], sg$sample_name)] + +summ <- function(row, cols) { + v <- as.numeric(row[cols]); c(mean = mean(v), sem = ifelse(length(v) > 1, sd(v) / sqrt(length(v)), 0)) +} +n_ok <- 0 +for (g in genes) { + cg <- deg[deg$gene_id == g, ] + if (nrow(cg) == 0) next + sc <- summ(cg[1, ], ctrl); st <- summ(cg[1, ], trt) + pd <- data.frame(group = factor(c(ctrl_name, trt_name), levels = c(ctrl_name, trt_name)), + mean_exp = c(sc["mean"], st["mean"]), sem = c(sc["sem"], st["sem"]), + difference = c("", cg$difference[1])) + pl <- ggplot2::ggplot(pd, ggplot2::aes(x = group, y = mean_exp, fill = group)) + + ggplot2::geom_col(width = 0.6) + + ggplot2::scale_fill_manual(values = setNames(c("#1f77b4", "#ff7f0e"), c(ctrl_name, trt_name)), name = "Group") + + ggplot2::geom_errorbar(ggplot2::aes(ymin = mean_exp - sem, ymax = mean_exp + sem), width = 0.2) + + ggplot2::geom_text(ggplot2::aes(y = mean_exp + sem, label = difference), vjust = 0, size = 5) + + ggplot2::labs(title = g, y = "Normalized Count", x = NULL) + + theme_rgraph() + + ggplot2::scale_y_continuous(expand = ggplot2::expansion(mult = c(0, 0.1))) + rgraph_save(pl, p$outdir, paste0(g, "_Expression"), + width = rgraph_opt(p, "width", 4), height = rgraph_opt(p, "height", 5), + dpi = rgraph_opt(p, "dpi", 300)) + n_ok <- n_ok + 1 +} +cat("RGRAPH_DONE: gene_bar genes=", n_ok, "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/gene_correlation.R b/plugins/Presisitence/rgraph/rscripts/gene_correlation.R new file mode 100644 index 0000000..85f2f35 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/gene_correlation.R @@ -0,0 +1,44 @@ +# gene_correlation.R —— 基因-基因相关性热图(气泡热图,圆点大小=|r|) +# params: fpkm(gene_fpkm.csv), sample_group, outdir, +# genes1(向量) 或 gene_list1(CSV), genes2(向量) 或 gene_list2(CSV, 缺省=genes1), +# dpi(300), width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "ggplot2", "reshape2", "scales")) + +fpkm <- rgraph_read_matrix(p$fpkm) +sg <- rgraph_read_sample_group(p$sample_group) +g1 <- if (!is.null(p$genes1)) p$genes1 else read.csv(p$gene_list1, check.names = FALSE)$gene_id +g2 <- if (!is.null(p$genes2)) p$genes2 else if (!is.null(p$gene_list2)) + read.csv(p$gene_list2, check.names = FALSE)$gene_id else g1 + +sub <- fpkm[fpkm$gene_id %in% unique(c(g1, g2)), c("gene_id", sg$sample_name)] +rownames(sub) <- sub$gene_id; sub$gene_id <- NULL +sub[is.na(sub)] <- 0 +logm <- log2(sub + 1) +cm <- cor(t(logm), use = "pairwise.complete.obs") # 基因 x 基因 +cm <- cm[rownames(cm) %in% g1, colnames(cm) %in% g2, drop = FALSE] + +rgraph_write_csv(as.data.frame(cbind(gene_id = rownames(cm), cm)), p$outdir, "correlation") + +df <- reshape2::melt(cm); df$size <- abs(df$value) +pl <- ggplot2::ggplot(df, ggplot2::aes(x = Var2, y = Var1)) + + ggplot2::geom_tile(color = "black", fill = "white", linewidth = 0.4) + + ggplot2::geom_point(ggplot2::aes(fill = value, size = size, color = value), shape = 21, + stroke = 0.4, alpha = 0.85) + + ggplot2::scale_color_gradient2(low = "#2167ad", mid = "#d8d4d4", high = "red", guide = "none") + + ggplot2::scale_size(range = c(4, 12), guide = "none") + + ggplot2::scale_fill_gradientn(colors = grDevices::colorRampPalette(c("#2167ad", "white", "red"))(100), + limits = c(-1, 1), name = "r") + + ggplot2::labs(title = "Pearson correlation between genes", x = NULL, y = NULL) + + theme_rgraph() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + + ggplot2::coord_fixed() + +nr <- nrow(cm); nc <- ncol(cm) +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "gene_correlation"), + width = rgraph_opt(p, "width", max(6, nc * 0.6)), + height = rgraph_opt(p, "height", max(6, nr * 0.6)), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: gene_correlation rows=", nr, "cols=", nc, "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/go_plot.R b/plugins/Presisitence/rgraph/rscripts/go_plot.R new file mode 100644 index 0000000..a43f3bb --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/go_plot.R @@ -0,0 +1,63 @@ +# go_plot.R —— GO 富集条形图/气泡图(BP/CC/MF 三面板,纵向拼接) +# 输入: GO_enrich.csv (列: ONTOLOGY, Description, GeneRatio, pvalue, Count) +# params: enrich, outdir, kind("bar"|"dot"), top_n(每类10), desc_len(120), +# dpi(300), width(11), height(10), name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "ggplot2", "stringr", "patchwork")) + +parse_ratio <- function(x) vapply(x, function(s) { + s <- as.character(s) + if (grepl("/", s)) { a <- as.numeric(strsplit(s, "/")[[1]]); a[1] / a[2] } else as.numeric(s) +}, numeric(1)) + +kind <- rgraph_opt(p, "kind", "bar") +top_n <- rgraph_opt(p, "top_n", 10) +desc_len <- rgraph_opt(p, "desc_len", 120) + +d <- read.csv(p$enrich, stringsAsFactors = FALSE, check.names = FALSE) +if (!"ONTOLOGY" %in% names(d)) stop("GO_enrich 缺少 ONTOLOGY 列") +d$neg_log10_pvalue <- -log10(d$pvalue) +d$GeneRatio <- parse_ratio(d$GeneRatio) +d$Description <- substr(d$Description, 1, desc_len) + +d <- d %>% dplyr::group_by(ONTOLOGY) %>% dplyr::arrange(pvalue) %>% + dplyr::slice(seq_len(min(top_n, dplyr::n()))) %>% dplyr::ungroup() + +# 三类各自的低/高配色(沿用课程) +pal <- list(BP = c("#f2c7b9", "#e74716"), CC = c("#90dee7", "#0a8a99"), MF = c("#b2eecf", "#41ae76")) +onts <- intersect(c("BP", "CC", "MF"), unique(d$ONTOLOGY)) + +panel <- function(ont, show_x) { + dd <- d %>% dplyr::filter(ONTOLOGY == ont) %>% + dplyr::arrange(if (kind == "dot") desc(GeneRatio) else desc(Count)) %>% + dplyr::mutate(Description = factor(Description, levels = rev(Description))) + cols <- pal[[ont]] + base <- ggplot2::ggplot(dd) + + ggplot2::facet_grid(ONTOLOGY ~ ., scales = "free_y", space = "free_y") + + theme_rgraph() + + ggplot2::theme(strip.text = ggplot2::element_text(face = "bold", size = 12, color = "white"), + strip.background = ggplot2::element_rect(fill = cols[2], color = "black"), + axis.text.y = ggplot2::element_text(size = 11, color = cols[2]), + axis.text.x = if (show_x) ggplot2::element_text(size = 10) else ggplot2::element_blank()) + + ggplot2::scale_y_discrete(labels = function(y) str_wrap(y, width = 60)) + + ggplot2::labs(x = if (show_x) (if (kind == "dot") "GeneRatio" else "-log10(pvalue)") else NULL, y = NULL) + if (kind == "dot") { + base + ggplot2::geom_point(ggplot2::aes(x = GeneRatio, y = Description, + fill = neg_log10_pvalue, size = Count), shape = 21, colour = cols[2]) + + ggplot2::scale_fill_gradient(low = cols[1], high = cols[2], name = "-log10(pvalue)") + + ggplot2::scale_size_continuous(range = c(3, 6), name = "Count") + } else { + base + ggplot2::geom_col(ggplot2::aes(x = neg_log10_pvalue, y = Description, fill = Count), width = 0.7) + + ggplot2::scale_fill_gradient(low = cols[1], high = cols[2], name = "Count") + } +} + +plots <- lapply(seq_along(onts), function(i) panel(onts[i], show_x = (i == length(onts)))) +pl <- Reduce(`/`, plots) + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", paste0("GO_", kind)), + width = rgraph_opt(p, "width", 11), height = rgraph_opt(p, "height", 10), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: go_plot", kind, "ontologies=", paste(onts, collapse = "/"), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/gsea.R b/plugins/Presisitence/rgraph/rscripts/gsea.R new file mode 100644 index 0000000..495042b --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/gsea.R @@ -0,0 +1,49 @@ +# gsea.R —— GSEA 分析(基于差异结果 log2FoldChange 排序 + gmt 基因集) +# params: result(Dse2_result.csv), gmt(路径), outdir, minGSSize(5), maxGSSize(1000), +# pvalueCutoff(1), desc_map(可选 CSV: 列 ID,Description), top_n(每方向出图数,默认5), +# draw(bool,默认T), seed(123) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("clusterProfiler", "enrichplot", "ggplot2", "dplyr")) + +deg <- read.csv(p$result, check.names = FALSE) +deg <- deg %>% dplyr::select(gene_id, log2FoldChange) %>% dplyr::filter(!is.na(log2FoldChange)) +gene_rank <- sort(setNames(deg$log2FoldChange, deg$gene_id), decreasing = TRUE) + +t2g <- clusterProfiler::read.gmt(p$gmt) +set.seed(rgraph_opt(p, "seed", 123)) +gsea <- clusterProfiler::GSEA(geneList = gene_rank, TERM2GENE = t2g, + pvalueCutoff = rgraph_opt(p, "pvalueCutoff", 1), + pAdjustMethod = "BH", + minGSSize = rgraph_opt(p, "minGSSize", 5), + maxGSSize = rgraph_opt(p, "maxGSSize", 1000)) +res <- gsea@result + +if (!is.null(p$desc_map)) { + dm <- read.csv(p$desc_map, stringsAsFactors = FALSE) + if (all(c("ID", "Description") %in% names(dm))) { + res$Description <- dm$Description[match(res$ID, dm$ID)] + } +} +rgraph_write_csv(res, p$outdir, rgraph_opt(p, "name", "GSEA_result")) + +# ---- 各方向 top-N 通路 ES 图 ----------------------------------------------- +if (isTRUE(rgraph_opt(p, "draw", TRUE)) && nrow(res) > 0) { + top_n <- rgraph_opt(p, "top_n", 5) + draw_dir <- function(sub, sign_dir) { + if (nrow(sub) == 0) return(invisible()) + sub <- sub[order(sub$pvalue), ] + sub <- utils::head(sub, top_n) + for (id in sub$ID) { + idx <- which(gsea@result$ID == id) + ttl <- paste0(res$Description[match(id, res$ID)], " [", id, "]") + pl <- enrichplot::gseaplot2(gsea, geneSetID = idx, ES_geom = "line", title = ttl) + rgraph_save(pl, file.path(p$outdir, sign_dir), gsub("[^A-Za-z0-9]", "_", id), + width = 8, height = 8, dpi = rgraph_opt(p, "dpi", 300)) + } + } + draw_dir(res[res$NES > 0, ], "UP") + draw_dir(res[res$NES < 0, ], "DOWN") +} +cat("RGRAPH_DONE: gsea genesets=", nrow(res), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/gsea_mountain.R b/plugins/Presisitence/rgraph/rscripts/gsea_mountain.R new file mode 100644 index 0000000..8192941 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/gsea_mountain.R @@ -0,0 +1,46 @@ +# gsea_mountain.R —— GSEA 山峦图(多通路 core_enrichment 基因的 log2FC 分布 + NES 点) +# params: gsea_result(gsea.R 产出的 CSV, 含 ID,Description,NES,pvalue,core_enrichment), +# result(Dse2_result.csv, 取 gene log2FC), outdir, +# geneset_ids(可选,要画的通路 ID 向量; 缺省取 pvalue 最小的前 n_top 个), +# n_top(15), dpi(300), width(15),height(13), name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2", "ggridges", "ggnewscale")) + +deg <- read.csv(p$result, check.names = FALSE) +gene_rank <- setNames(deg$log2FoldChange, deg$gene_id) + +gr <- read.csv(p$gsea_result, stringsAsFactors = FALSE, check.names = FALSE) +if (!is.null(p$geneset_ids)) { + gr <- gr[gr$ID %in% p$geneset_ids, ] +} else { + gr <- gr[order(gr$pvalue), ] + gr <- utils::head(gr, rgraph_opt(p, "n_top", 15)) +} +if (nrow(gr) == 0) stop("没有可用于山峦图的通路") + +df <- gr %>% + dplyr::select(Description, NES, pvalue, core_enrichment) %>% + tidyr::separate_rows(core_enrichment, sep = "/") %>% + dplyr::mutate(log2FoldChange = gene_rank[core_enrichment], + neg_log10_pvalue = -log10(pvalue)) %>% + dplyr::filter(!is.na(log2FoldChange)) + +min_fc <- min(df$log2FoldChange) - 1 +pl <- ggplot2::ggplot(df) + + ggridges::geom_density_ridges(ggplot2::aes(x = log2FoldChange, y = Description, + fill = neg_log10_pvalue), scale = 0.9) + + ggplot2::scale_fill_gradient(low = "#e6e2f9", high = "#998ec3", name = "-log10(pvalue)") + + ggnewscale::new_scale_fill() + + ggplot2::geom_point(ggplot2::aes(x = min_fc, y = Description, size = abs(NES), fill = NES), + shape = 21) + + ggplot2::scale_fill_gradient2(low = "#4575b4", mid = "#f7f7f7", high = "#d73027", name = "NES") + + ggplot2::scale_size(range = c(5, 10), guide = "none") + + theme_rgraph() + + ggplot2::labs(x = "log2FoldChange", y = NULL) + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "Mountain_plot"), + width = rgraph_opt(p, "width", 15), height = rgraph_opt(p, "height", 13), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: gsea_mountain genesets=", nrow(gr), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/heatmap.R b/plugins/Presisitence/rgraph/rscripts/heatmap.R new file mode 100644 index 0000000..34ceeca --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/heatmap.R @@ -0,0 +1,73 @@ +# heatmap.R —— 差异基因/指定基因表达聚类热图 +# 引擎: pheatmap(默认) 或 complexheatmap(按分组加顶部色块标签、列不聚类) +# params: fpkm, sample_group, outdir, deg(可选 Deg_all.csv 路径), genes(可选基因向量), +# engine("pheatmap"|"complexheatmap"), scale("row"|"column"|"none"), +# cluster_rows(T), cluster_cols(T), show_rownames(默认自动), +# low/mid/high 颜色, dpi(300), width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr")) + +fpkm <- rgraph_read_matrix(p$fpkm) +sg <- rgraph_read_sample_group(p$sample_group) +rgraph_check_samples(fpkm, sg$sample_name) + +# 选基因 +if (!is.null(p$genes)) { + genes <- p$genes +} else if (!is.null(p$deg)) { + deg <- read.csv(p$deg, check.names = FALSE) + genes <- deg$gene_id +} else { + genes <- fpkm$gene_id +} +sub <- fpkm[fpkm$gene_id %in% genes, c("gene_id", sg$sample_name)] +rownames(sub) <- sub$gene_id; sub$gene_id <- NULL +sub <- as.matrix(sub); sub[is.na(sub)] <- 0 +sub <- log2(sub + 1) +sub <- sub[rowSums(sub) != 0, , drop = FALSE] +sub <- sub[, colSums(sub) != 0, drop = FALSE] +if (nrow(sub) < 2) stop("可用于热图的基因不足 2 个") + +scale_by <- rgraph_opt(p, "scale", "row") +show_rn <- rgraph_opt(p, "show_rownames", nrow(sub) <= 50) +low <- rgraph_opt(p, "low", "blue"); mid <- rgraph_opt(p, "mid", "white"); high <- rgraph_opt(p, "high", "red") +nm <- rgraph_opt(p, "name", "heatmap") +w <- rgraph_opt(p, "width", 8); h <- rgraph_opt(p, "height", 8); dpi <- rgraph_opt(p, "dpi", 300) +engine <- rgraph_opt(p, "engine", "pheatmap") + +if (engine == "complexheatmap") { + rgraph_library(c("ComplexHeatmap", "circlize")) + sc <- t(scale(t(sub))) + col_fun <- circlize::colorRamp2(c(min(sc), 0, max(sc)), c(low, mid, high)) + grp <- factor(sg$group_name[match(colnames(sub), sg$sample_name)], + levels = unique(sg$group_name)) + gcol <- rgraph_palette(nlevels(grp), rgraph_opt(p, "palette", "course")) + names(gcol) <- levels(grp) + ha <- ComplexHeatmap::HeatmapAnnotation( + Group = ComplexHeatmap::anno_block(gp = grid::gpar(fill = gcol, col = NA), + labels = levels(grp), labels_gp = grid::gpar(col = "white", fontsize = 8, fontface = "bold")), + show_annotation_name = FALSE, height = grid::unit(0.3, "cm")) + ht <- ComplexHeatmap::Heatmap(sc, col = col_fun, cluster_rows = TRUE, cluster_columns = FALSE, + column_split = grp, show_row_names = show_rn, top_annotation = ha, + column_names_rot = 45, name = "Z-score") + rgraph_save_base(p$outdir, nm, function() ComplexHeatmap::draw(ht), + width = w, height = h, dpi = dpi) +} else { + rgraph_library(c("pheatmap")) + pal <- grDevices::colorRampPalette(c(low, mid, high))(100) + rgraph_ensure_dir(p$outdir) + for (fmt in c("png", "pdf")) { + f <- file.path(p$outdir, paste0(nm, ".", fmt)) + pheatmap::pheatmap(sub, scale = scale_by, + cluster_rows = rgraph_opt(p, "cluster_rows", TRUE), + cluster_cols = rgraph_opt(p, "cluster_cols", TRUE), + color = pal, border_color = NA, + show_rownames = show_rn, show_colnames = TRUE, + fontsize_col = 10, angle_col = 45, + width = w, height = h, filename = f) + rgraph_emit(f) + } +} +cat("RGRAPH_DONE: heatmap engine=", engine, "genes=", nrow(sub), "samples=", ncol(sub), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/interaction_diff.R b/plugins/Presisitence/rgraph/rscripts/interaction_diff.R new file mode 100644 index 0000000..722c0a5 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/interaction_diff.R @@ -0,0 +1,55 @@ +# interaction_diff.R —— 两因子交互作用差异分析(DESeq2 ~g1+g2+g1:g2,提取交互项) +# 输入: count, sample_group(需 group_name1,group_name2,TvsC1,TvsC2) +# params: count, sample_group, outdir, sig_metric("padj"|"pvalue"), pcut(0.05), log2fc(1), +# min_mean_count(2) +# 输出: 01.DEseq2_result.csv + 02/03/04 Deg_all/down/up(可直接喂 rgraph_volcano) +# 注: 需 DESeq2(本机默认缺,按提示安装)。 +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "DESeq2")) + +sig_metric <- rgraph_opt(p, "sig_metric", "padj") +pcut <- rgraph_opt(p, "pcut", 0.05) +log2fc <- rgraph_opt(p, "log2fc", 1) +min_mean <- rgraph_opt(p, "min_mean_count", 2) + +count <- rgraph_read_matrix(p$count) +sg <- read.csv(p$sample_group, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) +sg <- sg[!is.na(sg$sample_name) & sg$sample_name != "", , drop = FALSE] +for (col in c("group_name1", "group_name2", "TvsC1", "TvsC2")) + if (!col %in% names(sg)) stop(paste0("sample_group 缺少列: ", col, "(两因子交互需要)")) +rgraph_check_samples(count, sg$sample_name) + +t1 <- unique(sg$group_name1[sg$TvsC1 == "treatment"]); c1 <- unique(sg$group_name1[sg$TvsC1 == "control"]) +t2 <- unique(sg$group_name2[sg$TvsC2 == "treatment"]); c2 <- unique(sg$group_name2[sg$TvsC2 == "control"]) + +mat <- count[, c("gene_id", sg$sample_name)] +rownames(mat) <- mat$gene_id; mat$gene_id <- NULL; mat[is.na(mat)] <- 0 +cd <- data.frame(row.names = sg$sample_name, + g1 = factor(sg$group_name1, levels = c(c1, t1)), + g2 = factor(sg$group_name2, levels = c(c2, t2))) +dds <- DESeq2::DESeqDataSetFromMatrix(round(as.matrix(mat)), cd, design = ~ g1 + g2 + g1:g2) +dds <- DESeq2::DESeq(dds) +norm <- DESeq2::counts(dds, normalized = TRUE) +norm <- data.frame(gene_id = rownames(norm), norm, row.names = NULL, check.names = FALSE) + +# 自动定位交互项系数(形如 g1XXX.g2YYY) +inter <- grep("g1.*g2", DESeq2::resultsNames(dds), value = TRUE) +if (length(inter) == 0) stop("未找到交互项系数,请检查分组设计") +res <- as.data.frame(DESeq2::results(dds, name = inter[length(inter)])) +res <- data.frame(gene_id = rownames(res), res[, c("log2FoldChange", "pvalue", "padj")], row.names = NULL) +result <- merge(norm, res, by = "gene_id", all = TRUE) + +ns <- length(sg$sample_name); sc <- 2:(ns + 1) +result <- result[rowSums(result[, sc] == 0) < ns * 0.5, ] +result <- result[rowMeans(result[, sc]) >= min_mean, ] +result <- result[!is.na(result$pvalue) & !is.na(result$padj), ] +metric <- result[[sig_metric]] +deg_all <- result[metric < pcut & abs(result$log2FoldChange) > log2fc, ] + +rgraph_write_csv(result, p$outdir, "01.DEseq2_result") +rgraph_write_csv(deg_all, p$outdir, "02.Deg_all") +rgraph_write_csv(deg_all[deg_all$log2FoldChange < -log2fc, ], p$outdir, "03.Deg_down") +rgraph_write_csv(deg_all[deg_all$log2FoldChange > log2fc, ], p$outdir, "04.Deg_up") +cat("RGRAPH_DONE: interaction_diff term=", inter[length(inter)], "DEG=", nrow(deg_all), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/interactive_volcano.R b/plugins/Presisitence/rgraph/rscripts/interactive_volcano.R new file mode 100644 index 0000000..0db00d3 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/interactive_volcano.R @@ -0,0 +1,42 @@ +# interactive_volcano.R —— 交互式火山图(plotly HTML,悬停显示基因/log2FC/-log10p) +# params: result(Degs.csv 或 Dse2_result.csv), outdir, sig_metric("pvalue"|"padj"), +# pcut(0.05), log2fc(1), color_scheme("rb"|"rg"), include_no(bool,默认T), name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("ggplot2", "dplyr", "plotly", "htmlwidgets")) + +sig <- rgraph_opt(p, "sig_metric", "pvalue") +pcut <- rgraph_opt(p, "pcut", 0.05) +log2fc <- rgraph_opt(p, "log2fc", 1) + +d <- read.csv(p$result, header = TRUE, check.names = FALSE) +d <- d[!is.na(d[[sig]]) & !is.na(d$log2FoldChange), ] +m <- d[[sig]] +d$color <- ifelse(m < pcut & d$log2FoldChange > log2fc, "Up", + ifelse(m < pcut & d$log2FoldChange < -log2fc, "Down", "No")) +if (!isTRUE(rgraph_opt(p, "include_no", TRUE))) d <- d[d$color != "No", ] +cnt <- table(factor(d$color, levels = c("Up", "Down", "No"))) +labs3 <- c(Up = sprintf("Up (%d)", cnt["Up"]), Down = sprintf("Down (%d)", cnt["Down"]), + No = sprintf("No (%d)", cnt["No"])) +cols <- rgraph_updown_colors(rgraph_opt(p, "color_scheme", "rb")) + +gcol <- if ("gene_id" %in% names(d)) d$gene_id else rownames(d) +pl <- ggplot2::ggplot(d, ggplot2::aes(x = log2FoldChange, y = -log10(.data[[sig]]), color = color, + text = paste0("Gene: ", gcol, "
log2FC: ", round(log2FoldChange, 3), + "
-log10(", sig, "): ", round(-log10(.data[[sig]]), 3)))) + + ggplot2::geom_point(alpha = 0.5, size = 1) + + ggplot2::scale_color_manual(values = cols, breaks = c("Up", "Down", "No"), labels = labs3) + + ggplot2::geom_hline(yintercept = -log10(pcut), linetype = "dashed", color = "grey") + + ggplot2::geom_vline(xintercept = c(-log2fc, log2fc), linetype = "dashed", color = "grey") + + ggplot2::labs(x = "log2 (Fold Change)", y = sprintf("-log10(%s)", sig), color = "Trend") + + theme_rgraph() + +ip <- plotly::ggplotly(pl, tooltip = "text", width = 850, height = 620) +rgraph_ensure_dir(p$outdir) +f <- file.path(p$outdir, paste0(rgraph_opt(p, "name", "volcano_interactive"), ".html")) +# selfcontained=TRUE 需 pandoc;默认 FALSE(会伴生一个 _files 依赖目录),无需 pandoc +htmlwidgets::saveWidget(ip, f, selfcontained = isTRUE(rgraph_opt(p, "selfcontained", FALSE)), + title = "Interactive Volcano Plot") +rgraph_emit(f) +cat("RGRAPH_DONE: interactive_volcano up=", cnt["Up"], "down=", cnt["Down"], "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/kegg_plot.R b/plugins/Presisitence/rgraph/rscripts/kegg_plot.R new file mode 100644 index 0000000..0b6e8a5 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/kegg_plot.R @@ -0,0 +1,44 @@ +# kegg_plot.R —— KEGG 富集条形图/气泡图 +# 输入: KEGG_enrich.csv (列: Description, GeneRatio, pvalue, Count) +# params: enrich, outdir, kind("bar"|"dot"), top_n(30), dpi(300), width(9),height(7), name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "ggplot2", "stringr")) + +parse_ratio <- function(x) vapply(x, function(s) { + s <- as.character(s) + if (grepl("/", s)) { a <- as.numeric(strsplit(s, "/")[[1]]); a[1] / a[2] } else as.numeric(s) +}, numeric(1)) + +kind <- rgraph_opt(p, "kind", "dot") +top_n <- rgraph_opt(p, "top_n", 30) + +d <- read.csv(p$enrich, stringsAsFactors = FALSE, check.names = FALSE) +d$neg_log10_pvalue <- -log10(d$pvalue) +d$GeneRatio <- parse_ratio(d$GeneRatio) +d <- d %>% dplyr::arrange(pvalue) %>% dplyr::slice(seq_len(min(top_n, dplyr::n()))) +d$Description <- str_remove(d$Description, "\\(.*") # 去掉括号后缀 +ord <- if (kind == "dot") order(d$neg_log10_pvalue, decreasing = TRUE) else order(d$GeneRatio, decreasing = TRUE) +d$Description <- factor(d$Description, levels = rev(d$Description[ord])) + +if (kind == "dot") { + pl <- ggplot2::ggplot(d) + + ggplot2::geom_point(ggplot2::aes(x = GeneRatio, y = Description, fill = neg_log10_pvalue, + size = Count), shape = 21, colour = "grey20") + + ggplot2::scale_fill_gradient(low = "yellow", high = "red", name = "-log10(pvalue)") + + ggplot2::scale_size(range = c(2, 6), name = "Count") + + ggplot2::labs(x = "GeneRatio", y = NULL) +} else { + pl <- ggplot2::ggplot(d) + + ggplot2::geom_col(ggplot2::aes(x = neg_log10_pvalue, y = Description, fill = GeneRatio)) + + ggplot2::scale_fill_gradient(low = "yellow", high = "red", name = "GeneRatio") + + ggplot2::labs(x = "-log10(pvalue)", y = NULL) +} +pl <- pl + theme_rgraph() + + ggplot2::scale_y_discrete(labels = function(y) str_wrap(y, width = 55)) + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", paste0("KEGG_", kind)), + width = rgraph_opt(p, "width", 9), height = rgraph_opt(p, "height", 7), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: kegg_plot", kind, "pathways=", nrow(d), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/kmeans.R b/plugins/Presisitence/rgraph/rscripts/kmeans.R new file mode 100644 index 0000000..364bb99 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/kmeans.R @@ -0,0 +1,86 @@ +# kmeans.R —— 目标基因 K-means 聚类:肘部图 + 聚类热图 + 各簇表达趋势 + 成员表 +# params: fpkm, sample_group, outdir, gene_list(CSV, gene_id) 或 genes(向量), +# k(簇数,默认4), kmax(肘部图最大k,默认10), seed(123), dpi(300) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2", "pheatmap")) + +k <- rgraph_opt(p, "k", 4) +fpkm <- rgraph_read_matrix(p$fpkm) +sg <- rgraph_read_sample_group(p$sample_group) +genes <- if (!is.null(p$genes)) p$genes else read.csv(p$gene_list, check.names = FALSE)$gene_id + +sub <- fpkm[fpkm$gene_id %in% genes, c("gene_id", sg$sample_name), drop = FALSE] +rownames(sub) <- sub$gene_id; sub$gene_id <- NULL +m <- as.matrix(sub); m[is.na(m)] <- 0 +logm <- log2(m + 1) +logm <- logm[rowSums(logm) != 0, , drop = FALSE] +logm <- logm[, colSums(logm) != 0, drop = FALSE] +z <- t(scale(t(logm))) +z[is.na(z)] <- 0 +if (nrow(z) < k) stop("基因数少于簇数 k") + +# ---- 肘部图(base wss,不依赖 factoextra)--------------------------------- +set.seed(rgraph_opt(p, "seed", 123)) +kmax <- min(rgraph_opt(p, "kmax", 10), nrow(z) - 1) +wss <- sapply(1:kmax, function(kk) kmeans(z, centers = kk, nstart = 10)$tot.withinss) +elbow <- data.frame(k = 1:kmax, wss = wss) +pl_e <- ggplot2::ggplot(elbow, ggplot2::aes(k, wss)) + + ggplot2::geom_line(color = "#4575b4") + ggplot2::geom_point(size = 2, color = "#4575b4") + + ggplot2::scale_x_continuous(breaks = 1:kmax) + + ggplot2::labs(x = "Number of clusters k", y = "Total within-cluster SS", + title = "Elbow method") + theme_rgraph() +rgraph_save(pl_e, p$outdir, "elbow", width = 7, height = 5, dpi = rgraph_opt(p, "dpi", 300)) + +# ---- K-means + 成员表 ------------------------------------------------------- +set.seed(rgraph_opt(p, "seed", 123)) +km <- kmeans(z, centers = k, nstart = 25) +for (cn in 1:k) + rgraph_write_csv(data.frame(gene_id = names(km$cluster[km$cluster == cn])), + p$outdir, paste0("cluster_", cn, "_genes")) + +# ---- 聚类热图(pheatmap: 行按簇排序 + 行/列注释)-------------------------- +ord <- order(km$cluster) +zo <- z[ord, , drop = FALSE] +cl <- km$cluster[ord] +ann_row <- data.frame(Cluster = factor(cl)); rownames(ann_row) <- rownames(zo) +ann_col <- data.frame(Group = factor(sg$group_name[match(colnames(zo), sg$sample_name)], + levels = unique(sg$group_name))) +rownames(ann_col) <- colnames(zo) +gcol <- rgraph_palette(nlevels(ann_col$Group), rgraph_opt(p, "palette", "course")) +names(gcol) <- levels(ann_col$Group) +ccol <- rgraph_palette(k, "colorblind")[1:k]; names(ccol) <- as.character(1:k) +pal <- grDevices::colorRampPalette(c("blue", "white", "red"))(100) +gaps <- cumsum(as.numeric(table(cl))) +for (fmt in c("png", "pdf")) { + f <- file.path(p$outdir, paste0("kmeans_heatmap.", fmt)) + rgraph_ensure_dir(p$outdir) + pheatmap::pheatmap(zo, scale = "none", cluster_rows = FALSE, cluster_cols = FALSE, + color = pal, border_color = NA, show_rownames = nrow(zo) <= 60, + annotation_row = ann_row, annotation_col = ann_col, + annotation_colors = list(Cluster = ccol, Group = gcol), + gaps_row = gaps, fontsize_col = 9, angle_col = 45, + width = 8, height = max(5, nrow(zo) * 0.15 + 2), filename = f) + rgraph_emit(f) +} + +# ---- 各簇表达趋势线 -------------------------------------------------------- +long <- as.data.frame(z) %>% dplyr::mutate(gene_id = rownames(z), Cluster = factor(km$cluster)) %>% + tidyr::pivot_longer(cols = -c(gene_id, Cluster), names_to = "Sample", values_to = "Expression") +long$Sample <- factor(long$Sample, levels = colnames(z)) +for (cn in 1:k) { + cc <- long %>% dplyr::filter(Cluster == cn) + sm <- cc %>% dplyr::group_by(Sample) %>% + dplyr::summarise(mean = mean(Expression), lo = min(Expression), hi = max(Expression), .groups = "drop") + pl <- ggplot2::ggplot() + + ggplot2::geom_point(data = cc, ggplot2::aes(Sample, Expression), color = "grey70", size = 1.5, alpha = 0.5) + + ggplot2::geom_ribbon(data = sm, ggplot2::aes(Sample, ymin = lo, ymax = hi, group = 1), fill = "red", alpha = 0.2) + + ggplot2::geom_line(data = sm, ggplot2::aes(Sample, mean, group = 1), color = "red", linewidth = 1) + + ggplot2::geom_point(data = sm, ggplot2::aes(Sample, mean), fill = "red", shape = 21, size = 3, color = "white") + + ggplot2::labs(title = paste0("Cluster ", cn, " (n=", sum(km$cluster == cn), ")"), + x = NULL, y = "Z-score") + theme_rgraph() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + rgraph_save(pl, p$outdir, paste0("cluster_", cn, "_trend"), width = 8, height = 5, dpi = rgraph_opt(p, "dpi", 300)) +} +cat("RGRAPH_DONE: kmeans k=", k, "genes=", nrow(z), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/multilevel_scatter.R b/plugins/Presisitence/rgraph/rscripts/multilevel_scatter.R new file mode 100644 index 0000000..6cdb7d2 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/multilevel_scatter.R @@ -0,0 +1,43 @@ +# multilevel_scatter.R —— 多级分组基因表达散点图(逐基因;主分组 x,次分组着色) +# 输入: expr(gene_expression.csv: gene_id+各样本), sample_group(多级分组表) +# params: expr, sample_group, gene_list(CSV) 或 genes(向量), outdir, +# sample_col(默认自动 sample_name/sample_name1), x_col(默认自动 group_name/group_name1), +# color_col(默认 group_name2, 无则同 x_col), dpi(300), width(8),height(8) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2")) + +sg <- read.csv(p$sample_group, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) +pick <- function(cands, override) { + if (!is.null(override) && override %in% names(sg)) return(override) + hit <- intersect(cands, names(sg)); if (length(hit)) hit[1] else stop(paste0("缺少列: ", cands[1])) +} +sample_col <- pick(c("sample_name", "sample_name1"), rgraph_opt(p, "sample_col", NULL)) +x_col <- pick(c("group_name", "group_name1", "group"), rgraph_opt(p, "x_col", NULL)) +color_col <- rgraph_opt(p, "color_col", if ("group_name2" %in% names(sg)) "group_name2" else x_col) + +expr <- rgraph_read_matrix(p$expr) +genes <- if (!is.null(p$genes)) p$genes else read.csv(p$gene_list, check.names = FALSE)$gene_id +cols <- rgraph_palette(length(unique(sg[[color_col]])), rgraph_opt(p, "palette", "course")) + +n_ok <- 0 +for (g in genes) { + row <- expr[expr$gene_id == g, , drop = FALSE] + if (nrow(row) == 0) next + vals <- as.numeric(row[1, intersect(sg[[sample_col]], names(expr))]) + df <- data.frame(sample = intersect(sg[[sample_col]], names(expr)), expression = vals) + df <- merge(df, sg, by.x = "sample", by.y = sample_col) + df$.x <- factor(df[[x_col]], levels = unique(sg[[x_col]])) + df$.col <- factor(df[[color_col]], levels = unique(sg[[color_col]])) + pl <- ggplot2::ggplot(df, ggplot2::aes(x = .x, y = expression, color = .col)) + + ggplot2::geom_point(position = ggplot2::position_dodge(0.6), size = 4, alpha = 0.7) + + ggplot2::scale_color_manual(values = cols, name = color_col) + + ggplot2::labs(title = g, x = NULL, y = "Expression") + theme_rgraph() + rgraph_save(pl, p$outdir, paste0(gsub("[^A-Za-z0-9_.-]", "_", g), "_expression"), + width = rgraph_opt(p, "width", 8), height = rgraph_opt(p, "height", 8), + dpi = rgraph_opt(p, "dpi", 300)) + n_ok <- n_ok + 1 +} +if (n_ok == 0) stop("gene_list 中没有基因命中表达矩阵") +cat("RGRAPH_DONE: multilevel_scatter genes=", n_ok, "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/network.R b/plugins/Presisitence/rgraph/rscripts/network.R new file mode 100644 index 0000000..25aa8f8 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/network.R @@ -0,0 +1,87 @@ +# network.R —— ggNetView 网络图(相关/共表达/宏基因组共现网络,按模块着色) +# 两种输入模式: +# mode="matrix": expr(gene_id/OTU + 样本列) → build_graph_from_mat 构相关网络 +# (method: WGCNA/cor/SPARCC/SpiecEasi/Hmisc);适合共表达/宏基因组共现 +# mode="edge": edges(前两列 from,to [+ weight]) [+ nodes(首列=id,附注释)] +# → build_graph_from_node_edge;可直接吃 rgraph_wgcna 的 Cytoscape 边表、rgraph_ppi 的 Target_PPi +# params(通用): mode, outdir, name, layout("gephi"), layout_module("random"), +# group_by("Modularity"), fill_by("Modularity"), color_by, label(F 或前缀字串), +# pointlabel(如 "top5"), jitter(F), add_outer(F), dpi(300), width(9),height(8) +# params(matrix): expr, sample_group(可选), genes(可选子集), method("WGCNA"), cor_method("pearson"), +# r_threshold(0.7), p_threshold(0.05), transform("none"), node_annotation(CSV,首列=id), max_features(500) +# params(edge): edges, nodes(可选), from_col, to_col, weight_col, directed(F) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("ggNetView", "ggplot2", "dplyr")) + +mode <- rgraph_opt(p, "mode", "matrix") + +if (mode == "edge") { + edges <- read.csv(p$edges, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE, + sep = rgraph_opt(p, "sep", ",")) + if (ncol(edges) < 2) { + # 容错:可能是制表符/空格分隔 + edges <- read.table(p$edges, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + } + fc <- rgraph_opt(p, "from_col", names(edges)[1]) + tc <- rgraph_opt(p, "to_col", names(edges)[2]) + wc <- rgraph_opt(p, "weight_col", { wi <- which(tolower(names(edges)) == "weight"); if (length(wi)) names(edges)[wi[1]] else NA }) + edge_df <- data.frame(from = as.character(edges[[fc]]), to = as.character(edges[[tc]]), + stringsAsFactors = FALSE) + if (!is.na(wc) && wc %in% names(edges)) edge_df$weight <- as.numeric(edges[[wc]]) + if (!is.null(p$nodes)) { + node_df <- read.csv(p$nodes, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + } else { + node_df <- data.frame(name = unique(c(edge_df$from, edge_df$to)), stringsAsFactors = FALSE) + } + obj <- build_graph_from_node_edge(node = node_df, edge = edge_df, + directed = isTRUE(rgraph_opt(p, "directed", FALSE)), + top_modules = rgraph_opt(p, "top_modules", 15)) +} else { + expr <- rgraph_read_matrix(p$expr) + if (!is.null(p$sample_group)) { + sg <- rgraph_read_sample_group(p$sample_group); samples <- sg$sample_name + } else { + samples <- setdiff(names(expr), c("gene_id", "Length")) + } + if (!is.null(p$genes)) expr <- expr[expr$gene_id %in% p$genes, , drop = FALSE] + rownames(expr) <- expr$gene_id + mat <- as.matrix(expr[, samples, drop = FALSE]); mode(mat) <- "numeric" + mat[!is.finite(mat)] <- 0 + # 无基因子集时按方差取 top max_features,控制网络规模/耗时/可读性 + if (is.null(p$genes)) { + mf <- rgraph_opt(p, "max_features", 500) + if (nrow(mat) > mf) mat <- mat[order(apply(mat, 1, var), decreasing = TRUE)[seq_len(mf)], , drop = FALSE] + } + ann <- NULL + if (!is.null(p$node_annotation)) ann <- read.csv(p$node_annotation, check.names = FALSE, stringsAsFactors = FALSE) + obj <- build_graph_from_mat(mat = mat, + transfrom.method = rgraph_opt(p, "transform", "none"), + method = rgraph_opt(p, "method", "WGCNA"), + cor.method = rgraph_opt(p, "cor_method", "pearson"), + proc = rgraph_opt(p, "proc", "BH"), + r.threshold = rgraph_opt(p, "r_threshold", 0.7), + p.threshold = rgraph_opt(p, "p_threshold", 0.05), + node_annotation = ann, + top_modules = rgraph_opt(p, "top_modules", 15)) +} + +pl <- ggNetView(graph_obj = obj, + layout = rgraph_opt(p, "layout", "gephi"), + layout.module = rgraph_opt(p, "layout_module", "random"), + group.by = rgraph_opt(p, "group_by", "Modularity"), + fill.by = rgraph_opt(p, "fill_by", "Modularity"), + color.by = rgraph_opt(p, "color_by", NULL), + pointsize = rgraph_opt(p, "pointsize", c(1, 6)), + jitter = isTRUE(rgraph_opt(p, "jitter", FALSE)), + add_outer = isTRUE(rgraph_opt(p, "add_outer", FALSE)), + label = rgraph_opt(p, "label", FALSE), + pointlabel = rgraph_opt(p, "pointlabel", NULL), + linealpha = rgraph_opt(p, "linealpha", 0.2), + linecolor = rgraph_opt(p, "linecolor", "#d9d9d9")) + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "network"), + width = rgraph_opt(p, "width", 9), height = rgraph_opt(p, "height", 8), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: network mode=", mode, "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/norm_matrix.R b/plugins/Presisitence/rgraph/rscripts/norm_matrix.R new file mode 100644 index 0000000..286197d --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/norm_matrix.R @@ -0,0 +1,36 @@ +# norm_matrix.R —— 表达矩阵标准化(DESeq2 中位数比值 / limma voom logCPM) +# params: count(gene_count.csv), sample_group, outdir, method("deseq2"|"limma"), +# min_mean(limma 低表达过滤, 默认1), name(默认 normalized_count) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr")) + +method <- rgraph_opt(p, "method", "deseq2") +count <- rgraph_read_matrix(p$count) +sg <- rgraph_read_sample_group(p$sample_group) +rgraph_check_samples(count, sg$sample_name) + +mat <- count[, c("gene_id", sg$sample_name)] +rownames(mat) <- mat$gene_id; mat$gene_id <- NULL +mat[is.na(mat)] <- 0 +grp <- factor(sg$group[match(colnames(mat), sg$sample_name)]) + +if (method == "limma") { + rgraph_library(c("limma")) + ns <- length(sg$sample_name) + keep <- rowSums(mat == 0) < ns * 0.5 & rowMeans(mat) >= rgraph_opt(p, "min_mean", 1) + mat <- mat[keep, , drop = FALSE] + design <- model.matrix(~0 + grp) + v <- limma::voom(as.matrix(mat), design) + norm <- 2^v$E +} else { + rgraph_library(c("DESeq2")) + coldata <- data.frame(row.names = colnames(mat), group = grp) + dds <- DESeq2::DESeqDataSetFromMatrix(round(as.matrix(mat)), coldata, design = ~group) + dds <- DESeq2::estimateSizeFactors(dds) + norm <- DESeq2::counts(dds, normalized = TRUE) +} +out <- data.frame(gene_id = rownames(norm), norm, row.names = NULL, check.names = FALSE) +rgraph_write_csv(out, p$outdir, rgraph_opt(p, "name", "normalized_count")) +cat("RGRAPH_DONE: norm_matrix", method, "genes=", nrow(out), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/normalize.R b/plugins/Presisitence/rgraph/rscripts/normalize.R new file mode 100644 index 0000000..c2573d7 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/normalize.R @@ -0,0 +1,37 @@ +# normalize.R —— 由 count 计算 FPKM 或 TPM +# params: count(gene_count.csv 路径), sample_group(路径,可选), outdir, +# method("fpkm"|"tpm"), out_name(默认 gene_fpkm / gene_tpm) +# 输入列: gene_id, Length, <各样本 count 列> +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) # _common.R +p <- rgraph_load_params() # 读取 .args[[1]] 指向的参数文件 +rgraph_library(c("dplyr")) + +method <- rgraph_opt(p, "method", "fpkm") +gene_count <- rgraph_read_matrix(p$count) +if (!"Length" %in% names(gene_count)) stop("gene_count 缺少 Length 列 (基因长度),无法计算 FPKM/TPM") + +# 样本列:优先按 sample_group,否则用除 gene_id/Length 外的所有数值列 +if (!is.null(p$sample_group)) { + sg <- rgraph_read_sample_group(p$sample_group) + sample_list <- sg$sample_name + rgraph_check_samples(gene_count, sample_list) +} else { + sample_list <- setdiff(names(gene_count), c("gene_id", "Length")) +} + +out <- gene_count[, "gene_id", drop = FALSE] +for (s in sample_list) { + x <- gene_count[[s]] + x[is.na(x)] <- 0 + if (identical(method, "tpm")) { + rpk <- x / gene_count$Length + out[[s]] <- rpk / sum(rpk, na.rm = TRUE) * 1e6 + } else { + out[[s]] <- (x * 1e9) / gene_count$Length / sum(x, na.rm = TRUE) + } +} + +out_name <- rgraph_opt(p, "out_name", if (identical(method, "tpm")) "gene_tpm" else "gene_fpkm") +rgraph_write_csv(out, p$outdir, out_name) +cat("RGRAPH_DONE: normalize", method, "genes=", nrow(out), "samples=", length(sample_list), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/pathway_gene_heatmap.R b/plugins/Presisitence/rgraph/rscripts/pathway_gene_heatmap.R new file mode 100644 index 0000000..7db6b35 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/pathway_gene_heatmap.R @@ -0,0 +1,70 @@ +# pathway_gene_heatmap.R —— 通路富集基因的相关性热图 / 聚类热图(逐通路出图) +# 输入: enrich(GO/KEGG 富集表, 需 Description 及基因列 geneID/Gene_symbol), fpkm, sample_group +# params: enrich, fpkm, sample_group, outdir, kind("correlation"|"cluster"), +# id_col(富集表里与 fpkm gene_id 同类型的基因列, 默认 geneID), +# min_count(通路基因数下限,默认3), max_pathways(最多出图通路数,默认6), +# pathway_col(通路命名列, 默认 GOID/KEGGID/Description), dpi(300) +# 注: id_col 的基因需与 fpkm 的 gene_id 同一 ID 类型(本流程 rgraph_enrich→rgraph_diff 天然一致)。 +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2", "reshape2")) + +kind <- rgraph_opt(p, "kind", "correlation") +id_col <- rgraph_opt(p, "id_col", "geneID") +et <- read.csv(p$enrich, stringsAsFactors = FALSE, check.names = FALSE) +fpkm <- rgraph_read_matrix(p$fpkm) +sg <- rgraph_read_sample_group(p$sample_group) +if (!id_col %in% names(et)) stop(paste0("富集表缺少基因列: ", id_col)) + +# 通路命名列 +pcol <- rgraph_opt(p, "pathway_col", + intersect(c("GOID", "KEGGID", "Description"), names(et))[1]) +et <- et[et$Count >= rgraph_opt(p, "min_count", 3), , drop = FALSE] +et <- et[order(-et$Count), , drop = FALSE] +et <- utils::head(et, rgraph_opt(p, "max_pathways", 6)) + +sel_fpkm <- fpkm[, c("gene_id", sg$sample_name), drop = FALSE] +n_ok <- 0 +for (i in seq_len(nrow(et))) { + pid <- gsub("[^A-Za-z0-9_.-]", "_", as.character(et[[pcol]][i])) + genes <- unique(trimws(strsplit(as.character(et[[id_col]][i]), "/")[[1]])) + sub <- sel_fpkm[sel_fpkm$gene_id %in% genes, , drop = FALSE] + if (nrow(sub) < 2) next + rownames(sub) <- sub$gene_id; sub$gene_id <- NULL + m <- as.matrix(sub); m[is.na(m)] <- 0 + logm <- log2(m + 1) + logm <- logm[rowSums(logm) != 0, , drop = FALSE] + if (nrow(logm) < 2) next + + if (kind == "cluster") { + rgraph_library(c("pheatmap")) + pal <- grDevices::colorRampPalette(c("blue", "white", "red"))(100) + for (fmt in c("png", "pdf")) { + f <- file.path(p$outdir, paste0(pid, ".", fmt)) + rgraph_ensure_dir(p$outdir) + pheatmap::pheatmap(logm, scale = "row", cluster_rows = TRUE, cluster_cols = TRUE, + color = pal, border_color = "grey80", show_rownames = nrow(logm) <= 60, + fontsize_col = 10, angle_col = 45, + width = 8, height = max(4, nrow(logm) * 0.2 + 2), filename = f) + rgraph_emit(f) + } + } else { + cm <- cor(t(logm), use = "pairwise.complete.obs") + df <- reshape2::melt(cm); df$size <- abs(df$value) + pl <- ggplot2::ggplot(df, ggplot2::aes(x = Var2, y = Var1, fill = value)) + + ggplot2::geom_tile(color = "white", linewidth = 0.4) + + ggplot2::geom_text(ggplot2::aes(label = round(value, 2)), size = 2.6) + + ggplot2::scale_fill_gradientn(colors = grDevices::colorRampPalette(c("#f3ec98", "#f7440e"))(100), + name = expression(R^2)) + + ggplot2::labs(x = NULL, y = NULL, title = as.character(et[[pcol]][i])) + + theme_rgraph() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + + ggplot2::coord_fixed() + sz <- max(5, nrow(logm) * 0.5) + rgraph_save(pl, p$outdir, pid, width = sz, height = sz, dpi = rgraph_opt(p, "dpi", 300)) + } + n_ok <- n_ok + 1 +} +if (n_ok == 0) stop("无通路满足条件(可能 id_col 基因与 fpkm 的 gene_id 类型不一致)") +cat("RGRAPH_DONE: pathway_gene_heatmap", kind, "pathways=", n_ok, "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/pathway_trend.R b/plugins/Presisitence/rgraph/rscripts/pathway_trend.R new file mode 100644 index 0000000..b5cd7fe --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/pathway_trend.R @@ -0,0 +1,56 @@ +# pathway_trend.R —— 通路上下调基因统计图(发散条形/堆叠条形) +# 输入: enrich(富集表, 需 Description 与 geneID["a/b/c"]), deg(Dse2_result, 需 gene_id,log2FoldChange) +# params: enrich, deg, outdir, kind("diverge"|"stack"), top_n(默认全部), dpi(300), width(10),height(8), name +# 说明: 默认 enrich 的 geneID 与 deg 的 gene_id 为同一 ID 类型(直接匹配)。 +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2")) + +kind <- rgraph_opt(p, "kind", "diverge") +et <- read.csv(p$enrich, stringsAsFactors = FALSE, check.names = FALSE) +deg <- read.csv(p$deg, stringsAsFactors = FALSE, check.names = FALSE) +if (!"geneID" %in% names(et)) stop("富集表缺少 geneID 列") + +pg <- et %>% dplyr::rowwise() %>% + dplyr::mutate(gene_id = strsplit(as.character(geneID), "/")) %>% + tidyr::unnest(gene_id) %>% dplyr::ungroup() %>% + dplyr::select(Description, gene_id) %>% dplyr::distinct() +pg <- merge(pg, deg[, c("gene_id", "log2FoldChange")], by = "gene_id") +pg$Trend <- ifelse(pg$log2FoldChange > 0, "Up", "Down") + +stat <- pg %>% dplyr::group_by(Description, Trend) %>% + dplyr::summarise(Count = dplyr::n(), .groups = "drop") %>% + tidyr::pivot_wider(names_from = Trend, values_from = Count, values_fill = list(Count = 0)) +if (!"Up" %in% names(stat)) stat$Up <- 0 +if (!"Down" %in% names(stat)) stat$Down <- 0 + +# top_n 通路(按总数) +if (!is.null(p$top_n)) { + stat <- stat %>% dplyr::mutate(tot = Up + Down) %>% dplyr::arrange(dplyr::desc(tot)) %>% + dplyr::slice(seq_len(min(p$top_n, dplyr::n()))) %>% dplyr::select(-tot) +} + +if (kind == "stack") { + long <- stat %>% tidyr::pivot_longer(cols = c(Down, Up), names_to = "Trend", values_to = "Count") + pl <- ggplot2::ggplot(long, ggplot2::aes(x = Description, y = Count, fill = Trend)) + + ggplot2::geom_col(position = "stack", width = 0.7) + + ggplot2::geom_text(ggplot2::aes(label = Count), position = ggplot2::position_stack(vjust = 0.5), size = 3) + + ggplot2::scale_fill_manual(values = c(Down = "#157d9f", Up = "#f4736c")) + + ggplot2::labs(x = NULL, y = "Number of DEGs") + theme_rgraph() + ggplot2::coord_flip() +} else { + stat$Down <- -stat$Down + long <- stat %>% tidyr::pivot_longer(cols = c(Down, Up), names_to = "Trend", values_to = "Count") + xmax <- max(abs(long$Count)) + pl <- ggplot2::ggplot(long, ggplot2::aes(x = Count, y = Description, fill = Trend)) + + ggplot2::geom_col() + + ggplot2::geom_vline(xintercept = 0, lty = 4, lwd = 0.6, alpha = 0.8) + + ggplot2::scale_fill_manual(values = c(Down = "#157d9f", Up = "#f4736c")) + + ggplot2::scale_x_continuous(labels = function(x) abs(x), limits = c(-xmax, xmax)) + + ggplot2::labs(x = "Number of DEGs", y = NULL) + theme_rgraph() +} + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "pathway_trend"), + width = rgraph_opt(p, "width", 10), height = rgraph_opt(p, "height", 8), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: pathway_trend", kind, "pathways=", nrow(stat), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/pca.R b/plugins/Presisitence/rgraph/rscripts/pca.R new file mode 100644 index 0000000..076b974 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/pca.R @@ -0,0 +1,56 @@ +# pca.R —— 样本 PCA (log2(fpkm+1) + prcomp),支持置信椭圆与 ggrepel 标签 +# params: fpkm, sample_group, outdir, ellipse(bool,默认F), label(bool,默认T), +# palette("course"|"colorblind"), dpi(300), width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "ggplot2")) +use_repel <- requireNamespace("ggrepel", quietly = TRUE) + +fpkm <- rgraph_read_matrix(p$fpkm) +sg <- rgraph_read_sample_group(p$sample_group) +rgraph_check_samples(fpkm, sg$sample_name) + +mat <- fpkm[, sg$sample_name, drop = FALSE] +mat[is.na(mat)] <- 0 +mat <- mat[rowSums(mat != 0) > 0, , drop = FALSE] # 去全 0 行 +mat <- mat[apply(mat, 1, function(x) length(unique(x)) > 1), , drop = FALSE] # 去常数行 +logm <- log2(mat + 1) + +pca <- prcomp(t(logm), scale. = TRUE) +varexp <- (pca$sdev^2) / sum(pca$sdev^2) +scores <- as.data.frame(pca$x) +scores$Sample <- rownames(scores) +scores$Group <- sg$group_name[match(scores$Sample, sg$sample_name)] +scores$Group <- factor(scores$Group, levels = unique(sg$group_name)) + +cols <- rgraph_palette(length(levels(scores$Group)), rgraph_opt(p, "palette", "course")) + +pl <- ggplot2::ggplot(scores, ggplot2::aes(x = PC1, y = PC2, color = Group)) + + ggplot2::geom_point(size = 4, alpha = 0.85) + + ggplot2::scale_color_manual(values = cols) + + ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "#bebebe") + + ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "#bebebe") + + ggplot2::labs(title = "Principal Component Analysis (PCA)", + x = sprintf("PC1 (%.2f%%)", varexp[1] * 100), + y = sprintf("PC2 (%.2f%%)", varexp[2] * 100)) + + theme_rgraph() + +if (isTRUE(rgraph_opt(p, "ellipse", FALSE)) && nlevels(scores$Group) > 1) + pl <- pl + ggplot2::stat_ellipse(level = 0.95, linewidth = 0.3) + +if (isTRUE(rgraph_opt(p, "label", TRUE))) { + if (use_repel) { + pl <- pl + ggrepel::geom_text_repel(ggplot2::aes(label = Sample), size = 3, + color = "black", max.overlaps = Inf) + } else { + pl <- pl + ggplot2::geom_text(ggplot2::aes(label = Sample), vjust = -0.6, size = 2.5, + color = "black") + } +} + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "pca"), + width = rgraph_opt(p, "width", 6), height = rgraph_opt(p, "height", 5.5), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: pca genes_used=", nrow(logm), "PC1=", round(varexp[1] * 100, 2), + "PC2=", round(varexp[2] * 100, 2), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/ppi.R b/plugins/Presisitence/rgraph/rscripts/ppi.R new file mode 100644 index 0000000..d59e109 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/ppi.R @@ -0,0 +1,50 @@ +# ppi.R —— 基于 STRING 的差异基因蛋白互作网络边表 +# params: deg(Deg_all.csv, 列 gene_id), info(9606.protein.info.txt), +# links(9606.protein.links.txt), outdir, score(400), +# id_is_symbol(bool,默认F), orgdb(id 非 symbol 时用于 ENSEMBL->SYMBOL) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr")) + +deg <- read.csv(p$deg, stringsAsFactors = FALSE) +info <- read.csv(p$info, sep = "\t", stringsAsFactors = FALSE) +links <- read.csv(p$links, sep = " ", stringsAsFactors = FALSE) + +if (isTRUE(rgraph_opt(p, "id_is_symbol", FALSE))) { + genes <- data.frame(gene_id = unique(deg$gene_id)) +} else { + rgraph_library(c("clusterProfiler", p$orgdb)) + genes <- suppressWarnings(clusterProfiler::bitr(deg$gene_id, fromType = "ENSEMBL", + toType = "SYMBOL", OrgDb = p$orgdb)) + genes <- data.frame(gene_id = unique(genes$SYMBOL)) +} + +# symbol -> string_protein_id +genes <- merge(genes, info[, c("string_protein_id", "preferred_name")], + by.x = "gene_id", by.y = "preferred_name") + +dl <- merge(links, genes, by.x = "protein1", by.y = "string_protein_id") +names(dl)[names(dl) == "gene_id"] <- "from_symbol" +dl <- merge(dl, genes, by.x = "protein2", by.y = "string_protein_id") +names(dl)[names(dl) == "gene_id"] <- "to_symbol" +names(dl)[names(dl) == "protein1"] <- "from_protein" +names(dl)[names(dl) == "protein2"] <- "to_protein" + +dl <- dl %>% + dplyr::select(from_protein, to_protein, from_symbol, to_symbol, combined_score) %>% + dplyr::rowwise() %>% + dplyr::mutate(a = pmin(from_symbol, to_symbol), b = pmax(from_symbol, to_symbol)) %>% + dplyr::ungroup() %>% + dplyr::distinct(a, b, .keep_all = TRUE) %>% + dplyr::select(-a, -b) %>% + dplyr::filter(combined_score >= rgraph_opt(p, "score", 400)) + +rgraph_ensure_dir(p$outdir) +f1 <- file.path(p$outdir, "Target_PPi.tsv") +write.table(dl, f1, row.names = FALSE, sep = "\t"); rgraph_emit(f1) +nodes <- length(unique(c(dl$from_symbol, dl$to_symbol))) +f2 <- file.path(p$outdir, "Edge_Node_count.tsv") +write.table(data.frame(Node_count = nodes, Edge_count = nrow(dl)), f2, row.names = FALSE, sep = "\t") +rgraph_emit(f2) +cat("RGRAPH_DONE: ppi nodes=", nodes, "edges=", nrow(dl), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/quadrant.R b/plugins/Presisitence/rgraph/rscripts/quadrant.R new file mode 100644 index 0000000..910c23e --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/quadrant.R @@ -0,0 +1,66 @@ +# quadrant.R —— 九象限图 / 二象限图(两组差异比较的 log2FC 散点) +# params: x(比较1的DEG CSV,列 gene_id,log2FoldChange,pvalue), y(比较2的DEG CSV), +# x_name, y_name, mode("nine"|"four"), pcut(0.05), log2fc(1), outdir, dpi,width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "ggplot2")) + +mode <- rgraph_opt(p, "mode", "nine") +pcut <- rgraph_opt(p, "pcut", 0.05) +log2fc <- rgraph_opt(p, "log2fc", 1) +xn <- rgraph_opt(p, "x_name", "compare1") +yn <- rgraph_opt(p, "y_name", "compare2") + +read_cmp <- function(path, nm) { + d <- read.csv(path, check.names = FALSE) + d <- d[!is.na(d$pvalue) & d$pvalue < pcut, ] + if (mode == "four") d <- d[abs(d$log2FoldChange) > log2fc, ] + d <- d[, c("gene_id", "log2FoldChange")] + names(d)[2] <- nm + d +} +dx <- read_cmp(p$x, xn); dy <- read_cmp(p$y, yn) +m <- merge(dx, dy, by = "gene_id") +X <- m[[xn]]; Y <- m[[yn]] + +if (mode == "four") { + m <- m[(X > 0 & Y > 0) | (X < 0 & Y < 0), ] + X <- m[[xn]]; Y <- m[[yn]] + m$cls <- ifelse(abs(X) >= abs(Y), "High", "Low") + cnt <- table(m$cls) + cmap <- c(High = "#FF3333", Low = "#00CCCC") + labs2 <- c(High = sprintf("High (%d)", cnt["High"]), Low = sprintf("Low (%d)", cnt["Low"])) + brk <- c("High", "Low") +} else { + m$cls <- ifelse(X >= log2fc & Y >= log2fc, "Q1", + ifelse(abs(X) < log2fc & Y >= log2fc, "Q2", + ifelse(X <= -log2fc & Y >= log2fc, "Q3", + ifelse(X >= log2fc & abs(Y) < log2fc, "Q4", + ifelse(abs(X) < log2fc & abs(Y) < log2fc, "Q5", + ifelse(X <= -log2fc & abs(Y) < log2fc, "Q6", + ifelse(X >= log2fc & Y <= -log2fc, "Q7", + ifelse(abs(X) < log2fc & Y <= -log2fc, "Q8", + ifelse(X <= -log2fc & Y <= -log2fc, "Q9", "No"))))))))) + cnt <- table(factor(m$cls, levels = paste0("Q", 1:9))) + cmap <- setNames(c("#ede746", "#FF9933", "#FF3333", "#FFCCFF", "grey80", + "#FF66FF", "#CCFFFF", "#00FF99", "#00CCCC"), paste0("Q", 1:9)) + labs2 <- setNames(sprintf("quadrant.%d (%d)", 1:9, cnt), paste0("Q", 1:9)) + brk <- paste0("Q", 1:9) +} +rgraph_write_csv(m, p$outdir, if (mode == "four") "Four_quadrant_table" else "Nine_quadrant_table") + +pl <- ggplot2::ggplot(m, ggplot2::aes(x = .data[[xn]], y = .data[[yn]], color = cls)) + + ggplot2::geom_point(alpha = 0.8, size = 0.7) + + ggplot2::scale_color_manual(values = cmap, breaks = brk, labels = labs2) + + ggplot2::geom_hline(yintercept = c(-log2fc, log2fc), linetype = "dashed", color = "grey") + + ggplot2::geom_vline(xintercept = c(-log2fc, log2fc), linetype = "dashed", color = "grey") + + ggplot2::labs(x = sprintf("log2FoldChange(%s)", xn), y = sprintf("log2FoldChange(%s)", yn), + color = sprintf("pvalue<%s", pcut)) + + theme_rgraph() + + ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = 3))) + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", if (mode == "four") "Four_quadrant" else "Nine_quadrant"), + width = rgraph_opt(p, "width", 7.3), height = rgraph_opt(p, "height", 6), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: quadrant", mode, "genes=", nrow(m), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/sankey.R b/plugins/Presisitence/rgraph/rscripts/sankey.R new file mode 100644 index 0000000..6cfefdd --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/sankey.R @@ -0,0 +1,43 @@ +# sankey.R —— 富集通路-基因 桑基图(用 ggalluvial 重写;基因→通路流向) +# 输入: enrich(富集表, 需 Description/通路名 与 geneID["a/b/c"]) +# params: enrich, outdir, id_col(基因列,默认geneID), pathway_col(通路名列,默认Description), +# top_n(通路数,默认8), dpi(300), width(12),height(8), name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("dplyr", "tidyr", "ggplot2", "ggalluvial")) + +id_col <- rgraph_opt(p, "id_col", "geneID") +pcol <- rgraph_opt(p, "pathway_col", "Description") +et <- read.csv(p$enrich, stringsAsFactors = FALSE, check.names = FALSE) +if (!id_col %in% names(et)) stop(paste0("富集表缺少基因列: ", id_col)) +if ("Count" %in% names(et)) et <- et[order(-et$Count), , drop = FALSE] +et <- utils::head(et, rgraph_opt(p, "top_n", 8)) + +pairs <- et %>% dplyr::rowwise() %>% + dplyr::mutate(Gene = strsplit(as.character(.data[[id_col]]), "/")) %>% + tidyr::unnest(Gene) %>% dplyr::ungroup() %>% + dplyr::transmute(Gene = trimws(Gene), Pathway = .data[[pcol]]) %>% + dplyr::distinct() +pairs$freq <- 1 +pairs$id <- seq_len(nrow(pairs)) + +long <- ggalluvial::to_lodes_form(pairs, key = "axis", value = "stratum", + axes = c("Gene", "Pathway"), id = "id") +ncol_need <- length(unique(long$stratum)) +fill <- rgraph_palette(ncol_need, rgraph_opt(p, "palette", "course")) + +pl <- ggplot2::ggplot(long, ggplot2::aes(x = axis, stratum = stratum, alluvium = id, label = stratum)) + + ggalluvial::geom_flow(ggplot2::aes(fill = stratum), width = 0.3, alpha = 0.5) + + ggalluvial::geom_stratum(ggplot2::aes(fill = stratum), width = 0.3, color = "grey40") + + ggplot2::geom_text(stat = ggalluvial::StatStratum, size = 2.6) + + ggplot2::scale_fill_manual(values = fill, guide = "none") + + ggplot2::scale_x_discrete(limits = c("Gene", "Pathway"), expand = c(0.1, 0.1)) + + ggplot2::labs(x = NULL, y = NULL) + + ggplot2::theme_void() + + ggplot2::theme(axis.text.x = ggplot2::element_text(face = "bold", size = 12)) + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "sankey"), + width = rgraph_opt(p, "width", 12), height = rgraph_opt(p, "height", 8), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: sankey pathways=", nrow(et), "edges=", nrow(pairs), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/segmented_volcano.R b/plugins/Presisitence/rgraph/rscripts/segmented_volcano.R new file mode 100644 index 0000000..8d40445 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/segmented_volcano.R @@ -0,0 +1,46 @@ +# segmented_volcano.R —— 分段式火山图(y 轴断裂,适合个别通路 p 值极端时) +# params: result(Dse2_result.csv), outdir, sig_metric("pvalue"|"padj"), pcut(0.05), log2fc(1), +# color_scheme("rb"|"rg"), break_lower, break_upper(y 轴断裂区间), dpi,width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("ggplot2", "dplyr", "ggbreak")) + +sig <- rgraph_opt(p, "sig_metric", "pvalue") +pcut <- rgraph_opt(p, "pcut", 0.05) +log2fc <- rgraph_opt(p, "log2fc", 1) + +d <- read.csv(p$result, header = TRUE, check.names = FALSE) +d <- d[!is.na(d[[sig]]) & !is.na(d$log2FoldChange), ] +m <- d[[sig]] +d$color <- ifelse(m < pcut & d$log2FoldChange > log2fc, "Up", + ifelse(m < pcut & d$log2FoldChange < -log2fc, "Down", "No")) +cnt <- table(factor(d$color, levels = c("Up", "Down", "No"))) +labs3 <- c(Up = sprintf("Up (%d)", cnt["Up"]), Down = sprintf("Down (%d)", cnt["Down"]), + No = sprintf("No (%d)", cnt["No"])) +cols <- rgraph_updown_colors(rgraph_opt(p, "color_scheme", "rb")) + +pl <- ggplot2::ggplot(d, ggplot2::aes(x = log2FoldChange, y = -log10(.data[[sig]]), color = color)) + + ggplot2::geom_point(alpha = 0.5, size = 1) + + ggplot2::scale_color_manual(values = cols, breaks = c("Up", "Down", "No"), labels = labs3) + + ggplot2::geom_hline(yintercept = -log10(pcut), linetype = "dashed", color = "grey") + + ggplot2::geom_vline(xintercept = c(-log2fc, log2fc), linetype = "dashed", color = "grey") + + ggplot2::labs(x = "log2 (Fold Change)", y = sprintf("-log10(%s)", sig), + color = sprintf("%s<%s\n|log2FC|>%s", sig, pcut, log2fc)) + + theme_rgraph() + + ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = 3))) + +# y 轴断裂:默认自动在数据主体与极端值之间断开 +ymax <- max(-log10(d[[sig]]), na.rm = TRUE) +bl <- rgraph_opt(p, "break_lower", NULL); bu <- rgraph_opt(p, "break_upper", NULL) +if (is.null(bl) || is.null(bu)) { + q <- stats::quantile(-log10(d[[sig]]), 0.995, na.rm = TRUE) + if (ymax > q * 1.5) { bl <- ceiling(q); bu <- floor(ymax * 0.9) } +} +if (!is.null(bl) && !is.null(bu) && bu > bl) + pl <- pl + ggbreak::scale_y_break(c(bl, bu), space = 0.1, scales = c(0.5, 1)) + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "volcano_segmented"), + width = rgraph_opt(p, "width", 8), height = rgraph_opt(p, "height", 6.2), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: segmented_volcano up=", cnt["Up"], "down=", cnt["Down"], "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/ssgsea.R b/plugins/Presisitence/rgraph/rscripts/ssgsea.R new file mode 100644 index 0000000..a9f7922 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/ssgsea.R @@ -0,0 +1,40 @@ +# ssgsea.R —— ssGSEA 单样本基因集富集打分 + 各基因集分组箱线图 +# params: expr(gene_expression.csv: gene_id+各样本), gmt(geneset.gmt), sample_group, outdir, +# normalize(bool,默认T), dpi(300) +# 需 GSVA + GSEABase(本机默认缺,按提示 BiocManager::install(c("GSVA","GSEABase"))) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("GSVA", "GSEABase", "dplyr", "ggplot2")) + +sg <- rgraph_read_sample_group(p$sample_group) +expr <- rgraph_read_matrix(p$expr) +expr <- expr[, c("gene_id", sg$sample_name), drop = FALSE] +ns <- length(sg$sample_name) +keep <- rowSums(expr[, -1, drop = FALSE] > 0) / ns >= 0.5 & + rowMeans(expr[, -1, drop = FALSE]) >= 1 +expr <- expr[keep, , drop = FALSE] +rownames(expr) <- expr$gene_id; expr$gene_id <- NULL +mat <- as.matrix(expr) + +gs <- GSEABase::getGmt(p$gmt, geneIdType = GSEABase::SymbolIdentifier()) +param <- GSVA::ssgseaParam(exprData = mat, geneSets = gs, + normalize = isTRUE(rgraph_opt(p, "normalize", TRUE))) +score <- GSVA::gsva(param, verbose = FALSE) +rgraph_write_csv(data.frame(geneset = rownames(score), score, check.names = FALSE), + p$outdir, "ssGSEA_score", row.names = FALSE) + +cols <- rgraph_palette(length(unique(sg$group_name)), rgraph_opt(p, "palette", "course")) +for (g in rownames(score)) { + df <- data.frame(Score = as.numeric(score[g, ]), sample_name = colnames(score)) + df$group_name <- factor(sg$group_name[match(df$sample_name, sg$sample_name)], + levels = unique(sg$group_name)) + pl <- ggplot2::ggplot(df, ggplot2::aes(x = group_name, y = Score, fill = group_name, colour = group_name)) + + ggplot2::geom_boxplot(alpha = 0.5, linewidth = 0.6, width = 0.7, outlier.alpha = 0) + + ggplot2::scale_fill_manual(values = cols) + ggplot2::scale_color_manual(values = cols) + + ggplot2::labs(title = g, x = NULL, y = "ssGSEA Score") + theme_rgraph() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + rgraph_save(pl, p$outdir, paste0(gsub("[^A-Za-z0-9_.-]", "_", g), "_box"), + width = 6, height = 6, dpi = rgraph_opt(p, "dpi", 300)) +} +cat("RGRAPH_DONE: ssgsea genesets=", nrow(score), "samples=", ncol(score), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/stacked_volcano.R b/plugins/Presisitence/rgraph/rscripts/stacked_volcano.R new file mode 100644 index 0000000..0521010 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/stacked_volcano.R @@ -0,0 +1,62 @@ +# stacked_volcano.R —— 堆叠火山图(多比较组 DEG 的 log2FC 抖动散点) +# 输入: DEG_data.csv (列 gene_id, log2FoldChange, pvalue[, TvsC 分组]) +# params: deg_data, outdir, sig_metric("pvalue"|"padj"), pcut(0.05), log2fc(1), +# group_col("TvsC"), top_n(每组每方向标注数,5), dpi(300), width(10),height(9), name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("ggplot2", "dplyr")) +use_repel <- requireNamespace("ggrepel", quietly = TRUE) + +sig <- rgraph_opt(p, "sig_metric", "pvalue") +pcut <- rgraph_opt(p, "pcut", 0.05) +log2fc <- rgraph_opt(p, "log2fc", 1) +gcol <- rgraph_opt(p, "group_col", "TvsC") +top_n <- rgraph_opt(p, "top_n", 5) + +d <- read.csv(p$deg_data, check.names = FALSE) +if (!gcol %in% names(d)) d[[gcol]] <- "DEG" # 无分组列时单组 +d$regulate <- ifelse(d[[sig]] < pcut & d$log2FoldChange > log2fc, "Up", + ifelse(d[[sig]] < pcut & d$log2FoldChange < -log2fc, "Down", "No")) +d <- d[d$regulate != "No", ] +d[[gcol]] <- factor(d[[gcol]], levels = unique(d[[gcol]])) +set.seed(123) +d$xj <- as.numeric(d[[gcol]]) + runif(nrow(d), -0.18, 0.18) + +# 每组背景柱(min/max*1.2) +bar <- d %>% dplyr::group_by(.g = .data[[gcol]]) %>% + dplyr::summarise(minv = min(log2FoldChange) * 1.2, maxv = max(log2FoldChange) * 1.2, + center = mean(as.numeric(.data[[gcol]])), .groups = "drop") +bar$xc <- as.numeric(factor(bar$.g, levels = levels(d[[gcol]]))) +bcols <- rgraph_palette(nrow(bar), rgraph_opt(p, "palette", "course")) + +# 标注 top_n +lab <- d %>% dplyr::group_by(.data[[gcol]]) %>% + dplyr::slice_max(log2FoldChange, n = top_n) %>% dplyr::ungroup() +lab <- rbind(lab, d %>% dplyr::group_by(.data[[gcol]]) %>% + dplyr::slice_min(log2FoldChange, n = top_n) %>% dplyr::ungroup()) +lab <- lab[!duplicated(lab$gene_id), ] + +pl <- ggplot2::ggplot(d, ggplot2::aes(x = xj, y = log2FoldChange)) + + ggplot2::geom_col(data = bar, ggplot2::aes(x = xc, y = minv), fill = rep(bcols, 1), alpha = 0.15, width = 0.8) + + ggplot2::geom_col(data = bar, ggplot2::aes(x = xc, y = maxv), fill = rep(bcols, 1), alpha = 0.15, width = 0.8) + + ggplot2::geom_point(ggplot2::aes(size = abs(log2FoldChange), alpha = abs(log2FoldChange), color = regulate)) + + ggplot2::scale_color_manual(values = c(Down = "#3288bd", Up = "#d92523")) + + ggplot2::scale_size_continuous(range = c(0.5, 4)) + + ggplot2::scale_alpha_continuous(range = c(0.2, 0.85)) + + ggplot2::geom_tile(data = bar, ggplot2::aes(x = xc, y = 0), height = 0.8, width = 0.9, + fill = bcols, alpha = 0.85) + + ggplot2::geom_text(data = bar, ggplot2::aes(x = xc, y = 0, label = .g), size = 5) + + ggplot2::scale_x_continuous(breaks = bar$xc, labels = bar$.g) + + ggplot2::labs(x = NULL, y = "log2FoldChange", color = "Trend", size = "|log2FoldChange|") + + theme_rgraph() + + ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = 3)), alpha = "none") + +if (use_repel && nrow(lab) > 0) + pl <- pl + ggrepel::geom_text_repel(data = lab, ggplot2::aes(x = xj, y = log2FoldChange, label = gene_id), + size = 3.5, max.overlaps = Inf, segment.size = 0.3, segment.color = "gray40") + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "multi_volcano"), + width = rgraph_opt(p, "width", 10), height = rgraph_opt(p, "height", 9), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: stacked_volcano groups=", nrow(bar), "DEG=", nrow(d), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/venn.R b/plugins/Presisitence/rgraph/rscripts/venn.R new file mode 100644 index 0000000..33d70a3 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/venn.R @@ -0,0 +1,49 @@ +# venn.R —— 韦恩图(2-5 组)+ 交集成员表 +# params: sets(CSV 路径向量) 或 folder(装 CSV 的目录), col(元素列名,默认自动: Element_name/gene_id/首列), +# names(可选,各集合名), outdir, name("venn"), dpi(300), width(7),height(7) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("VennDiagram")) +suppressWarnings(try(futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger"), + silent = TRUE)) + +pick_col <- function(df, col) { + if (!is.null(col) && col %in% names(df)) return(col) + for (c in c("Element_name", "gene_id", "Gene_symbol")) if (c %in% names(df)) return(c) + names(df)[1] +} + +# 收集集合 +if (!is.null(p$folder)) { + files <- list.files(p$folder, pattern = "\\.csv$", full.names = TRUE) + set_names <- tools::file_path_sans_ext(basename(files)) +} else { + files <- p$sets + set_names <- if (!is.null(p$names)) p$names else tools::file_path_sans_ext(basename(files)) +} +if (length(files) < 2 || length(files) > 5) stop("韦恩图支持 2-5 个集合") + +elem <- list() +for (i in seq_along(files)) { + df <- read.csv(files[i], header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + v <- unique(df[[pick_col(df, p$col)]]) + elem[[set_names[i]]] <- v[!is.na(v) & v != ""] +} + +fill <- rgraph_palette(length(elem), rgraph_opt(p, "palette", "course")) +grob <- VennDiagram::venn.diagram(elem, category.names = names(elem), filename = NULL, + fill = fill, col = NA, cat.col = "black", cat.cex = 0.9, + cex = 1.1, margin = 0.08) + +rgraph_save_base(p$outdir, rgraph_opt(p, "name", "venn"), + function() { grid::grid.newpage(); grid::grid.draw(grob) }, + width = rgraph_opt(p, "width", 7), height = rgraph_opt(p, "height", 7), + dpi = rgraph_opt(p, "dpi", 300)) + +# 成员表 +all_g <- unique(unlist(elem)) +tab <- data.frame(Element = all_g) +for (nm in names(elem)) tab[[nm]] <- all_g %in% elem[[nm]] +rgraph_write_csv(tab, p$outdir, "veen_table") +cat("RGRAPH_DONE: venn sets=", length(elem), "elements=", length(all_g), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/virtual_knockout.R b/plugins/Presisitence/rgraph/rscripts/virtual_knockout.R new file mode 100644 index 0000000..b153da2 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/virtual_knockout.R @@ -0,0 +1,69 @@ +# virtual_knockout.R —— bulk RNA 虚拟敲除(GENIE3 调控网络 + bulkKnk 扰动 + 双轴火山图) +# 输入: count(gene_count.csv), sample_group(需 group_name, pheno_data) +# params: count, sample_group, outdir, ko_gene(敲除基因), ko_group(敲除样本所在分组), +# logfc(0.5), hv_num(高变基因数,5000), dpi(300) +# 依赖: bulkKnk(课程作者本地源码包, install.packages("bulkKnk-main", repos=NULL, type="source")) +# + GENIE3(Bioconductor)。二者缺失将返回缺包提示。 +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("bulkKnk", "GENIE3", "ggplot2", "dplyr")) +use_repel <- requireNamespace("ggrepel", quietly = TRUE) + +ko_gene <- p$ko_gene; ko_group <- p$ko_group +logfc <- rgraph_opt(p, "logfc", 0.5) +hv_num <- rgraph_opt(p, "hv_num", 5000) + +sg <- read.csv(p$sample_group, check.names = FALSE, stringsAsFactors = FALSE) +expr <- read.csv(p$count, check.names = FALSE, stringsAsFactors = FALSE) +if (length(unique(expr$gene_id)) < nrow(expr)) stop("gene_id 存在重复,请先去重") +rownames(expr) <- expr$gene_id +mat <- t(expr[, sg$sample_name]) + +# 高变基因 + 确保 ko_gene 在内 +sub_expr <- bulkKnk::select_advanced_features(mat, n_features = hv_num) +genes <- unique(colnames(sub_expr)) +if (!(ko_gene %in% genes)) sub_expr <- mat[, c(genes, ko_gene)] +sub_expr <- as.matrix(sub_expr); mode(sub_expr) <- "numeric" + +# WGCNA 找相关模块 → 调控网络(GENIE3)→ DPI 剪枝 +wg <- bulkKnk::identify_hub_modules(sub_expr, trait_vec = sg$pheno_data) +best <- unique(wg$best_genes); if (!(ko_gene %in% best)) best <- c(best, ko_gene) +focus <- as.matrix(sub_expr[, best]) +netc <- bulkKnk::prune_network_dpi_fast(bulkKnk::infer_causal_network(focus)$clean_matrix) +W <- sweep(t(netc), 2, colSums(t(netc)) + 1e-9, FUN = "/") + +pheno <- setNames(sg$group_name, sg$sample_name) +common <- intersect(rownames(focus), names(pheno)) +ko <- bulkKnk::batch_virtual_knockout(expr_matrix = focus[common, ], adj_matrix = W, + target_genes = ko_gene, target_pheno = ko_group, pheno_vec = pheno[common]) + +# 汇总 +gs <- data.frame(gene_id = colnames(ko$WT), + WT_Mean = colMeans(ko$WT, na.rm = TRUE), KO_Mean = colMeans(ko$KO, na.rm = TRUE)) +gs$Diff <- gs$KO_Mean - gs$WT_Mean +gs$log2FoldChange <- log2((gs$KO_Mean + 1e-9) / (gs$WT_Mean + 1e-9)) +rgraph_write_csv(gs, p$outdir, "1.virtual_knockout_summary") +sig <- gs[abs(gs$log2FoldChange) > logfc & gs$gene_id != ko_gene, ] +rgraph_write_csv(sig, p$outdir, "2.significant_perturbed_genes") + +# 双轴火山图(x=虚拟 log2FC, y=WT 基线表达) +gs$color <- ifelse(gs$log2FoldChange > logfc, "Up regulated", + ifelse(gs$log2FoldChange < -logfc, "Down regulated", "Stable")) +gs$color[gs$gene_id == ko_gene] <- "KO gene" +cmap <- c("KO gene" = "#2ca02c", "Up regulated" = "#c0392b", "Down regulated" = "#2a80b9", "Stable" = "grey80") +pl <- ggplot2::ggplot(gs, ggplot2::aes(x = log2FoldChange, y = WT_Mean, color = color)) + + ggplot2::geom_point(alpha = 0.7, size = 1) + + ggplot2::scale_color_manual(values = cmap, breaks = names(cmap)) + + ggplot2::geom_vline(xintercept = c(-logfc, logfc), linetype = "dashed", color = "grey") + + ggplot2::labs(title = "Virtual Knockout Effect", x = "Virtual Log2FC (KO vs WT)", + y = "Baseline Expression (WT)", color = NULL) + + theme_rgraph() + + ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = 3))) +lab <- gs[gs$gene_id == ko_gene, ] +if (use_repel && nrow(lab) > 0) + pl <- pl + ggrepel::geom_label_repel(data = lab, ggplot2::aes(label = gene_id), color = "black", + size = 3, min.segment.length = 0, max.overlaps = Inf) +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "dual_volcano"), + width = 7, height = 5, dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: virtual_knockout ko=", ko_gene, "perturbed=", nrow(sig), "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/volcano.R b/plugins/Presisitence/rgraph/rscripts/volcano.R new file mode 100644 index 0000000..352d878 --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/volcano.R @@ -0,0 +1,61 @@ +# volcano.R —— 火山图(可标注 top-N 或自定义基因;色盲友好红蓝配色) +# params: result(Dse2_result.csv 路径), outdir, sig_metric("pvalue"|"padj",默认padj), +# pcut(0.05), log2fc(1), color_scheme("rb"|"rg"), label(bool), +# label_n(默认0), genes(可选,自定义标注基因向量), gene_col(默认gene_id), +# dpi(300), width,height, name +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("ggplot2", "dplyr")) +use_repel <- requireNamespace("ggrepel", quietly = TRUE) + +sig_metric <- rgraph_opt(p, "sig_metric", "padj") +pcut <- rgraph_opt(p, "pcut", 0.05) +log2fc <- rgraph_opt(p, "log2fc", 1) +gene_col <- rgraph_opt(p, "gene_col", "gene_id") + +d <- read.csv(p$result, header = TRUE, check.names = FALSE) +if (!sig_metric %in% names(d)) stop(paste0("结果表缺少列: ", sig_metric)) +d <- d[!is.na(d[[sig_metric]]) & !is.na(d$log2FoldChange), ] +m <- d[[sig_metric]] +d$color <- ifelse(m < pcut & d$log2FoldChange > log2fc, "Up", + ifelse(m < pcut & d$log2FoldChange < -log2fc, "Down", "No")) + +cnt <- table(factor(d$color, levels = c("Up", "Down", "No"))) +labs3 <- c(Up = sprintf("Up (%d)", cnt["Up"]), Down = sprintf("Down (%d)", cnt["Down"]), + No = sprintf("No (%d)", cnt["No"])) +cols <- rgraph_updown_colors(rgraph_opt(p, "color_scheme", "rb")) + +pl <- ggplot2::ggplot(d, ggplot2::aes(x = log2FoldChange, y = -log10(.data[[sig_metric]]), + color = color)) + + ggplot2::geom_point(alpha = 0.5, size = 1) + + ggplot2::scale_color_manual(values = cols, breaks = c("Up", "Down", "No"), labels = labs3) + + ggplot2::geom_hline(yintercept = -log10(pcut), linetype = "dashed", color = "grey") + + ggplot2::geom_vline(xintercept = c(-log2fc, log2fc), linetype = "dashed", color = "grey") + + ggplot2::labs(x = "log2 (Fold Change)", y = sprintf("-log10(%s)", sig_metric), + color = sprintf("%s<%s\n|log2FC|>%s", sig_metric, pcut, log2fc)) + + theme_rgraph() + + ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = 3))) + +# ---- 标注基因:自定义清单优先,其次 top-N ---------------------------------- +sel <- NULL +if (!is.null(p$genes)) { + sel <- d[d[[gene_col]] %in% p$genes, ] +} else if (rgraph_opt(p, "label_n", 0) > 0) { + sel <- d[d$color != "No", ] + sel <- sel[order(sel[[sig_metric]]), ] + sel <- utils::head(sel, rgraph_opt(p, "label_n", 10)) +} +if (!is.null(sel) && nrow(sel) > 0) { + aes_lab <- ggplot2::aes(x = log2FoldChange, y = -log10(.data[[sig_metric]]), + label = .data[[gene_col]]) + if (use_repel) pl <- pl + ggrepel::geom_text_repel(data = sel, mapping = aes_lab, + size = 3, color = "black", max.overlaps = Inf) + else pl <- pl + ggplot2::geom_text(data = sel, mapping = aes_lab, size = 3, color = "black", + vjust = -0.6) +} + +rgraph_save(pl, p$outdir, rgraph_opt(p, "name", "volcano"), + width = rgraph_opt(p, "width", 8), height = rgraph_opt(p, "height", 6.2), + dpi = rgraph_opt(p, "dpi", 300)) +cat("RGRAPH_DONE: volcano up=", cnt["Up"], "down=", cnt["Down"], "no=", cnt["No"], "\n") diff --git a/plugins/Presisitence/rgraph/rscripts/wgcna.R b/plugins/Presisitence/rgraph/rscripts/wgcna.R new file mode 100644 index 0000000..9312a3c --- /dev/null +++ b/plugins/Presisitence/rgraph/rscripts/wgcna.R @@ -0,0 +1,125 @@ +# wgcna.R —— WGCNA 加权共表达网络分析(软阈值/模块/模块-分组相关/特征向量网络/Cytoscape 导出) +# 输入: expr(gene_expression.csv: gene_id + 各样本), sample_group(需 group_name) +# params: expr, sample_group, outdir, power(0=自动, 否则手动), min_module_size(30), +# merge_cut(0.25), network_type("unsigned"|"signed"), min_mean(1), +# max_genes(方差过滤上限,默认3000), export_module(可选:导出该颜色模块的 Cytoscape 边表), +# tomplot(bool,默认F,较慢), n_threads(0=全部), dpi(300) +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- rgraph_load_params() +rgraph_library(c("WGCNA", "dplyr")) +options(stringsAsFactors = FALSE) + +expr <- rgraph_read_matrix(p$expr) +sg <- rgraph_read_sample_group(p$sample_group) +rgraph_check_samples(expr, sg$sample_name) +rownames(expr) <- expr$gene_id +ge <- expr[, sg$sample_name, drop = FALSE] +ge[is.na(ge)] <- 0 +# 低表达过滤 + 方差 top 过滤(控制规模/耗时) +ns <- ncol(ge) +ge <- ge[rowSums(ge == 0) < ns * 0.5 & rowMeans(ge) >= rgraph_opt(p, "min_mean", 1), , drop = FALSE] +maxg <- rgraph_opt(p, "max_genes", 3000) +if (nrow(ge) > maxg) ge <- ge[order(apply(ge, 1, var), decreasing = TRUE)[seq_len(maxg)], , drop = FALSE] +datExpr <- as.matrix(t(ge)); mode(datExpr) <- "numeric" +nGenes <- ncol(datExpr); nSamples <- nrow(datExpr) + +nt <- rgraph_opt(p, "n_threads", 0) +if (is.numeric(nt) && nt >= 2) WGCNA::enableWGCNAThreads(nThreads = nt) else try(WGCNA::allowWGCNAThreads(), silent = TRUE) +net_type <- rgraph_opt(p, "network_type", "unsigned") + +# ---- 软阈值 ---------------------------------------------------------------- +powers <- 1:20 +sft <- WGCNA::pickSoftThreshold(datExpr, powerVector = powers, verbose = 0, networkType = net_type) +rgraph_save_base(p$outdir, "01.power", function() { + par(mfrow = c(1, 2)); cex1 <- 0.9 + plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], + xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Model Fit, signed R^2", + type = "n", main = "Scale independence") + text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels = powers, cex = cex1, col = "red") + abline(h = 0.85, col = "red") + plot(sft$fitIndices[, 1], sft$fitIndices[, 5], xlab = "Soft Threshold (power)", + ylab = "Mean Connectivity", type = "n", main = "Mean connectivity") + text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers, cex = cex1, col = "red") +}, width = 12, height = 7, dpi = rgraph_opt(p, "dpi", 300)) +power <- rgraph_opt(p, "power", 0) +if (is.null(power) || power <= 0) power <- if (!is.na(sft$powerEstimate)) sft$powerEstimate else 6 +cat("RGRAPH_INFO: power=", power, "\n") + +# ---- 构建网络 -------------------------------------------------------------- +net <- WGCNA::blockwiseModules(datExpr, corType = "pearson", power = power, + networkType = net_type, TOMType = net_type, + deepSplit = 2, minModuleSize = rgraph_opt(p, "min_module_size", 30), + mergeCutHeight = rgraph_opt(p, "merge_cut", 0.25), numericLabels = TRUE, + maxBlockSize = nGenes + 1, saveTOMs = FALSE, nThreads = 0) +moduleColors <- WGCNA::labels2colors(net$colors) + +rgraph_save_base(p$outdir, "02.dendrogram", function() { + WGCNA::plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], + "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05) +}, width = 9, height = 5, dpi = rgraph_opt(p, "dpi", 300)) + +rgraph_write_csv(data.frame(gene_id = colnames(datExpr), Module = moduleColors), + p$outdir, "gene_module") + +# ---- 模块-分组相关性 ------------------------------------------------------- +MEs <- WGCNA::orderMEs(WGCNA::moduleEigengenes(datExpr, moduleColors)$eigengenes) +gfac <- as.factor(sg$group_name[match(rownames(datExpr), sg$sample_name)]) +gdesign <- model.matrix(~ gfac - 1); colnames(gdesign) <- levels(gfac) +mgc <- cor(MEs, gdesign, use = "p") +mgp <- WGCNA::corPvalueStudent(mgc, nSamples) +txt <- paste(signif(mgc, 2), "\n(", signif(mgp, 1), ")", sep = ""); dim(txt) <- dim(mgc) +rgraph_save_base(p$outdir, "05.module-group", function() { + par(mar = c(8, 10, 4, 4)) + WGCNA::labeledHeatmap(Matrix = mgc, xLabels = colnames(gdesign), yLabels = names(MEs), + ySymbols = names(MEs), colorLabels = FALSE, colors = WGCNA::blueWhiteRed(50), + textMatrix = txt, setStdMargins = FALSE, cex.text = 0.7, zlim = c(-1, 1), + main = "Module-group relationships") +}, width = 10, height = 8, dpi = rgraph_opt(p, "dpi", 300)) + +# ---- MM(基因-模块隶属度) + geneInfo --------------------------------------- +modNames <- substring(names(MEs), 3) +MM <- as.data.frame(cor(datExpr, MEs, use = "p")); names(MM) <- paste0("MM", modNames) +MMp <- as.data.frame(WGCNA::corPvalueStudent(as.matrix(MM), nSamples)); names(MMp) <- paste0("p.MM", modNames) +geneInfo <- cbind(gene_id = colnames(datExpr), Module = moduleColors, MM, MMp) +rgraph_write_csv(geneInfo, p$outdir, "geneInfo") + +# ---- 特征向量邻接热图 ------------------------------------------------------ +mecor <- cor(MEs, use = "p") +rgraph_save_base(p$outdir, "08.eigengene_adjacency", function() { + par(mar = c(8, 10, 4, 3)) + WGCNA::labeledHeatmap(Matrix = mecor, xLabels = names(MEs), yLabels = names(MEs), + ySymbols = names(MEs), colorLabels = FALSE, colors = WGCNA::blueWhiteRed(50), + setStdMargins = FALSE, cex.text = 0.7, zlim = c(-1, 1), main = "Eigengene adjacency") +}, width = 8, height = 7, dpi = rgraph_opt(p, "dpi", 300)) + +# ---- 可选: TOM 网络热图(较慢)-------------------------------------------- +if (isTRUE(rgraph_opt(p, "tomplot", FALSE))) { + TOM <- WGCNA::TOMsimilarityFromExpr(datExpr, power = power, networkType = net_type) + nSelect <- min(400, nGenes); set.seed(10); sel <- sample(nGenes, nSelect) + selTOM <- (1 - TOM)[sel, sel]; selTree <- hclust(as.dist(selTOM), method = "average") + plotDiss <- selTOM^7; diag(plotDiss) <- NA + rgraph_save_base(p$outdir, "06.network_heatmap", + function() WGCNA::TOMplot(plotDiss, selTree, moduleColors[sel], main = "Network heatmap (selected genes)"), + width = 8, height = 8, dpi = rgraph_opt(p, "dpi", 300)) + # 可选 Cytoscape 导出(与 rgraph_network edge 模式串联) + em <- rgraph_opt(p, "export_module", NULL) + if (is.null(em)) { # 未指定时自动取最大非-grey 模块 + tb <- sort(table(moduleColors[moduleColors != "grey"]), decreasing = TRUE) + if (length(tb) > 0) em <- names(tb)[1] + } + if (!is.null(em)) { + inMod <- moduleColors == em; probes <- colnames(datExpr)[inMod] + modTOM <- TOM[inMod, inMod]; dimnames(modTOM) <- list(probes, probes) + ef <- file.path(p$outdir, paste0("Cytoscape_edges_", em, ".txt")) + nf <- file.path(p$outdir, paste0("Cytoscape_nodes_", em, ".txt")) + WGCNA::exportNetworkToCytoscape(modTOM, weighted = TRUE, + threshold = rgraph_opt(p, "export_threshold", 0.1), + nodeNames = probes, nodeAttr = moduleColors[inMod], + edgeFile = ef, nodeFile = nf) + rgraph_emit(ef); rgraph_emit(nf) + cat("RGRAPH_INFO: exported module=", em, "nodes=", length(probes), "\n") + } +} +cat("RGRAPH_DONE: wgcna genes=", nGenes, "samples=", nSamples, "modules=", + length(unique(moduleColors)), "\n") diff --git a/plugins/Presisitence/rgraph/server.py b/plugins/Presisitence/rgraph/server.py new file mode 100644 index 0000000..bea2bfa --- /dev/null +++ b/plugins/Presisitence/rgraph/server.py @@ -0,0 +1,421 @@ +"""rgraph —— RNAseq 下游可视化/分析 MCP Server (FastMCP)。 + +把一套清洗过、参数化的 R 出图脚本(rscripts/*.R)暴露为 MCP 工具,由本机 Rscript 引擎渲染, +产出与课程一致风格的 png+pdf。定位不到 R 时返回可手动运行的脚本与命令;R 端缺包时返回缺失 +包名与安装建议。 + +运行: + uv run --directory server.py +""" +from __future__ import annotations + +from typing import Any + +from mcp.server.fastmcp import FastMCP + +from rgraph_toolkit import runner + +mcp = FastMCP("rgraph") + + +def _clean(**kw: Any) -> dict[str, Any]: + """丢弃值为 None 的参数,避免污染 R 端 params。""" + return {k: v for k, v in kw.items() if v is not None} + + +def _run(script: str, params: dict[str, Any], outdir: str) -> dict[str, Any]: + return runner.run_script(script, params, outdir=outdir) + + +# ============================ 环境自检 ============================ +@mcp.tool() +def rgraph_env() -> dict[str, Any]: + """检查 R 引擎与关键 R 包是否就绪(首次使用/报错排查先跑这个)。返回 Rscript 路径、 + 各包安装状态与缺包安装建议。""" + rscript = runner.find_rscript() + pkgs = [ + "ggplot2", "dplyr", "tidyr", "stringr", "pheatmap", "ggrepel", "reshape2", + "DESeq2", "edgeR", "limma", "clusterProfiler", "enrichplot", "pathview", + "circlize", "ComplexHeatmap", "WGCNA", "VennDiagram", "GSVA", "GSEABase", "fgsea", + "ggpubr", "patchwork", "ggridges", "ggbreak", "ggalluvial", "plotly", "factoextra", + "igraph", "ggraph", "tidygraph", "ggforce", "ggNetView", "GENIE3", + ] + status = runner.check_packages(pkgs, rscript=rscript) if rscript else {p: False for p in pkgs} + missing = [p for p, ok in status.items() if not ok] + return { + "rscript": rscript or "未找到 (请设置环境变量 RGRAPH_RSCRIPT 指向 Rscript.exe)", + "rscripts_dir": str(runner.rscripts_dir()), + "packages": status, + "missing": missing, + "install_hint": runner._install_hint(missing) if missing else "", + } + + +# ============================ 01 表达量归一化 ============================ +@mcp.tool() +def rgraph_normalize(count: str, outdir: str, method: str = "fpkm", + sample_group: str | None = None) -> dict[str, Any]: + """由 count 计算 FPKM/TPM。count 为 gene_count.csv(列: gene_id,Length,<各样本>)。 + method: fpkm/tpm。sample_group 可选(限定样本列)。产出 gene_fpkm.csv 或 gene_tpm.csv。""" + return _run("normalize.R", _clean(count=count, method=method, sample_group=sample_group), outdir) + + +# ============================ 02 相关性 / PCA ============================ +@mcp.tool() +def rgraph_correlation(fpkm: str, sample_group: str, outdir: str, + dpi: int = 300, width: float = 6, height: float = 6) -> dict[str, Any]: + """样本间 Pearson 相关性热图(基于 log2(fpkm+1))。fpkm: gene_fpkm.csv;sample_group 含 sample_name。""" + return _run("correlation_heatmap.R", + _clean(fpkm=fpkm, sample_group=sample_group, dpi=dpi, width=width, height=height), outdir) + + +@mcp.tool() +def rgraph_pca(fpkm: str, sample_group: str, outdir: str, ellipse: bool = False, + label: bool = True, palette: str = "course", dpi: int = 300) -> dict[str, Any]: + """样本 PCA 散点图(log2(fpkm+1)+prcomp)。ellipse: 加 95% 置信椭圆;label: ggrepel 样本标签; + palette: course/colorblind。sample_group 需含 group_name。""" + return _run("pca.R", _clean(fpkm=fpkm, sample_group=sample_group, ellipse=ellipse, + label=label, palette=palette, dpi=dpi), outdir) + + +# ============================ 03 表达分布 ============================ +@mcp.tool() +def rgraph_distribution(fpkm: str, sample_group: str, outdir: str, kind: str = "box", + palette: str = "course", dpi: int = 300) -> dict[str, Any]: + """表达分布图。kind: box(箱线)/violin(小提琴)/density(密度)。按 group_name 着色。""" + return _run("distribution.R", _clean(fpkm=fpkm, sample_group=sample_group, kind=kind, + palette=palette, dpi=dpi), outdir) + + +# ============================ 04 差异分析 ============================ +@mcp.tool() +def rgraph_diff(count: str, sample_group: str, outdir: str, method: str = "deseq2", + sig_metric: str = "padj", pcut: float = 0.05, log2fc: float = 1.0, + min_mean_count: float = 1.0) -> dict[str, Any]: + """差异表达分析。method: deseq2/edger/edger_norep(无生物学重复)/limma。 + sig_metric: padj(推荐,默认)/pvalue;pcut、log2fc 为阈值。分组由 sample_group 的 TvsC 列 + (treatment/control)决定。产出 01.Dse2_result.csv 及 Deg_all/up/down。 + 注:DESeq2 需已安装(本机默认缺,可改 method=edger/limma 或按提示安装)。""" + return _run("diff.R", _clean(count=count, sample_group=sample_group, method=method, + sig_metric=sig_metric, pcut=pcut, log2fc=log2fc, + min_mean_count=min_mean_count), outdir) + + +@mcp.tool() +def rgraph_volcano(result: str, outdir: str, sig_metric: str = "padj", pcut: float = 0.05, + log2fc: float = 1.0, color_scheme: str = "rb", label_n: int = 0, + genes: list[str] | None = None) -> dict[str, Any]: + """火山图。result 为差异结果表(含 gene_id,log2FoldChange,pvalue,padj)。 + color_scheme: rb(色盲友好红蓝,默认)/rg(课程红绿);label_n: 标注显著性最高的前 N 个; + genes: 自定义标注基因清单(优先于 label_n)。""" + return _run("volcano.R", _clean(result=result, sig_metric=sig_metric, pcut=pcut, log2fc=log2fc, + color_scheme=color_scheme, label_n=label_n, genes=genes), outdir) + + +@mcp.tool() +def rgraph_heatmap(fpkm: str, sample_group: str, outdir: str, deg: str | None = None, + genes: list[str] | None = None, engine: str = "pheatmap", + scale: str = "row", cluster_cols: bool = True, + show_rownames: bool | None = None) -> dict[str, Any]: + """差异/指定基因表达聚类热图(log2(fpkm+1))。基因来源: genes 优先,其次 deg(Deg_all.csv), + 否则全部。engine: pheatmap(默认) 或 complexheatmap(按 group_name 加顶部分组色块、列不聚类)。""" + return _run("heatmap.R", _clean(fpkm=fpkm, sample_group=sample_group, deg=deg, genes=genes, + engine=engine, scale=scale, cluster_cols=cluster_cols, + show_rownames=show_rownames), outdir) + + +# ============================ 05 富集分析 (GO/KEGG) ============================ +@mcp.tool() +def rgraph_enrich(gene_list: str, outdir: str, type: str = "go", orgdb: str = "org.Hs.eg.db", + id_type: str = "ENSEMBL", ont: str = "all", go_source: str = "orgdb", + kegg_source: str = "online", kegg_species: str = "hsa", + gmt: str | None = None, kegg_info: str | None = None, + pvalueCutoff: float = 1.0, qvalueCutoff: float = 1.0) -> dict[str, Any]: + """GO/KEGG 富集分析。gene_list 为含 gene_id 列的 CSV。type: go/kegg。 + orgdb: 物种注释包(模式物种如 org.Hs.eg.db;非模式需先自建 org.Xxx.eg.db)。 + go_source: orgdb(enrichGO) 或 custom(从 OrgDb 的 GO 列自建, 非模式物种)。 + kegg_source: online(rest.kegg.jp, 需 kegg_species) 或 gmt(用自建 gmt)。 + 产出 GO_enrich.csv / KEGG_enrich.csv。注:OrgDb 需已安装,否则按提示安装。""" + return _run("enrich.R", _clean(gene_list=gene_list, type=type, orgdb=orgdb, id_type=id_type, + ont=ont, go_source=go_source, kegg_source=kegg_source, + kegg_species=kegg_species, gmt=gmt, kegg_info=kegg_info, + pvalueCutoff=pvalueCutoff, qvalueCutoff=qvalueCutoff), outdir) + + +@mcp.tool() +def rgraph_go_plot(enrich: str, outdir: str, kind: str = "bar", top_n: int = 10) -> dict[str, Any]: + """GO 富集图(BP/CC/MF 三面板纵向拼接)。enrich 为 GO_enrich.csv(需 ONTOLOGY/Description/ + GeneRatio/pvalue/Count)。kind: bar(x=-log10p,填充Count) / dot(x=GeneRatio,填充-log10p,大小Count)。""" + return _run("go_plot.R", _clean(enrich=enrich, kind=kind, top_n=top_n), outdir) + + +@mcp.tool() +def rgraph_kegg_plot(enrich: str, outdir: str, kind: str = "dot", top_n: int = 30) -> dict[str, Any]: + """KEGG 富集图。enrich 为 KEGG_enrich.csv(需 Description/GeneRatio/pvalue/Count)。 + kind: dot(气泡) / bar(柱形)。top_n: 取 pvalue 最小的前 N 条通路。""" + return _run("kegg_plot.R", _clean(enrich=enrich, kind=kind, top_n=top_n), outdir) + + +# ============================ 06 PPI ============================ +@mcp.tool() +def rgraph_ppi(deg: str, info: str, links: str, outdir: str, score: int = 400, + id_is_symbol: bool = False, orgdb: str = "org.Hs.eg.db") -> dict[str, Any]: + """由差异基因 + STRING 库构建蛋白互作边表。deg: Deg_all.csv(gene_id);info/links 为 STRING 的 + protein.info 与 protein.links 文件;score: combined_score 阈值。id_is_symbol=False 时用 orgdb 做 + ENSEMBL→SYMBOL。产出 Target_PPi.tsv + Edge_Node_count.tsv(可导入 Cytoscape)。""" + return _run("ppi.R", _clean(deg=deg, info=info, links=links, score=score, + id_is_symbol=id_is_symbol, orgdb=orgdb), outdir) + + +# ============================ 07 GSEA ============================ +@mcp.tool() +def rgraph_build_gmt(outdir: str, orgdb: str, type: str = "go", source: str = "orgdb", + keytype: str = "GID", kegg_species: str = "hsa", + min_count: int = 5) -> dict[str, Any]: + """为 GSEA 构建 gmt 基因集。type: go/kegg。source: orgdb(从注释包的 GO/Pathway 列) 或 + online(KEGG 从 rest.kegg.jp + OrgDb 映射)。产出 GSEA_GO.gmt / GSEA_KEGG.gmt。需 orgdb 已安装。""" + return _run("build_gmt.R", _clean(type=type, orgdb=orgdb, source=source, keytype=keytype, + kegg_species=kegg_species, min_count=min_count), outdir) + + +@mcp.tool() +def rgraph_gsea(result: str, gmt: str, outdir: str, minGSSize: int = 5, maxGSSize: int = 1000, + pvalueCutoff: float = 1.0, top_n: int = 5, desc_map: str | None = None, + draw: bool = True) -> dict[str, Any]: + """GSEA 分析。result 为差异结果表(按 log2FoldChange 排序),gmt 为基因集文件。 + 产出 GSEA_result.csv 及各方向 top-N 通路的 ES 图(UP/DOWN 子目录)。 + desc_map(可选 CSV: ID,Description) 用于补充通路描述。""" + return _run("gsea.R", _clean(result=result, gmt=gmt, minGSSize=minGSSize, maxGSSize=maxGSSize, + pvalueCutoff=pvalueCutoff, top_n=top_n, desc_map=desc_map, + draw=draw), outdir) + + +@mcp.tool() +def rgraph_gsea_mountain(gsea_result: str, result: str, outdir: str, + geneset_ids: list[str] | None = None, n_top: int = 15) -> dict[str, Any]: + """GSEA 山峦图(多通路 core_enrichment 基因的 log2FC 分布 + NES 点)。gsea_result 为 rgraph_gsea + 产出的 CSV;result 为差异结果表(取基因 log2FC)。geneset_ids 指定通路,缺省取 pvalue 最小的前 n_top。""" + return _run("gsea_mountain.R", _clean(gsea_result=gsea_result, result=result, + geneset_ids=geneset_ids, n_top=n_top), outdir) + + +# ============================ 非流程化绘图 (高级) ============================ +@mcp.tool() +def rgraph_venn(outdir: str, sets: list[str] | None = None, folder: str | None = None, + col: str | None = None, names: list[str] | None = None) -> dict[str, Any]: + """韦恩图(2-5 组) + 交集成员表。sets: 各集合 CSV 路径列表;folder: 装 CSV 的目录(二选一)。 + col: 元素列名(默认自动识别 Element_name/gene_id/首列)。产出 venn.png/pdf + veen_table.csv。""" + return _run("venn.R", _clean(sets=sets, folder=folder, col=col, names=names), outdir) + + +@mcp.tool() +def rgraph_quadrant(x: str, y: str, outdir: str, mode: str = "nine", x_name: str = "compare1", + y_name: str = "compare2", pcut: float = 0.05, log2fc: float = 1.0) -> dict[str, Any]: + """九象限图/二象限图。x、y 为两个比较的差异结果 CSV(各含 gene_id,log2FoldChange,pvalue)。 + mode: nine(九象限) / four(二象限,仅同向基因 High/Low)。产出散点图 + 象限表。""" + return _run("quadrant.R", _clean(x=x, y=y, mode=mode, x_name=x_name, y_name=y_name, + pcut=pcut, log2fc=log2fc), outdir) + + +@mcp.tool() +def rgraph_stacked_volcano(deg_data: str, outdir: str, sig_metric: str = "pvalue", + pcut: float = 0.05, log2fc: float = 1.0, group_col: str = "TvsC", + top_n: int = 5) -> dict[str, Any]: + """堆叠火山图(多比较组)。deg_data 含 gene_id,log2FoldChange,pvalue 及分组列(group_col,默认 TvsC)。 + 按组抽动散点,点大小/深浅=|log2FC|,红上蓝下,每组每方向标注 top_n 个基因。""" + return _run("stacked_volcano.R", _clean(deg_data=deg_data, sig_metric=sig_metric, pcut=pcut, + log2fc=log2fc, group_col=group_col, top_n=top_n), outdir) + + +@mcp.tool() +def rgraph_gene_bar(deg: str, sample_group: str, outdir: str, gene_list: str | None = None, + genes: list[str] | None = None) -> dict[str, Any]: + """单基因表达柱形图(对照 vs 处理,均值±SEM,显著性星号)。deg 为 Deg_all.csv(含 pvalue 与各 + 样本归一化列);sample_group 需 TvsC。基因来源: genes 优先,否则 gene_list(CSV)。每基因一张图。""" + return _run("gene_bar.R", _clean(deg=deg, sample_group=sample_group, gene_list=gene_list, + genes=genes), outdir) + + +@mcp.tool() +def rgraph_gene_correlation(fpkm: str, sample_group: str, outdir: str, + gene_list1: str | None = None, gene_list2: str | None = None, + genes1: list[str] | None = None, + genes2: list[str] | None = None) -> dict[str, Any]: + """基因-基因相关性气泡热图(Pearson, 基于 log2(fpkm+1))。行基因=genes1/gene_list1, + 列基因=genes2/gene_list2(缺省=genes1)。圆点大小=|r|,蓝-白-红。产出热图 + correlation.csv。""" + return _run("gene_correlation.R", _clean(fpkm=fpkm, sample_group=sample_group, + gene_list1=gene_list1, gene_list2=gene_list2, + genes1=genes1, genes2=genes2), outdir) + + +@mcp.tool() +def rgraph_pathway_trend(enrich: str, deg: str, outdir: str, kind: str = "diverge", + top_n: int | None = None) -> dict[str, Any]: + """通路上下调基因统计图。enrich 为富集表(需 Description 与 geneID 列, 形如 a/b/c); + deg 为差异结果表(需 gene_id,log2FoldChange)。kind: diverge(左右发散条形) / stack(堆叠条形)。 + 注: 默认 enrich 的 geneID 与 deg 的 gene_id 为同一 ID 类型(直接匹配)。""" + return _run("pathway_trend.R", _clean(enrich=enrich, deg=deg, kind=kind, top_n=top_n), outdir) + + +@mcp.tool() +def rgraph_segmented_volcano(result: str, outdir: str, sig_metric: str = "pvalue", + pcut: float = 0.05, log2fc: float = 1.0, color_scheme: str = "rb", + break_lower: float | None = None, + break_upper: float | None = None) -> dict[str, Any]: + """分段式火山图(y 轴断裂,适合个别基因 p 值极端时)。break_lower/upper 手动指定断裂区间, + 缺省自动推断。需 ggbreak 包(本机默认缺,按提示 install.packages('ggbreak'))。""" + return _run("segmented_volcano.R", _clean(result=result, sig_metric=sig_metric, pcut=pcut, + log2fc=log2fc, color_scheme=color_scheme, + break_lower=break_lower, break_upper=break_upper), outdir) + + +@mcp.tool() +def rgraph_pathway_gene_heatmap(enrich: str, fpkm: str, sample_group: str, outdir: str, + kind: str = "correlation", id_col: str = "geneID", + min_count: int = 3, max_pathways: int = 6) -> dict[str, Any]: + """通路富集基因的相关性热图/聚类热图(逐通路出图)。enrich 为 GO/KEGG 富集表 + (需 Description 与基因列 id_col, 默认 geneID);fpkm 为表达矩阵。kind: correlation(基因-基因 + Pearson) / cluster(pheatmap 聚类热图)。注: id_col 基因需与 fpkm 的 gene_id 同 ID 类型。""" + return _run("pathway_gene_heatmap.R", _clean(enrich=enrich, fpkm=fpkm, sample_group=sample_group, + kind=kind, id_col=id_col, min_count=min_count, + max_pathways=max_pathways), outdir) + + +@mcp.tool() +def rgraph_kmeans(fpkm: str, sample_group: str, outdir: str, gene_list: str | None = None, + genes: list[str] | None = None, k: int = 4, kmax: int = 10) -> dict[str, Any]: + """目标基因 K-means 聚类:肘部图 + 聚类热图(行按簇排序、行/列注释) + 各簇表达趋势线 + + 各簇基因成员表。基因来源: genes 优先,否则 gene_list(CSV)。k: 簇数(可先看肘部图定 k)。""" + return _run("kmeans.R", _clean(fpkm=fpkm, sample_group=sample_group, gene_list=gene_list, + genes=genes, k=k, kmax=kmax), outdir) + + +@mcp.tool() +def rgraph_interactive_volcano(result: str, outdir: str, sig_metric: str = "pvalue", + pcut: float = 0.05, log2fc: float = 1.0, + color_scheme: str = "rb", include_no: bool = True) -> dict[str, Any]: + """交互式火山图(plotly HTML,悬停显示基因/log2FC/-log10p)。result 为差异结果表。 + 产出 .html(默认非 selfcontained,伴生 _files 目录,无需 pandoc)。""" + return _run("interactive_volcano.R", _clean(result=result, sig_metric=sig_metric, pcut=pcut, + log2fc=log2fc, color_scheme=color_scheme, + include_no=include_no), outdir) + + +@mcp.tool() +def rgraph_sankey(enrich: str, outdir: str, id_col: str = "geneID", + pathway_col: str = "Description", top_n: int = 8) -> dict[str, Any]: + """富集通路-基因桑基图(ggalluvial; 基因→通路流向)。enrich 为富集表(需通路名列 + pathway_col 与基因列 id_col, 形如 a/b/c)。top_n 限制展示通路数。""" + return _run("sankey.R", _clean(enrich=enrich, id_col=id_col, pathway_col=pathway_col, + top_n=top_n), outdir) + + +@mcp.tool() +def rgraph_ssgsea(expr: str, gmt: str, sample_group: str, outdir: str, + normalize: bool = True) -> dict[str, Any]: + """ssGSEA 单样本基因集富集打分 + 各基因集分组箱线图。expr 为表达矩阵,gmt 为基因集。 + 需 GSVA + GSEABase(本机默认缺,按提示 BiocManager::install(c('GSVA','GSEABase')))。""" + return _run("ssgsea.R", _clean(expr=expr, gmt=gmt, sample_group=sample_group, + normalize=normalize), outdir) + + +@mcp.tool() +def rgraph_multilevel_scatter(expr: str, sample_group: str, outdir: str, + gene_list: str | None = None, genes: list[str] | None = None, + x_col: str | None = None, color_col: str | None = None, + sample_col: str | None = None) -> dict[str, Any]: + """多级分组基因表达散点图(逐基因)。主分组为 x、次分组着色。sample_group 可含 + sample_name(1)/group_name(1)/group_name2 等多级列(缺省自动识别)。基因来源: genes 或 gene_list。""" + return _run("multilevel_scatter.R", _clean(expr=expr, sample_group=sample_group, + gene_list=gene_list, genes=genes, x_col=x_col, + color_col=color_col, sample_col=sample_col), outdir) + + +@mcp.tool() +def rgraph_norm_matrix(count: str, sample_group: str, outdir: str, method: str = "deseq2", + min_mean: float = 1.0) -> dict[str, Any]: + """输出标准化表达矩阵。method: deseq2(中位数比值标准化 count) / limma(voom logCPM)。 + 产出 normalized_count.csv。注: deseq2 需安装 DESeq2;limma 本机已装。""" + return _run("norm_matrix.R", _clean(count=count, sample_group=sample_group, method=method, + min_mean=min_mean), outdir) + + +@mcp.tool() +def rgraph_interaction_diff(count: str, sample_group: str, outdir: str, sig_metric: str = "padj", + pcut: float = 0.05, log2fc: float = 1.0, + min_mean_count: float = 2.0) -> dict[str, Any]: + """两因子交互作用差异分析(DESeq2 ~g1+g2+g1:g2,提取交互项)。sample_group 需含 + group_name1,group_name2,TvsC1,TvsC2。产出 01.DEseq2_result + Deg_all/down/up(可喂 rgraph_volcano)。 + 需 DESeq2(本机默认缺,按提示安装)。""" + return _run("interaction_diff.R", _clean(count=count, sample_group=sample_group, + sig_metric=sig_metric, pcut=pcut, log2fc=log2fc, + min_mean_count=min_mean_count), outdir) + + +@mcp.tool() +def rgraph_wgcna(expr: str, sample_group: str, outdir: str, power: int = 0, + min_module_size: int = 30, merge_cut: float = 0.25, + network_type: str = "unsigned", max_genes: int = 3000, + tomplot: bool = False, export_module: str | None = None, + export_threshold: float = 0.1) -> dict[str, Any]: + """WGCNA 加权共表达网络分析。产出:软阈值图、模块树状图、模块-分组相关热图、 + 特征向量邻接热图、gene_module/geneInfo 表。power=0 自动选最佳软阈值;max_genes 按方差限制 + 基因数(控制耗时)。tomplot=True 时额外绘网络热图并导出 Cytoscape 边/点表(用于与 rgraph_network + 串联);export_module 指定导出模块(缺省自动取最大非-grey 模块),export_threshold 控制边数。""" + return _run("wgcna.R", _clean(expr=expr, sample_group=sample_group, power=power, + min_module_size=min_module_size, merge_cut=merge_cut, + network_type=network_type, max_genes=max_genes, + tomplot=tomplot, export_module=export_module, + export_threshold=export_threshold), outdir) + + +@mcp.tool() +def rgraph_enrich_circos(enrich: str, outdir: str, type: str = "go", topN: int = 10, + split_count: bool | None = None, show_axis: bool = True) -> dict[str, Any]: + """富集分析圈图(circlize 四轨道: 分类/背景基因数/前景基因(上下调)/RichFactor)。 + enrich 需含 pvalue,RichFactor,Bg_gene,Count[,Up,Down] 及 ONTOLOGY(GO)/Category(KEGG)。 + type: go/kegg。图例需 ComplexHeatmap(缺失则略去图例,圈图主体照常输出)。""" + return _run("enrich_circos.R", _clean(enrich=enrich, type=type, topN=topN, + split_count=split_count, show_axis=show_axis), outdir) + + +@mcp.tool() +def rgraph_virtual_knockout(count: str, sample_group: str, outdir: str, ko_gene: str, + ko_group: str, logfc: float = 0.5, hv_num: int = 5000) -> dict[str, Any]: + """bulk RNA 虚拟敲除(GENIE3 推断调控网络 + bulkKnk 扰动)。sample_group 需含 group_name + 与 pheno_data。ko_gene/ko_group 为要敲除的基因与样本分组。产出扰动汇总/显著扰动基因 + 表 + 双轴火山图。依赖 bulkKnk(课程本地源码包,install.packages('bulkKnk-main',repos=NULL,type='source')) + 与 GENIE3(Bioconductor)——本机默认缺,会返回缺包提示。""" + return _run("virtual_knockout.R", _clean(count=count, sample_group=sample_group, + ko_gene=ko_gene, ko_group=ko_group, logfc=logfc, + hv_num=hv_num), outdir) + + +@mcp.tool() +def rgraph_network(outdir: str, mode: str = "matrix", expr: str | None = None, + sample_group: str | None = None, genes: list[str] | None = None, + method: str = "WGCNA", cor_method: str = "pearson", r_threshold: float = 0.7, + p_threshold: float = 0.05, transform: str = "none", max_features: int = 500, + node_annotation: str | None = None, edges: str | None = None, + nodes: str | None = None, from_col: str | None = None, to_col: str | None = None, + weight_col: str | None = None, directed: bool = False, layout: str = "gephi", + layout_module: str = "random", group_by: str = "Modularity", + fill_by: str = "Modularity", label: bool = False, pointlabel: str | None = None, + add_outer: bool = False) -> dict[str, Any]: + """网络图(ggNetView):共表达/相关/宏基因组共现网络,按模块着色。 + mode="matrix": expr(gene_id/OTU+样本) 构相关网络(method: WGCNA/cor/SPARCC/SpiecEasi/Hmisc; + genes 限定子集否则按方差取 max_features)。mode="edge": edges(前两列 from,to[+weight]) + [+nodes] —— 可直接吃 rgraph_wgcna 的 Cytoscape 边表、rgraph_ppi 的 Target_PPi。 + layout: gephi 等;label/pointlabel(如 'top5')/add_outer 控制标签与模块边界。需 ggNetView 包。""" + return _run("network.R", _clean(mode=mode, expr=expr, sample_group=sample_group, genes=genes, + method=method, cor_method=cor_method, r_threshold=r_threshold, + p_threshold=p_threshold, transform=transform, + max_features=max_features, node_annotation=node_annotation, + edges=edges, nodes=nodes, from_col=from_col, to_col=to_col, + weight_col=weight_col, directed=directed, layout=layout, + layout_module=layout_module, group_by=group_by, fill_by=fill_by, + label=label, pointlabel=pointlabel, add_outer=add_outer), outdir) + + +if __name__ == "__main__": + mcp.run() diff --git a/plugins/Presisitence/rgraph/skills/rgraph/SKILL.md b/plugins/Presisitence/rgraph/skills/rgraph/SKILL.md new file mode 100644 index 0000000..d6005e7 --- /dev/null +++ b/plugins/Presisitence/rgraph/skills/rgraph/SKILL.md @@ -0,0 +1,19 @@ +--- +name: rgraph +description: Use when the user already has an RNA-seq count or expression matrix and needs downstream analysis or figures (normalize, PCA, DESeq2/edgeR/limma, volcano, heatmap, GO/KEGG, GSEA, WGCNA). Drive the rgraph MCP tools with local Rscript; do not redraw in Python. +--- + +# rgraph + +Use the `rgraph` MCP tools. Prefer R-rendered png+pdf over matplotlib copies. + +## Typical order + +1. `rgraph_env` — confirm Rscript and packages +2. `rgraph_normalize` / `rgraph_pca` / `rgraph_correlation` as needed +3. `rgraph_diff` — default significance metric is **padj**, not raw p-value +4. `rgraph_volcano`, `rgraph_heatmap`, `rgraph_enrich`, `rgraph_gsea`, `rgraph_wgcna` as requested + +Required table columns are in the Plugin README (`sample_name`/`group`, `gene_id`, counts). + +If Rscript is missing, return the generated `.R` script and the command to run it. If a package is missing, return the install hint from the tool. Do not pretend the figure was drawn. diff --git a/plugins/Presisitence/rgraph/tests/data/gene_count.csv b/plugins/Presisitence/rgraph/tests/data/gene_count.csv new file mode 100644 index 0000000..54b478e --- /dev/null +++ b/plugins/Presisitence/rgraph/tests/data/gene_count.csv @@ -0,0 +1,11 @@ +gene_id,Length,A1,A2,A3,B1,B2,B3 +g1,1000,100,120,110,10,12,9 +g2,2000,50,40,45,200,220,210 +g3,1500,0,0,0,0,0,0 +g4,500,30,35,28,33,31,36 +g5,800,5,7,6,100,90,95 +g6,1200,80,88,84,20,18,22 +g7,900,300,320,290,15,17,14 +g8,1100,12,10,11,140,150,160 +g9,700,60,64,58,62,66,59 +g10,1600,2,3,1,80,85,90 diff --git a/plugins/Presisitence/rgraph/tests/data/sample_group.csv b/plugins/Presisitence/rgraph/tests/data/sample_group.csv new file mode 100644 index 0000000..4c68281 --- /dev/null +++ b/plugins/Presisitence/rgraph/tests/data/sample_group.csv @@ -0,0 +1,7 @@ +sample_name,group_name,TvsC +A1,A,treatment +A2,A,treatment +A3,A,treatment +B1,B,control +B2,B,control +B3,B,control From f4cce672f13ff5d2d91aa9d9ab191706408d5e8a Mon Sep 17 00:00:00 2001 From: Presisitence Date: Wed, 19 Aug 2026 22:12:05 +0800 Subject: [PATCH 2/2] Rename plugin rgraph to rnaseq-plot Catalog name now states RNA-seq downstream plotting. MCP tool functions stay rgraph_*. --- plugins/Presisitence/rgraph/plugin.json | 22 ----------------- .../{rgraph => rnaseq-plot}/LICENSE | 0 .../{rgraph => rnaseq-plot}/README.md | 12 ++++++---- .../{rgraph => rnaseq-plot}/mcp.json | 2 +- plugins/Presisitence/rnaseq-plot/plugin.json | 24 +++++++++++++++++++ .../{rgraph => rnaseq-plot}/pyproject.toml | 6 ++--- .../rgraph_toolkit/__init__.py | 0 .../rgraph_toolkit/cli.py | 0 .../rgraph_toolkit/runner.py | 0 .../rscripts/_common.R | 0 .../rscripts/build_gmt.R | 0 .../rscripts/correlation_heatmap.R | 0 .../{rgraph => rnaseq-plot}/rscripts/diff.R | 0 .../rscripts/distribution.R | 0 .../{rgraph => rnaseq-plot}/rscripts/enrich.R | 0 .../rscripts/enrich_circos.R | 0 .../rscripts/gene_bar.R | 0 .../rscripts/gene_correlation.R | 0 .../rscripts/go_plot.R | 0 .../{rgraph => rnaseq-plot}/rscripts/gsea.R | 0 .../rscripts/gsea_mountain.R | 0 .../rscripts/heatmap.R | 0 .../rscripts/interaction_diff.R | 0 .../rscripts/interactive_volcano.R | 0 .../rscripts/kegg_plot.R | 0 .../{rgraph => rnaseq-plot}/rscripts/kmeans.R | 0 .../rscripts/multilevel_scatter.R | 0 .../rscripts/network.R | 0 .../rscripts/norm_matrix.R | 0 .../rscripts/normalize.R | 0 .../rscripts/pathway_gene_heatmap.R | 0 .../rscripts/pathway_trend.R | 0 .../{rgraph => rnaseq-plot}/rscripts/pca.R | 0 .../{rgraph => rnaseq-plot}/rscripts/ppi.R | 0 .../rscripts/quadrant.R | 0 .../{rgraph => rnaseq-plot}/rscripts/sankey.R | 0 .../rscripts/segmented_volcano.R | 0 .../{rgraph => rnaseq-plot}/rscripts/ssgsea.R | 0 .../rscripts/stacked_volcano.R | 0 .../{rgraph => rnaseq-plot}/rscripts/venn.R | 0 .../rscripts/virtual_knockout.R | 0 .../rscripts/volcano.R | 0 .../{rgraph => rnaseq-plot}/rscripts/wgcna.R | 0 .../{rgraph => rnaseq-plot}/server.py | 9 +++---- .../skills/rnaseq-plot}/SKILL.md | 8 +++---- .../tests/data/gene_count.csv | 0 .../tests/data/sample_group.csv | 0 47 files changed, 44 insertions(+), 39 deletions(-) delete mode 100644 plugins/Presisitence/rgraph/plugin.json rename plugins/Presisitence/{rgraph => rnaseq-plot}/LICENSE (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/README.md (73%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/mcp.json (91%) create mode 100644 plugins/Presisitence/rnaseq-plot/plugin.json rename plugins/Presisitence/{rgraph => rnaseq-plot}/pyproject.toml (64%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rgraph_toolkit/__init__.py (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rgraph_toolkit/cli.py (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rgraph_toolkit/runner.py (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/_common.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/build_gmt.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/correlation_heatmap.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/diff.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/distribution.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/enrich.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/enrich_circos.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/gene_bar.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/gene_correlation.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/go_plot.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/gsea.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/gsea_mountain.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/heatmap.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/interaction_diff.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/interactive_volcano.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/kegg_plot.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/kmeans.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/multilevel_scatter.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/network.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/norm_matrix.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/normalize.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/pathway_gene_heatmap.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/pathway_trend.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/pca.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/ppi.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/quadrant.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/sankey.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/segmented_volcano.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/ssgsea.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/stacked_volcano.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/venn.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/virtual_knockout.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/volcano.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/rscripts/wgcna.R (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/server.py (98%) rename plugins/Presisitence/{rgraph/skills/rgraph => rnaseq-plot/skills/rnaseq-plot}/SKILL.md (64%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/tests/data/gene_count.csv (100%) rename plugins/Presisitence/{rgraph => rnaseq-plot}/tests/data/sample_group.csv (100%) diff --git a/plugins/Presisitence/rgraph/plugin.json b/plugins/Presisitence/rgraph/plugin.json deleted file mode 100644 index 2db2c7d..0000000 --- a/plugins/Presisitence/rgraph/plugin.json +++ /dev/null @@ -1,22 +0,0 @@ -{ - "$schema": "https://agent-plugins.org/schemas/1.0.0/plugin.schema.json", - "name": "rgraph", - "version": "0.1.0", - "description": "RNA-seq downstream analysis and publication figures via local Rscript (DESeq2/edgeR/limma, volcano, heatmap, GSEA, WGCNA).", - "author": { - "name": "Presisitence", - "url": "https://github.com/Presisitence" - }, - "homepage": "https://github.com/Presisitence/rgraph-mcp", - "repository": "https://github.com/Presisitence/rgraph-mcp", - "license": "MIT", - "keywords": [ - "minimax-code", - "plugin", - "mcp", - "rnaseq", - "deseq2", - "ggplot2", - "wgcna" - ] -} diff --git a/plugins/Presisitence/rgraph/LICENSE b/plugins/Presisitence/rnaseq-plot/LICENSE similarity index 100% rename from plugins/Presisitence/rgraph/LICENSE rename to plugins/Presisitence/rnaseq-plot/LICENSE diff --git a/plugins/Presisitence/rgraph/README.md b/plugins/Presisitence/rnaseq-plot/README.md similarity index 73% rename from plugins/Presisitence/rgraph/README.md rename to plugins/Presisitence/rnaseq-plot/README.md index d89561d..c8a0ee6 100644 --- a/plugins/Presisitence/rgraph/README.md +++ b/plugins/Presisitence/rnaseq-plot/README.md @@ -1,12 +1,14 @@ -# rgraph +# rnaseq-plot -RNA-seq **downstream** analysis and figures for MiniMax Code. Parameterized R scripts -(ggplot2, pheatmap, clusterProfiler, edgeR, limma, WGCNA, …) are exposed as MCP tools and -rendered by the user's local **Rscript** to png + pdf. +RNA-seq **downstream plotting** for MiniMax Code (formerly published as `rgraph`). +Parameterized R scripts (ggplot2, pheatmap, clusterProfiler, edgeR, limma, WGCNA, …) +are exposed as MCP tools and rendered by the user's local **Rscript** to png + pdf. The Plugin consumes the user's own count/FPKM tables. It does not ship genomes or experimental matrices. `tests/data/` is a tiny synthetic `g1`–`g10` table for smoke tests. +MCP **tool names remain `rgraph_*`** (for example `rgraph_volcano`) so existing prompts keep working. + ## Try it ```text @@ -34,4 +36,4 @@ by default. If R or a package is missing, the tool returns install commands inst ## License -MIT. See [LICENSE](LICENSE). +MIT. See [LICENSE](LICENSE). Source: https://github.com/Presisitence/rnaseq-plot-mcp diff --git a/plugins/Presisitence/rgraph/mcp.json b/plugins/Presisitence/rnaseq-plot/mcp.json similarity index 91% rename from plugins/Presisitence/rgraph/mcp.json rename to plugins/Presisitence/rnaseq-plot/mcp.json index 1cdbe22..0e13111 100644 --- a/plugins/Presisitence/rgraph/mcp.json +++ b/plugins/Presisitence/rnaseq-plot/mcp.json @@ -1,7 +1,7 @@ { "$schema": "https://agent-plugins.org/schemas/1.0.0/mcp.schema.json", "mcpServers": { - "rgraph": { + "rnaseq-plot": { "type": "stdio", "command": "uv", "args": ["run", "server.py"], diff --git a/plugins/Presisitence/rnaseq-plot/plugin.json b/plugins/Presisitence/rnaseq-plot/plugin.json new file mode 100644 index 0000000..ab3aac5 --- /dev/null +++ b/plugins/Presisitence/rnaseq-plot/plugin.json @@ -0,0 +1,24 @@ +{ + "$schema": "https://agent-plugins.org/schemas/1.0.0/plugin.schema.json", + "name": "rnaseq-plot", + "version": "0.1.1", + "description": "RNA-seq downstream plotting via local Rscript: DESeq2/edgeR/limma, volcano, heatmap, GSEA, WGCNA (png+pdf).", + "author": { + "name": "Presisitence", + "url": "https://github.com/Presisitence" + }, + "homepage": "https://github.com/Presisitence/rnaseq-plot-mcp", + "repository": "https://github.com/Presisitence/rnaseq-plot-mcp", + "license": "MIT", + "keywords": [ + "minimax-code", + "plugin", + "mcp", + "rnaseq", + "deseq2", + "volcano", + "heatmap", + "ggplot2", + "wgcna" + ] +} diff --git a/plugins/Presisitence/rgraph/pyproject.toml b/plugins/Presisitence/rnaseq-plot/pyproject.toml similarity index 64% rename from plugins/Presisitence/rgraph/pyproject.toml rename to plugins/Presisitence/rnaseq-plot/pyproject.toml index 3808ce4..17c9184 100644 --- a/plugins/Presisitence/rgraph/pyproject.toml +++ b/plugins/Presisitence/rnaseq-plot/pyproject.toml @@ -1,7 +1,7 @@ [project] -name = "rgraph-mcp" -version = "0.1.0" -description = "RNAseq 下游可视化/分析 MCP 服务:把一套 R (ggplot2/pheatmap/clusterProfiler/WGCNA...) 出图脚本参数化封装为 MCP 工具,由 Rscript 引擎渲染" +name = "rnaseq-plot-mcp" +version = "0.1.1" +description = "RNA-seq 下游出图 MCP:把一套 R (ggplot2/pheatmap/clusterProfiler/WGCNA...) 脚本封装为工具,由本机 Rscript 渲染 png+pdf" readme = "README.md" requires-python = ">=3.10" license = "MIT" diff --git a/plugins/Presisitence/rgraph/rgraph_toolkit/__init__.py b/plugins/Presisitence/rnaseq-plot/rgraph_toolkit/__init__.py similarity index 100% rename from plugins/Presisitence/rgraph/rgraph_toolkit/__init__.py rename to plugins/Presisitence/rnaseq-plot/rgraph_toolkit/__init__.py diff --git a/plugins/Presisitence/rgraph/rgraph_toolkit/cli.py b/plugins/Presisitence/rnaseq-plot/rgraph_toolkit/cli.py similarity index 100% rename from plugins/Presisitence/rgraph/rgraph_toolkit/cli.py rename to plugins/Presisitence/rnaseq-plot/rgraph_toolkit/cli.py diff --git a/plugins/Presisitence/rgraph/rgraph_toolkit/runner.py b/plugins/Presisitence/rnaseq-plot/rgraph_toolkit/runner.py similarity index 100% rename from plugins/Presisitence/rgraph/rgraph_toolkit/runner.py rename to plugins/Presisitence/rnaseq-plot/rgraph_toolkit/runner.py diff --git a/plugins/Presisitence/rgraph/rscripts/_common.R b/plugins/Presisitence/rnaseq-plot/rscripts/_common.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/_common.R rename to plugins/Presisitence/rnaseq-plot/rscripts/_common.R diff --git a/plugins/Presisitence/rgraph/rscripts/build_gmt.R b/plugins/Presisitence/rnaseq-plot/rscripts/build_gmt.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/build_gmt.R rename to plugins/Presisitence/rnaseq-plot/rscripts/build_gmt.R diff --git a/plugins/Presisitence/rgraph/rscripts/correlation_heatmap.R b/plugins/Presisitence/rnaseq-plot/rscripts/correlation_heatmap.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/correlation_heatmap.R rename to plugins/Presisitence/rnaseq-plot/rscripts/correlation_heatmap.R diff --git a/plugins/Presisitence/rgraph/rscripts/diff.R b/plugins/Presisitence/rnaseq-plot/rscripts/diff.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/diff.R rename to plugins/Presisitence/rnaseq-plot/rscripts/diff.R diff --git a/plugins/Presisitence/rgraph/rscripts/distribution.R b/plugins/Presisitence/rnaseq-plot/rscripts/distribution.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/distribution.R rename to plugins/Presisitence/rnaseq-plot/rscripts/distribution.R diff --git a/plugins/Presisitence/rgraph/rscripts/enrich.R b/plugins/Presisitence/rnaseq-plot/rscripts/enrich.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/enrich.R rename to plugins/Presisitence/rnaseq-plot/rscripts/enrich.R diff --git a/plugins/Presisitence/rgraph/rscripts/enrich_circos.R b/plugins/Presisitence/rnaseq-plot/rscripts/enrich_circos.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/enrich_circos.R rename to plugins/Presisitence/rnaseq-plot/rscripts/enrich_circos.R diff --git a/plugins/Presisitence/rgraph/rscripts/gene_bar.R b/plugins/Presisitence/rnaseq-plot/rscripts/gene_bar.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/gene_bar.R rename to plugins/Presisitence/rnaseq-plot/rscripts/gene_bar.R diff --git a/plugins/Presisitence/rgraph/rscripts/gene_correlation.R b/plugins/Presisitence/rnaseq-plot/rscripts/gene_correlation.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/gene_correlation.R rename to plugins/Presisitence/rnaseq-plot/rscripts/gene_correlation.R diff --git a/plugins/Presisitence/rgraph/rscripts/go_plot.R b/plugins/Presisitence/rnaseq-plot/rscripts/go_plot.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/go_plot.R rename to plugins/Presisitence/rnaseq-plot/rscripts/go_plot.R diff --git a/plugins/Presisitence/rgraph/rscripts/gsea.R b/plugins/Presisitence/rnaseq-plot/rscripts/gsea.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/gsea.R rename to plugins/Presisitence/rnaseq-plot/rscripts/gsea.R diff --git a/plugins/Presisitence/rgraph/rscripts/gsea_mountain.R b/plugins/Presisitence/rnaseq-plot/rscripts/gsea_mountain.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/gsea_mountain.R rename to plugins/Presisitence/rnaseq-plot/rscripts/gsea_mountain.R diff --git a/plugins/Presisitence/rgraph/rscripts/heatmap.R b/plugins/Presisitence/rnaseq-plot/rscripts/heatmap.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/heatmap.R rename to plugins/Presisitence/rnaseq-plot/rscripts/heatmap.R diff --git a/plugins/Presisitence/rgraph/rscripts/interaction_diff.R b/plugins/Presisitence/rnaseq-plot/rscripts/interaction_diff.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/interaction_diff.R rename to plugins/Presisitence/rnaseq-plot/rscripts/interaction_diff.R diff --git a/plugins/Presisitence/rgraph/rscripts/interactive_volcano.R b/plugins/Presisitence/rnaseq-plot/rscripts/interactive_volcano.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/interactive_volcano.R rename to plugins/Presisitence/rnaseq-plot/rscripts/interactive_volcano.R diff --git a/plugins/Presisitence/rgraph/rscripts/kegg_plot.R b/plugins/Presisitence/rnaseq-plot/rscripts/kegg_plot.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/kegg_plot.R rename to plugins/Presisitence/rnaseq-plot/rscripts/kegg_plot.R diff --git a/plugins/Presisitence/rgraph/rscripts/kmeans.R b/plugins/Presisitence/rnaseq-plot/rscripts/kmeans.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/kmeans.R rename to plugins/Presisitence/rnaseq-plot/rscripts/kmeans.R diff --git a/plugins/Presisitence/rgraph/rscripts/multilevel_scatter.R b/plugins/Presisitence/rnaseq-plot/rscripts/multilevel_scatter.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/multilevel_scatter.R rename to plugins/Presisitence/rnaseq-plot/rscripts/multilevel_scatter.R diff --git a/plugins/Presisitence/rgraph/rscripts/network.R b/plugins/Presisitence/rnaseq-plot/rscripts/network.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/network.R rename to plugins/Presisitence/rnaseq-plot/rscripts/network.R diff --git a/plugins/Presisitence/rgraph/rscripts/norm_matrix.R b/plugins/Presisitence/rnaseq-plot/rscripts/norm_matrix.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/norm_matrix.R rename to plugins/Presisitence/rnaseq-plot/rscripts/norm_matrix.R diff --git a/plugins/Presisitence/rgraph/rscripts/normalize.R b/plugins/Presisitence/rnaseq-plot/rscripts/normalize.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/normalize.R rename to plugins/Presisitence/rnaseq-plot/rscripts/normalize.R diff --git a/plugins/Presisitence/rgraph/rscripts/pathway_gene_heatmap.R b/plugins/Presisitence/rnaseq-plot/rscripts/pathway_gene_heatmap.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/pathway_gene_heatmap.R rename to plugins/Presisitence/rnaseq-plot/rscripts/pathway_gene_heatmap.R diff --git a/plugins/Presisitence/rgraph/rscripts/pathway_trend.R b/plugins/Presisitence/rnaseq-plot/rscripts/pathway_trend.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/pathway_trend.R rename to plugins/Presisitence/rnaseq-plot/rscripts/pathway_trend.R diff --git a/plugins/Presisitence/rgraph/rscripts/pca.R b/plugins/Presisitence/rnaseq-plot/rscripts/pca.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/pca.R rename to plugins/Presisitence/rnaseq-plot/rscripts/pca.R diff --git a/plugins/Presisitence/rgraph/rscripts/ppi.R b/plugins/Presisitence/rnaseq-plot/rscripts/ppi.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/ppi.R rename to plugins/Presisitence/rnaseq-plot/rscripts/ppi.R diff --git a/plugins/Presisitence/rgraph/rscripts/quadrant.R b/plugins/Presisitence/rnaseq-plot/rscripts/quadrant.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/quadrant.R rename to plugins/Presisitence/rnaseq-plot/rscripts/quadrant.R diff --git a/plugins/Presisitence/rgraph/rscripts/sankey.R b/plugins/Presisitence/rnaseq-plot/rscripts/sankey.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/sankey.R rename to plugins/Presisitence/rnaseq-plot/rscripts/sankey.R diff --git a/plugins/Presisitence/rgraph/rscripts/segmented_volcano.R b/plugins/Presisitence/rnaseq-plot/rscripts/segmented_volcano.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/segmented_volcano.R rename to plugins/Presisitence/rnaseq-plot/rscripts/segmented_volcano.R diff --git a/plugins/Presisitence/rgraph/rscripts/ssgsea.R b/plugins/Presisitence/rnaseq-plot/rscripts/ssgsea.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/ssgsea.R rename to plugins/Presisitence/rnaseq-plot/rscripts/ssgsea.R diff --git a/plugins/Presisitence/rgraph/rscripts/stacked_volcano.R b/plugins/Presisitence/rnaseq-plot/rscripts/stacked_volcano.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/stacked_volcano.R rename to plugins/Presisitence/rnaseq-plot/rscripts/stacked_volcano.R diff --git a/plugins/Presisitence/rgraph/rscripts/venn.R b/plugins/Presisitence/rnaseq-plot/rscripts/venn.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/venn.R rename to plugins/Presisitence/rnaseq-plot/rscripts/venn.R diff --git a/plugins/Presisitence/rgraph/rscripts/virtual_knockout.R b/plugins/Presisitence/rnaseq-plot/rscripts/virtual_knockout.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/virtual_knockout.R rename to plugins/Presisitence/rnaseq-plot/rscripts/virtual_knockout.R diff --git a/plugins/Presisitence/rgraph/rscripts/volcano.R b/plugins/Presisitence/rnaseq-plot/rscripts/volcano.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/volcano.R rename to plugins/Presisitence/rnaseq-plot/rscripts/volcano.R diff --git a/plugins/Presisitence/rgraph/rscripts/wgcna.R b/plugins/Presisitence/rnaseq-plot/rscripts/wgcna.R similarity index 100% rename from plugins/Presisitence/rgraph/rscripts/wgcna.R rename to plugins/Presisitence/rnaseq-plot/rscripts/wgcna.R diff --git a/plugins/Presisitence/rgraph/server.py b/plugins/Presisitence/rnaseq-plot/server.py similarity index 98% rename from plugins/Presisitence/rgraph/server.py rename to plugins/Presisitence/rnaseq-plot/server.py index bea2bfa..b4cbced 100644 --- a/plugins/Presisitence/rgraph/server.py +++ b/plugins/Presisitence/rnaseq-plot/server.py @@ -1,8 +1,9 @@ -"""rgraph —— RNAseq 下游可视化/分析 MCP Server (FastMCP)。 +"""rnaseq-plot —— RNA-seq 下游出图 MCP Server (FastMCP)。 把一套清洗过、参数化的 R 出图脚本(rscripts/*.R)暴露为 MCP 工具,由本机 Rscript 引擎渲染, -产出与课程一致风格的 png+pdf。定位不到 R 时返回可手动运行的脚本与命令;R 端缺包时返回缺失 -包名与安装建议。 +产出 png+pdf。工具函数名保持 `rgraph_*`,与既有客户端配置兼容。 + +定位不到 R 时返回可手动运行的脚本与命令;R 端缺包时返回缺失包名与安装建议。 运行: uv run --directory server.py @@ -15,7 +16,7 @@ from rgraph_toolkit import runner -mcp = FastMCP("rgraph") +mcp = FastMCP("rnaseq-plot") def _clean(**kw: Any) -> dict[str, Any]: diff --git a/plugins/Presisitence/rgraph/skills/rgraph/SKILL.md b/plugins/Presisitence/rnaseq-plot/skills/rnaseq-plot/SKILL.md similarity index 64% rename from plugins/Presisitence/rgraph/skills/rgraph/SKILL.md rename to plugins/Presisitence/rnaseq-plot/skills/rnaseq-plot/SKILL.md index d6005e7..ce057cf 100644 --- a/plugins/Presisitence/rgraph/skills/rgraph/SKILL.md +++ b/plugins/Presisitence/rnaseq-plot/skills/rnaseq-plot/SKILL.md @@ -1,11 +1,11 @@ --- -name: rgraph -description: Use when the user already has an RNA-seq count or expression matrix and needs downstream analysis or figures (normalize, PCA, DESeq2/edgeR/limma, volcano, heatmap, GO/KEGG, GSEA, WGCNA). Drive the rgraph MCP tools with local Rscript; do not redraw in Python. +name: rnaseq-plot +description: Use when the user already has an RNA-seq count or expression matrix and needs downstream plots or analysis (normalize, PCA, DESeq2/edgeR/limma, volcano, heatmap, GO/KEGG, GSEA, WGCNA). Call the rnaseq-plot MCP tools named rgraph_*; do not redraw in Python. --- -# rgraph +# rnaseq-plot -Use the `rgraph` MCP tools. Prefer R-rendered png+pdf over matplotlib copies. +Use the `rnaseq-plot` MCP server. Tool functions are still named `rgraph_*`. Prefer R-rendered png+pdf over matplotlib copies. ## Typical order diff --git a/plugins/Presisitence/rgraph/tests/data/gene_count.csv b/plugins/Presisitence/rnaseq-plot/tests/data/gene_count.csv similarity index 100% rename from plugins/Presisitence/rgraph/tests/data/gene_count.csv rename to plugins/Presisitence/rnaseq-plot/tests/data/gene_count.csv diff --git a/plugins/Presisitence/rgraph/tests/data/sample_group.csv b/plugins/Presisitence/rnaseq-plot/tests/data/sample_group.csv similarity index 100% rename from plugins/Presisitence/rgraph/tests/data/sample_group.csv rename to plugins/Presisitence/rnaseq-plot/tests/data/sample_group.csv