From abdd486ed8cf8ae56c3edaa0fc602388e4638434 Mon Sep 17 00:00:00 2001 From: potaser Date: Sun, 30 Aug 2026 00:17:25 +0800 Subject: [PATCH] =?UTF-8?q?WGCNA=20=E6=A8=A1=E5=9D=97-=E6=80=A7=E7=8A=B6?= =?UTF-8?q?=E7=83=AD=E5=9B=BE=E6=94=B9=E7=94=A8=20BH=20=E6=A0=A1=E6=AD=A3?= =?UTF-8?q?=20p=20=E5=80=BC=EF=BC=8C=E5=B9=B6=E6=94=AF=E6=8C=81=20job=20?= =?UTF-8?q?=E7=BD=AE=E9=A1=B6=E4=B8=8E=E4=BC=98=E5=85=88=E7=BA=A7=E6=8E=92?= =?UTF-8?q?=E5=BA=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- apps/cli/run_rnaseq_wgcna.py | 19 +++++++++++--- apps/web/vue3/main.py | 25 +++++++++++++++++-- .../stages/rnaseq_wgcna_contract_v1.md | 9 +++++++ .../test_rnaseq_wgcna_adjusted_pvalues.py | 16 ++++++++++++ 4 files changed, 63 insertions(+), 6 deletions(-) create mode 100644 tests/unit/test_rnaseq_wgcna_adjusted_pvalues.py diff --git a/apps/cli/run_rnaseq_wgcna.py b/apps/cli/run_rnaseq_wgcna.py index 8a38c7d36..62be2eb1c 100755 --- a/apps/cli/run_rnaseq_wgcna.py +++ b/apps/cli/run_rnaseq_wgcna.py @@ -360,23 +360,29 @@ colnames(condition_mm) <- sub("^condition", "", colnames(condition_mm)) traitData <- as.data.frame(condition_mm) moduleTraitCor <- cor(MEs, traitData, use = "p") moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples = nrow(datExpr)) +moduleTraitPadj <- matrix( + p.adjust(as.vector(moduleTraitPvalue), method = "BH"), + nrow = nrow(moduleTraitPvalue), ncol = ncol(moduleTraitPvalue), + dimnames = dimnames(moduleTraitPvalue) +) write.table(moduleTraitCor, file = paste0(kv$out_prefix, ".module_trait_correlation.tsv"), sep = "\t", quote = FALSE, col.names = NA) write.table(moduleTraitPvalue, file = paste0(kv$out_prefix, ".module_trait_pvalue.tsv"), sep = "\t", quote = FALSE, col.names = NA) +write.table(moduleTraitPadj, file = paste0(kv$out_prefix, ".module_trait_padj.tsv"), sep = "\t", quote = FALSE, col.names = NA) write.table(MEs, file = paste0(kv$out_prefix, ".module_eigengenes.tsv"), sep = "\t", quote = FALSE, col.names = NA) -textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "") +textMatrix <- paste(signif(moduleTraitCor, 2), "\n(FDR ", signif(moduleTraitPadj, 2), ")", sep = "") dim(textMatrix) <- dim(moduleTraitCor) png(paste0(kv$out_prefix, ".module_trait_heatmap.png"), width = 1800, height = 1600, res = 220) labeledHeatmap(Matrix = moduleTraitCor, xLabels = colnames(traitData), yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.7, - zlim = c(-1, 1), main = "Module-trait relationships") + zlim = c(-1, 1), main = "Module-trait relationships (BH-adjusted P values)") dev.off() pdf(paste0(kv$out_prefix, ".module_trait_heatmap.pdf"), width = 8, height = 7) labeledHeatmap(Matrix = moduleTraitCor, xLabels = colnames(traitData), yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.7, - zlim = c(-1, 1), main = "Module-trait relationships") + zlim = c(-1, 1), main = "Module-trait relationships (BH-adjusted P values)") dev.off() geneModule <- data.frame(gene_id = colnames(datExpr), module = moduleColors) @@ -394,7 +400,8 @@ colnames(moduleSizes) <- c("module", "genes") write.table(moduleSizes, file = paste0(kv$out_prefix, ".module_sizes.tsv"), sep = "\t", quote = FALSE, row.names = FALSE) saveRDS(list(datExpr = datExpr, moduleColors = moduleColors, MEs = MEs, softPower = softPower, - moduleTraitCor = moduleTraitCor, moduleTraitPvalue = moduleTraitPvalue), + moduleTraitCor = moduleTraitCor, moduleTraitPvalue = moduleTraitPvalue, + moduleTraitPadj = moduleTraitPadj), file = paste0(kv$out_prefix, ".wgcna_result.rds")) summary <- data.frame( @@ -475,6 +482,7 @@ def write_report(path: Path, manifest: dict[str, Any]) -> None: "", "- WGCNA is most reliable with larger sample sizes; interpret very small RNA-seq designs cautiously.", "- Genes are filtered by expression and variance before TOM construction to keep the analysis tractable.", + "- Module-trait heatmap labels use Benjamini-Hochberg adjusted P values (FDR) across all displayed module-trait tests; raw P values are retained as a separate artifact.", "", "## Artifacts", "", @@ -581,6 +589,7 @@ def main(argv: list[str] | None = None) -> int: "gene_dendrogram_modules_pdf": (Path(f"{out_prefix}.gene_dendrogram_modules.pdf"), "pdf"), "module_trait_correlation": (Path(f"{out_prefix}.module_trait_correlation.tsv"), "tsv"), "module_trait_pvalue": (Path(f"{out_prefix}.module_trait_pvalue.tsv"), "tsv"), + "module_trait_padj": (Path(f"{out_prefix}.module_trait_padj.tsv"), "tsv"), "module_trait_heatmap_png": (Path(f"{out_prefix}.module_trait_heatmap.png"), "image"), "module_trait_heatmap_pdf": (Path(f"{out_prefix}.module_trait_heatmap.pdf"), "pdf"), "module_eigengenes": (Path(f"{out_prefix}.module_eigengenes.tsv"), "tsv"), @@ -607,6 +616,7 @@ def main(argv: list[str] | None = None) -> int: metrics = read_summary(Path(f"{out_prefix}.wgcna_summary.tsv")) metrics["conditions_detail"] = {sample: groups[sample] for sample in samples} metrics["matrix_type"] = matrix_type + metrics["module_trait_pvalue_adjustment"] = "BH_global" manifest = { "manifest_version": "1.0", "stage": "rnaseq_wgcna", @@ -637,6 +647,7 @@ def main(argv: list[str] | None = None) -> int: "min_module_size": args.min_module_size, "merge_cut_height": args.merge_cut_height, "powers": args.powers, + "module_trait_pvalue_adjustment": "BH_global", "cpu": threads, "threads": threads, }, diff --git a/apps/web/vue3/main.py b/apps/web/vue3/main.py index e76beb7f5..d88b9f061 100644 --- a/apps/web/vue3/main.py +++ b/apps/web/vue3/main.py @@ -8970,9 +8970,28 @@ async def list_existing_jobs(): except OSError as exc: skipped.append({"job_id": job_dir.name, "reason": str(exc)}) continue - job_entries.append((job_dir, stat_result.st_mtime)) + pinned = False + priority = 0 + priority_file = job_dir / ".plantai_job_priority.json" + if priority_file.is_file(): + try: + priority_data = json.loads(priority_file.read_text(encoding="utf-8")) + pinned = priority_data.get("pinned") is True + priority = int(priority_data.get("priority", 0)) + except (OSError, ValueError, TypeError, json.JSONDecodeError) as exc: + skipped.append({ + "job_id": job_dir.name, + "reason": f"invalid_job_priority_metadata: {exc}", + }) + job_entries.append((job_dir, stat_result.st_mtime, pinned, priority)) - for job_dir, st_mtime in sorted(job_entries, key=lambda item: item[1], reverse=True): + # Pinned jobs come first, then explicit priority, while mtime remains + # the stable default ordering for jobs without metadata. + for job_dir, st_mtime, pinned, priority in sorted( + job_entries, + key=lambda item: (item[2], item[3], item[1]), + reverse=True, + ): if job_dir.is_dir() and job_dir.name.startswith("job_"): # 检查目录内容 final_dir = job_dir / "final" @@ -8990,6 +9009,8 @@ async def list_existing_jobs(): jobs.append({ "job_id": job_dir.name, "created_at": mtime.strftime("%Y-%m-%d %H:%M:%S"), + "pinned": pinned, + "priority": priority, "completed_stages": completed_stages, "has_input": (INPUT_ROOT / job_dir.name).exists(), "workspace_dir": str(job_dir) diff --git a/gfa_core/schemas/stages/rnaseq_wgcna_contract_v1.md b/gfa_core/schemas/stages/rnaseq_wgcna_contract_v1.md index 88c3824e3..5d058f9a6 100755 --- a/gfa_core/schemas/stages/rnaseq_wgcna_contract_v1.md +++ b/gfa_core/schemas/stages/rnaseq_wgcna_contract_v1.md @@ -46,6 +46,7 @@ - `{prefix}.rnaseq_wgcna.module_eigengenes.tsv` - `{prefix}.rnaseq_wgcna.module_trait_correlation.tsv` - `{prefix}.rnaseq_wgcna.module_trait_pvalue.tsv` +- `{prefix}.rnaseq_wgcna.module_trait_padj.tsv` - `{prefix}.rnaseq_wgcna.wgcna_summary.tsv` ### final/ 图形产物 @@ -103,8 +104,16 @@ moduleColors <- labels2colors(dynamicMods) merge <- mergeCloseModules(datExpr, moduleColors, cutHeight = mergeCutHeight) MEs <- merge$newMEs moduleTraitCor <- cor(MEs, traits, use = "p") +moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples = nrow(datExpr)) +moduleTraitPadj <- matrix( + p.adjust(as.vector(moduleTraitPvalue), method = "BH"), + nrow = nrow(moduleTraitPvalue), ncol = ncol(moduleTraitPvalue), + dimnames = dimnames(moduleTraitPvalue) +) ``` +模块-性状热图单元格显示相关系数与 `moduleTraitPadj`。校正范围为当前热图中的全部模块×性状检验，方法为 Benjamini-Hochberg；原始 P 值继续写入 `module_trait_pvalue.tsv` 供追溯和兼容旧下游。 + ## 7. 解释限制 WGCNA 对样本数量敏感。少样本数据可以完成计算，但模块稳定性和模块-性状相关性需要谨慎解释，建议后续进行模块稳定性或重采样验证。 diff --git a/tests/unit/test_rnaseq_wgcna_adjusted_pvalues.py b/tests/unit/test_rnaseq_wgcna_adjusted_pvalues.py new file mode 100644 index 000000000..26a04e3c0 --- /dev/null +++ b/tests/unit/test_rnaseq_wgcna_adjusted_pvalues.py @@ -0,0 +1,16 @@ +from pathlib import Path + +from apps.cli.run_rnaseq_wgcna import write_r_script + + +def test_wgcna_heatmap_uses_bh_adjusted_pvalues(tmp_path: Path) -> None: + script_path = tmp_path / "run_wgcna.R" + write_r_script(script_path, threads=2) + script = script_path.read_text(encoding="utf-8") + + assert 'p.adjust(as.vector(moduleTraitPvalue), method = "BH")' in script + assert '".module_trait_pvalue.tsv"' in script + assert '".module_trait_padj.tsv"' in script + assert '"\\n(FDR ", signif(moduleTraitPadj, 2)' in script + assert "Module-trait relationships (BH-adjusted P values)" in script + assert "moduleTraitPadj = moduleTraitPadj" in script