diff --git a/plugins/Presisitence/microbe/LICENSE b/plugins/Presisitence/microbe/LICENSE new file mode 100644 index 0000000..9637979 --- /dev/null +++ b/plugins/Presisitence/microbe/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/microbe/README.md b/plugins/Presisitence/microbe/README.md new file mode 100644 index 0000000..564e051 --- /dev/null +++ b/plugins/Presisitence/microbe/README.md @@ -0,0 +1,39 @@ +# microbe + +Downstream **16S / ITS** microbiome analysis for MiniMax Code. vegan, DESeq2, edgeR, igraph, +Hmisc, randomForest, and ggtree scripts are wrapped as MCP tools and rendered by local **Rscript**. + +Figure types follow two public papers for layout (the Plugin does not include their data): + +- Liu et al. 2023, *Nature Microbiology* +- Zhou et al. 2022, *Nature Communications* + +Does not run DADA2/QIIME2. Start from a feature table. + +## Try it + +```text +I have feature_table.csv, taxonomy.csv, and metadata.csv (groups DP vs DSP). Draw alpha and beta +diversity, a genus stacked bar, DESeq2 differential abundance, and a co-occurrence network. +``` + +Expected result: the agent calls `microbe_env`, then `microbe_alpha`, `microbe_beta`, +`microbe_composition`, `microbe_diff`, and `microbe_network`. png+pdf paths and key statistics +are returned. Synthetic tables can be generated with `python tests/prep_test.py` for a dry run. + +## Requirements + +- Python 3.10+ and [uv](https://docs.astral.sh/uv/) on PATH. +- R with `Rscript` on PATH, or set `MICROBE_RSCRIPT`. +- R packages such as vegan, ggplot2, igraph; DESeq2/edgeR optional (tools degrade with install hints). +- Windows, macOS, and Linux. + +## Data and network + +- Analyses are local on user CSVs. No telemetry. +- No credentials in the package. +- `tests/prep_test.py` writes synthetic OTU-like tables only. + +## License + +MIT. See [LICENSE](LICENSE). diff --git a/plugins/Presisitence/microbe/mcp.json b/plugins/Presisitence/microbe/mcp.json new file mode 100644 index 0000000..8ae6a27 --- /dev/null +++ b/plugins/Presisitence/microbe/mcp.json @@ -0,0 +1,11 @@ +{ + "$schema": "https://agent-plugins.org/schemas/1.0.0/mcp.schema.json", + "mcpServers": { + "microbe": { + "type": "stdio", + "command": "uv", + "args": ["run", "server.py"], + "cwd": "${PLUGIN_ROOT}" + } + } +} diff --git a/plugins/Presisitence/microbe/microbe_toolkit/__init__.py b/plugins/Presisitence/microbe/microbe_toolkit/__init__.py new file mode 100644 index 0000000..10c4b28 --- /dev/null +++ b/plugins/Presisitence/microbe/microbe_toolkit/__init__.py @@ -0,0 +1,6 @@ +"""microbe_toolkit: 把一套下游微生物组/扩增子分析 R 脚本参数化封装为 MCP 工具。""" +from __future__ import annotations + +from . import runner + +__all__ = ["runner"] diff --git a/plugins/Presisitence/microbe/microbe_toolkit/cli.py b/plugins/Presisitence/microbe/microbe_toolkit/cli.py new file mode 100644 index 0000000..b1c0063 --- /dev/null +++ b/plugins/Presisitence/microbe/microbe_toolkit/cli.py @@ -0,0 +1,33 @@ +"""命令行自检入口：`microbe-cli` 检查 R 引擎与关键微生物组分析包是否就绪。""" +from __future__ import annotations + +import sys + +from . import runner + +_KEY_PKGS = [ + "vegan", "ggplot2", "ggpubr", "dplyr", "tidyr", "reshape2", "pheatmap", + "ggrepel", "igraph", "Hmisc", "randomForest", "ape", "scales", "circlize", + "ggalluvial", "edgeR", "DESeq2", "ggtree", "phyloseq", "SpiecEasi", +] + + +def main() -> int: + rscript = runner.find_rscript() + print(f"Rscript: {rscript or '未找到 (设置 MICROBE_RSCRIPT 或安装 R)'}") + print(f"mscripts 目录: {runner.mscripts_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/microbe/microbe_toolkit/runner.py b/plugins/Presisitence/microbe/microbe_toolkit/runner.py new file mode 100644 index 0000000..d0e76ab --- /dev/null +++ b/plugins/Presisitence/microbe/microbe_toolkit/runner.py @@ -0,0 +1,293 @@ +"""Rscript 执行引擎：定位 R、序列化参数、运行 mscripts/ 下的微生物组 R 模板并解析结果。 + +设计要点（与 rgraph-mcp 一致的架构，便于统一维护） +-------------------------------------------------------- +* **忠实 R 出图**：不在 Python 里重画，而是驱动一套清洗过的参数化 R 脚本（vegan/DESeq2/ + edgeR/igraph/Hmisc/randomForest/ggtree/ggplot2 等），产出出版级 png/pdf。 +* **零额外 R 依赖传参**：参数以「R 源文件」形式落盘（`params <- list(...)`），R 端 `source()` + 即可，不依赖 jsonlite 等包。 +* **优雅降级**：定位不到 Rscript → 返回可直接手动运行的脚本与命令；R 端缺包 → 返回缺失包名与 + 安装建议（区分 CRAN/Bioconductor），而不是抛出晦涩的 R 报错。 +* **可复现**：以 `--no-init-file --no-site-file` 运行，屏蔽用户 .Rprofile/Rprofile.site + 的自动加载横幅，避免污染 stdout 标记解析。 +""" +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：环境变量 MICROBE_RSCRIPT → PATH → 常见安装目录/注册表。""" + env = os.environ.get("MICROBE_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 + + 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) + + 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 mscripts_dir() -> Path: + """定位 mscripts 目录（源码布局为包的同级；wheel 布局为包内）。""" + here = Path(__file__).resolve().parent + for cand in (here.parent / "mscripts", here / "mscripts"): + if cand.is_dir(): + return cand + return here.parent / "mscripts" + + +# --------------------------------------------------------------------------- +# 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 = mscripts_dir() / script_name + common = mscripts_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 = 900.0, +) -> dict[str, Any]: + """运行一个 R 模板脚本，返回结构化结果。 + + 返回字段： + status: ok | missing_packages | r_not_found | error + outputs: 产物文件绝对路径列表（R 端通过 `MICROBE_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() + param_dir = Path(outdir) if outdir else Path(tempfile.mkdtemp(prefix="microbe_")) + params_file = param_dir / f".microbe_params_{script_name.replace('.R', '')}.R" + write_params_file(params, params_file) + + script_path = mscripts_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 = mscripts_dir() / "_common.R" + manual = f'"" {" ".join(_START_FLAGS)} "{script_path}" "{params_file}" "{common}"' + return { + "status": "r_not_found", + "message": ("未找到 Rscript。请安装 R 或设置环境变量 MICROBE_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] = [] + metrics: list[str] = [] + for line in combined.splitlines(): + line = line.strip() + if line.startswith("MICROBE_OUTPUT:"): + outputs.append(line[len("MICROBE_OUTPUT:"):].strip()) + elif line.startswith("MICROBE_MISSING_PACKAGES:"): + pk = line[len("MICROBE_MISSING_PACKAGES:"):].strip() + missing = [x.strip() for x in pk.replace(";", ",").split(",") if x.strip()] + elif line.startswith("MICROBE_METRIC:"): + metrics.append(line[len("MICROBE_METRIC:"):].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, + } + + result: dict[str, Any] = { + "status": "ok", + "outputs": outputs, + "n_outputs": len(outputs), + "script": str(script_path), + "params_file": str(params_file), + "log": log_tail, + } + if metrics: + result["metrics"] = metrics + return result + + +# --------------------------------------------------------------------------- +# 环境自检 / 安装建议 +# --------------------------------------------------------------------------- + +# CRAN vs Bioconductor 归属，用于给出正确的安装命令 +_BIOC = { + "DESeq2", "edgeR", "limma", "phyloseq", "microbiome", "ALDEx2", "metagenomeSeq", + "ANCOMBC", "ggtree", "treeio", "ComplexHeatmap", "DirichletMultinomial", + "BiocManager", "SummarizedExperiment", "Biobase", +} + + +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/microbe/mscripts/_common.R b/plugins/Presisitence/microbe/mscripts/_common.R new file mode 100644 index 0000000..1e470a9 --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/_common.R @@ -0,0 +1,318 @@ +# ============================================================================= +# _common.R — microbe-mcp 共享 R 库 +# 所有下游微生物组/扩增子分析模板都 source() 本文件，统一：参数读取、包预检、 +# 元数据/特征表/分类表读入与对齐、按分类层级汇总、相对丰度、抽平、出版级主题、 +# 配色、png+pdf 导出、产物/指标回传标记。 +# +# 与 rgraph-mcp 的 _common.R 保持同构（便于统一维护），仅面向微生物组数据。 +# * 不 setwd()/不 rm(list=ls())；一切走参数与绝对/相对入参。 +# * 缺包时输出 MICROBE_MISSING_PACKAGES 标记并以状态码 3 退出，交由上层友好提示。 +# * 产物统一通过 MICROBE_OUTPUT 标记回传；关键统计量通过 MICROBE_METRIC 回传。 +# ============================================================================= + +options(stringsAsFactors = FALSE, warn = 1) + +# ---- 参数读取 --------------------------------------------------------------- +microbe_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) +} + +microbe_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 + +# ---- 包预检 ----------------------------------------------------------------- +microbe_need <- function(pkgs) { + inst <- rownames(installed.packages()) + miss <- setdiff(pkgs, inst) + if (length(miss) > 0) { + cat("MICROBE_MISSING_PACKAGES:", paste(miss, collapse = ","), "\n") + quit(save = "no", status = 3) + } + invisible(TRUE) +} + +microbe_library <- function(pkgs) { + microbe_need(pkgs) + for (p in pkgs) suppressMessages(suppressWarnings(library(p, character.only = TRUE))) + invisible(TRUE) +} + +microbe_has <- function(pkg) requireNamespace(pkg, quietly = TRUE) + +# ---- 产物 / 指标回传 / 目录 ------------------------------------------------- +microbe_emit <- function(path) { + cat("MICROBE_OUTPUT:", normalizePath(path, winslash = "/", mustWork = FALSE), "\n") +} + +microbe_metric <- function(...) { + cat("MICROBE_METRIC:", paste0(...), "\n") +} + +microbe_ensure_dir <- function(path) { + if (!dir.exists(path)) dir.create(path, recursive = TRUE, showWarnings = FALSE) + invisible(path) +} + +# ---- 读入：样本元数据 (sample_name + group[+group_name][+环境/代谢物列]) ----- +microbe_read_metadata <- function(path) { + md <- read.csv(path, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + # 首列容错为 sample_name + cand <- c("sample_name", "SampleID", "sample", "Sample", "sample_id", "#SampleID") + hit <- intersect(cand, names(md)) + if (length(hit) == 0) names(md)[1] <- "sample_name" else if (hit[1] != "sample_name") { + names(md)[which(names(md) == hit[1])] <- "sample_name" + } + md <- md[!is.na(md$sample_name) & md$sample_name != "", , drop = FALSE] + has_g <- "group" %in% names(md) + has_gn <- "group_name" %in% names(md) + if (has_g && !has_gn) md$group_name <- md$group + if (has_gn && !has_g) md$group <- md$group_name + if (!has_g && !has_gn) { md$group <- md$sample_name; md$group_name <- md$sample_name } + rownames(md) <- md$sample_name + md +} + +# ---- 读入：特征表 (feature_id + 各样本计数列) ------------------------------ +microbe_read_feature <- function(path) { + df <- read.csv(path, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + fid <- c("feature_id", "OTU_ID", "OTUID", "#OTU ID", "ASV_ID", "ASV", "OTU", "ID", "taxon") + hit <- intersect(fid, names(df)) + if (length(hit) == 0) names(df)[1] <- "feature_id" else if (hit[1] != "feature_id") { + names(df)[which(names(df) == hit[1])] <- "feature_id" + } + df <- df[!duplicated(df$feature_id), , drop = FALSE] + df +} + +# 转为数值矩阵：行=feature_id，列=样本 +microbe_feature_matrix <- function(feature, samples = NULL) { + m <- feature + rownames(m) <- m$feature_id + m$feature_id <- NULL + if (!is.null(samples)) { + miss <- setdiff(samples, colnames(m)) + if (length(miss) > 0) stop(paste0("以下样本不在特征表中: ", paste(miss, collapse = ", "))) + m <- m[, samples, drop = FALSE] + } + m <- as.matrix(m) + suppressWarnings(storage.mode(m) <- "double") + m[is.na(m)] <- 0 + m +} + +# ---- 读入：分类表 (feature_id + 界门纲目科属种 或 单列 QIIME 分号串) --------- +.MICROBE_RANKS <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species") +.MICROBE_RANK_ALIAS <- list( + Kingdom = c("Kingdom", "Domain", "kingdom", "domain", "k"), + Phylum = c("Phylum", "phylum", "p"), + Class = c("Class", "class", "c"), + Order = c("Order", "order", "o"), + Family = c("Family", "family", "f"), + Genus = c("Genus", "genus", "g"), + Species = c("Species", "species", "s") +) + +microbe_read_taxonomy <- function(path) { + tx <- read.csv(path, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) + fid <- c("feature_id", "OTU_ID", "OTUID", "#OTU ID", "ASV_ID", "ASV", "OTU", "ID", "taxon") + hit <- intersect(fid, names(tx)) + if (length(hit) == 0) names(tx)[1] <- "feature_id" else if (hit[1] != "feature_id") { + names(tx)[which(names(tx) == hit[1])] <- "feature_id" + } + # 情况一：已含分列的分类等级 + present <- character(0) + for (r in .MICROBE_RANKS) { + al <- .MICROBE_RANK_ALIAS[[r]] + h <- intersect(al, names(tx)) + if (length(h) > 0 && h[1] != r) names(tx)[which(names(tx) == h[1])] <- r + if (r %in% names(tx)) present <- c(present, r) + } + if (length(present) >= 2) { + out <- tx[, c("feature_id", present), drop = FALSE] + return(.microbe_clean_tax(out, present)) + } + # 情况二：单列分号串 (k__..;p__..;c__..;o__..;f__..;g__..;s__..) + taxcol <- setdiff(names(tx), "feature_id") + taxcol <- taxcol[which.max(vapply(taxcol, function(cc) mean(grepl(";", tx[[cc]])), 0))] + parts <- strsplit(as.character(tx[[taxcol]]), "\\s*;\\s*") + mat <- matrix("", nrow = nrow(tx), ncol = 7, + dimnames = list(NULL, .MICROBE_RANKS)) + for (i in seq_along(parts)) { + v <- parts[[i]] + v <- sub("^[a-zA-Z]__", "", v) # 去 k__/p__ 前缀 + v <- trimws(v) + n <- min(length(v), 7) + if (n > 0) mat[i, seq_len(n)] <- v[seq_len(n)] + } + out <- data.frame(feature_id = tx$feature_id, mat, check.names = FALSE, stringsAsFactors = FALSE) + .microbe_clean_tax(out, .MICROBE_RANKS) +} + +.microbe_clean_tax <- function(out, ranks) { + for (r in ranks) { + x <- as.character(out[[r]]) + x[is.na(x) | x %in% c("", "NA", "unclassified", "Unclassified", "unidentified", "__")] <- NA + out[[r]] <- x + } + out +} + +# ---- 对齐特征表与元数据（取样本交集，保持元数据顺序） ---------------------- +microbe_align <- function(feature, meta) { + fsamp <- setdiff(names(feature), "feature_id") + common <- intersect(meta$sample_name, fsamp) + if (length(common) < 2) stop("特征表与元数据的公共样本不足 2 个，请检查样本名是否一致") + meta2 <- meta[meta$sample_name %in% common, , drop = FALSE] + list(samples = meta2$sample_name, meta = meta2) +} + +# ---- 按分类层级汇总：返回 行=taxa 列=样本 的计数矩阵 ----------------------- +microbe_collapse <- function(feature, taxonomy, level, samples) { + if (!level %in% names(taxonomy)) { + stop(paste0("分类表不含层级 '", level, "'，可用: ", + paste(setdiff(names(taxonomy), "feature_id"), collapse = ", "))) + } + m <- microbe_feature_matrix(feature, samples) + key <- taxonomy[[level]][match(rownames(m), taxonomy$feature_id)] + key[is.na(key) | key == ""] <- "Unclassified" + agg <- rowsum(m, group = key, reorder = TRUE) + agg +} + +# ---- 相对丰度（按列/样本归一到比例） --------------------------------------- +microbe_relabund <- function(mat) { + cs <- colSums(mat) + cs[cs == 0] <- 1 + sweep(mat, 2, cs, "/") +} + +# ---- 取 top-N 类群，其余合并为 Others（基于平均相对丰度排序） -------------- +microbe_top_taxa <- function(mat, top_n = 10, others_label = "Others") { + ra <- microbe_relabund(mat) + ord <- order(rowMeans(ra), decreasing = TRUE) + mat <- mat[ord, , drop = FALSE] + if (nrow(mat) <= top_n) return(mat) + top <- mat[seq_len(top_n), , drop = FALSE] + others <- colSums(mat[(top_n + 1):nrow(mat), , drop = FALSE]) + out <- rbind(top, others) + rownames(out)[nrow(out)] <- others_label + out +} + +# ---- 抽平（vegan::rrarefy；depth=NULL 时取最小样本深度） ------------------- +microbe_rarefy <- function(mat, depth = NULL, seed = 123) { + microbe_library(c("vegan")) + set.seed(seed) + cm <- t(round(mat)) + d <- if (is.null(depth)) min(rowSums(cm)) else depth + keep <- rowSums(cm) >= d + cm <- cm[keep, , drop = FALSE] + rr <- vegan::rrarefy(cm, sample = d) + list(mat = t(rr), depth = d, dropped = sum(!keep)) +} + +# ---- 出版级主题 ------------------------------------------------------------- +theme_microbe <- function(base_size = 12) { + ggplot2::theme_bw(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), + legend.key = ggplot2::element_blank(), + plot.title = ggplot2::element_text(hjust = 0.5, colour = "black") + ) +} + +# ---- 配色 ------------------------------------------------------------------- +# 分组配色（Okabe-Ito 色盲友好 + 课程红橙蓝） +.MICROBE_PAL_GROUP <- c("#E64B35", "#4DBBD5", "#00A087", "#3C5488", "#F39B7F", + "#8491B4", "#91D1C2", "#DC0000", "#7E6148", "#B09C85") +# 类群堆叠柱配色（20 色，末位灰给 Others/Unclassified） +.MICROBE_PAL_TAXA <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", + "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf", + "#aec7e8", "#ffbb78", "#98df8a", "#ff9896", "#c5b0d5", + "#c49c94", "#f7b6d2", "#dbdb8d", "#9edae5", "#393b79") + +microbe_palette <- function(n, name = "group") { + base <- switch(name, + taxa = .MICROBE_PAL_TAXA, + group = .MICROBE_PAL_GROUP, + .MICROBE_PAL_GROUP) + if (n <= length(base)) base[seq_len(n)] else grDevices::colorRampPalette(base)(n) +} + +# 给分组返回命名颜色向量 +microbe_group_colors <- function(groups, name = "group") { + lv <- levels(as.factor(groups)) + cols <- microbe_palette(length(lv), name) + stats::setNames(cols, lv) +} + +# 给类群返回命名颜色（Others/Unclassified 固定灰） +microbe_taxa_colors <- function(taxa) { + grey_lab <- intersect(c("Others", "Unclassified"), taxa) + main <- setdiff(taxa, grey_lab) + cols <- microbe_palette(length(main), "taxa") + cols <- stats::setNames(cols, main) + for (g in grey_lab) cols[g] <- "#bdbdbd" + cols[taxa] +} + +# 两两比较列表（供 ggpubr::stat_compare_means comparisons） +microbe_comparisons <- function(groups) { + lv <- levels(as.factor(groups)) + if (length(lv) < 2) return(list()) + cb <- utils::combn(lv, 2, simplify = FALSE) + cb +} + +# ---- 保存：ggplot 对象 → png+pdf 并回传 ------------------------------------ +microbe_save <- function(plot, outdir, base, width = 7, height = 6, dpi = 300, + formats = c("png", "pdf")) { + microbe_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", + limitsize = FALSE) + microbe_emit(f) + } + invisible(TRUE) +} + +# ---- 保存：base graphics (pheatmap/plot 等) → png+pdf ---------------------- +microbe_save_base <- function(outdir, base, draw_fn, width = 7, height = 6, dpi = 300, + formats = c("png", "pdf")) { + microbe_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() + microbe_emit(f) + } + invisible(TRUE) +} + +# ---- 写出 CSV 并回传 -------------------------------------------------------- +microbe_write_csv <- function(df, outdir, base, row.names = FALSE) { + microbe_ensure_dir(outdir) + f <- file.path(outdir, paste0(base, ".csv")) + write.csv(df, f, row.names = row.names) + microbe_emit(f) + invisible(f) +} diff --git a/plugins/Presisitence/microbe/mscripts/alpha.R b/plugins/Presisitence/microbe/mscripts/alpha.R new file mode 100644 index 0000000..79a74ec --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/alpha.R @@ -0,0 +1,83 @@ +# alpha.R —— α多样性指数箱线图 + 组间检验（对齐 Liu2023 Fig1a,b,d,e / Zhou2022 Fig1e,f） +# params: feature_table, metadata, outdir, +# metrics(默认 c("Observed","Chao1","Shannon","Simpson")) +# test("wilcox"|"t.test"|"anova"|"kruskal"，默认 wilcox) +# rarefy(bool,默认F), depth(抽平深度,缺省取最小样本), palette, +# pairwise(bool,默认T 画两两比较括号), dpi,width,height,name +# 输出: alpha_diversity.csv(各样本各指数) + alpha_boxplot.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("vegan", "ggplot2", "dplyr", "tidyr")) +use_pubr <- microbe_has("ggpubr") + +feature <- microbe_read_feature(p$feature_table) +meta <- microbe_read_metadata(p$metadata) +al <- microbe_align(feature, meta) +samples <- al$samples; meta <- al$meta + +m <- microbe_feature_matrix(feature, samples) # 行=feature 列=样本 +if (isTRUE(microbe_opt(p, "rarefy", FALSE))) { + rr <- microbe_rarefy(m, depth = microbe_opt(p, "depth", NULL)) + m <- rr$mat + samples <- colnames(m) + meta <- meta[meta$sample_name %in% samples, , drop = FALSE] + microbe_metric("rarefy depth=", rr$depth, " dropped_samples=", rr$dropped) +} + +comm <- t(round(m)) # vegan 要求 行=样本 +metrics <- microbe_opt(p, "metrics", c("Observed", "Chao1", "Shannon", "Simpson")) + +est <- vegan::estimateR(comm) # S.obs, S.chao1, se.chao1, S.ACE, se.ACE +idx <- data.frame( + sample_name = rownames(comm), + Observed = est["S.obs", ], + Chao1 = est["S.chao1", ], + ACE = est["S.ACE", ], + Shannon = vegan::diversity(comm, index = "shannon"), + Simpson = vegan::diversity(comm, index = "simpson"), + InvSimpson = vegan::diversity(comm, index = "invsimpson"), + Pielou = vegan::diversity(comm, index = "shannon") / log(vegan::specnumber(comm)), + check.names = FALSE, row.names = NULL +) +idx$group <- meta$group_name[match(idx$sample_name, meta$sample_name)] +microbe_write_csv(idx, p$outdir, "alpha_diversity") + +metrics <- intersect(metrics, names(idx)) +long <- tidyr::pivot_longer(idx[, c("sample_name", "group", metrics)], + cols = all_of(metrics), names_to = "metric", values_to = "value") +long$metric <- factor(long$metric, levels = metrics) +long$group <- factor(long$group, levels = unique(meta$group_name)) + +cols <- microbe_group_colors(long$group) +pl <- ggplot2::ggplot(long, ggplot2::aes(x = group, y = value, fill = group)) + + ggplot2::geom_boxplot(width = 0.6, outlier.shape = NA, alpha = 0.85) + + ggplot2::geom_jitter(width = 0.15, size = 1.2, alpha = 0.6, color = "black") + + ggplot2::scale_fill_manual(values = cols) + + ggplot2::facet_wrap(~metric, scales = "free_y", nrow = 1) + + ggplot2::labs(x = NULL, y = "Alpha diversity index", fill = "Group") + + theme_microbe() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 30, hjust = 1)) + +test <- microbe_opt(p, "test", "wilcox") +if (use_pubr && nlevels(long$group) >= 2) { + if (nlevels(long$group) == 2 || isTRUE(microbe_opt(p, "pairwise", TRUE))) { + pw_method <- if (identical(test, "t.test")) "t.test" else "wilcox.test" + pl <- pl + ggpubr::stat_compare_means(comparisons = microbe_comparisons(long$group), + method = pw_method, size = 2.8) + } + if (nlevels(long$group) > 2) { + glob <- if (identical(test, "anova")) "anova" else "kruskal.test" + pl <- pl + ggpubr::stat_compare_means(method = glob, label = "p.format", + size = 2.8, label.y.npc = "bottom") + } +} else if (!use_pubr) { + microbe_metric("ggpubr 未安装：箱线图正常输出，但未叠加显著性标注") +} + +n_metric <- length(metrics) +microbe_save(pl, p$outdir, microbe_opt(p, "name", "alpha_boxplot"), + width = microbe_opt(p, "width", max(3, n_metric * 2.2)), + height = microbe_opt(p, "height", 4.5), dpi = microbe_opt(p, "dpi", 300)) +cat("MICROBE_DONE: alpha samples=", length(samples), "groups=", nlevels(long$group), + "metrics=", paste(metrics, collapse = "/"), "\n") diff --git a/plugins/Presisitence/microbe/mscripts/beta.R b/plugins/Presisitence/microbe/mscripts/beta.R new file mode 100644 index 0000000..54ae7b6 --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/beta.R @@ -0,0 +1,97 @@ +# beta.R —— β多样性排序 (PCoA/NMDS) + PERMANOVA/ANOSIM（对齐 Liu2023 Fig1c,f / Zhou2022 Fig1c,d） +# params: feature_table, metadata, outdir, +# method("pcoa"|"nmds"，默认 pcoa), distance("bray"|"jaccard"|"euclidean"，默认 bray), +# test("permanova"|"anosim"|"both"，默认 permanova), permutations(默认 999), +# relabund(bool,默认T 先转相对丰度), ellipse(bool,默认T), label(bool,默认F), +# palette, dpi,width,height,name +# 输出: beta_scores.csv, beta_stats.csv, beta_.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("vegan", "ggplot2")) +use_repel <- microbe_has("ggrepel") + +feature <- microbe_read_feature(p$feature_table) +meta <- microbe_read_metadata(p$metadata) +al <- microbe_align(feature, meta) +samples <- al$samples; meta <- al$meta + +m <- microbe_feature_matrix(feature, samples) +if (isTRUE(microbe_opt(p, "relabund", TRUE))) m <- microbe_relabund(m) +comm <- t(m) # 行=样本 列=feature + +distance <- microbe_opt(p, "distance", "bray") +d <- vegan::vegdist(comm, method = distance) +grp <- factor(meta$group_name[match(rownames(comm), meta$sample_name)], + levels = unique(meta$group_name)) + +method <- microbe_opt(p, "method", "pcoa") +if (identical(method, "nmds")) { + set.seed(123) + ord <- vegan::metaMDS(d, k = 2, trymax = 50, trace = 0) + pts <- as.data.frame(ord$points) + names(pts)[1:2] <- c("Axis1", "Axis2") + xlab <- "NMDS1"; ylab <- "NMDS2" + subtitle <- sprintf("NMDS (%s) stress = %.3f", distance, ord$stress) +} else { + pc <- stats::cmdscale(d, k = 2, eig = TRUE) + pts <- as.data.frame(pc$points) + names(pts) <- c("Axis1", "Axis2") + eig <- pc$eig; pos <- eig[eig > 0] + ve <- 100 * eig[1:2] / sum(pos) + xlab <- sprintf("PCoA1 (%.2f%%)", ve[1]); ylab <- sprintf("PCoA2 (%.2f%%)", ve[2]) + subtitle <- sprintf("PCoA on %s distance", distance) +} +pts$sample_name <- rownames(pts) +pts$group <- grp + +# ---- 组间差异检验 ---------------------------------------------------------- +perm <- microbe_opt(p, "permutations", 999) +test <- microbe_opt(p, "test", "permanova") +stat_rows <- list() +ann <- character(0) +if (test %in% c("permanova", "both")) { + ad <- vegan::adonis2(d ~ grp, permutations = perm) + r2 <- ad$R2[1]; pv <- ad$`Pr(>F)`[1]; fval <- ad$F[1] + stat_rows[["PERMANOVA"]] <- data.frame(test = "PERMANOVA", statistic = fval, R2 = r2, p_value = pv) + ann <- c(ann, sprintf("PERMANOVA: R2 = %.3f, P = %.3f", r2, pv)) +} +if (test %in% c("anosim", "both")) { + an <- vegan::anosim(d, grp, permutations = perm) + stat_rows[["ANOSIM"]] <- data.frame(test = "ANOSIM", statistic = an$statistic, R2 = NA, p_value = an$signif) + ann <- c(ann, sprintf("ANOSIM: R = %.3f, P = %.3f", an$statistic, an$signif)) +} +stats_df <- do.call(rbind, stat_rows) +microbe_write_csv(stats_df, p$outdir, "beta_stats") +microbe_write_csv(pts[, c("sample_name", "group", "Axis1", "Axis2")], p$outdir, "beta_scores") +for (a in ann) microbe_metric(a) + +cols <- microbe_group_colors(grp) +pl <- ggplot2::ggplot(pts, ggplot2::aes(x = Axis1, y = Axis2, color = group)) + + ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "grey80") + + ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "grey80") + + ggplot2::geom_point(size = 3, alpha = 0.85) + + ggplot2::scale_color_manual(values = cols) + + ggplot2::labs(x = xlab, y = ylab, color = "Group", + title = subtitle, + subtitle = paste(ann, collapse = " ")) + + theme_microbe() + + ggplot2::theme(plot.subtitle = ggplot2::element_text(hjust = 0.5, size = 9)) + +if (isTRUE(microbe_opt(p, "ellipse", TRUE)) && nlevels(grp) > 1 && + all(table(grp) >= 3)) { + pl <- pl + ggplot2::stat_ellipse(ggplot2::aes(fill = group), geom = "polygon", + alpha = 0.12, level = 0.95, show.legend = FALSE) + + ggplot2::scale_fill_manual(values = cols) +} +if (isTRUE(microbe_opt(p, "label", FALSE))) { + if (use_repel) pl <- pl + ggrepel::geom_text_repel(ggplot2::aes(label = sample_name), + size = 2.6, color = "black", max.overlaps = Inf) + else pl <- pl + ggplot2::geom_text(ggplot2::aes(label = sample_name), size = 2.4, + color = "black", vjust = -0.7) +} + +microbe_save(pl, p$outdir, microbe_opt(p, "name", paste0("beta_", method)), + width = microbe_opt(p, "width", 6.5), height = microbe_opt(p, "height", 5.5), + dpi = microbe_opt(p, "dpi", 300)) +cat("MICROBE_DONE: beta", method, distance, "|", paste(ann, collapse = " | "), "\n") diff --git a/plugins/Presisitence/microbe/mscripts/composition.R b/plugins/Presisitence/microbe/mscripts/composition.R new file mode 100644 index 0000000..829d838 --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/composition.R @@ -0,0 +1,67 @@ +# composition.R —— 物种组成堆叠柱状图（对齐 Liu2023 Fig1g / Zhou2022 Fig4k,l, Fig5） +# params: feature_table, taxonomy, metadata, outdir, +# level("Phylum"|"Class"|"Order"|"Family"|"Genus"|"Species"，默认 Phylum), +# top_n(默认 10 其余合并 Others), mode("group"|"sample"，默认 group 组内均值), +# dpi,width,height,name +# 输出: composition_.csv(相对丰度表) + composition_.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("ggplot2", "reshape2")) + +feature <- microbe_read_feature(p$feature_table) +taxonomy <- microbe_read_taxonomy(p$taxonomy) +meta <- microbe_read_metadata(p$metadata) +al <- microbe_align(feature, meta) +samples <- al$samples; meta <- al$meta + +level <- microbe_opt(p, "level", "Phylum") +top_n <- microbe_opt(p, "top_n", 10) +mode <- microbe_opt(p, "mode", "group") + +agg <- microbe_collapse(feature, taxonomy, level, samples) # counts taxa x samples +agg <- microbe_top_taxa(agg, top_n) # top_n + Others +ra <- microbe_relabund(agg) # 每样本相对丰度 + +if (identical(mode, "group")) { + grp <- meta$group_name[match(colnames(ra), meta$sample_name)] + glv <- unique(meta$group_name) + gm <- sapply(glv, function(g) rowMeans(ra[, grp == g, drop = FALSE])) + colnames(gm) <- glv + ra <- microbe_relabund(gm) # 组内均值后再归一 + xlab <- "Group" +} else { + ra <- ra[, meta$sample_name[order(meta$group_name)], drop = FALSE] + xlab <- "Sample" +} + +tab <- data.frame(Taxon = rownames(ra), ra, check.names = FALSE) +microbe_write_csv(tab, p$outdir, paste0("composition_", level)) + +# 类群顺序：Others/Unclassified 置底，其余按总丰度降序 +taxa_order <- rownames(ra) +special <- intersect(c("Others", "Unclassified"), taxa_order) +main <- setdiff(taxa_order, special) +main <- main[order(rowSums(ra[main, , drop = FALSE]), decreasing = FALSE)] # 底->顶 +taxa_levels <- c(special, main) + +long <- reshape2::melt(as.matrix(ra), varnames = c("Taxon", "Unit"), value.name = "RelAbund") +long$Taxon <- factor(long$Taxon, levels = taxa_levels) +long$Unit <- factor(long$Unit, levels = colnames(ra)) + +cols <- microbe_taxa_colors(taxa_levels) +pl <- ggplot2::ggplot(long, ggplot2::aes(x = Unit, y = RelAbund, fill = Taxon)) + + ggplot2::geom_bar(stat = "identity", width = 0.75, color = "white", linewidth = 0.1) + + ggplot2::scale_fill_manual(values = cols, breaks = rev(taxa_levels)) + + ggplot2::scale_y_continuous(expand = c(0, 0), labels = scales::percent_format(accuracy = 1)) + + ggplot2::labs(x = xlab, y = "Relative abundance", fill = level) + + theme_microbe() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.key.size = ggplot2::unit(0.4, "cm")) + +n_unit <- ncol(ra) +microbe_save(pl, p$outdir, microbe_opt(p, "name", paste0("composition_", level)), + width = microbe_opt(p, "width", max(4.5, n_unit * 0.5 + 2.5)), + height = microbe_opt(p, "height", 5.5), dpi = microbe_opt(p, "dpi", 300)) +cat("MICROBE_DONE: composition level=", level, "units=", n_unit, "taxa=", nrow(ra), + "mode=", mode, "\n") diff --git a/plugins/Presisitence/microbe/mscripts/corr.R b/plugins/Presisitence/microbe/mscripts/corr.R new file mode 100644 index 0000000..3c883ef --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/corr.R @@ -0,0 +1,83 @@ +# corr.R —— 微生物-代谢物/环境因子 相关热图 (Spearman/Pearson, r 填色 + 显著性星号) +# 对齐 Liu2023 Fig2f (氨基酸/毒素与丰度的相关矩阵, r + P) +# params: feature_table, metadata, outdir, taxonomy(可选,配合 level 汇总), +# level(可选,如 "Genus"), variables(metadata 中的数值列名向量=代谢物/环境因子; +# 缺省自动取除分组外的全部数值列), features(可选,限定微生物类群向量), +# method("spearman"|"pearson"，默认 spearman), top_n(微生物类群上限,默认 30), +# padjust(bool,默认T BH), star(bool,默认T 显著性星号), dpi,width,height,name +# 输出: corr_r.csv, corr_p.csv, corr_heatmap.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("Hmisc", "ggplot2", "reshape2")) + +feature <- microbe_read_feature(p$feature_table) +meta <- microbe_read_metadata(p$metadata) +al <- microbe_align(feature, meta) +samples <- al$samples; meta <- al$meta + +# 微生物矩阵（层级汇总或 feature），相对丰度 +level <- microbe_opt(p, "level", NULL) +if (!is.null(level) && !is.null(p$taxonomy)) { + taxonomy <- microbe_read_taxonomy(p$taxonomy) + mat <- microbe_collapse(feature, taxonomy, level, samples) +} else { + mat <- microbe_feature_matrix(feature, samples) +} +if (!is.null(p$features)) mat <- mat[rownames(mat) %in% p$features, , drop = FALSE] +ra <- microbe_relabund(mat) +top_n <- microbe_opt(p, "top_n", 30) +if (nrow(ra) > top_n) ra <- ra[order(rowMeans(ra), decreasing = TRUE)[seq_len(top_n)], , drop = FALSE] + +# 环境/代谢物变量：指定列或自动取数值列 +reserved <- c("sample_name", "group", "group_name", "TvsC") +vars <- microbe_opt(p, "variables", NULL) +if (is.null(vars)) { + numcols <- names(meta)[vapply(meta, function(x) is.numeric(x) || !any(is.na(suppressWarnings(as.numeric(as.character(x))))), logical(1))] + vars <- setdiff(numcols, reserved) +} +vars <- intersect(vars, names(meta)) +if (length(vars) < 1) stop("未找到可用的数值型变量(代谢物/环境因子)，请通过 variables 指定 metadata 列名") +env <- sapply(vars, function(v) as.numeric(as.character(meta[[v]][match(colnames(ra), meta$sample_name)]))) +env <- as.matrix(env); rownames(env) <- colnames(ra) +env <- env[, apply(env, 2, function(x) stats::sd(x, na.rm = TRUE) > 0), drop = FALSE] +if (ncol(env) < 1) stop("数值变量方差为 0，无法计算相关") + +# ---- 相关：微生物(行) × 变量(列) ------------------------------------------ +method <- microbe_opt(p, "method", "spearman") +combined <- cbind(t(ra), env) +rc <- Hmisc::rcorr(as.matrix(combined), type = method) +taxa <- rownames(ra); vv <- colnames(env) +R <- rc$r[taxa, vv, drop = FALSE] +P <- rc$P[taxa, vv, drop = FALSE] +if (isTRUE(microbe_opt(p, "padjust", TRUE))) { + P[] <- p.adjust(as.vector(P), method = "BH") +} +microbe_write_csv(data.frame(taxon = rownames(R), R, check.names = FALSE), p$outdir, "corr_r") +microbe_write_csv(data.frame(taxon = rownames(P), P, check.names = FALSE), p$outdir, "corr_p") + +long <- reshape2::melt(R, varnames = c("taxon", "variable"), value.name = "r") +longp <- reshape2::melt(P, varnames = c("taxon", "variable"), value.name = "p") +long$p <- longp$p +long$star <- ifelse(long$p < 0.001, "***", ifelse(long$p < 0.01, "**", ifelse(long$p < 0.05, "*", ""))) +long$taxon <- factor(long$taxon, levels = rev(taxa)) +long$variable <- factor(long$variable, levels = vv) + +pl <- ggplot2::ggplot(long, ggplot2::aes(x = variable, y = taxon, fill = r)) + + ggplot2::geom_tile(color = "white", linewidth = 0.4) + + ggplot2::scale_fill_gradient2(low = "#4575b4", mid = "white", high = "#d73027", + midpoint = 0, limits = c(-1, 1), name = "r") + + ggplot2::labs(x = NULL, y = NULL, + title = sprintf("%s correlation", tools::toTitleCase(method))) + + theme_microbe() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + panel.border = ggplot2::element_blank(), axis.ticks = ggplot2::element_blank()) +if (isTRUE(microbe_opt(p, "star", TRUE))) { + pl <- pl + ggplot2::geom_text(ggplot2::aes(label = star), color = "black", size = 3) +} + +microbe_save(pl, p$outdir, microbe_opt(p, "name", "corr_heatmap"), + width = microbe_opt(p, "width", max(4, ncol(env) * 0.5 + 2.5)), + height = microbe_opt(p, "height", max(3.5, nrow(R) * 0.28 + 1.5)), + dpi = microbe_opt(p, "dpi", 300)) +cat("MICROBE_DONE: corr taxa=", nrow(R), "variables=", ncol(env), "method=", method, "\n") diff --git a/plugins/Presisitence/microbe/mscripts/diff.R b/plugins/Presisitence/microbe/mscripts/diff.R new file mode 100644 index 0000000..afdd51b --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/diff.R @@ -0,0 +1,140 @@ +# diff.R —— 差异丰度分析（DESeq2 / edgeR / Wilcoxon）（对齐 Liu2023 Fig3a,b / Zhou2022 EdgeR） +# params: feature_table, metadata, outdir, taxonomy(可选，配合 level 汇总), +# level(可选，如 "Genus"；给出则按该层级汇总后做差异，否则按 feature), +# method("auto"|"deseq2"|"edger"|"wilcox"，默认 auto: 有 DESeq2 用 deseq2 否则 edger), +# group_test/group_ref(指定两组；缺省用 metadata 的 TvsC=treatment/control，再缺省取前两组), +# padj(0.05), log2fc(1), min_count(平均计数下限,默认1), min_prev(最小检出样本数,默认2), +# kind("bar"|"volcano"，默认 bar), top_n(bar 展示上限,默认 30), dpi,width,height,name +# 输出: diff_result.csv + diff_.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("ggplot2", "dplyr")) + +feature <- microbe_read_feature(p$feature_table) +meta <- microbe_read_metadata(p$metadata) +al <- microbe_align(feature, meta) +samples <- al$samples; meta <- al$meta + +# ---- 计数矩阵：按层级汇总或按 feature ------------------------------------- +level <- microbe_opt(p, "level", NULL) +if (!is.null(level) && !is.null(p$taxonomy)) { + taxonomy <- microbe_read_taxonomy(p$taxonomy) + mat <- microbe_collapse(feature, taxonomy, level, samples) + unit_label <- level +} else { + mat <- microbe_feature_matrix(feature, samples) + unit_label <- "feature" +} +mat <- round(mat) + +# ---- 选两组：group_test/group_ref → TvsC → 前两组 -------------------------- +gname <- meta$group_name[match(colnames(mat), meta$sample_name)] +if (!is.null(p$group_test) && !is.null(p$group_ref)) { + g_test <- p$group_test; g_ref <- p$group_ref +} else if ("TvsC" %in% names(meta) && all(c("treatment", "control") %in% meta$TvsC)) { + g_test <- unique(meta$group_name[meta$TvsC == "treatment"])[1] + g_ref <- unique(meta$group_name[meta$TvsC == "control"])[1] +} else { + gs <- unique(gname) + if (length(gs) < 2) stop("分组不足两组，无法差异分析") + g_ref <- gs[1]; g_test <- gs[2] +} +keep_s <- gname %in% c(g_test, g_ref) +mat <- mat[, keep_s, drop = FALSE] +grp <- factor(gname[keep_s], levels = c(g_ref, g_test)) + +# ---- 过滤低丰度 ------------------------------------------------------------ +min_count <- microbe_opt(p, "min_count", 1) +min_prev <- microbe_opt(p, "min_prev", 2) +keep_f <- rowMeans(mat) >= min_count & rowSums(mat > 0) >= min_prev +mat <- mat[keep_f, , drop = FALSE] +if (nrow(mat) < 2) stop("过滤后可分析的类群不足 2 个，请放宽 min_count/min_prev") + +# ---- 方法选择 -------------------------------------------------------------- +method <- microbe_opt(p, "method", "auto") +if (identical(method, "auto")) method <- if (microbe_has("DESeq2")) "deseq2" else "edger" + +if (method == "deseq2") { + microbe_library(c("DESeq2")) + coldata <- data.frame(row.names = colnames(mat), group = grp) + dds <- DESeq2::DESeqDataSetFromMatrix(mat, coldata, design = ~group) + dds <- DESeq2::estimateSizeFactors(dds, type = "poscounts") # 微生物组多零，用 poscounts + dds <- DESeq2::DESeq(dds, fitType = "local", quiet = TRUE) + res <- as.data.frame(DESeq2::results(dds, contrast = c("group", g_test, g_ref))) + stats <- data.frame(taxon = rownames(res), baseMean = res$baseMean, + log2FoldChange = res$log2FoldChange, pvalue = res$pvalue, padj = res$padj) +} else if (method == "edger") { + microbe_library(c("edgeR")) + dge <- edgeR::DGEList(counts = mat, group = grp) + dge <- edgeR::calcNormFactors(dge, method = "TMM") + dge <- edgeR::estimateDisp(dge) + et <- edgeR::exactTest(dge, pair = c(g_ref, g_test)) # logFC = test/ref + tt <- edgeR::topTags(et, n = Inf, adjust.method = "BH", sort.by = "none")$table + stats <- data.frame(taxon = rownames(tt), baseMean = 2^tt$logCPM, + log2FoldChange = tt$logFC, pvalue = tt$PValue, padj = tt$FDR) +} else if (method == "wilcox") { + ra <- microbe_relabund(mat) + is_t <- grp == g_test + pv <- apply(ra, 1, function(x) tryCatch(stats::wilcox.test(x[is_t], x[!is_t])$p.value, error = function(e) NA)) + mt <- rowMeans(ra[, is_t, drop = FALSE]); mr <- rowMeans(ra[, !is_t, drop = FALSE]) + lfc <- log2((mt + 1e-6) / (mr + 1e-6)) + stats <- data.frame(taxon = rownames(ra), baseMean = rowMeans(ra), + log2FoldChange = lfc, pvalue = pv, padj = p.adjust(pv, "BH")) +} else { + stop(paste0("未知 method: ", method)) +} + +padj_cut <- microbe_opt(p, "padj", 0.05) +log2fc <- microbe_opt(p, "log2fc", 1) +stats <- stats[!is.na(stats$log2FoldChange), ] +stats$padj[is.na(stats$padj)] <- 1 +stats$sig <- stats$padj < padj_cut & abs(stats$log2FoldChange) > log2fc +stats$enriched_in <- ifelse(stats$log2FoldChange > 0, g_test, g_ref) +stats$enriched_in[!stats$sig] <- "n.s." +stats <- stats[order(stats$padj, -abs(stats$log2FoldChange)), ] +microbe_write_csv(stats, p$outdir, "diff_result") +microbe_metric(sprintf("diff method=%s test=%s ref=%s sig=%d (up_in_%s=%d, up_in_%s=%d)", + method, g_test, g_ref, sum(stats$sig), g_test, + sum(stats$sig & stats$log2FoldChange > 0), g_ref, + sum(stats$sig & stats$log2FoldChange < 0))) + +kind <- microbe_opt(p, "kind", "bar") +gcols <- stats::setNames(microbe_palette(2, "group"), c(g_ref, g_test)) +gcols["n.s."] <- "grey75" + +if (identical(kind, "volcano")) { + stats$neglogp <- -log10(pmax(stats$padj, 1e-300)) + pl <- ggplot2::ggplot(stats, ggplot2::aes(x = log2FoldChange, y = neglogp, color = enriched_in)) + + ggplot2::geom_point(alpha = 0.75, size = 1.8) + + ggplot2::scale_color_manual(values = gcols) + + ggplot2::geom_hline(yintercept = -log10(padj_cut), linetype = "dashed", color = "grey60") + + ggplot2::geom_vline(xintercept = c(-log2fc, log2fc), linetype = "dashed", color = "grey60") + + ggplot2::labs(x = sprintf("log2 fold change (%s / %s)", g_test, g_ref), + y = "-log10(FDR)", color = "Enriched in") + + theme_microbe() + w <- microbe_opt(p, "width", 7); h <- microbe_opt(p, "height", 6) +} else { + sig <- stats[stats$sig, , drop = FALSE] + if (nrow(sig) == 0) { + sig <- utils::head(stats[order(stats$padj), ], min(10, nrow(stats))) + microbe_metric("无显著差异类群，bar 图展示 padj 最小的若干条以供检视") + } + sig <- utils::head(sig[order(-abs(sig$log2FoldChange)), ], microbe_opt(p, "top_n", 30)) + sig <- sig[order(sig$log2FoldChange), ] + sig$taxon <- factor(sig$taxon, levels = sig$taxon) + pl <- ggplot2::ggplot(sig, ggplot2::aes(x = log2FoldChange, y = taxon, fill = enriched_in)) + + ggplot2::geom_col(width = 0.7) + + ggplot2::scale_fill_manual(values = gcols) + + ggplot2::geom_vline(xintercept = 0, color = "black", linewidth = 0.3) + + ggplot2::labs(x = sprintf("log2 fold change (%s / %s)", g_test, g_ref), + y = unit_label, fill = "Enriched in") + + theme_microbe() + w <- microbe_opt(p, "width", 7) + h <- microbe_opt(p, "height", max(3, nrow(sig) * 0.22 + 1.5)) +} + +microbe_save(pl, p$outdir, microbe_opt(p, "name", paste0("diff_", kind)), + width = w, height = h, dpi = microbe_opt(p, "dpi", 300)) +cat("MICROBE_DONE: diff method=", method, "unit=", unit_label, "n=", nrow(stats), + "sig=", sum(stats$sig), "\n") diff --git a/plugins/Presisitence/microbe/mscripts/network.R b/plugins/Presisitence/microbe/mscripts/network.R new file mode 100644 index 0000000..63cb026 --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/network.R @@ -0,0 +1,145 @@ +# network.R —— 微生物共现网络 (Spearman/Hmisc + igraph 拓扑 + keystone) +# 对齐 Liu2023 Fig1i (Spearman ρ>0.7 & P<0.05, BH校正; 正/负边着色; 点大小=degree) +# 与 Zhou2022 Fig2a-f (网络拓扑: 节点/边/度/平均路径长度/模块度; keystone 驱动菌) +# params: feature_table, metadata(可选,配合 group 只取某组样本), taxonomy(可选,按门着色), +# outdir, level(可选,如 "Genus" 先汇总), group(可选,只用该组样本建网), +# method("spearman"|"pearson"，默认 spearman), r_threshold(0.7), p_threshold(0.05), +# padjust(bool,默认T BH), min_prev(最小检出样本数,默认样本数的1/3), top_n(限最丰 N 类群,默认150), +# layout("fr"|"circle"，默认 fr), label(bool,默认F 仅标 keystone), dpi,width,height,name +# 输出: network_edges.csv, network_nodes.csv, network_topology.csv, network.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("Hmisc", "igraph")) + +feature <- microbe_read_feature(p$feature_table) +fsamp <- setdiff(names(feature), "feature_id") + +# 可选：按 metadata 的某组取样本子集 +if (!is.null(p$metadata)) { + meta <- microbe_read_metadata(p$metadata) + fsamp <- intersect(meta$sample_name, fsamp) + if (!is.null(p$group)) { + sel <- meta$sample_name[meta$group_name == p$group] + fsamp <- intersect(sel, fsamp) + } +} +if (length(fsamp) < 4) stop("建网样本数不足（<4），Spearman 相关不稳健") + +# 汇总层级或按 feature +level <- microbe_opt(p, "level", NULL) +if (!is.null(level) && !is.null(p$taxonomy)) { + taxonomy <- microbe_read_taxonomy(p$taxonomy) + mat <- microbe_collapse(feature, taxonomy, level, fsamp) +} else { + taxonomy <- if (!is.null(p$taxonomy)) microbe_read_taxonomy(p$taxonomy) else NULL + mat <- microbe_feature_matrix(feature, fsamp) +} + +# 过滤：检出率 + 取最丰 top_n +min_prev <- microbe_opt(p, "min_prev", max(2, floor(ncol(mat) / 3))) +mat <- mat[rowSums(mat > 0) >= min_prev, , drop = FALSE] +top_n <- microbe_opt(p, "top_n", 150) +if (nrow(mat) > top_n) mat <- mat[order(rowMeans(microbe_relabund(mat)), decreasing = TRUE)[seq_len(top_n)], , drop = FALSE] +if (nrow(mat) < 3) stop("过滤后类群不足 3 个，无法建网") + +# ---- 相关 (Hmisc::rcorr) + 阈值筛边 --------------------------------------- +method <- microbe_opt(p, "method", "spearman") +rc <- Hmisc::rcorr(t(mat), type = method) +r <- rc$r; P <- rc$P +diag(P) <- 1 +if (isTRUE(microbe_opt(p, "padjust", TRUE))) { + ut <- upper.tri(P) + padj <- P; padj[ut] <- p.adjust(P[ut], method = "BH") + padj[lower.tri(padj)] <- t(padj)[lower.tri(padj)] + P <- padj +} +r_thr <- microbe_opt(p, "r_threshold", 0.7) +p_thr <- microbe_opt(p, "p_threshold", 0.05) + +ut <- which(upper.tri(r) & abs(r) >= r_thr & P < p_thr, arr.ind = TRUE) +taxa <- rownames(mat) +edges <- data.frame( + from = taxa[ut[, 1]], to = taxa[ut[, 2]], + r = r[ut], p = P[ut], + sign = ifelse(r[ut] > 0, "positive", "negative"), + stringsAsFactors = FALSE +) +microbe_write_csv(edges, p$outdir, "network_edges") + +# ---- 建图 + 拓扑 ----------------------------------------------------------- +g <- igraph::graph_from_data_frame(edges[, c("from", "to")], directed = FALSE, + vertices = data.frame(name = taxa)) +igraph::E(g)$weight <- abs(edges$r) +igraph::E(g)$sign <- edges$sign + +# 社区/模块度（Louvain，权重=|r|） +comm <- tryCatch(igraph::cluster_louvain(g, weights = igraph::E(g)$weight), + error = function(e) igraph::cluster_fast_greedy(igraph::simplify(g))) +mem <- igraph::membership(comm) +modu <- igraph::modularity(comm) +deg <- igraph::degree(g) +clo <- suppressWarnings(igraph::closeness(g, normalized = TRUE)); clo[is.na(clo)] <- 0 +betw <- igraph::betweenness(g, normalized = TRUE) + +# keystone：度与接近中心性均居前的枢纽节点 +hub <- (deg >= stats::quantile(deg, 0.90)) & (clo >= stats::median(clo)) & (deg > 0) +nodes <- data.frame(taxon = taxa, module = as.integer(mem[taxa]), + degree = deg[taxa], closeness = round(clo[taxa], 4), + betweenness = round(betw[taxa], 4), keystone = hub[taxa], + stringsAsFactors = FALSE) +if (!is.null(taxonomy)) { + ann_level <- if (!is.null(level)) NULL else "Phylum" + if (!is.null(ann_level) && ann_level %in% names(taxonomy)) { + nodes$Phylum <- taxonomy$Phylum[match(nodes$taxon, taxonomy$feature_id)] + } +} +nodes <- nodes[order(-nodes$degree), ] +microbe_write_csv(nodes, p$outdir, "network_nodes") + +n_pos <- sum(edges$sign == "positive"); n_neg <- sum(edges$sign == "negative") +topo <- data.frame( + metric = c("nodes", "edges", "positive_edges", "negative_edges", "positive_pct", + "avg_degree", "avg_path_length", "clustering_coef", "modularity", + "n_modules", "n_keystone"), + value = c(igraph::vcount(g), igraph::ecount(g), n_pos, n_neg, + round(100 * n_pos / max(1, igraph::ecount(g)), 2), + round(mean(deg), 3), + round(suppressWarnings(igraph::mean_distance(g)), 3), + round(igraph::transitivity(g, type = "global"), 3), + round(modu, 3), length(unique(mem)), sum(hub)) +) +microbe_write_csv(topo, p$outdir, "network_topology") +microbe_metric(sprintf("network nodes=%d edges=%d pos%%=%.1f avg_degree=%.2f modularity=%.3f keystone=%d", + igraph::vcount(g), igraph::ecount(g), + 100 * n_pos / max(1, igraph::ecount(g)), mean(deg), modu, sum(hub))) + +# ---- 绘图 (igraph base graphics) ------------------------------------------- +lay_name <- microbe_opt(p, "layout", "fr") +set.seed(123) +lay <- if (identical(lay_name, "circle")) igraph::layout_in_circle(g) else igraph::layout_with_fr(g) +mcols <- microbe_palette(length(unique(mem)), "taxa") +vcol <- mcols[as.integer(factor(mem[taxa]))] +ecol <- ifelse(edges$sign == "positive", grDevices::adjustcolor("#d6604d", 0.55), + grDevices::adjustcolor("#4393c3", 0.55)) +vsize <- 3 + 7 * (deg[taxa] / max(1, max(deg))) +show_label <- isTRUE(microbe_opt(p, "label", FALSE)) +vlabel <- if (show_label) taxa else ifelse(hub[taxa], taxa, NA) + +microbe_save_base(p$outdir, microbe_opt(p, "name", "network"), function() { + graphics::par(mar = c(1, 1, 2, 1)) + igraph::plot.igraph(g, layout = lay, vertex.color = vcol, vertex.frame.color = "grey40", + vertex.size = vsize, vertex.label = vlabel, vertex.label.cex = 0.6, + vertex.label.color = "black", vertex.label.dist = 0.6, + edge.color = ecol, edge.width = 1.2 * abs(edges$r)) + graphics::title(sprintf("Co-occurrence network (%s |r|>=%.2f, P<%.2g)\nnodes=%d edges=%d +%.0f%%/-%.0f%% modularity=%.2f", + method, r_thr, p_thr, igraph::vcount(g), igraph::ecount(g), + 100 * n_pos / max(1, igraph::ecount(g)), + 100 * n_neg / max(1, igraph::ecount(g)), modu), cex.main = 0.8) + graphics::legend("bottomleft", legend = c("positive", "negative"), + col = c("#d6604d", "#4393c3"), lty = 1, lwd = 2, bty = "n", cex = 0.7) +}, width = microbe_opt(p, "width", 7.5), height = microbe_opt(p, "height", 7), + dpi = microbe_opt(p, "dpi", 300)) + +cat("MICROBE_DONE: network nodes=", igraph::vcount(g), "edges=", igraph::ecount(g), + "modularity=", round(modu, 3), "keystone=", sum(hub), "\n") diff --git a/plugins/Presisitence/microbe/mscripts/rf.R b/plugins/Presisitence/microbe/mscripts/rf.R new file mode 100644 index 0000000..4312cd4 --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/rf.R @@ -0,0 +1,94 @@ +# rf.R —— 随机森林分类 biomarker + 十折交叉验证误差曲线（对齐 Zhou2022 Fig2g-j） +# params: feature_table, metadata, outdir, taxonomy(可选,配合 level 汇总), +# level(可选,如 "Genus"), top_n(重要性展示条数,默认 20), ntree(默认 1000), +# cv(bool,默认T 做 10 折 rfcv), cv_fold(默认10), seed(默认123), dpi,width,height,name +# 输出: rf_importance.csv, rf_cv.csv(有CV时), rf_importance.png/pdf, rf_cv_error.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("randomForest", "ggplot2")) + +feature <- microbe_read_feature(p$feature_table) +meta <- microbe_read_metadata(p$metadata) +al <- microbe_align(feature, meta) +samples <- al$samples; meta <- al$meta + +level <- microbe_opt(p, "level", NULL) +if (!is.null(level) && !is.null(p$taxonomy)) { + taxonomy <- microbe_read_taxonomy(p$taxonomy) + mat <- microbe_collapse(feature, taxonomy, level, samples) + unit_label <- level +} else { + mat <- microbe_feature_matrix(feature, samples) + unit_label <- "feature" +} +ra <- microbe_relabund(mat) # 相对丰度，行=taxa 列=样本 +X <- t(ra) # 样本 x 特征 +grp <- factor(meta$group_name[match(rownames(X), meta$sample_name)], + levels = unique(meta$group_name)) + +seed <- microbe_opt(p, "seed", 123) +ntree <- microbe_opt(p, "ntree", 1000) +set.seed(seed) +rf <- randomForest::randomForest(x = X, y = grp, ntree = ntree, importance = TRUE) +oob <- rf$err.rate[nrow(rf$err.rate), "OOB"] +microbe_metric(sprintf("random_forest classes=%d features=%d ntree=%d OOB_error=%.3f", + nlevels(grp), ncol(X), ntree, oob)) + +imp <- randomForest::importance(rf) +gini <- if ("MeanDecreaseGini" %in% colnames(imp)) imp[, "MeanDecreaseGini"] else imp[, ncol(imp)] +# 每个特征富集于哪一组（均值相对丰度最大的组） +enr <- apply(ra, 1, function(x) { + mv <- tapply(x, grp[match(colnames(ra), rownames(X))], mean) + names(mv)[which.max(mv)] +}) +impdf <- data.frame(taxon = rownames(imp), MeanDecreaseGini = as.numeric(gini), + enriched_in = enr[rownames(imp)], stringsAsFactors = FALSE) +impdf <- impdf[order(-impdf$MeanDecreaseGini), ] +microbe_write_csv(impdf, p$outdir, "rf_importance") + +top_n <- min(microbe_opt(p, "top_n", 20), nrow(impdf)) +show <- impdf[seq_len(top_n), ] +show$taxon <- factor(show$taxon, levels = rev(show$taxon)) +gcols <- microbe_group_colors(grp) +pl <- ggplot2::ggplot(show, ggplot2::aes(x = MeanDecreaseGini, y = taxon, color = enriched_in)) + + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = MeanDecreaseGini, yend = taxon), + color = "grey70", linewidth = 0.5) + + ggplot2::geom_point(size = 3) + + ggplot2::scale_color_manual(values = gcols) + + ggplot2::labs(x = "Mean Decrease Gini (importance)", y = unit_label, + color = "Enriched in", + title = sprintf("Random-forest biomarkers (OOB error = %.1f%%)", 100 * oob)) + + theme_microbe() +microbe_save(pl, p$outdir, microbe_opt(p, "name", "rf_importance"), + width = microbe_opt(p, "width", 7), + height = microbe_opt(p, "height", max(3, top_n * 0.28 + 1.5)), + dpi = microbe_opt(p, "dpi", 300)) + +# ---- 十折交叉验证误差曲线（确定最优 biomarker 数）------------------------- +if (isTRUE(microbe_opt(p, "cv", TRUE)) && ncol(X) >= 4) { + set.seed(seed) + cvres <- randomForest::rfcv(trainx = X, trainy = grp, + cv.fold = microbe_opt(p, "cv_fold", 10)) + cvdf <- data.frame(n_features = cvres$n.var, cv_error = cvres$error.cv) + cvdf <- cvdf[order(cvdf$n_features), ] + microbe_write_csv(cvdf, p$outdir, "rf_cv") + optn <- cvdf$n_features[which.min(cvdf$cv_error)] + microbe_metric(sprintf("rfcv optimal_n_features=%d min_cv_error=%.3f", optn, min(cvdf$cv_error))) + plcv <- ggplot2::ggplot(cvdf, ggplot2::aes(x = n_features, y = cv_error)) + + ggplot2::geom_line(color = "#3C5488", linewidth = 0.7) + + ggplot2::geom_point(size = 2, color = "#3C5488") + + ggplot2::geom_vline(xintercept = optn, linetype = "dashed", color = "#E64B35") + + ggplot2::scale_x_log10() + + ggplot2::annotate("text", x = optn, y = max(cvdf$cv_error), + label = sprintf("optimal = %d", optn), color = "#E64B35", + hjust = -0.1, size = 3) + + ggplot2::labs(x = "Number of features (log10)", y = "Cross-validation error rate", + title = "Ten-fold cross-validation") + + theme_microbe() + microbe_save(plcv, p$outdir, "rf_cv_error", + width = microbe_opt(p, "width", 6), height = 4.5, + dpi = microbe_opt(p, "dpi", 300)) +} + +cat("MICROBE_DONE: rf unit=", unit_label, "features=", ncol(X), "OOB=", round(oob, 3), "\n") diff --git a/plugins/Presisitence/microbe/mscripts/tree.R b/plugins/Presisitence/microbe/mscripts/tree.R new file mode 100644 index 0000000..700624a --- /dev/null +++ b/plugins/Presisitence/microbe/mscripts/tree.R @@ -0,0 +1,99 @@ +# tree.R —— 分类树/系统发育树 + 丰度环（对齐 Zhou2022 Fig3a,b: 树 + 相对丰度热环 + 类群着色） +# 无测序数据时以「分类学层级」构树(ape::as.phylo formula)；失败则退化为按丰度谱聚类构树。 +# params: feature_table, taxonomy, metadata, outdir, +# level(树尖层级,默认 "Genus"), color_by(按此分类等级给树尖着色,默认 "Phylum"), +# layout("circular"|"fan"|"rectangular"，默认 circular), ring(bool,默认T 加组均值丰度环), +# top_n(树尖上限,默认 60), dpi,width,height,name +# 输出: tree_tip_abundance.csv + tree.png/pdf +.args <- commandArgs(trailingOnly = TRUE) +source(.args[[2]]) +p <- microbe_load_params() +microbe_library(c("ape", "ggtree", "ggplot2")) + +feature <- microbe_read_feature(p$feature_table) +taxonomy <- microbe_read_taxonomy(p$taxonomy) +meta <- microbe_read_metadata(p$metadata) +al <- microbe_align(feature, meta) +samples <- al$samples; meta <- al$meta + +level <- microbe_opt(p, "level", "Genus") +color_by <- microbe_opt(p, "color_by", "Phylum") +top_n <- microbe_opt(p, "top_n", 60) + +agg <- microbe_collapse(feature, taxonomy, level, samples) # counts taxa x samples +ra <- microbe_relabund(agg) +tips <- rownames(ra)[order(rowMeans(ra), decreasing = TRUE)] +tips <- setdiff(tips, c("Unclassified")) +tips <- utils::head(tips, top_n) +if (length(tips) < 3) stop("可用树尖类群不足 3 个") +ra <- ra[tips, , drop = FALSE] + +# 树尖对应的高阶分类（每个 tip 取首个匹配 feature 的分类行） +ranks_all <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species") +idx <- match(level, ranks_all) +use_ranks <- intersect(ranks_all[seq_len(idx)], names(taxonomy)) +taxrow <- function(tip, r) { + hit <- which(taxonomy[[level]] == tip) + if (length(hit) == 0) return(NA) + v <- taxonomy[[r]][hit[1]] + if (is.na(v) || v == "") paste0("unclassified_", r) else v +} +taxcols <- lapply(use_ranks, function(r) unname(vapply(tips, taxrow, character(1), r = r))) +names(taxcols) <- use_ranks +taxdf <- as.data.frame(taxcols, stringsAsFactors = FALSE, check.names = FALSE) +taxdf[[level]] <- make.unique(as.character(tips)) # 去重 tip 标签 +tips <- taxdf[[level]] +rownames(ra) <- tips +for (r in names(taxdf)) taxdf[[r]] <- factor(taxdf[[r]]) + +# ---- 构树：分类学 formula 树，失败则退化为丰度谱聚类树 --------------------- +tr <- NULL +if (length(use_ranks) >= 2) { + frm <- stats::as.formula(paste("~", paste(use_ranks, collapse = "/"))) + tr <- tryCatch(ape::as.phylo(frm, data = taxdf, collapse.singles = FALSE), + error = function(e) NULL) +} +if (is.null(tr)) { + microbe_metric("分类学构树失败，退化为按相对丰度谱 (1-Pearson) 层次聚类构树") + dd <- stats::as.dist(1 - stats::cor(t(ra))) + tr <- ape::as.phylo(stats::hclust(dd, method = "average")) +} + +# 树尖着色注释 +tipmeta <- data.frame(label = tips, + color_group = if (color_by %in% names(taxonomy)) { + vapply(rownames(ra), function(t) { + hit <- which(taxonomy[[level]] == sub("\\.[0-9]+$", "", t)) + v <- if (length(hit)) taxonomy[[color_by]][hit[1]] else NA + if (is.na(v) || v == "") "Unclassified" else v + }, character(1)) + } else "NA", + stringsAsFactors = FALSE) + +lay <- microbe_opt(p, "layout", "circular") +gt <- ggtree::ggtree(tr, layout = lay, linewidth = 0.4) %<+% tipmeta + + ggtree::geom_tippoint(ggplot2::aes(color = color_group), size = 1.8) + + ggplot2::scale_color_manual(values = microbe_palette(length(unique(tipmeta$color_group)), "taxa"), + name = color_by) + + ggtree::geom_tiplab(size = 1.8, offset = 0.02, align = TRUE) + +# ---- 组均值相对丰度环 ------------------------------------------------------ +if (isTRUE(microbe_opt(p, "ring", TRUE))) { + grp <- meta$group_name[match(colnames(ra), meta$sample_name)] + glv <- unique(meta$group_name) + gm <- sapply(glv, function(g) rowMeans(ra[, grp == g, drop = FALSE])) + gm <- as.data.frame(gm); rownames(gm) <- tips; colnames(gm) <- glv + microbe_write_csv(data.frame(tip = rownames(gm), gm, check.names = FALSE), + p$outdir, "tree_tip_abundance") + gt <- ggtree::gheatmap(gt, gm, offset = 0.15, width = 0.25, font.size = 2, + colnames_angle = 90, colnames_offset_y = 0) + + ggplot2::scale_fill_gradient(low = "white", high = "#238b45", name = "Rel. abund.") +} else { + microbe_write_csv(data.frame(tip = rownames(ra), ra, check.names = FALSE), + p$outdir, "tree_tip_abundance") +} + +microbe_save(gt, p$outdir, microbe_opt(p, "name", "tree"), + width = microbe_opt(p, "width", 8), height = microbe_opt(p, "height", 8), + dpi = microbe_opt(p, "dpi", 300)) +cat("MICROBE_DONE: tree tips=", length(tips), "level=", level, "layout=", lay, "\n") diff --git a/plugins/Presisitence/microbe/plugin.json b/plugins/Presisitence/microbe/plugin.json new file mode 100644 index 0000000..9bd1b5c --- /dev/null +++ b/plugins/Presisitence/microbe/plugin.json @@ -0,0 +1,22 @@ +{ + "$schema": "https://agent-plugins.org/schemas/1.0.0/plugin.schema.json", + "name": "microbe", + "version": "0.1.0", + "description": "16S/ITS downstream microbiome analysis via local Rscript (alpha/beta diversity, composition, differential abundance, co-occurrence, random forest).", + "author": { + "name": "Presisitence", + "url": "https://github.com/Presisitence" + }, + "homepage": "https://github.com/Presisitence/microbe-mcp", + "repository": "https://github.com/Presisitence/microbe-mcp", + "license": "MIT", + "keywords": [ + "minimax-code", + "plugin", + "mcp", + "16s", + "its", + "microbiome", + "vegan" + ] +} diff --git a/plugins/Presisitence/microbe/pyproject.toml b/plugins/Presisitence/microbe/pyproject.toml new file mode 100644 index 0000000..98b8353 --- /dev/null +++ b/plugins/Presisitence/microbe/pyproject.toml @@ -0,0 +1,23 @@ +[project] +name = "microbe-mcp" +version = "0.1.0" +description = "下游微生物组/扩增子(16S·ITS) 分析 MCP 服务：把 vegan/DESeq2/edgeR/igraph/Hmisc/randomForest/ggtree 等一套微生物组分析出图脚本参数化封装为 MCP 工具，由本机 Rscript 引擎渲染 png+pdf" +readme = "README.md" +requires-python = ">=3.10" +license = "MIT" +dependencies = [ + "mcp[cli]>=1.2.0", +] + +[project.scripts] +microbe-cli = "microbe_toolkit.cli:main" + +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[tool.hatch.build.targets.wheel] +packages = ["microbe_toolkit"] + +[tool.hatch.build.targets.wheel.force-include] +"mscripts" = "microbe_toolkit/mscripts" diff --git a/plugins/Presisitence/microbe/server.py b/plugins/Presisitence/microbe/server.py new file mode 100644 index 0000000..015b1a6 --- /dev/null +++ b/plugins/Presisitence/microbe/server.py @@ -0,0 +1,195 @@ +"""microbe —— 下游微生物组/扩增子(16S·ITS) 分析 MCP Server (FastMCP)。 + +把一套清洗、参数化的微生物组分析 R 脚本（mscripts/*.R；vegan/DESeq2/edgeR/igraph/Hmisc/ +randomForest/ggtree/ggplot2）暴露为 MCP 工具，由本机 Rscript 引擎渲染，产出出版级 png+pdf。 +定位不到 R 时返回可手动运行的脚本与命令；R 端缺包时返回缺失包名与安装建议。 + +工具映射两篇经典微生物组论文的核心图： + Liu 2023 (Nat Microbiol, 稻曲病叶际菌群) 与 Zhou 2022 (Nat Commun, 番茄枯萎病跨界合成菌群)。 + +数据契约（3 张 CSV）： + feature_table.csv : feature_id + 各样本计数列（OTU/ASV/属丰度表） + taxonomy.csv : feature_id + 界门纲目科属种(分列) 或 单列 QIIME 分号串(k__;p__;...) + metadata.csv : sample_name + group[/group_name][/TvsC][+ 环境/代谢物数值列] + +运行: + uv run --directory server.py +""" +from __future__ import annotations + +from typing import Any + +from mcp.server.fastmcp import FastMCP + +from microbe_toolkit import runner + +mcp = FastMCP("microbe") + + +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 microbe_env() -> dict[str, Any]: + """检查 R 引擎与关键微生物组分析包是否就绪（首次使用/报错排查先跑这个）。返回 Rscript + 路径、各包安装状态与缺包安装建议（区分 CRAN / Bioconductor）。""" + rscript = runner.find_rscript() + pkgs = [ + "vegan", "ggplot2", "ggpubr", "dplyr", "tidyr", "reshape2", "pheatmap", "ggrepel", + "scales", "igraph", "Hmisc", "randomForest", "ape", "circlize", "ggalluvial", + "edgeR", "DESeq2", "ggtree", "phyloseq", "SpiecEasi", + ] + 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 "未找到 (请设置环境变量 MICROBE_RSCRIPT 指向 Rscript.exe)", + "mscripts_dir": str(runner.mscripts_dir()), + "packages": status, + "missing": missing, + "install_hint": runner._install_hint(missing) if missing else "", + } + + +# ============================ α 多样性 ============================ +@mcp.tool() +def microbe_alpha(feature_table: str, metadata: str, outdir: str, + metrics: list[str] | None = None, test: str = "wilcox", + rarefy: bool = False, depth: int | None = None, + pairwise: bool = True, dpi: int = 300) -> dict[str, Any]: + """α多样性指数箱线图 + 组间检验（Liu2023 Fig1a,b,d,e / Zhou2022 Fig1e,f）。 + metrics 可选 Observed/Chao1/ACE/Shannon/Simpson/InvSimpson/Pielou(默认前四)。 + test: wilcox(默认)/t.test/anova/kruskal；rarefy=True 先抽平(depth 缺省取最小样本深度)。 + 产出 alpha_diversity.csv(各样本指数) + alpha_boxplot.png/pdf。需 vegan(已装)； + 显著性标注需 ggpubr(已装)。""" + return _run("alpha.R", _clean(feature_table=feature_table, metadata=metadata, metrics=metrics, + test=test, rarefy=rarefy, depth=depth, pairwise=pairwise, + dpi=dpi), outdir) + + +# ============================ β 多样性 ============================ +@mcp.tool() +def microbe_beta(feature_table: str, metadata: str, outdir: str, method: str = "pcoa", + distance: str = "bray", test: str = "permanova", permutations: int = 999, + relabund: bool = True, ellipse: bool = True, label: bool = False, + dpi: int = 300) -> dict[str, Any]: + """β多样性排序 (PCoA/NMDS) + 组间差异检验（Liu2023 Fig1c,f / Zhou2022 Fig1c,d）。 + method: pcoa(默认)/nmds；distance: bray(默认)/jaccard/euclidean； + test: permanova(adonis2,默认)/anosim/both；relabund=True 先转相对丰度；ellipse 加 95% 椭圆。 + 产出 beta_scores.csv、beta_stats.csv(R2/P)、beta_.png/pdf(副标题含 PERMANOVA R2/P)。需 vegan。""" + return _run("beta.R", _clean(feature_table=feature_table, metadata=metadata, method=method, + distance=distance, test=test, permutations=permutations, + relabund=relabund, ellipse=ellipse, label=label, dpi=dpi), outdir) + + +# ============================ 物种组成 ============================ +@mcp.tool() +def microbe_composition(feature_table: str, taxonomy: str, metadata: str, outdir: str, + level: str = "Phylum", top_n: int = 10, mode: str = "group", + dpi: int = 300) -> dict[str, Any]: + """物种组成堆叠柱状图（Liu2023 Fig1g / Zhou2022 Fig4k,l, Fig5）。 + level: Phylum(默认)/Class/Order/Family/Genus/Species；top_n 其余合并为 Others； + mode: group(组内均值,默认)/sample(逐样本)。产出 composition_.csv(相对丰度表) + + composition_.png/pdf。taxonomy 支持分列或 QIIME 分号串。""" + return _run("composition.R", _clean(feature_table=feature_table, taxonomy=taxonomy, + metadata=metadata, level=level, top_n=top_n, mode=mode, + dpi=dpi), outdir) + + +# ============================ 差异丰度 ============================ +@mcp.tool() +def microbe_diff(feature_table: str, metadata: str, outdir: str, taxonomy: str | None = None, + level: str | None = None, method: str = "auto", + group_test: str | None = None, group_ref: str | None = None, + padj: float = 0.05, log2fc: float = 1.0, min_count: float = 1.0, + min_prev: int = 2, kind: str = "bar", top_n: int = 30, + dpi: int = 300) -> dict[str, Any]: + """差异丰度分析（Liu2023 Fig3a,b DESeq2 / Zhou2022 EdgeR）。 + method: auto(有 DESeq2 用 deseq2 否则 edger,默认)/deseq2/edger/wilcox。 + 给 taxonomy+level(如 Genus) 则按该层级汇总后做差异，否则按 feature。 + 两组来源: group_test/group_ref → metadata 的 TvsC(treatment/control) → 前两组。 + kind: bar(显著类群 log2FC 条形,按富集组着色,默认)/volcano。产出 diff_result.csv + diff_.png/pdf。 + 注: DESeq2 用 poscounts 处理微生物组多零；本机已装 DESeq2/edgeR。""" + return _run("diff.R", _clean(feature_table=feature_table, metadata=metadata, taxonomy=taxonomy, + level=level, method=method, group_test=group_test, + group_ref=group_ref, padj=padj, log2fc=log2fc, + min_count=min_count, min_prev=min_prev, kind=kind, top_n=top_n, + dpi=dpi), outdir) + + +# ============================ 共现网络 ============================ +@mcp.tool() +def microbe_network(feature_table: str, outdir: str, metadata: str | None = None, + taxonomy: str | None = None, level: str | None = None, + group: str | None = None, method: str = "spearman", + r_threshold: float = 0.7, p_threshold: float = 0.05, padjust: bool = True, + min_prev: int | None = None, top_n: int = 150, layout: str = "fr", + label: bool = False, dpi: int = 300) -> dict[str, Any]: + """微生物共现网络（Liu2023 Fig1i / Zhou2022 Fig2a-f）。Hmisc::rcorr 算 Spearman ρ， + |r|>=r_threshold(0.7) 且 P dict[str, Any]: + """随机森林分类 biomarker + 十折交叉验证误差曲线（Zhou2022 Fig2g-j）。 + 以相对丰度训练 randomForest，输出 MeanDecreaseGini 重要性(按富集组着色)与 OOB 误差； + cv=True 做 rfcv 十折交叉验证并画误差-特征数曲线以定最优 biomarker 数。 + 给 taxonomy+level(如 Genus) 则按该层级建模，否则按 feature。产出 rf_importance.csv/png/pdf、 + rf_cv.csv、rf_cv_error.png/pdf。需 randomForest(已装)。""" + return _run("rf.R", _clean(feature_table=feature_table, metadata=metadata, taxonomy=taxonomy, + level=level, top_n=top_n, ntree=ntree, cv=cv, cv_fold=cv_fold, + seed=seed, dpi=dpi), outdir) + + +# ============================ 微生物-代谢物/环境 相关热图 ============================ +@mcp.tool() +def microbe_corr(feature_table: str, metadata: str, outdir: str, taxonomy: str | None = None, + level: str | None = None, variables: list[str] | None = None, + features: list[str] | None = None, method: str = "spearman", top_n: int = 30, + padjust: bool = True, star: bool = True, dpi: int = 300) -> dict[str, Any]: + """微生物-代谢物/环境因子 相关热图（Liu2023 Fig2f: 丰度与氨基酸/毒素相关, r + P）。 + variables: metadata 中的数值列名向量(代谢物/环境因子；缺省自动取除分组外的数值列)； + features: 限定微生物类群(缺省取最丰 top_n)。method: spearman(默认)/pearson； + star=True 叠加显著性星号(*<0.05 **<0.01 ***<0.001，默认 BH 校正)。 + 产出 corr_r.csv、corr_p.csv、corr_heatmap.png/pdf(蓝-白-红填 r)。需 Hmisc(已装)。""" + return _run("corr.R", _clean(feature_table=feature_table, metadata=metadata, taxonomy=taxonomy, + level=level, variables=variables, features=features, method=method, + top_n=top_n, padjust=padjust, star=star, dpi=dpi), outdir) + + +# ============================ 分类树/发育树 + 丰度环 ============================ +@mcp.tool() +def microbe_tree(feature_table: str, taxonomy: str, metadata: str, outdir: str, + level: str = "Genus", color_by: str = "Phylum", layout: str = "circular", + ring: bool = True, top_n: int = 60, dpi: int = 300) -> dict[str, Any]: + """分类树/系统发育树 + 组均值丰度环（Zhou2022 Fig3a,b）。无测序数据时以分类学层级 + (ape::as.phylo formula)构树；失败自动退化为按丰度谱(1-Pearson)聚类构树。 + level: 树尖层级(默认 Genus)；color_by: 树尖着色的分类等级(默认 Phylum)； + layout: circular(默认)/fan/rectangular；ring=True 加组均值相对丰度热环。 + 产出 tree_tip_abundance.csv + tree.png/pdf。需 ape+ggtree(已装)。""" + return _run("tree.R", _clean(feature_table=feature_table, taxonomy=taxonomy, metadata=metadata, + level=level, color_by=color_by, layout=layout, ring=ring, + top_n=top_n, dpi=dpi), outdir) + + +if __name__ == "__main__": + mcp.run() diff --git a/plugins/Presisitence/microbe/skills/microbe/SKILL.md b/plugins/Presisitence/microbe/skills/microbe/SKILL.md new file mode 100644 index 0000000..73a8093 --- /dev/null +++ b/plugins/Presisitence/microbe/skills/microbe/SKILL.md @@ -0,0 +1,19 @@ +--- +name: microbe +description: Use when the user has a 16S or ITS feature table plus taxonomy and metadata and needs alpha/beta diversity, composition bars, differential abundance, co-occurrence networks, random-forest biomarkers, or a taxonomy tree. Drive the microbe MCP tools with local Rscript. +--- + +# microbe + +Use the `microbe` MCP tools. Do not reimplement vegan/DESeq2 plots in Python. + +## Typical order + +1. `microbe_env` +2. `microbe_alpha`, `microbe_beta` +3. `microbe_composition`, `microbe_diff` +4. `microbe_network`, `microbe_rf`, `microbe_corr`, `microbe_tree` as requested + +Expected columns: `feature_id` + sample counts; taxonomy ranks or a QIIME semicolon string; `sample_name` + `group`. + +Missing R or packages: return the tool's repair command. This Plugin does not go from raw reads to ASVs. diff --git a/plugins/Presisitence/microbe/tests/list_tools.py b/plugins/Presisitence/microbe/tests/list_tools.py new file mode 100644 index 0000000..fd30e31 --- /dev/null +++ b/plugins/Presisitence/microbe/tests/list_tools.py @@ -0,0 +1,7 @@ +import asyncio +import server + +tools = asyncio.run(server.mcp.list_tools()) +print("TOOL_COUNT:", len(tools)) +for t in tools: + print(" -", t.name) diff --git a/plugins/Presisitence/microbe/tests/prep_test.py b/plugins/Presisitence/microbe/tests/prep_test.py new file mode 100644 index 0000000..ba409d7 --- /dev/null +++ b/plugins/Presisitence/microbe/tests/prep_test.py @@ -0,0 +1,137 @@ +"""生成合成扩增子测试数据（stdlib-only，无需 numpy）到 tests/data/： + feature_table.csv : feature_id + 12 样本计数 (DP1-6 病株, DSP1-6 抑病株) + taxonomy.csv : feature_id + Kingdom..Genus + metadata.csv : sample_name, group, group_name, TvsC, 代谢物列(Leucine/Isoleucine/Valine/Ustiloxin) + +数据结构（刻意复刻两篇论文的信号，便于冒烟验证）： + * L1 = 抑病轴：DSP≈+1 / DP≈-1。L1+ 类群富集于 DSP、L1- 富集于 DP（供 diff/rf/composition）。 + * L2/L3 = 与分组无关的共变模块（供 network 出现满足 |ρ|>0.7 的边）。 + * BCAA(Leucine/Isoleucine/Valine)=30+10*L1（DSP 高）；Ustiloxin=100-40*L1（DP 高，与 BCAA 负相关，复刻 Liu Fig2f）。 +""" +import csv +import math +import random +from pathlib import Path + +random.seed(42) +OUT = Path(__file__).resolve().parent / "data" +OUT.mkdir(parents=True, exist_ok=True) + +# ---- 样本与潜变量 ---------------------------------------------------------- +dp = [f"DP{i}" for i in range(1, 7)] +dsp = [f"DSP{i}" for i in range(1, 7)] +samples = dp + dsp +L1 = {s: (random.gauss(-1.0, 0.3) if s in dp else random.gauss(1.0, 0.3)) for s in samples} +L2 = {s: random.gauss(0, 1) for s in samples} +L3 = {s: random.gauss(0, 1) for s in samples} + +# ---- 分类谱系 (Kingdom, Phylum, Class, Order, Family, Genus) --------------- +LINEAGES = [ + ("Bacteria", "Proteobacteria", "Gammaproteobacteria", "Pseudomonadales", "Pseudomonadaceae", "Pseudomonas"), + ("Bacteria", "Proteobacteria", "Gammaproteobacteria", "Enterobacterales", "Enterobacteriaceae", "Enterobacter"), + ("Bacteria", "Proteobacteria", "Gammaproteobacteria", "Xanthomonadales", "Xanthomonadaceae", "Stenotrophomonas"), + ("Bacteria", "Proteobacteria", "Alphaproteobacteria", "Sphingomonadales", "Sphingomonadaceae", "Sphingomonas"), + ("Bacteria", "Proteobacteria", "Alphaproteobacteria", "Rhizobiales", "Rhizobiaceae", "Rhizobium"), + ("Bacteria", "Proteobacteria", "Alphaproteobacteria", "Rhizobiales", "Methylobacteriaceae", "Methylobacterium"), + ("Bacteria", "Proteobacteria", "Betaproteobacteria", "Burkholderiales", "Comamonadaceae", "Variovorax"), + ("Bacteria", "Proteobacteria", "Betaproteobacteria", "Burkholderiales", "Burkholderiaceae", "Ralstonia"), + ("Bacteria", "Firmicutes", "Bacilli", "Lactobacillales", "Lactobacillaceae", "Lactobacillus"), + ("Bacteria", "Firmicutes", "Bacilli", "Bacillales", "Bacillaceae", "Bacillus"), + ("Bacteria", "Firmicutes", "Bacilli", "Bacillales", "Paenibacillaceae", "Paenibacillus"), + ("Bacteria", "Firmicutes", "Clostridia", "Clostridiales", "Clostridiaceae", "Clostridium"), + ("Bacteria", "Actinobacteria", "Actinomycetia", "Micrococcales", "Microbacteriaceae", "Curtobacterium"), + ("Bacteria", "Actinobacteria", "Actinomycetia", "Micrococcales", "Microbacteriaceae", "Leifsonia"), + ("Bacteria", "Actinobacteria", "Actinomycetia", "Streptomycetales", "Streptomycetaceae", "Streptomyces"), + ("Bacteria", "Bacteroidetes", "Sphingobacteriia", "Sphingobacteriales", "Sphingobacteriaceae", "Pedobacter"), + ("Bacteria", "Bacteroidetes", "Flavobacteriia", "Flavobacteriales", "Flavobacteriaceae", "Flavobacterium"), + ("Bacteria", "Bacteroidetes", "Cytophagia", "Cytophagales", "Cytophagaceae", "Dyadobacter"), + ("Fungi", "Ascomycota", "Eurotiomycetes", "Eurotiales", "Aspergillaceae", "Aspergillus"), + ("Fungi", "Ascomycota", "Eurotiomycetes", "Eurotiales", "Aspergillaceae", "Penicillium"), + ("Fungi", "Ascomycota", "Sordariomycetes", "Hypocreales", "Nectriaceae", "Fusarium"), + ("Fungi", "Ascomycota", "Sordariomycetes", "Hypocreales", "Clavicipitaceae", "Ustilaginoidea"), + ("Fungi", "Ascomycota", "Sordariomycetes", "Hypocreales", "Hypocreaceae", "Trichoderma"), + ("Fungi", "Ascomycota", "Dothideomycetes", "Pleosporales", "Pleosporaceae", "Alternaria"), + ("Fungi", "Ascomycota", "Dothideomycetes", "Capnodiales", "Cladosporiaceae", "Cladosporium"), + ("Fungi", "Ascomycota", "Sordariomycetes", "Hypocreales", "Bionectriaceae", "Acremonium"), + ("Fungi", "Basidiomycota", "Ustilaginomycetes", "Ustilaginales", "Ustilaginaceae", "Moesziomyces"), + ("Fungi", "Basidiomycota", "Tremellomycetes", "Tremellales", "Trimorphomycetaceae", "Nigrospora"), + ("Fungi", "Mortierellomycota", "Mortierellomycetes", "Mortierellales", "Mortierellaceae", "Mortierella"), + ("Fungi", "Ascomycota", "Sordariomycetes", "Sordariales", "Chaetomiaceae", "Chaetomium"), +] + +# 富集/共变模块指派：genus -> ("dsp"|"dp"|"m2"|"m3"|"none") +# 复刻论文：Lactobacillus/Aspergillus 等富集于抑病株；Ustilaginoidea/Fusarium 富集于病株。 +MODULE = { + "Lactobacillus": "dsp", "Aspergillus": "dsp", "Sphingomonas": "dsp", "Bacillus": "dsp", + "Pseudomonas": "dsp", "Curtobacterium": "dsp", "Paenibacillus": "dsp", "Nigrospora": "dsp", + "Ustilaginoidea": "dp", "Fusarium": "dp", "Ralstonia": "dp", "Alternaria": "dp", + "Enterobacter": "dp", "Cladosporium": "dp", + "Streptomyces": "m2", "Rhizobium": "m2", "Methylobacterium": "m2", "Variovorax": "m2", + "Pedobacter": "m2", "Flavobacterium": "m2", + "Trichoderma": "m3", "Penicillium": "m3", "Mortierella": "m3", "Chaetomium": "m3", + "Acremonium": "m3", "Dyadobacter": "m3", +} + +RANKS = ["Kingdom", "Phylum", "Class", "Order", "Family", "Genus"] +features = [] # (feature_id, lineage, module, base_mean, loading) +oid = 0 +for lin in LINEAGES: + genus = lin[5] + mod = MODULE.get(genus, "none") + n_otu = random.choice([2, 3, 4]) + for _ in range(n_otu): + oid += 1 + base = math.exp(random.gauss(3.2, 1.1)) # 基础丰度 + load = random.uniform(0.7, 1.0) # 模块载荷 + features.append((f"OTU{oid:03d}", lin, mod, base, load)) + + +def scale(mod, load, s): + if mod == "dsp": + return math.exp(load * L1[s] + random.gauss(0, 0.25)) + if mod == "dp": + return math.exp(-load * L1[s] + random.gauss(0, 0.25)) + if mod == "m2": + return math.exp(load * L2[s] + random.gauss(0, 0.25)) + if mod == "m3": + return math.exp(load * L3[s] + random.gauss(0, 0.25)) + return math.exp(random.gauss(0, 0.5)) + + +# ---- 特征表 ---------------------------------------------------------------- +with open(OUT / "feature_table.csv", "w", newline="", encoding="utf-8") as f: + w = csv.writer(f) + w.writerow(["feature_id"] + samples) + for fid, lin, mod, base, load in features: + row = [fid] + for s in samples: + val = base * scale(mod, load, s) + # 稀有类群偶发 0 + cnt = 0 if random.random() < 0.05 else int(round(val)) + row.append(cnt) + w.writerow(row) + +# ---- 分类表 ---------------------------------------------------------------- +with open(OUT / "taxonomy.csv", "w", newline="", encoding="utf-8") as f: + w = csv.writer(f) + w.writerow(["feature_id"] + RANKS) + for fid, lin, mod, base, load in features: + w.writerow([fid] + list(lin)) + +# ---- 元数据（含代谢物，复刻 BCAA 高于 DSP、Ustiloxin 高于 DP 且负相关）------ +with open(OUT / "metadata.csv", "w", newline="", encoding="utf-8") as f: + w = csv.writer(f) + w.writerow(["sample_name", "group", "group_name", "TvsC", + "Leucine", "Isoleucine", "Valine", "Ustiloxin"]) + for s in samples: + tvsc = "treatment" if s in dsp else "control" + leu = 30 + 10 * L1[s] + random.gauss(0, 1.5) + ile = 18 + 6 * L1[s] + random.gauss(0, 1.2) + val = 22 + 7 * L1[s] + random.gauss(0, 1.2) + ust = max(0.0, 100 - 40 * L1[s] + random.gauss(0, 6)) + w.writerow([s, s[:-1] if s[-1].isdigit() else s, + "DSP" if s in dsp else "DP", tvsc, + round(leu, 2), round(ile, 2), round(val, 2), round(ust, 2)]) + +print(f"OK: {len(features)} OTUs x {len(samples)} samples -> {OUT}") +print(" feature_table.csv / taxonomy.csv / metadata.csv") diff --git a/plugins/Presisitence/microbe/tests/smoke.py b/plugins/Presisitence/microbe/tests/smoke.py new file mode 100644 index 0000000..be6b51f --- /dev/null +++ b/plugins/Presisitence/microbe/tests/smoke.py @@ -0,0 +1,41 @@ +"""端到端冒烟测试：Python runner -> Rscript -> _common.R -> 各 mscripts/*.R。 +先运行 prep_test.py 生成 tests/data/ 的三张 CSV，再依次跑 9 个工具。""" +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from microbe_toolkit import runner # noqa: E402 + +DATA = Path(__file__).resolve().parent / "data" +OUT = Path(__file__).resolve().parent / "out" +OUT.mkdir(parents=True, exist_ok=True) +FT = str(DATA / "feature_table.csv") +TAX = str(DATA / "taxonomy.csv") +MD = str(DATA / "metadata.csv") + + +def show(tag, res): + st = res["status"] + if st == "ok": + extra = f"outputs={res['n_outputs']}" + for mtr in res.get("metrics", []): + extra += f"\n · {mtr}" + elif st == "missing_packages": + extra = f"missing={res['packages']}" + else: + extra = (res.get("error", "") + " | " + (res.get("log", "")[-500:])) + print(f"[{st:16}] {tag:22} {extra}") + + +print("Rscript:", runner.find_rscript()) +show("alpha", runner.run_script("alpha.R", {"feature_table": FT, "metadata": MD, "test": "wilcox"}, outdir=str(OUT))) +show("beta.pcoa", runner.run_script("beta.R", {"feature_table": FT, "metadata": MD, "method": "pcoa", "test": "both"}, outdir=str(OUT))) +show("beta.nmds", runner.run_script("beta.R", {"feature_table": FT, "metadata": MD, "method": "nmds"}, outdir=str(OUT / "nmds"))) +show("composition.Phylum", runner.run_script("composition.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "level": "Phylum"}, outdir=str(OUT))) +show("composition.Genus", runner.run_script("composition.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "level": "Genus", "top_n": 12, "mode": "sample"}, outdir=str(OUT / "genus"))) +show("diff.auto(deseq2)", runner.run_script("diff.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "level": "Genus", "group_test": "DSP", "group_ref": "DP"}, outdir=str(OUT))) +show("diff.edger", runner.run_script("diff.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "level": "Genus", "method": "edger", "group_test": "DSP", "group_ref": "DP", "kind": "volcano"}, outdir=str(OUT / "edger"))) +show("network", runner.run_script("network.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "r_threshold": 0.7, "label": True}, outdir=str(OUT))) +show("rf", runner.run_script("rf.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "level": "Genus", "top_n": 15}, outdir=str(OUT))) +show("corr", runner.run_script("corr.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "level": "Genus", "variables": ["Leucine", "Isoleucine", "Valine", "Ustiloxin"]}, outdir=str(OUT))) +show("tree", runner.run_script("tree.R", {"feature_table": FT, "taxonomy": TAX, "metadata": MD, "level": "Genus", "layout": "circular"}, outdir=str(OUT)))