diff --git a/docs/bioflow_corpus/PlantAI_现有分析全流程_基因家族_RNAseq_WGS_GWAS.md b/docs/bioflow_corpus/PlantAI_现有分析全流程_基因家族_RNAseq_WGS_GWAS.md new file mode 100644 index 0000000..7df1ca3 --- /dev/null +++ b/docs/bioflow_corpus/PlantAI_现有分析全流程_基因家族_RNAseq_WGS_GWAS.md @@ -0,0 +1,2359 @@ +# PlantAI 现有分析全流程说明 + +> 覆盖范围：基因家族分析、RNA-seq、WGS 全基因组重测序、群体遗传分析与 GWAS +> 核对依据：服务器当前项目 `/home/xzq/5090/teamWorkspace/sj/PlantAI-dev` 的阶段契约、运行入口和技能说明 +> 核对版本：Git 提交 `47a40ee41`，分支 `fix/rnaseq-explicit-m4-interaction-20260828` +> 整理日期：2026-09-04 +> 说明：本文描述的是当前项目已经具备的分析模块及其标准衔接。不同物种、数据规模和研究设计可能需要调整参数。文中“输出文件名”以 `{prefix}` 表示用户设置的项目名称前缀。 +> 关联规划：基于本文的「对话式生信流程编排」实施方案见 [[PlantAI-对话式生信流程编排-实施计划]]（2026-09-05，规划中）。 + +--- + +## 1. 文档阅读说明 + +本文把每个分析阶段拆成六个问题： + +1. **目的**：这一步解决什么科学或数据问题。 +2. **前置条件**：必须先完成哪些阶段。 +3. **输入**：用户需要提供什么，或系统从上一阶段读取什么。 +4. **主要处理**：系统实际执行的核心分析。 +5. **输出**：供用户查看、下载或供下一阶段使用的结果。 +6. **质控与失败条件**：怎样判定结果可信，什么情况下应停止而不是继续传递错误。 + +### 1.1 四条主流程的关系 + +```mermaid +flowchart LR + A[参考基因组/注释/蛋白] --> GF[基因家族分析] + A --> RNA[RNA-seq] + A --> WGS[WGS 重测序] + WGS --> VCF[高质量群体 SNP VCF] + VCF --> POP[群体遗传分析] + VCF --> GWAS[GWAS] + POP --> GWAS + RNA --> CAND[候选基因与家族] + GF --> CAND + GWAS --> CAND +``` + +- **基因家族分析**回答“目标物种中哪些基因属于这个家族，它们的结构、进化、调控和表达有什么特点”。 +- **RNA-seq**回答“不同处理或组织中哪些基因发生表达变化，哪些通路和共表达模块与性状有关”。 +- **WGS**把原始重测序 reads 处理成经过严格过滤的群体 SNP 数据。 +- **群体遗传分析**描述样本亲缘、群体结构、连锁不平衡和 SNP 分布。 +- **GWAS**把基因型与表型联系起来，寻找与性状显著相关的遗传位点。 + +### 1.2 项目统一的目录与结果契约 + +每个任务使用独立工作区： + +```text +workspace// +├─ logs/ # 各阶段标准日志 +├─ /analysis/ # 工具原始结果和中间文件 +├─ /reproducibility/ # 命令、环境、校验和、来源记录 +└─ final/ # 稳定的正式结果，供下载和下游读取 + ├─ {prefix}._report.md + └─ {prefix}..manifest.json +``` + +统一规则： + +- 下游阶段优先读取上游成功的 Manifest，而不是随意在目录中搜索同名文件。 +- 正式文件统一使用 `{prefix}` 前缀，避免多个任务的结果混淆。 +- Manifest 记录阶段状态、输入、参数、指标、结果文件、文件大小、SHA256、开始结束时间和运行命令。 +- 成功阶段必须有非空、经过验证的科学结果；失败阶段的 `artifacts` 为空，避免旧结果被误认为本次成功结果。 +- WGS、群体遗传和 GWAS 使用分析身份与样本集合校验，阻止不同参考基因组、不同样本批次或不同参数的文件被错误混用。 +- 关键结果先在本次运行目录完成校验，再原子发布到 `final/`；较旧任务不能覆盖较新的同身份结果。 + +--- + +# 第一部分：基因家族分析 + +## 2. 基因家族完整流程总览 + +```mermaid +flowchart TD + P0[genome_prepare
整理基因组、GFF、CDS、蛋白] --> P1[identify
HMM/参考蛋白检索] + P1 --> P2[identify_refine
成员复核与精炼] + P2 --> P3[phylogeny
多序列比对与系统树] + P3 --> P4[reference_classification
亚家族分类与确认] + P4 --> P5[motif
保守基序] + P4 --> P6[domain
蛋白结构域] + P4 --> P7[structure
基因结构] + P4 --> P8[promoter/tfbs
启动子与顺式元件] + P4 --> P9[synteny/karyotype
染色体定位与共线性] + P4 --> P10[ppi/kaks/pep
互作、选择压力、蛋白性质] + P4 --> P11[expression
家族表达模式] + P5 --> P12[combo
综合图] + P6 --> P12 + P7 --> P12 + P3 --> P12 + P12 --> P13[gene_family_report
正式报告] +``` + +### 2.1 建议的最小输入包 + +| 文件 | 是否必需 | 内容 | 主要用途 | +|---|---:|---|---| +| 目标物种蛋白 FASTA | 是 | 全蛋白组，ID 唯一 | 家族成员鉴定、系统树、motif、domain、蛋白性质 | +| 目标物种 CDS FASTA | 建议 | 编码序列，ID 可与蛋白对应 | Ka/Ks、序列导出、ID 关联 | +| 目标物种基因组 FASTA | 多数后续阶段需要 | 染色体或 scaffold 序列 | 启动子提取、TFBS、染色体信息 | +| 目标物种 GFF3/GTF | 多数后续阶段需要 | 基因坐标、转录本和外显子/CDS | 基因结构、启动子、染色体定位、共线性 | +| Pfam HMM 模型 | 鉴定必需 | 一个或多个目标结构域模型 | HMM 家族成员搜索和结构域复核 | +| 已知家族参考蛋白 | 可选但推荐 | 已报道成员，需明确声明为参考序列 | BLAST 辅助鉴定、树上参照、亚家族命名 | +| 比较物种蛋白/GFF/基因组 | 视分析而定 | 一个或多个相关物种 | 跨物种系统树、种间共线性、Ka/Ks | +| 表达矩阵与样本信息 | 表达分析需要 | TPM/FPKM 或 RNA-seq 自动流程产物 | 家族表达热图、组织特异性和相关性 | + +## 3. `genome_prepare`：基因组数据整理 + +| 项目 | 说明 | +|---|---| +| 目的 | 把不同来源、不同命名习惯、可能压缩的基因组与注释文件整理成后续分析可稳定读取的标准文件。 | +| 前置条件 | 无，是基因家族流程的推荐起点。 | +| 输入 | 基因组 FASTA、GFF3/GTF、CDS FASTA、蛋白 FASTA；可为压缩或未压缩文件；可附参考/比较物种数据。 | +| 主要处理 | 识别文件类型；检查 FASTA 和注释格式；规范物种与文件命名；核对序列 ID；建立必要的 ID 别名关系；为参考家族序列生成组合 FASTA。 | +| 主要输出 | `{species}.genome.fasta`、`{species}.gff3`、`{species}.cds.fasta`、`{species}.pep.fasta`；必要时 `{species}.id_alias.tsv`；参考序列可形成 `Ref_combined.{family}.pep.fasta`/CDS；报告和 Manifest。 | +| 质控 | 文件必须非空且格式可解析；FASTA ID 不应重复；GFF 的染色体名称应能在基因组 FASTA 中找到；蛋白/CDS/注释 ID 应可建立对应关系。 | +| 下一步 | 标准蛋白组进入 `identify`；基因组和 GFF 留给 structure、promoter、TFBS、karyotype、synteny。 | + +## 4. `identify`：家族成员初筛与结构域确认 + +| 项目 | 说明 | +|---|---| +| 目的 | 从目标物种全蛋白组中找出可能属于指定基因家族的成员，并用目标结构域证据进行复核。 | +| 前置条件 | 已准备目标蛋白组；已明确一个或多个 Pfam/HMM 模型；如使用参考蛋白，必须明确声明其角色。 | +| 输入 | 目标物种全蛋白 FASTA；Pfam HMM；可选的已知家族参考蛋白 FASTA；E-value、覆盖度、CPU、证据合并方式等参数。 | +| 主要处理 | ① HMMER 搜索目标结构域；② 可选地以已知家族蛋白为 query、目标蛋白组为 database 做 BLAST；③ 合并 HMM 与 BLAST 候选；④ 再次检查结构域完整性、覆盖度和显著性；⑤ 形成最终成员集合。 | +| 主要输出 | `final/{prefix}.final_members.id`；`final/{prefix}.final_members.pep.fa`；`final/{prefix}.identify_report.md`；`final/{prefix}.identify.manifest.json`；中间目录保留 HMM 命中表、`domtblout`、BLAST 原始结果、候选 ID 和来源证据矩阵。 | +| 关键指标 | 请求的 Pfam 模型数、成功加载模型数、缺失模型、每个模型命中数、HMM 候选数、BLAST 候选数、并集/交集成员数、最终成员数。 | +| 质控 | 请求的 Pfam 模型缺失时应失败；BLAST 命中 ID 必须来自目标蛋白组的 subject ID；不能仅凭文件名猜测某 FASTA 是参考家族序列；空成员集不能伪装为成功。 | +| 下一步 | 成员 ID 和蛋白序列进入 `identify_refine`、`phylogeny`、motif、domain、PEP 等阶段。 | + +### 4.1 HMM 与参考蛋白证据的含义 + +- **HMM 证据**：从保守结构域出发，适合识别序列相似度较低但仍保留核心域的成员。 +- **参考蛋白 BLAST 证据**：从已知家族成员出发，适合补充漏检候选，但不能代替结构域复核。 +- 常用合并方式是取并集后再次做结构域筛选；若研究要求更严格，也可要求多类证据同时支持。 + +## 5. `identify_refine`：成员精炼与人工复核 + +| 项目 | 说明 | +|---|---| +| 目的 | 在自动初筛后，处理重复转录本、可疑短序列、结构域残缺成员、ID 混乱和边界成员，形成可用于正式下游分析的成员集。 | +| 输入 | `identify` 的成员 ID、蛋白序列和证据表；必要时原始蛋白组、结构域结果、参考家族序列和用户确认意见。 | +| 主要处理 | 核对每个候选的证据来源；按基因/转录本规则去冗余；检查序列长度、终止符、非法字符和结构域覆盖；标记保留、排除、待确认；重新导出精炼集合。 | +| 输出 | 精炼成员 ID、精炼蛋白 FASTA、成员决策/证据表、阶段报告和 Manifest。若形成经人工确认的亚家族集合，可供后续分类与综合图使用。 | +| 质控 | 所有排除或保留操作应可追溯；不得静默改变序列；精炼后的 FASTA ID 必须与成员表完全一致。 | +| 下一步 | 作为 `phylogeny` 的权威成员集合，并供所有家族特征分析复用。 | + +## 6. `phylogeny`：多序列比对与系统发育树 + +| 项目 | 说明 | +|---|---| +| 目的 | 比较家族蛋白的序列关系，判断成员间的亲缘远近、分支和可能的亚家族结构。 | +| 前置条件 | 已有确认或精炼后的家族蛋白 FASTA。 | +| 输入 | 默认使用 `{prefix}.final_members.pep.fa`；若走亚家族证据路线，可使用 curated members；可附参考物种或比较物种已知家族蛋白；分析范围可设为仅目标物种、加入参考、加入比较物种或全部。 | +| 主要处理 | 复制到阶段目录；仅在派生序列中清理末端 `*`；MAFFT 多序列比对；IQ-TREE 建树与模型选择/支持度计算；检查叶名称与输入成员的一致性。 | +| 主要输出 | `final/{prefix}.pep.msa.fa`；`final/{prefix}.tree.nwk`；系统树报告与 Manifest；中间目录保留 MAFFT、IQ-TREE 原始结果。 | +| 质控 | 树的叶节点集合必须与本次分析成员集合一致；不得混入其他任务或其他 scope 的序列；过短、全缺口或不可比对序列应明确报告。 | +| 下一步 | 树用于亚家族分类、基因结构/motif/domain 顺序排列、启动子图、综合图和报告。 | + +## 7. `phylogeny_msa`、`msa_region_plot` 与 `tree_plot_revise`：比对和树图辅助展示 + +### 7.1 `phylogeny_msa` + +- **目的**：把系统树使用的多序列比对结果整理成可阅读、可复用的比对数据或图形。 +- **输入**：`phylogeny` 产生的 MSA，必要时附成员分组信息。 +- **处理**：检查序列顺序、缺口比例、保守位点；生成适合展示的比对视图。 +- **输出**：整理后的 MSA、比对图或展示数据、报告和 Manifest。 +- **质控**：展示的序列和原始树叶节点必须同源，不得因排序或截取丢失成员。 + +### 7.2 `msa_region_plot` + +- **目的**：放大展示 MSA 中某一段关键保守区域，便于说明结构域或功能位点。 +- **输入**：MSA、起止位置、要展示的成员或分组。 +- **处理**：按对齐坐标截取区域，计算/展示保守性和残基差异。 +- **输出**：指定区域比对图及其参数记录。 +- **质控**：需说明使用的是对齐坐标还是原序列坐标，避免位置解释错误。 + +### 7.3 `tree_plot_revise` + +- **目的**：在不改变树拓扑证据的前提下，调整树的显示方式、标签、颜色、分组和版式。 +- **输入**：Newick 树、成员分类表、标签或配色要求。 +- **处理**：矩形/圆形等布局、分支或亚家族着色、标签替换、图例和尺寸调整。 +- **输出**：修订版树图和绘图参数记录。 +- **质控**：仅改变展示，不应悄悄重建或修改树拓扑；若重建必须回到 `phylogeny` 阶段。 + +## 8. `reference_classification`：参考成员辅助的亚家族分类 + +| 项目 | 说明 | +|---|---| +| 目的 | 利用已知参考成员在系统树上的位置，给目标成员分配候选亚家族标签。 | +| 输入 | 系统树、目标成员、带可靠标签的参考家族成员、必要的结构域或 motif 证据。 | +| 主要处理 | 将参考与目标成员放在同一树中；依据分支聚类、参考标签和辅助证据生成亚家族预测；输出待确认清单。 | +| 输出 | 亚家族预测表、证据表、确认后的 `confirmed_subfamily.tsv` 或同类正式表、报告和 Manifest。 | +| 质控 | 自动标签属于“预测”；存在参考不足、分支支持度低或证据冲突时必须标记；正式亚家族命名应经过用户确认。 | +| 下一步 | 分类结果用于树图着色、motif/domain/structure 分组、表达模式比较和综合图。 | + +## 9. `motif`：蛋白保守基序分析 + +| 项目 | 说明 | +|---|---| +| 目的 | 发现家族蛋白中较短、重复出现的保守序列模式，比较不同亚家族的 motif 组成。 | +| 输入 | 确认后的家族蛋白 FASTA；可选树顺序和亚家族分类；motif 数量、宽度等 MEME 参数。 | +| 主要处理 | MEME motif 搜索；解析 motif 坐标、宽度、显著性；按系统树或亚家族顺序绘图。 | +| 输出 | motif 表、MEME XML/文本结果、motif 图 PNG/PDF、报告和 Manifest。 | +| 质控 | motif 是统计模式，不等于已知功能域；成员数太少或序列长度差异过大时应谨慎解释。 | +| 下一步 | motif 数据进入 `combo` 综合图，也可与亚家族功能差异联合解释。 | + +## 10. `domain`：蛋白结构域注释 + +| 项目 | 说明 | +|---|---| +| 目的 | 确认家族核心结构域，并比较成员间结构域的数量、位置和组合差异。 | +| 输入 | 家族蛋白 FASTA、Pfam HMM 数据库/指定模型、显著性和覆盖度阈值。 | +| 主要处理 | HMMER 搜索；解析 `domtblout`；计算结构域位置和覆盖；生成成员结构域图。 | +| 输出 | `domtblout`、`domains.tsv`、结构域复核成员表、domain PNG/PDF、报告和 Manifest。 | +| 质控 | 应区分完整域、部分域和低显著性命中；结构域坐标不得超出蛋白长度；关键域缺失成员需单独标记。 | +| 下一步 | 用于成员确认、亚家族比较和 `combo` 综合图。 | + +## 11. `structure`：基因外显子—内含子结构 + +> 此阶段分析的是**基因结构**，不是蛋白三维结构预测。 + +| 项目 | 说明 | +|---|---| +| 目的 | 比较家族基因的外显子、内含子、CDS 和 UTR 组织方式，辅助判断家族分化和结构保守性。 | +| 输入 | 目标物种 GFF3/GTF、最终家族成员 ID、系统树；必要时 ID 别名表。 | +| 主要处理 | 从注释中筛选家族基因及其代表转录本；解析 gene/transcript/exon/CDS/UTR；按树叶顺序生成结构数据并绘图。 | +| 输出 | `final/{prefix}.structure.png`、`.structure.pdf`、报告和 Manifest；中间文件包括筛选后的 GFF3 和结构坐标表。 | +| 质控 | 家族 ID 必须能映射到注释；同一基因多转录本的选择规则要明确；染色体坐标、链方向和区段长度必须合法。 | +| 下一步 | 进入 `combo`，并与 motif/domain/亚家族结构联合解释。 | + +## 12. `promoter`：启动子顺式作用元件分析 + +| 项目 | 说明 | +|---|---| +| 目的 | 提取家族基因上游启动子区域，分析潜在顺式作用元件，推测激素、胁迫和发育调控线索。 | +| 输入 | 一个目标物种的基因组 FASTA、GFF3、家族成员 ID；标准绘图通常需要系统树；可选亚家族表；上游长度默认 2000 bp。 | +| 主要处理 | 按链方向提取转录起始位点上游序列；处理染色体边界；提交 PlantCARE 或解析已有注释；汇总元件类型、数量和位置；按树顺序绘图。 | +| 输出 | 上游序列 FASTA、PlantCARE 原始表、element counts、BED/坐标表、启动子元件 PNG/PDF、报告和 Manifest。 | +| 外部服务约束 | PlantCARE 单批建议不超过 30 条序列且约 95 KB；系统可分批提交并轮询结果。外部服务不可用时应明确失败或使用用户提供的结果，不能伪造注释。 | +| 质控 | 上游区域必须考虑正负链和染色体边界；未匹配 ID、空序列和重复元件要报告；顺式元件表示潜在调控线索，不代表已经实验证实。 | +| 下一步 | 可进入启动子修图、与表达/亚家族联合解释，或进一步做 TFBS 扫描。 | + +## 13. `promoter_plot_revise`：启动子图修订 + +- **目的**：按用户关注的元件类别、顺序、颜色和标签，重新组织启动子展示图。 +- **输入**：promoter 元件坐标与计数、树、亚家族表、展示参数。 +- **处理**：筛选/合并元件类别、调整颜色与图例、改变基因顺序和画布尺寸。 +- **输出**：修订版 PNG/PDF、绘图参数和 Manifest。 +- **质控**：不得通过修图改变原始元件注释；被隐藏或合并的类别应在参数中可追溯。 + +## 14. `tfbs`：转录因子结合位点扫描 + +| 项目 | 说明 | +|---|---| +| 目的 | 在家族基因启动子中扫描已知植物转录因子 motif，构建候选 TF—靶基因调控关系。 | +| 输入 | GFF3/GTF、基因组 FASTA、家族成员 ID、植物 motif 数据库；优先使用匹配物种的 PlantTFDB，必要时使用近缘物种或 JASPAR；默认上游 2000 bp。 | +| 主要处理 | 提取启动子；准备 motif；FIMO 扫描；按 p/q 值过滤；汇总每个基因和 motif 的命中；做富集、位置密度和网络展示。默认 FIMO `p=1e-4`、`q=0.05`。 | +| 主要输出 | `{prefix}.tfbs_sites.tsv`、motif 汇总、gene×motif 矩阵、富集表、motif 来源决策表；系统树热图、启动子结构图、富集气泡图、位置密度图、TF—gene 网络、来源汇总图和全 motif 热图；报告与 Manifest。 | +| 质控 | 必须记录 motif 来源与物种匹配程度；背景模型若未显式提供不应擅自使用全基因组背景；位点预测和网络边都是候选关系，需要实验或独立数据验证。 | +| 下一步 | 与表达、WGCNA、候选基因和启动子元件一起筛选可能的上游调控因子。 | + +## 15. `synteny`、`synteny_plot_revise` 与 `interspecies_synteny`：共线性分析 + +### 15.1 `synteny`：物种内共线性 + +| 项目 | 说明 | +|---|---| +| 目的 | 识别同一基因组内由大片段复制、全基因组复制等产生的家族同源基因对。 | +| 输入 | 目标物种蛋白 FASTA、GFF3、家族成员 ID；必要时基因组染色体信息。 | +| 主要处理 | 全蛋白相似性搜索；构建基因坐标；识别共线性区块和家族成员对；生成连线图或 Circos 类图。 | +| 输出 | 共线性 block、同源基因对表、家族共线性对、PNG/PDF、报告和 Manifest。 | +| 质控 | 蛋白 ID 与 GFF ID 必须一致；区块需满足锚定基因数等标准；相似性不等同于共线性。 | + +### 15.2 `interspecies_synteny`：物种间共线性 + +| 项目 | 说明 | +|---|---| +| 目的 | 比较目标物种与参考/比较物种间保守的染色体区段和家族直系同源关系。 | +| 输入 | 两个或多个物种的蛋白 FASTA、GFF3、染色体长度和目标家族成员。 | +| 主要处理 | 跨物种蛋白相似性、坐标整合、共线区块识别、目标家族基因对筛选和可视化。 | +| 输出 | 跨物种共线性区块、同源基因对、家族连接数据、比较图、报告和 Manifest。 | +| 质控 | 物种名称和 ID 命名空间必须明确；需区分全基因组共线性背景与目标家族高亮。 | + +### 15.3 `synteny_plot_revise` + +- **目的**：调整共线性图的物种顺序、染色体显示、家族高亮、颜色、标签和布局。 +- **输入**：已确认的共线性区块/基因对和绘图参数。 +- **输出**：修订版图和参数记录。 +- **质控**：只改变展示，不改变已识别的同源对或区块。 + +## 16. `karyotype`：染色体定位 + +| 项目 | 说明 | +|---|---| +| 目的 | 展示家族基因在各染色体上的位置和分布密度。 | +| 输入 | 基因组/染色体长度、GFF3、家族成员 ID；可选亚家族标签。 | +| 主要处理 | 提取染色体长度和基因坐标；按染色体排序；映射家族基因；绘制染色体位置图。 | +| 输出 | 家族基因染色体坐标表、karyotype PNG/PDF、报告和 Manifest。 | +| 质控 | scaffold 是否纳入需明确；基因坐标必须在染色体范围内；未定位成员单独报告。 | + +## 17. `kaks`：Ka/Ks 与选择压力分析 + +| 项目 | 说明 | +|---|---| +| 目的 | 对同源基因对估算非同义替换率 Ka、同义替换率 Ks 和 Ka/Ks，辅助判断净化选择、正选择或分化时间线索。 | +| 输入 | 物种 1/物种 2 或同一物种的 CDS FASTA、经过确认的同源/直系同源基因对；必要时蛋白用于密码子对齐。 | +| 主要处理 | ID 配对；蛋白或密码子层面对齐；估算 Ka、Ks、Ka/Ks；汇总分布；可进行 codeml 等模型分析和分化时间估算。 | +| 输出 | `kaks_results.csv`、Ka/Ks 分布 PNG/PDF、分化时间结果、codeml 结果、报告和 Manifest。 | +| 质控 | CDS 必须为正确阅读框且可对应；饱和的 Ks、极短序列或不可估计值要标记；Ka/Ks>1 只是候选正选择线索，不能单独作为定论。 | + +## 18. `ppi`：蛋白互作网络 + +| 项目 | 说明 | +|---|---| +| 目的 | 利用 STRING 等已知数据库，为家族成员构建候选蛋白互作网络，寻找枢纽蛋白和功能模块。 | +| 输入 | 家族基因/蛋白 ID、物种信息；非模式植物通常需要到近缘模式物种的 ID 映射。 | +| 主要处理 | ID 解析和映射；查询数据库；过滤互作置信度；构建网络和节点/边属性。 | +| 输出 | `{prefix}.ppi_network.tsv`、ID 解析/映射表、网络图或 STRING 图片、报告和 Manifest。 | +| 质控 | 记录物种和映射来源；区分实验、数据库、文本挖掘等证据；跨物种映射产生的关系属于间接证据。 | + +## 19. `pep`：蛋白理化性质、跨膜区与亚细胞定位 + +| 项目 | 说明 | +|---|---| +| 目的 | 描述家族蛋白的长度、分子量、等电点、稳定性、疏水性、跨膜拓扑和可能的亚细胞定位。 | +| 输入 | 最终家族蛋白 FASTA。 | +| 主要处理 | ProtParam 计算理化性质；DeepTMHMM 或等效工具预测跨膜区；亚细胞定位工具预测细胞位置；汇总与绘图。 | +| 输出 | `{prefix}.pep_properties.pdf/png`、`pep/analysis/{prefix}_properties.tsv`、`protparam_full.tsv`；TMHMM 图和 summary；predicted topologies；`{prefix}.pep_subcellular.pdf/png`；`final/subcellular_data.tsv`；报告和 Manifest。 | +| 质控 | 输入蛋白不得含大面积非法字符；预测工具和版本需记录；跨膜与定位结果是计算预测，需避免写成实验结论。 | + +## 20. `expression`：家族表达模式 + +| 项目 | 说明 | +|---|---| +| 目的 | 查看家族成员在不同组织、时期、处理或基因型中的表达差异，筛选响应型或组织特异型成员。 | +| 两种模式 | ① standalone：直接上传 TPM/FPKM 表和样本表；② `rnaseq_auto`：复用本项目 RNA-seq quant、sample_qc、DE 等结果。 | +| 输入 | 家族成员 ID；表达矩阵；样本分组/组织/处理信息；可选 DE 结果。RNA-seq 自动模式优先读取成功 Manifest。 | +| 主要处理 | ID 交集；表达量变换与标准化；热图；基因相关性；表达分布；组织特异性；按树或亚家族排序。 | +| 输出 | 家族表达热图、基因相关性热图、表达分布/直方图、组织特异性结果、整理后的矩阵、报告和 Manifest。 | +| RNA-seq 矩阵选择 | 仅当 transcript-aware matrix bridge QC 通过且回退比例不超过 5% 时，优先使用 tximport/转录本感知 TPM；否则使用 featureCounts TPM。 | +| 质控 | 样本名必须与样本表一致；表达矩阵 ID 必须与家族 ID 可映射；FPKM/TPM 适合展示，不应直接替代 DESeq2 原始计数做差异检验。 | + +## 21. `combo`：系统树—motif—结构域—基因结构综合图 + +| 项目 | 说明 | +|---|---| +| 目的 | 将家族分类的主要证据放在同一张图中，直观比较“亲缘关系—motif—domain—基因结构”的一致性与差异。 | +| 输入 | 系统树 Newick、motif XML/表、结构域 `domtblout`/表、筛选后的 GFF3/基因结构表、确认的亚家族标签。 | +| 主要处理 | 统一成员 ID 和顺序；按树叶排列各面板；统一比例尺、颜色和标签；生成多面板图。 | +| 输出 | `combo.pdf`、`combo.png`、报告和 Manifest。 | +| 质控 | 每个面板的成员集合与顺序必须核对；缺失某类数据时应留空并说明，不能错位拼接。 | + +## 22. `gene_family_report`：基因家族正式报告 + +| 项目 | 说明 | +|---|---| +| 目的 | 将当前计划范围内成功的基因家族阶段整合成结构化报告，并给出结果覆盖度和证据链。 | +| 输入 | 同一 prefix、同一成员 scope 下各成功/部分成功阶段的 Manifest、报告、表格和图。 | +| 主要处理 | 校验成员范围和文件哈希；按流程组织方法、结果、图表与限制；区分已完成、未运行和失败阶段。 | +| 输出 | 正式 Markdown 报告、证据/模型数据、报告阶段 Manifest。 | +| 质控 | 报告只能引用本次有效 Manifest 的结果；不能把自动预测写成实验证实；覆盖不完整时明确标注 partial。 | + +--- + +# 第二部分：RNA-seq 分析 + +## 23. RNA-seq 完整流程总览 + +```mermaid +flowchart TD + R0[FASTQ + samples.yaml] --> R1[rnaseq_qc
fastp 质控] + R1 --> R2[rnaseq_align
HISAT2 比对] + R2 --> R3[rnaseq_quant
featureCounts/可选转录本定量] + R3 --> R4[rnaseq_sample_qc
PCA/相关性/聚类] + R3 --> R5[rnaseq_de
DESeq2 差异表达] + A0[蛋白 FASTA/已有注释] --> A1[annotation
GO/KEGG 注释] + R5 --> E1[GO 富集] + A1 --> E1 + R5 --> E2[KEGG 富集] + A1 --> E2 + R5 --> E3[FGSEA] + A1 --> E3 + R4 --> W[rnaseq_wgcna] + R5 --> C[rnaseq_candidate_family] + E1 --> C + E2 --> C + E3 --> C + W --> C + C --> RR[rnaseq_report] +``` + +## 24. RNA-seq 样本配置 + +典型字段如下，实际任务应通过界面或标准配置确认： + +```yaml +samples: + CK_1: + r1: /data/CK_1_R1.fastq.gz + r2: /data/CK_1_R2.fastq.gz + condition: CK + replicate: 1 + include: true + group_role: control + library_strandness: non_stranded + T_1: + r1: /data/T_1_R1.fastq.gz + r2: /data/T_1_R2.fastq.gz + condition: Treatment + replicate: 1 + include: true + group_role: compare + library_strandness: non_stranded +groups: + CK_1: CK + T_1: Treatment +experiment: + design: two_group + control_group: CK + treatment_groups: [Treatment] +compare: + - [Treatment, CK] +de: + strategy: auto + design: two_group + control_group: CK + pairwise_modes: [control] +``` + +配置原则： + +- `sample_id` 必须唯一，并贯穿 FASTQ、BAM、计数矩阵和样本表。 +- `include=false` 或 `group_role=ignore` 的样本保留来源记录，但不进入分组和差异分析。 +- `library_strandness` 必须明确为 `non_stranded`、`forward` 或 `reverse`，不能凭经验猜测。 +- 显式 `compare` 优先，格式为 `[比较组, 对照组]`；若未显式给出，系统才根据 `group_role` 构造比较。 +- 双组、多组或时间序列需要唯一且明确的 control/baseline；对照不明确时应停止确认。 +- 建议每组至少 3 个生物学重复；重复太少会明显降低离散度估计、差异分析和 WGCNA 的可靠性。 + +## 25. `rnaseq_qc`：原始数据质控与清洗 + +| 项目 | 说明 | +|---|---| +| 目的 | 去除接头、低质量碱基和过短 reads，检查测序质量，为比对提供 clean reads。 | +| 输入 | 双端 FASTQ/FASTQ.GZ、samples.yaml；质控工具和参数。 | +| 主要处理 | 默认使用 fastp；统计原始 reads；接头识别/去除；低质量过滤；长度过滤；输出清洗后成对 reads。只有用户明确要求时使用 cutadapt，不做静默切换。 | +| 输出 | 每个样本 clean R1/R2；fastp HTML、JSON 和日志；阶段报告和 Manifest。 | +| 关键指标 | raw reads、clean reads、保留率、Q20/Q30、GC、接头 reads、过滤 reads、平均读长。 | +| 质控 | R1/R2 必须成对且非空；样本集合应与配置一致；缺文件、工具失败或关键样本无有效 clean reads 时失败。部分样本失败只能标记 partial，不能把缺失样本继续当完整设计。 | +| 下一步 | clean FASTQ 进入 `rnaseq_align`。 | + +## 26. `rnaseq_align`：参考基因组比对 + +| 项目 | 说明 | +|---|---| +| 目的 | 将 clean reads 定位到参考基因组，生成排序并建立索引的 BAM。 | +| 输入 | QC 成功的 clean R1/R2、参考基因组 FASTA、完整 HISAT2 索引、samples.yaml 中明确的链特异性；可选 GTF/GFF 和已知剪接位点。 | +| 主要处理 | 检查参考和索引完整性；按样本运行 HISAT2；必要时使用已知剪接位点；通过 samtools 转换、排序和索引。 | +| 输出 | 每个样本 sorted BAM、BAI、HISAT2 比对日志；可选 `splice_sites.ss`；报告和 Manifest。 | +| 关键指标 | 总 reads、总体比对率、唯一比对率、多位置比对率、未比对率、成对比对情况。 | +| 质控 | 链特异性缺失或非法时失败；BAM 必须可读取且有索引；参考染色体名应与注释一致；明显偏低的比对率需排查物种、污染、参考版本或样本质量。 | +| 下一步 | sorted BAM 进入 `rnaseq_quant`。 | + +## 27. `rnaseq_quant`：基因/转录本表达定量 + +| 项目 | 说明 | +|---|---| +| 目的 | 从 BAM 和注释中统计每个基因的原始 reads 数，并生成适合展示的 FPKM/TPM；可选转录本层定量。 | +| 输入 | sorted BAM+BAI、GTF/GFF、每个样本的链特异性；可选 Salmon/Kallisto 结果或索引。 | +| 主要处理 | featureCounts 按基因计数；混合链特异性样本会分开运行再按基因合并；计算基因长度、FPKM 和 TPM；可选通过 Salmon/Kallisto 与 tximport 汇总转录本结果；执行 matrix bridge QC。 | +| 主要输出 | `{prefix}.count_matrix.txt`、`.fpkm_matrix.txt`、`.tpm_matrix.txt`、featureCounts summary 和 command；可选 gene-level tximport TPM/counts、length matrix、tx2gene；matrix bridge QC；报告和 Manifest。 | +| 质控 | 计数矩阵必须为非负整数；列名与样本集合严格一致；注释 feature/attribute 设置必须与文件匹配；分配率异常需检查注释版本和链特异性。 | +| 下游使用规则 | DESeq2 默认只使用 featureCounts **原始 count matrix**；TPM/FPKM用于展示。表达展示仅在桥接 QC 通过且转录本映射回退比例≤5%时优先使用转录本感知 TPM，否则使用 featureCounts TPM。 | + +## 28. `rnaseq_sample_qc`：样本层质量评估 + +| 项目 | 说明 | +|---|---| +| 目的 | 在做差异分析前判断重复样本是否聚在一起，发现离群样本、批次效应、标签错误或整体表达异常。 | +| 输入 | 原始 count matrix；可选样本分组表。 | +| 主要处理 | DESeq2 VST，`blind=TRUE`、`design=~1`；PCA；样本相关性；样本距离；层次聚类。 | +| 输出 | VST 矩阵、PCA 坐标和解释率、相关系数矩阵、距离矩阵、PCA/相关性/距离/聚类图、报告和 Manifest。 | +| 质控 | 至少 2 个样本；矩阵必须为数值；样本列名不能重复；缺包或变换失败时失败。样本离群只作为检查信号，是否排除必须结合原始质控和研究设计确认。 | +| 下一步 | VST 矩阵供 WGCNA；QC 结论用于确认是否进入 DE。 | + +## 29. `rnaseq_de`：差异表达分析 + +| 项目 | 说明 | +|---|---| +| 目的 | 比较指定处理/组织/时期与对照，识别显著上调和下调基因。 | +| 输入 | featureCounts 原始 count matrix、samples.yaml、显式比较关系；阈值参数。 | +| 主要处理 | 样本与分组核对；建立 DESeq2 设计；估计 size factor 和离散度；拟合模型；提取每个 comparison 的结果；进行多重检验校正；绘制 MA/火山图。 | +| 默认阈值 | `|log2FoldChange| ≥ 1` 且 `padj < 0.05`。可根据项目预先调整，但报告必须记录实际阈值。 | +| 输出 | 每个比较的完整 DE 结果表、显著 DEG 表、上调/下调列表、DESeq2 中间对象、MA/火山图、报告和 Manifest。 | +| 关键指标 | 计划比较、实际比较、参与样本数、有效对照组、可检验基因数、显著基因数、上调数、下调数。 | +| 质控 | 样本名与矩阵列不一致、比较非法、对照不唯一、某组无有效样本或 DESeq2 失败时停止；不得自动猜测对照方向，因为方向会直接改变 log2FC 正负。 | +| 下一步 | DEG 进入 GO/KEGG 富集和候选家族；完整 DE 排名进入 FGSEA。 | + +## 30. `annotation`：GO 与 KEGG 功能注释 + +| 项目 | 说明 | +|---|---| +| 目的 | 建立基因到 GO term 和 KEGG KO 的映射，为富集分析提供背景注释。 | +| 输入 | 全蛋白组或目标蛋白 FASTA；也可导入已有 gene2go、gene2ko、eggNOG 结果；本地数据库和工具路径。 | +| 主要处理 | GO 可通过 InterProScan；KEGG 可通过 KofamScan；解析原始输出；进行基因/蛋白 ID 映射；生成标准 gene2go/gene2ko。 | +| 输出 | gene2go、可选 gene2ko、注释 summary、ID map、原始结果、报告、Manifest 和重现信息。 | +| 质控 | GO 与 KEGG 是两个独立 scope：GO 成功不代表 KEGG 已完成；数据库和版本应记录；映射率过低需检查 ID 或蛋白集。 | +| 下一步 | gene2go 进入 GO 富集；gene2ko 进入 KEGG 富集和 FGSEA。 | + +## 31. `rnaseq_go_enrich`：GO 富集 + +| 项目 | 说明 | +|---|---| +| 目的 | 判断差异基因是否在某些生物过程、分子功能或细胞组分中出现过度代表。 | +| 输入 | 某个 comparison 的 DEG 前景集；gene2go/OrgDb/自动注释；背景基因集。 | +| 主要处理 | 构造前景和背景；GO 类别统计检验；多重校正；按 BP/CC/MF 汇总。 | +| 输出 | 每个比较的 `GO_enrich.txt`、富集图、报告、Manifest、命令/环境/校验和。 | +| 质控 | 前景和背景不能为空；前景应是背景子集；必须报告实际被注释的 DEG 数；无显著条目可以是合法结果，但空输入不是成功。 | + +## 32. `rnaseq_kegg_enrich`：KEGG 富集 + +| 项目 | 说明 | +|---|---| +| 目的 | 识别差异基因集中显著富集的代谢或信号通路。 | +| 输入 | DEG、gene2ko；若无 gene2ko，可在完整 Kofam 环境下自动生成。 | +| 主要处理 | 基因—KO—通路映射；构造背景；过度富集检验；多重校正；结果汇总和绘图。 | +| 输出 | 每个比较的 `KEGG_enrich.txt`、富集图、报告和 Manifest。 | +| 质控 | Kofam 环境不完整或 gene2ko 为空时失败；应报告 DEG 到 KO 的映射率；通路名称和物种来源要明确。 | + +## 33. `rnaseq_fgsea`：基于全基因排序的通路富集 + +| 项目 | 说明 | +|---|---| +| 目的 | 不只看超过阈值的 DEG，而是利用全部可检验基因的变化方向和强度，识别整体上调或下调的通路。 | +| 输入 | 完整 DE 结果，至少含 `gene_id` 和 `log2FoldChange`；gene2ko；本地 KEGG pathway cache；可选物种、comparison。 | +| 主要处理 | 按 log2FC 建立全基因排序；构建 KO/通路 gene set；植物相关通路过滤与物种/KO 回退；运行 fgsea；计算 NES 和 padj；绘制 top 通路。 | +| 默认参数 | `min_size=15`、`max_size=500`、`nperm=10000`、`seed=1`。 | +| 输出 | `fgsea_kegg.tsv`、summary、dotplot/top-enrichment PNG/PDF、报告和 Manifest。 | +| 质控 | 排名值需有限且基因 ID 唯一；gene set 过小/过大被过滤；缺少本地缓存、R 包或有效映射时失败。NES 正负代表通路在排序两端的方向。 | + +## 34. `rnaseq_wgcna`：加权基因共表达网络 + +| 项目 | 说明 | +|---|---| +| 目的 | 将表达模式相似的基因聚成模块，并寻找与处理、组织或连续性状相关的模块和枢纽候选基因。 | +| 输入 | 优先使用 `rnaseq_sample_qc` 的 VST 矩阵，也可使用合适的表达矩阵；可选样本性状/分组表。 | +| 主要处理 | 低表达/低变化过滤；选择 soft-threshold power；构建邻接和 TOM；动态剪切模块；相近模块合并；计算模块特征基因；与样本性状相关。 | +| 默认参数 | `min_mean=1`、`max_genes=8000`、`min_module_size=30`、`merge_cut_height=0.25`。 | +| 输出 | soft-threshold 结果、gene-module 表、模块大小、module eigengenes、trait correlation/p/padj、模块基因列表、诊断和相关性图、报告和 Manifest。 | +| 质控 | WGCNA 对样本量很敏感；样本少时只作为探索性结果；必须检查离群样本、soft power 拟合、模块规模和性状编码；相关不代表因果。 | + +## 35. `rnaseq_candidate_family`：候选基因与候选家族优先级 + +| 项目 | 说明 | +|---|---| +| 目的 | 综合差异表达、功能注释、富集、FGSEA、WGCNA 和研究问题，把大量基因缩小为可继续验证的候选基因/家族列表。 | +| 输入 | DEG、全蛋白组或可复用 InterPro 注释；可选 GO/KEGG 富集、FGSEA、WGCNA、研究背景和 RAG 文献证据。 | +| 主要处理 | 汇总每个基因的差异方向/显著性、功能域、通路、模块和研究上下文；构建证据矩阵；计算排序；划分推荐层级；生成实验验证建议。 | +| 输出 | candidate ranking、gene annotation、enrichment context、family enrichment、evidence matrix、recommendations、key genes、experiment context、RAG plan/references、报告和 Manifest。 | +| 质控 | 这是**候选优先级排序**，不是正式的全基因组家族鉴定；自动推荐必须保留证据来源；缺乏上下文时使用中性推荐，不能过度解释。 | +| 下一步 | 候选家族可送入第一部分的 `identify` 做完整家族鉴定；候选基因可与 GWAS 位点、启动子和 PPI 联合验证。 | + +## 36. `rnaseq_report`：RNA-seq 正式报告 + +| 项目 | 说明 | +|---|---| +| 目的 | 汇总同一任务中 RNA-seq、注释和功能分析阶段的有效证据，形成可交付报告。 | +| 输入 | 同一 prefix/scope 下成功或部分成功的 RNA-seq 与 annotation Manifest。 | +| 主要处理 | 验证输入哈希和阶段身份；提取指标、表格和图；按“数据质控—比对—定量—样本 QC—差异—功能—网络—候选”组织。 | +| 输出 | report Markdown、evidence JSON、model JSON、report Manifest。 | +| 质控 | 即使不调用模型也可生成确定性报告；“阶段成功”和“报告覆盖完整”分开记录；没有任何成功证据时不生成空报告。 | + +--- + +# 第三部分：WGS 全基因组重测序 + +## 37. WGS 完整流程总览 + +```mermaid +flowchart TD + W0[FASTQ + wgs_samples.yaml] --> W1[wgs_qc
fastp] + REF[参考基因组 FASTA] --> W2[wgs_index
FAI/DICT/BWA-MEM2] + W1 --> W3[wgs_align
BWA-MEM2] + W2 --> W3 + W3 --> W4[wgs_markdup
去重复] + W4 --> W5[wgs_bam_qc
比对与深度质控] + W5 --> W6[wgs_haplotype_calling
单样本 gVCF] + W2 --> W6 + W6 --> W7[wgs_gvcf_combine
合并 gVCF] + W7 --> W8[wgs_joint_genotyping
联合分型] + W8 --> W9[wgs_variant_select
SNP/INDEL 分离] + W9 --> W10[wgs_variant_filter
硬过滤] + W10 --> W11[wgs_impute
Beagle 填补/定相] + W11 --> W12[wgs_maf_missing_filter
MAF/缺失率过滤] + W12 --> V[最终群体 SNP VCF] +``` + +## 38. WGS 样本配置 + +当前标准 `wgs_samples.yaml` 只接受非空的 `samples` 映射；每个样本只允许 `r1` 和 `r2`： + +```yaml +samples: + S01: + r1: /data/S01_R1.fastq.gz + r2: /data/S01_R2.fastq.gz + S02: + r1: /data/S02_R1.fastq.gz + r2: /data/S02_R2.fastq.gz +``` + +严格规则： + +- 样本 ID 仅允许字母、数字、点、下划线和短横线，且首字符为字母或数字。 +- 每个样本必须恰好包含 `r1`、`r2`，不接受 `read1`、`read2`、`include`、`exclude`、`group_role`。 +- 路径可绝对或相对 YAML 所在目录；只接受 `.fastq`、`.fastq.gz`、`.fq`、`.fq.gz`。 +- 文件必须存在、非空，且 R1 与 R2 不能是同一文件。 +- 配置中的全部样本都是本次 WGS 必需样本；当前标准流程不允许中途静默丢样本。 + +## 39. `wgs_qc`：WGS 原始 reads 质控 + +| 项目 | 说明 | +|---|---| +| 目的 | 清理每个样本的双端重测序 reads，并建立队列级样本集合。 | +| 输入 | 确认后的 `wgs_samples.yaml`、所有 R1/R2 FASTQ。 | +| 主要处理 | fastp 成对质控、接头和低质量过滤、长度过滤；并发处理多个样本。 | +| 默认参数 | `qualified_quality_phred=5`、`unqualified_percent_limit=50`、`n_base_limit=10`、`length_required=50`；线程数和最大并发任务可配置。 | +| 输出 | 每样本 `{prefix}.wgs_qc.{sample_id}.clean.R1.fastq.gz`、`.clean.R2.fastq.gz`、`.fastp.html`、`.fastp.json`；报告和 Manifest。 | +| 质控 | 所有声明样本都必须成功，样本集合必须精确一致；任一样本缺失或输出无效时阶段失败，失败 Manifest 的 artifacts 为空。 | +| 下一步 | clean reads 与参考索引一起进入 `wgs_align`。 | + +## 40. `wgs_index`：参考基因组索引 + +| 项目 | 说明 | +|---|---| +| 目的 | 为比对、BAM 处理和 GATK 变异检测准备同一套参考基因组及其索引。 | +| 输入 | 显式指定的参考基因组 FASTA。 | +| 主要处理 | 复制为本任务标准参考；samtools faidx；GATK CreateSequenceDictionary；BWA-MEM2 index。 | +| 输出 | 标准 reference FASTA、`.fai`、`.dict`、BWA-MEM2 索引组件 `.0123/.amb/.ann/.bwt.2bit.64/.pac`；报告和 Manifest。 | +| 质控 | contig 名必须唯一且可解析；FASTA、FAI、DICT 和 BWA 索引必须属于同一参考文件；参考哈希形成 `wgs_reference_id`。 | +| 下一步 | 提供给 align、HaplotypeCaller、CombineGVCFs 和 GenotypeGVCFs。 | + +## 41. `wgs_align`：BWA-MEM2 比对 + +| 项目 | 说明 | +|---|---| +| 目的 | 将每个 WGS 样本的 clean reads 比对到参考基因组并生成坐标排序 BAM。 | +| 输入 | `wgs_qc` clean reads；`wgs_index` 参考 FASTA 和完整索引；样本 ID。 | +| 主要处理 | BWA-MEM2 比对；设置 read group；samtools 转换和坐标排序；建立 BAI。默认平台为 ILLUMINA。 | +| 输出 | 每样本 sorted BAM、BAI、比对日志；报告和 Manifest。 | +| 质控 | BAM read group 的 `SM` 必须等于样本 ID；输入和输出样本集合完全一致；BAM 可读取、有索引、参考序列一致。 | +| 下一步 | sorted BAM 进入 `wgs_markdup`。 | + +## 42. `wgs_markdup`：重复 reads 标记/移除 + +| 项目 | 说明 | +|---|---| +| 目的 | 降低 PCR/光学重复对深度和变异证据的偏倚。 | +| 输入 | 每个样本 sorted BAM+BAI。 | +| 主要处理 | GATK duplicate 处理；当前默认 `remove_duplicates=true`；创建索引和重复统计。 | +| 输出 | 每样本 markdup BAM、BAI、metrics、日志、报告和 Manifest。 | +| 质控 | 处理前后样本 ID 不变；输出 BAM 需可读并索引；重复比例异常高可能提示文库复杂度不足或测序重复。 | +| 下一步 | markdup BAM 同时进入 BAM QC 和变异检测。 | + +## 43. `wgs_bam_qc`：比对后质量评估 + +| 项目 | 说明 | +|---|---| +| 目的 | 评估每个样本的比对情况、覆盖深度和参考基因组覆盖比例，判断是否适合继续变异检测。 | +| 输入 | markdup BAM+BAI、参考 FAI。 | +| 主要处理 | samtools flagstat；逐位点深度统计；队列汇总；默认覆盖阈值为 ≥1×。 | +| 输出 | 每样本 flagstat、depth 文件、队列 summary、报告和 Manifest。 | +| 关键指标 | 总 reads、mapped reads/rate、properly paired、平均深度、中位深度、≥1×覆盖比例等。 | +| 质控 | 覆盖比例分母包含参考基因组的零覆盖位置；低覆盖样本不会被自动删除，而是报告并要求研究者决定；标准 WGS 链路仍保持完整样本集合。 | +| 下一步 | 成功的样本 BAM 进入 `wgs_haplotype_calling`。 | + +## 44. `wgs_haplotype_calling`：单样本 gVCF + +| 项目 | 说明 | +|---|---| +| 目的 | 对每个样本独立发现候选变异，并以 gVCF 形式记录变异和参考区段证据，为联合分型做准备。 | +| 输入 | markdup BAM+BAI、成功 BAM QC、reference FASTA/FAI/DICT。 | +| 主要处理 | GATK HaplotypeCaller，`ERC=GVCF`；按样本调用并索引。 | +| 输出 | 每样本 `.g.vcf.gz`、`.tbi`、日志、报告和 Manifest。 | +| 质控 | gVCF 内样本 ID 必须与任务样本 ID 一致；参考 contig 和顺序必须匹配；每个声明样本都要有可索引 gVCF。 | +| 下一步 | 全部 gVCF 进入 combine。 | + +## 45. `wgs_gvcf_combine`：队列 gVCF 合并 + +| 项目 | 说明 | +|---|---| +| 目的 | 将所有单样本 gVCF 合并为一个队列级 gVCF，供一次性联合分型。 | +| 输入 | 所有样本 gVCF+TBI、reference FASTA/FAI/DICT。 | +| 主要处理 | GATK CombineGVCFs；合并样本 header 和参考区段记录。 | +| 输出 | cohort combined gVCF、TBI、日志、报告和 Manifest。 | +| 质控 | 每个样本必须参与；header 中样本集合与 `wgs_samples.yaml` 完全一致；缺一不可静默继续。 | +| 下一步 | combined gVCF 进入 joint genotyping。 | + +## 46. `wgs_joint_genotyping`：联合分型 + +| 项目 | 说明 | +|---|---| +| 目的 | 在同一位点框架下联合判断所有样本的基因型，生成群体级原始 VCF。 | +| 输入 | combined gVCF、reference FASTA/FAI/DICT。 | +| 主要处理 | GATK GenotypeGVCFs；统一计算位点与样本基因型。 | +| 默认参数 | 最大替代等位基因数 6；calling confidence 10；默认不输出纯参考的非变异位点。 | +| 输出 | joint cohort VCF、TBI、summary、日志、报告和 Manifest。 | +| 质控 | 样本集合精确保存；VCF 必须可索引、contig 与参考一致；必须有可用变异记录。 | +| 下一步 | 原始 VCF 进入 SNP/INDEL 分离。 | + +## 47. `wgs_variant_select`：SNP 与 INDEL 分离 + +| 项目 | 说明 | +|---|---| +| 目的 | 将联合 VCF 中的 SNP 和 INDEL 分开，为不同硬过滤规则做准备。 | +| 输入 | joint VCF+TBI。 | +| 主要处理 | 默认 GATK SelectVariants 分别选择 SNP 和 INDEL；只有 GATK 失败且用户确认后才允许用 bcftools，不能自动静默回退。 | +| 输出 | selected SNP VCF+TBI、selected INDEL VCF+TBI、报告和 Manifest。 | +| 质控 | 两类 VCF 必须保持相同样本集合；记录实际使用的方法；变异类型必须与文件角色一致。 | +| 下一步 | 分别进入硬过滤。 | + +## 48. `wgs_variant_filter`：变异硬过滤 + +| 项目 | 说明 | +|---|---| +| 目的 | 根据位点质量、链偏倚、比对质量和位置偏倚去除低可信 SNP/INDEL。 | +| 输入 | selected SNP VCF、selected INDEL VCF 及索引。 | +| 主要处理 | 按 SNP 和 INDEL 的不同阈值标记过滤原因；提取 PASS 记录；创建索引。 | +| SNP 默认阈值 | QD<2、QUAL<30、SOR>3、FS>60、MQ<40、MQRankSum<-12.5、ReadPosRankSum<-8 时标记。 | +| INDEL 默认阈值 | QD<2、QUAL<30、FS>200、ReadPosRankSum<-20 时标记。 | +| 输出 | filtered SNP/INDEL VCF、PASS SNP/INDEL VCF、全部 TBI、统计/报告、Manifest。 | +| 质控 | 必须保留 FILTER 原因和过滤前后数量；报告实际阈值和方法；参考、样本和阶段身份不能改变。 | +| 下一步 | PASS SNP 进入 Beagle 填补；PASS INDEL 可作为独立结果保存。 | + +## 49. `wgs_impute`：基因型填补与定相 + +| 项目 | 说明 | +|---|---| +| 目的 | 利用群体单倍型信息填补缺失基因型，并对基因型进行定相，提高下游群体分析和 GWAS 的可用性。 | +| 输入 | PASS SNP VCF+TBI。 | +| 主要处理 | Beagle imputation/phasing；按窗口处理；索引输出。 | +| 默认参数 | memory 8 GB、window 40、overlap 2、`ne=100000`、iterations 12、phase_states 280、imp_states 1600、cluster 0.005。 | +| 输出 | imputed SNP VCF+TBI、Beagle 日志、统计、报告和 Manifest。 | +| 质控 | 输出样本集合必须保持一致；检查填补前后缺失率、位点数和基因型完整性；小群体或群体结构复杂时需谨慎解释填补质量。 | +| 下一步 | 进入 MAF/缺失率终筛。 | + +## 50. `wgs_maf_missing_filter`：MAF 与缺失率终筛 + +| 项目 | 说明 | +|---|---| +| 目的 | 去除过于稀有或缺失过多的 SNP，得到适合群体遗传和 GWAS 的最终基因型集。 | +| 输入 | imputed SNP VCF+TBI。 | +| 主要处理 | bcftools 按等位基因频率和位点缺失率过滤；逐记录复核；生成统计。 | +| 默认参数 | `MAF ≥ 0.05`；`max_missing=0.5`，即位点允许的最大缺失比例 50%。具体研究可预先调整。 | +| 输出 | final SNP VCF+TBI、位点统计、过滤前后统计、报告和 Manifest。 | +| 质控 | 输出至少包含一个可用于 GWAS 的变异；样本集合不变；每条记录满足阈值；最终 VCF 是连接 WGS 与群体/GWAS 的权威基因型输入。 | + +--- + +# 第四部分：群体遗传分析与 GWAS + +## 51. 群体遗传分析总览 + +```mermaid +flowchart TD + V[最终群体 SNP VCF] --> P1[pop_pca] + V --> P2[pop_phylogeny] + V --> P3[pop_structure] + V --> P4[pop_ld_decay] + V --> P5[pop_snp_density] + P1 --> G[GWAS 协变量/解释] + P3 --> G +``` + +群体遗传模块既可使用 WGS 最终 VCF，也可使用用户提供并经过验证的群体 VCF。五个阶段可以按研究目的并行执行。 + +## 52. `pop_pca`：群体主成分分析 + +| 项目 | 说明 | +|---|---| +| 目的 | 用少数主成分概括样本间主要遗传差异，发现群体分层、离群样本和潜在批次结构，并为 GWAS 提供 PCA 协变量。 | +| 输入 | population VCF；可选 sample group 表。 | +| 主要处理 | 保留双等位 SNP；PLINK LD pruning，默认窗口 50、步长 10、`r²=0.2`；对剪枝后的标记做 PCA。PC 数不超过 `min(样本数-1, 剪枝后标记数)`。 | +| 输出 | PC scores、eigenvalues、variance explained、剪枝标记列表、PCA 图、scree 图、TASSEL covariate 表、报告、环境和 Manifest。 | +| 质控 | 样本 ID 与 VCF 一致；剪枝后需有足够标记；解释率和样本数必须可核对；PCA 是结构描述，不是群体标签的唯一判据。 | + +## 53. `pop_phylogeny`：群体样本树 + +| 项目 | 说明 | +|---|---| +| 目的 | 根据群体 SNP 描述样本间遗传距离和聚类关系。 | +| 输入 | population VCF。 | +| 主要处理 | 提取双等位 SNP；TASSEL 生成 Phylip interleaved alignment；FastTree `-nt -gtr` 建树；R/ape 绘制矩形和圆形树。 | +| 输出 | SNP alignment、Newick 树、矩形树图、圆形树图、QC、报告、环境和 Manifest。 | +| 质控 | 树叶必须与 VCF 样本完全对应；缺失和标记数要报告；这是基于筛选 SNP 的探索性样本树，不等同于经过 ascertainment correction 的发表级物种系统树。 | + +## 54. `pop_structure`：群体结构/祖先成分 + +| 项目 | 说明 | +|---|---| +| 目的 | 估计每个样本在不同祖先群体中的成分比例，为群体分层解释和 GWAS 协变量提供依据。 | +| 输入 | population VCF。 | +| 主要处理 | 双等位 SNP 和 LD pruning；生成 PLINK BED/BIM/FAM；ADMIXTURE 对 `K=2..10` 做 5 折交叉验证，固定随机种子，且 K 必须小于样本数。 | +| 输出 | 每个 K 的 Q/P/log、CV error 表、best-K 结果和结构图、报告、环境和 Manifest。 | +| 质控 | 以最小 CV error 选择 best K；并列时选较小 K；best K 是数据驱动建议，不等于用户最终用于 GWAS 的 K，正式协变量 K 应明确确认。 | + +## 55. `pop_ld_decay`：连锁不平衡衰减 + +| 项目 | 说明 | +|---|---| +| 目的 | 估计 SNP 间连锁不平衡随物理距离增加的衰减速度，为群体历史、标记密度和候选区间范围提供参考。 | +| 输入 | population VCF；可选样本分组表。 | +| 主要处理 | PopLDdecay 计算成对 r²；按距离分箱；可对全部样本和各分组分别绘制衰减曲线。 | +| 默认参数 | `MAF=0.05`、`Miss=0.2`、最大距离 500 kb；近距离/远距离分箱 10/500；分组最少 3 个样本。 | +| 输出 | 全体与各组 LD statistics、衰减曲线、组合图、group QC、报告、环境和 Manifest。 | +| 质控 | 分组表非法时默认失败而非忽略；组内样本太少不应计算；需记录 r² 计算与衰减定义。 | + +## 56. `pop_snp_density`：SNP 密度 + +| 项目 | 说明 | +|---|---| +| 目的 | 描述变异在染色体上的空间分布，寻找高密度、低密度或可能受过滤/组装影响的区域。 | +| 输入 | population VCF；染色体选择策略。 | +| 主要处理 | 统计染色体坐标；按 1 Mb 和 500 kb 窗口计数；生成染色体密度图和 500 kb 热图。 | +| 染色体策略 | 可选 all、显式列表或 rice1-12；系统不会对非水稻数据隐式套用 1–12 号染色体过滤。 | +| 输出 | SNP coordinates、1 Mb/500 kb 密度表与图、500 kb heatmap、chromosome counts、报告、环境和 Manifest。 | +| 质控 | 所有窗口计数之和必须等于所选染色体的 SNP 总数；窗口边界和染色体长度需一致。 | + +## 57. GWAS 三条模型路线总览 + +```mermaid +flowchart TD + V[VCF] --> F[gwas_fmt
格式标准化] + PH[表型表] --> T[gwas_trait] + Q[ADMIXTURE Q + FAM] --> C[gwas_cov] + F --> K[gwas_kinship] + + F --> M[gwas_mlm
MLM + CMLM] + K --> M + T --> M + C --> M + M --> MP[gwas_plot] + + F --> E[gwas_emmax] + T --> E + E --> EP[gwas_emplot] + + F --> G[gwas_glm
GLM + Q-GLM] + T --> G + C --> G +``` + +三条路线的定位： + +- **MLM/CMLM**：显式使用亲缘矩阵 K 和群体结构 Q，控制混杂较完整，是当前标准主路线。 +- **EMMAX**：从 VCF 内部构建关系矩阵和 PCA 协变量，计算效率高，适合大批量标记。 +- **GLM/Q-GLM**：结构简单，便于比较有无 Q 协变量的影响，但不能像混合模型那样充分控制亲缘关系，假阳性风险通常更高。 + +## 58. GWAS 输入表格式建议 + +### 58.1 表型表 + +至少包括唯一样本 ID 和一个非恒定、有限数值性状，例如： + +```text +sample_id\theight_cm +S01\t102.4 +S02\t95.8 +S03\t110.1 +``` + +- 样本 ID 不得重复。 +- 性状不能全相同，不能全部缺失或含无法解析的字符。 +- 表型样本可为 VCF 样本的子集；真正建模时取表型、基因型、协变量/亲缘矩阵的严格交集并保持确定顺序。 +- 多环境或多重复表型建议在进入 GWAS 前明确使用原始值、均值、BLUE 或 BLUP，当前项目不会替研究者猜测表型统计口径。 + +### 58.2 ADMIXTURE Q 与 FAM + +- Q 文件每行对应一个样本，每列为一个祖先成分比例。 +- FAM 提供与 Q 文件行完全一致的样本 ID 顺序。 +- 至少需要两列 Q；建模时省略最后一列 Q，避免所有成分和为 1 导致完全共线。 + +## 59. `gwas_fmt`：基因型格式标准化 + +| 项目 | 说明 | +|---|---| +| 目的 | 将 VCF 规范化并转换为 TASSEL 可使用的 HapMap，同时验证样本和位点守恒。此阶段只做格式转换，不执行关联分析。 | +| 输入 | 可读取、含样本和位点的 VCF/VCF.GZ；TASSEL 环境；不支持标记的处理策略。 | +| 主要处理 | VCF 排序和校验；转换 HapMap；对不支持的标记默认 mask，或按显式策略失败。 | +| 输出 | sorted VCF、HapMap、环境记录、报告和 Manifest。 | +| 质控 | 样本数、样本顺序和变异数量需守恒；`genotype_scope_id` 由输入内容哈希建立；不能混入另一个 VCF 的下游文件。 | +| 下一步 | sorted VCF 进入 kinship、MLM、EMMAX、GLM；HapMap 供 TASSEL 模型。 | + +## 60. `gwas_kinship`：亲缘矩阵 + +| 项目 | 说明 | +|---|---| +| 目的 | 量化样本两两遗传相似度，为混合线性模型控制亲缘关系。 | +| 输入 | `gwas_fmt` sorted VCF。 | +| 主要处理 | TASSEL centered IBS 计算 kinship。 | +| 输出 | `{stem}.kinship.txt`、环境、报告和 Manifest。 | +| 质控 | 矩阵必须为方阵、对称、半正定；行列样本顺序必须与 VCF 精确一致。 | +| 下一步 | 进入 `gwas_mlm`。 | + +## 61. `gwas_trait`：表型标准化 + +| 项目 | 说明 | +|---|---| +| 目的 | 将用户表型表整理为 GWAS 模型可稳定读取的 trait table，并建立表型 scope。 | +| 输入 | 外部 `phenotype_table`。 | +| 主要处理 | 识别样本 ID 和性状列；检查重复、缺失、有限数值和变异；标准化表头/格式；记录输入哈希。 | +| 输出 | `{stem}.trait.txt`、环境、报告和 Manifest。 | +| 质控 | 性状必须非恒定且有足够有效样本；样本 ID 唯一；输入来源不明确或有多个候选时要求确认。 | +| 下一步 | trait table 可进入 MLM、EMMAX、GLM。 | + +## 62. `gwas_cov`：群体结构协变量 + +| 项目 | 说明 | +|---|---| +| 目的 | 将 ADMIXTURE Q 矩阵与样本 ID 对齐，生成 GWAS 可使用的群体结构协变量。 | +| 输入 | 显式 ADMIXTURE Q 文件、对应 sample FAM 文件。 | +| 主要处理 | 读取 FAM 样本顺序；为 Q 列命名；检查每行比例与数值；输出带样本 ID 的协变量表。 | +| 输出 | `{stem}.covariates.txt`、环境、报告和 Manifest。 | +| 质控 | Q 与 FAM 行数和顺序必须一致；至少两列 Q；模型中去掉最后一列以避免共线；所选 K 应由用户明确。 | +| 下一步 | covariate table 进入 MLM/CMLM 和 GLM/Q-GLM。 | + +## 63. `gwas_mlm`：TASSEL MLM 与 CMLM + +| 项目 | 说明 | +|---|---| +| 目的 | 同时控制群体结构 Q 和亲缘关系 K，检测 SNP 与性状的关联，并比较 MLM 与压缩混合模型 CMLM。 | +| 前置阶段 | `gwas_fmt`、`gwas_kinship`、`gwas_trait`、`gwas_cov` 全部成功且 scope 相容。 | +| 输入 | sorted VCF/HapMap、kinship、trait、covariate。 | +| 主要处理 | 样本交集与顺序校验；TASSEL P3D MLM compression=None；TASSEL CMLM compression=Optimum；提取可估计标记的效应和 p 值。两种模型均为标准成功要求。 | +| 输出 | `mlm_full`、`cmlm_full`、`mlm_pvalues`、`cmlm_pvalues`、日志、报告、环境和 Manifest。标准 p 文件表头为 `chr,pos,p`。 | +| 质控 | 四类输入身份和哈希一致；样本顺序严格对齐；剔除无法估计的标记并记录原因；两种模型必须正常结束且 p 值在 0–1。 | +| 下一步 | p-value 文件进入 `gwas_plot`。 | + +## 64. `gwas_plot`：MLM/CMLM Manhattan 与 QQ 图 + +| 项目 | 说明 | +|---|---| +| 目的 | 展示全基因组关联峰，并用 QQ 图检查统计量膨胀或模型校正情况。 | +| 输入 | MLM 和 CMLM 标准 p-value 文件、权威 tested marker 数。 | +| 主要处理 | 染色体/contig 确定性排序；计算累积坐标；Bonferroni 阈值；绘制 Manhattan 和 QQ。 | +| 输出 | MLM Manhattan、MLM QQ、CMLM Manhattan、CMLM QQ 共四张 PNG；报告、环境和 Manifest。 | +| 默认阈值 | Bonferroni `0.05 / 实际检验标记数`。 | +| 质控 | 阈值分母必须来自模型实际 tested marker 数，而不是原始 VCF 位点数；p 值必须有效；contig 映射规则可追溯。 | + +## 65. `gwas_emmax`：EMMAX 混合模型 + +| 项目 | 说明 | +|---|---| +| 目的 | 通过高效混合模型控制样本关系，执行大规模 SNP—性状关联分析。 | +| 前置阶段 | `gwas_fmt` 和 `gwas_trait`。 | +| 输入 | sorted VCF、trait table。 | +| 主要处理 | 从 VCF 内部生成关系矩阵和 PCA 协变量；对齐表型样本；运行 EMMAX；过滤无法估计标记；标准化结果。 | +| 输出 | `emmax_full`，表头 `SNPID,Chr,Pos,beta,SE,p`；`emmax_pvalues`，表头 `chr,pos,p`；报告、环境和 Manifest。 | +| 质控 | EMMAX 必须以成功状态退出；样本 ID/顺序一致；p 值合法；只在工具支持的三字段输出场景允许 SE 为字面量 `NA`，不能无依据补造标准误。 | +| 下一步 | `emmax_pvalues` 进入 `gwas_emplot`。 | + +## 66. `gwas_emplot`：EMMAX Manhattan 与 QQ 图 + +| 项目 | 说明 | +|---|---| +| 目的 | 展示 EMMAX 关联峰和 p 值分布，辅助判断显著位点及系统性膨胀/过度校正。 | +| 输入 | EMMAX 标准 p-value 文件、模型权威 tested count。 | +| 主要处理 | 染色体排序、累积坐标、Bonferroni 阈值、Manhattan/QQ 绘图。 | +| 输出 | EMMAX Manhattan PNG、QQ PNG、报告、环境和 Manifest。 | +| 质控 | 与模型的 source scope 精确一致；tested count 不能重新从其他 VCF 推算。 | + +## 67. `gwas_glm`：TASSEL GLM 与 Q-GLM + +| 项目 | 说明 | +|---|---| +| 目的 | 比较不加群体结构的普通 GLM 与加入 Q 协变量的结构校正 GLM，作为关联结果的补充或模型对照。 | +| 前置阶段 | `gwas_fmt`、`gwas_trait`、`gwas_cov`。 | +| 输入 | sorted VCF/HapMap、trait table、covariate table。 | +| 主要处理 | 严格对齐样本；运行 unstructured GLM；运行 Q-structured GLM；结构模型省略最后一列 Q；标准化完整结果和 p 值。 | +| 输出 | `glm_full`、`glm_structure_full`、`glm_pvalues`、`glm_structure_pvalues`、日志、报告、环境和 Manifest。 | +| 质控 | 两个模型的输入样本与标记 scope 一致；Q 列避免完全共线；不做未声明的 permutations；无 kinship 的 GLM 对亲缘混杂较敏感，需结合 QQ 图和混合模型判断。 | + +## 68. GWAS 结果解释与候选区域衔接 + +当前已有阶段完成到模型结果与标准图。后续解释建议按以下顺序进行： + +1. 核对 QQ 图：严重膨胀可能来自群体结构、亲缘、批次或表型分布问题。 +2. 明确显著性阈值：Bonferroni 严格阈值与建议性阈值应分开呈现。 +3. 将显著 SNP 映射到参考基因组坐标。 +4. 依据 LD 衰减距离或局部 LD 确定候选区间，不能机械地固定上下游长度。 +5. 用 GFF 注释区间内或邻近基因。 +6. 联合 RNA-seq DEG、WGCNA module、功能注释、基因家族、启动子和 PPI 证据排序。 +7. 在独立群体、不同环境或实验中验证候选位点和候选基因。 + +--- + +# 第五部分：跨流程输入输出衔接 + +## 69. WGS → 群体遗传 → GWAS + +| 上游输出 | 下游输入 | 作用 | +|---|---|---| +| `wgs_maf_missing_filter` final SNP VCF | pop_pca/pop_structure/pop_phylogeny/pop_ld_decay/pop_snp_density | 描述群体结构、亲缘、LD 和变异分布 | +| final SNP VCF | `gwas_fmt` | 生成 GWAS 标准基因型 scope | +| pop_structure 的 Q + FAM | `gwas_cov` | 生成群体结构协变量 | +| pop_pca 的 PC covariate | 研究者确认后用于适用模型 | 用主成分控制结构；当前标准三分支各有明确协变量策略，不应随意混用 | +| `gwas_fmt` sorted VCF | kinship/MLM/EMMAX/GLM | 保证三条模型路线使用同一基因型集合 | + +## 70. RNA-seq → 基因家族 + +| 上游输出 | 下游输入 | 作用 | +|---|---|---| +| quant TPM/FPKM | family `expression` | 展示家族成员表达模式 | +| DE gene list | candidate_family | 找到响应处理或性状的候选基因 | +| candidate family | gene family `identify` | 对候选家族做正式全蛋白组鉴定 | +| WGCNA modules | candidate_family / family expression | 寻找与性状相关的家族成员和枢纽基因 | +| GO/KEGG/FGSEA | candidate_family | 为候选基因提供功能和通路证据 | + +## 71. GWAS + RNA-seq + 基因家族联合分析 + +```mermaid +flowchart LR + G[GWAS 显著位点/候选区间] --> J[候选基因交集与排序] + D[RNA-seq DEG] --> J + W[WGCNA 性状相关模块] --> J + F[家族成员/系统树/结构域] --> J + T[TFBS/PPI/功能注释] --> J + J --> H[高优先级候选基因] + H --> V[独立群体验证/表达验证/功能实验] +``` + +联合证据的合理表述： + +- “位于显著关联区间且在处理下差异表达”属于多数据支持的候选。 +- “同时位于性状相关 WGCNA 模块并具有相关功能域”可进一步提高优先级。 +- 以上仍属于计算证据链，不能直接写成“该基因决定性状”；最终需要遗传、分子或生理实验验证。 + +--- + +# 第六部分：常见质控清单 + +## 72. 所有流程通用检查 + +- [ ] 项目 prefix、物种、参考版本和样本批次已明确。 +- [ ] 输入文件非空、格式可解析、压缩文件完整。 +- [ ] 样本/基因/蛋白 ID 唯一，跨文件映射规则明确。 +- [ ] 上游 Manifest 状态为 success 或被明确允许的 partial。 +- [ ] 输入路径、大小和 SHA256 与 Manifest 一致。 +- [ ] 软件、数据库版本和实际参数已记录。 +- [ ] 输出文件非空，核心表的行列数与指标一致。 +- [ ] 图形数据来自正式表，而不是手工修改后的临时文件。 +- [ ] 自动预测、数据库推断和实验结论在报告中明确区分。 + +## 73. 基因家族专项检查 + +- [ ] Pfam/HMM 模型名称和 accession 明确，模型未缺失。 +- [ ] 参考家族蛋白被显式声明，未仅凭文件名猜测。 +- [ ] 最终成员 ID、蛋白 FASTA、树叶和各分析面板成员一致。 +- [ ] 多转录本去冗余规则已说明。 +- [ ] 启动子提取考虑正负链和染色体边界。 +- [ ] 跨物种分析未混淆物种命名空间。 +- [ ] 亚家族分类和 TFBS/PPI 等预测标记为候选证据。 + +## 74. RNA-seq 专项检查 + +- [ ] R1/R2 配对正确，样本名和分组一致。 +- [ ] 链特异性明确且 featureCounts 参数对应。 +- [ ] 参考基因组和注释来自相同版本，染色体命名一致。 +- [ ] 样本 PCA、相关性和比对/分配率已检查。 +- [ ] 差异方向 compare vs control 已明确。 +- [ ] DE 使用原始 counts，不用 TPM/FPKM 替代。 +- [ ] 富集分析背景集是本次实际可检验/有注释的基因。 +- [ ] WGCNA 结合样本量限制做探索性解释。 + +## 75. WGS 专项检查 + +- [ ] `wgs_samples.yaml` 只包含已确认的完整样本集合。 +- [ ] 参考 FASTA、FAI、DICT、BWA 索引完全一致。 +- [ ] 每个阶段样本集合、样本顺序和 `SM` 标签守恒。 +- [ ] BAM 比对率、深度、覆盖率和重复比例已查看。 +- [ ] 每个单样本 gVCF 都成功后才做 cohort combine。 +- [ ] SNP/INDEL 使用各自过滤阈值，报告过滤前后数量。 +- [ ] 填补前后缺失率和位点数已核对。 +- [ ] 最终 GWAS VCF 中至少有一个有效位点，且样本不丢失。 + +## 76. GWAS 专项检查 + +- [ ] 基因型、表型、Q、K 的样本 ID 和顺序严格一致。 +- [ ] 表型非恒定，有足够有效样本，统计口径明确。 +- [ ] Q 的 K 值经过确认，最后一列在模型中正确省略。 +- [ ] 亲缘矩阵方正、对称、半正定。 +- [ ] tested marker 数来自实际模型输出。 +- [ ] Manhattan 阈值分母与实际 tested marker 数一致。 +- [ ] QQ 图用于判断膨胀/过度校正，不只看最小 p 值。 +- [ ] 候选区间根据 LD 或明确规则确定，并记录参考基因组版本。 + +--- + +# 第七部分：交付结果建议 + +## 77. 推荐的最终交付目录 + +```text +project_delivery/ +├─ 00_input_and_design/ +│ ├─ sample_information.tsv +│ ├─ comparisons.tsv +│ ├─ phenotype.tsv +│ └─ parameter_summary.md +├─ 01_gene_family/ +│ ├─ members/ +│ ├─ phylogeny/ +│ ├─ motif_domain_structure/ +│ ├─ promoter_tfbs/ +│ ├─ synteny_karyotype_kaks/ +│ ├─ ppi_pep_expression/ +│ └─ report/ +├─ 02_rnaseq/ +│ ├─ qc_alignment/ +│ ├─ quant_sample_qc/ +│ ├─ differential_expression/ +│ ├─ annotation_enrichment/ +│ ├─ fgsea_wgcna_candidates/ +│ └─ report/ +├─ 03_wgs/ +│ ├─ qc_alignment_bamqc/ +│ ├─ gvcf_joint_calling/ +│ ├─ filtered_variants/ +│ ├─ imputation_final_vcf/ +│ └─ report/ +├─ 04_population_gwas/ +│ ├─ pca_phylogeny_structure/ +│ ├─ ld_snp_density/ +│ ├─ mlm_cmlm/ +│ ├─ emmax/ +│ ├─ glm/ +│ └─ report/ +└─ 99_manifests_and_reproducibility/ + ├─ manifests/ + ├─ environments/ + ├─ commands/ + └─ checksums/ +``` + +## 78. 报告中建议固定包含的内容 + +1. 数据与研究设计：样本、分组、比较、表型、参考版本。 +2. 软件和数据库：名称、版本、来源和日期。 +3. 实际运行参数：不得只写“默认参数”。 +4. 每阶段输入、输出和关键指标。 +5. 失败、跳过、部分成功和人工确认项。 +6. 核心图表的结果解释与适用边界。 +7. 候选基因/家族/位点的证据矩阵。 +8. 结果文件索引、Manifest 和 SHA256。 + +--- + +# 79. 当前流程边界与解释注意事项 + +- 基因家族的 reference classification、TFBS、PPI、亚细胞定位等属于计算预测，不能替代实验验证。 +- RNA-seq 的候选家族阶段是“优先级筛选”，完整家族成员鉴定仍需进入基因家族 `identify`。 +- WGS 的硬过滤默认阈值是当前标准值，低深度、高倍性、特殊群体或非模式物种可能需要专门调整。 +- 群体样本树用于描述样本遗传聚类，不应直接当作物种系统发育结论。 +- GWAS 的显著关联表示统计相关，不等于因果；群体结构、亲缘、表型质量、环境和样本量都会影响结果。 +- GLM 不含亲缘随机效应，主要用于对照或特定设计；正式结论应优先结合 MLM/CMLM、EMMAX 和 QQ 图。 +- 当前 GWAS 标准模块重点覆盖基因型格式、协变量、三类模型和 Manhattan/QQ 展示；显著位点注释、局部 LD、单倍型和候选区间精细化可作为后续扩展分析。 + +--- + +## 80. 一页式流程索引 + +| 分析 | 起始输入 | 核心阶段 | 最终主要输出 | +|---|---|---|---| +| 基因家族 | 蛋白、HMM、基因组、GFF、CDS、参考家族 | prepare → identify/refine → phylogeny/classification → motif/domain/structure/promoter/TFBS/synteny/kaks/PPI/PEP/expression → combo/report | 最终成员、系统树、分类、特征表、系列图和正式报告 | +| RNA-seq | FASTQ、samples.yaml、参考基因组、注释 | QC → align → quant → sample QC → DE → annotation/enrichment/FGSEA/WGCNA/candidate → report | clean reads、BAM、counts/TPM、DEG、功能结果、模块、候选和报告 | +| WGS | paired FASTQ、wgs_samples.yaml、参考基因组 | QC/index → align/markdup/BAM QC → gVCF → joint call → select/filter → impute → MAF/missing | 高质量、可索引、样本守恒的最终 SNP VCF | +| 群体遗传 | 最终群体 SNP VCF | PCA、样本树、ADMIXTURE、LD decay、SNP density | 群体结构、亲缘聚类、LD 和变异分布结果 | +| GWAS | VCF、表型、Q/FAM；MLM 另需 K | fmt/trait/cov/kinship → MLM/CMLM 或 EMMAX 或 GLM → plot | 完整关联表、标准 p 表、Manhattan/QQ 图和可追溯 Manifest | + +--- + +# 附录 A：各阶段调用软件与关键代码简版 + +## A.1 代码示例的使用说明 + +当前项目采用“两层调用”方式： + +1. **项目阶段入口**：`python -m apps.cli.run_*`。负责读取上游 Manifest、检查样本/参考/分析身份、组织工作目录、写日志、校验产物并发布到 `final/`。 +2. **底层分析软件**：HMMER、BLAST、MAFFT、GATK、DESeq2 等。负责具体科学计算。 + +下面的代码是为了说明每一步的核心逻辑而整理的**简洁版骨架**，省略了绝对路径、重试、并发、日志、SHA256、Manifest 和原子发布等工程代码。正式运行应优先使用项目阶段入口，而不是单独复制底层命令。 + +通用占位符： + +```bash +WS=workspace/ +PREFIX= +CPU=16 +``` + +推荐环境： + +| 分析范围 | 项目环境 | 主要用途 | +|---|---|---| +| 基因家族、WGS、群体遗传、GWAS | `test4py311` | Python 阶段入口及多数命令行软件 | +| RNA-seq 主流程 | `rnaseq` | fastp、HISAT2、samtools、featureCounts、Salmon/Kallisto 等 | +| RNA-seq R 分析与绘图 | `rnaR` | DESeq2、clusterProfiler、fgsea、WGCNA、ggplot2、pheatmap | + +## A.2 基因家族阶段：软件与项目入口 + +| 阶段 | 主要软件/代码 | 当前项目入口 | +|---|---|---| +| `genome_prepare` | Python、Biopython、GFF/FASTA 解析 | `python -m apps.cli.run_genome_prepare` | +| `identify` | HMMER3：`hmmsearch/hmmfetch/hmmpress`；BLAST+：`makeblastdb/blastp` | `python -m apps.cli.run_identify` | +| `identify_refine` | Python/Biopython；HMM/BLAST 证据和成员决策表解析 | `python -m apps.cli.run_identify_refine` | +| `phylogeny` | MAFFT 或 MUSCLE；IQ-TREE 2 或 FastTree | `python -m apps.cli.run_phylogeny` | +| `phylogeny_msa` | MSA 结果解析与内部绘图代码 | 复用 `run_phylogeny` 的 MSA 产物 | +| `msa_region_plot` | Python 比对坐标截取与绘图 | `python -m apps.cli.run_msa_region_plot` | +| `tree_plot_revise` | Newick 解析、Python/R 树图绘制 | `python -m apps.cli.run_tree_plot` | +| `reference_classification` | Newick 树解析、参考标签传播和内部规则 | `python -m apps.cli.run_reference_classification` | +| 亚家族确认 | 人工确认表、Python 校验 | `python -m apps.cli.run_subfamily_confirm` | +| `motif` | MEME Suite：`meme` | `python -m apps.cli.run_motif` | +| `domain` | HMMER3：`hmmscan` | `python -m apps.cli.run_domain` | +| `structure` | Python GFF3 解析与绘图 | `python -m apps.cli.run_structure` | +| `promoter` | Python/Biopython 序列提取；PlantCARE 在线服务 | `python -m apps.cli.run_promoter` | +| `promoter_plot_revise` | Python/R 绘图 | `python -m apps.cli.run_promoter_plot_revise` | +| `tfbs` | MEME Suite：`fimo`；PlantTFDB/JASPAR motif | `python -m apps.cli.run_tfbs` | +| `synteny` | BLAST+、MCScanX、Python/R 绘图 | `python -m apps.cli.run_synteny` | +| `synteny_plot_revise` | 共线数据解析与内部绘图代码 | 复用 synteny 数据重新绘图 | +| `karyotype` | GFF/染色体长度解析、Python/R 绘图 | `python -m apps.cli.run_karyotype` | +| `interspecies_synteny` | BLAST+、MCScanX、跨物种绘图 | `python -m apps.cli.run_interspecies_synteny` | +| `kaks` | MUSCLE；PAML `yn00`、`codeml` | `python -m apps.cli.run_kaks_v2` | +| `ppi` | STRING REST API；Python 网络解析/绘图 | `python -m apps.cli.run_ppi_v2` | +| `pep` | Biopython `ProtParam`；DeepTMHMM；配置的亚细胞定位预测器 | `python -m apps.cli.run_pep` | +| `expression` | Python/R；表达矩阵处理、聚类热图、相关性 | `python -m apps.cli.run_expression` | +| `combo` | Python/R 多面板排版 | `python -m apps.cli.run_combo` | +| `gene_family_report` | Python Manifest 汇总与 Markdown 报告生成 | `python -m apps.cli.run_gene_family_report` | + +### A.2.1 数据整理 + +项目入口： + +```bash +python -m apps.cli.run_genome_prepare \ + --input input/ \ + --output "$WS/genome_prepare" \ + --workspace "$WS" --prefix "$PREFIX" +``` + +关键 Python 逻辑简版： + +```python +from Bio import SeqIO + +records = list(SeqIO.parse("target.pep.fa", "fasta")) +ids = [r.id for r in records] +assert records and len(ids) == len(set(ids)), "蛋白 FASTA 为空或 ID 重复" +SeqIO.write(records, "prepared/target.pep.fasta", "fasta") +``` + +### A.2.2 HMM + BLAST 家族成员鉴定 + +项目入口： + +```bash +python -m apps.cli.run_identify \ + --input-dir input/ \ + --workspace "$WS" --prefix "$PREFIX" \ + --required-pfams PFxxxxx \ + --ref-file reference_family.pep.fa \ + --dom-i-evalue-max 1e-5 --hmm-cov-min 0.5 --cpu "$CPU" +``` + +底层关键命令： + +```bash +# HMM 结构域搜索 +hmmsearch --cpu "$CPU" --domtblout family.domtblout \ + -E 1e-5 family.hmm target.pep.fasta > family.hmmsearch.txt + +# 固定方向 BLAST：参考家族蛋白是 query，目标物种全蛋白组是 database +makeblastdb -in target.pep.fasta -dbtype prot -out target_proteome_db +blastp -query reference_family.pep.fa -db target_proteome_db \ + -evalue 1e-5 -outfmt 6 -num_threads "$CPU" > family.blastp.tsv +``` + +候选合并的关键代码： + +```python +hmm_ids = parse_domtblout("family.domtblout", min_coverage=0.5) +blast_ids = {row.sseqid for row in read_blast6("family.blastp.tsv")} +candidate_ids = sorted(hmm_ids | blast_ids) +final_ids = domain_quality_control(candidate_ids, "family.domtblout") +``` + +### A.2.3 成员精炼 + +```bash +python -m apps.cli.run_identify_refine \ + --input-dir input/ \ + --workspace "$WS" --prefix "$PREFIX" \ + --id-file "$WS/final/$PREFIX.final_members.id" \ + --source-fasta "$WS/final/$PREFIX.final_members.pep.fa" +``` + +```python +# 简化示意：按基因座保留代表转录本，并记录每个成员的决定 +for gene_id, transcripts in group_transcripts(candidates).items(): + representative = max(transcripts, key=lambda x: (x.domain_coverage, x.length)) + write_decision(gene_id, representative.id, reason="best_domain_then_longest") +``` + +### A.2.4 多序列比对与系统树 + +```bash +python -m apps.cli.run_phylogeny \ + --workspace "$WS" --prefix "$PREFIX" \ + --proteins "$WS/final/$PREFIX.final_members.pep.fa" \ + --mafft-params='--auto' \ + --iqtree-params='-m MFP -bb 1000 -alrt 1000' \ + --cpu "$CPU" +``` + +```bash +mafft --auto --thread "$CPU" family.pep.fa > family.pep.msa.fa +iqtree2 -s family.pep.msa.fa -m MFP -bb 1000 -alrt 1000 -nt AUTO + +# 快速树备选，仅在任务明确选择 FastTree 时使用 +FastTree -wag family.pep.msa.fa > family.tree.nwk +``` + +### A.2.5 MSA 区域图与树图修订 + +```bash +python -m apps.cli.run_msa_region_plot \ + --workspace "$WS" --prefix "$PREFIX" \ + --msa-file "$WS/final/$PREFIX.pep.msa.fa" \ + --start-col 100 --end-col 180 + +python -m apps.cli.run_tree_plot \ + --workspace "$WS" --prefix "$PREFIX" \ + --manifest "$WS/final/$PREFIX.phylogeny.manifest.json" +``` + +```python +from Bio import AlignIO +alignment = AlignIO.read("family.pep.msa.fa", "fasta") +region = alignment[:, 99:180] # Python 为 0-based，展示区间须在报告中说明 +AlignIO.write(region, "family.msa.100_180.fa", "fasta") +``` + +### A.2.6 参考辅助亚家族分类 + +```bash +python -m apps.cli.run_reference_classification \ + --workspace "$WS" --prefix "$PREFIX" \ + --tree "$WS/final/$PREFIX.tree.nwk" \ + --reference-subfamily reference_subfamily.tsv + +python -m apps.cli.run_subfamily_confirm \ + --workspace "$WS" --prefix "$PREFIX" \ + --input-dir input/ \ + --edited-subfamily confirmed_subfamily.tsv +``` + +```python +# 只生成候选标签；正式标签仍需确认 +for target_leaf in target_leaves: + nearest_refs = nearest_reference_leaves(tree, target_leaf) + prediction[target_leaf] = consensus_label(nearest_refs) +``` + +### A.2.7 MEME motif 分析 + +```bash +python -m apps.cli.run_motif \ + --workspace "$WS" --prefix "$PREFIX" \ + --fasta "$WS/final/$PREFIX.final_members.pep.fa" \ + --nmotifs 10 --minw 6 --maxw 50 --cpu "$CPU" + +meme family.pep.fa -protein -oc meme_out \ + -mod zoops -nmotifs 10 -minw 6 -maxw 50 +``` + +### A.2.8 Pfam 结构域分析 + +```bash +python -m apps.cli.run_domain \ + --workspace "$WS" --prefix "$PREFIX" \ + --proteins "$WS/final/$PREFIX.final_members.pep.fa" \ + --pfam-base Pfam-A.hmm --evalue 1e-5 --coverage 0.5 --cpu "$CPU" + +hmmscan --cpu "$CPU" --domtblout family.domtblout \ + Pfam-A.hmm family.pep.fa > family.hmmscan.txt +``` + +### A.2.9 基因结构 + +```bash +python -m apps.cli.run_structure \ + --workspace "$WS" --prefix "$PREFIX" \ + --gff3 target.gff3 \ + --tree "$WS/final/$PREFIX.tree.nwk" +``` + +```python +features = parse_gff3("target.gff3") +family_features = [x for x in features if x.gene_id in family_ids] +assert {x.gene_id for x in family_features} <= family_ids +plot_exon_intron_structure(family_features, order=tree_leaf_order) +``` + +### A.2.10 启动子与顺式元件 + +```bash +python -m apps.cli.run_promoter \ + --input-dir input/ \ + --workspace "$WS" --prefix "$PREFIX" \ + --genome target.genome.fa --gff3 target.gff3 \ + --gene-id "$WS/final/$PREFIX.final_members.id" \ + --upstream-length 2000 +``` + +```python +# 正链：[TSS-upstream, TSS)；负链：[TSS, TSS+upstream)，之后反向互补 +if gene.strand == "+": + seq = genome[gene.chrom][max(0, gene.start-upstream):gene.start] +else: + seq = genome[gene.chrom][gene.end:min(chrom_len, gene.end+upstream)].reverse_complement() +``` + +提取后的 FASTA 由项目按 PlantCARE 的条数和大小限制分批提交；PlantCARE 是网页服务，不应写成伪造的本地命令。 + +### A.2.11 TFBS 扫描 + +```bash +python -m apps.cli.run_tfbs \ + --input-dir input/ \ + --workspace "$WS" --prefix "$PREFIX" \ + --genome target.genome.fa --gff3 target.gff3 \ + --gene-ids "$WS/final/$PREFIX.final_members.id" \ + --motif-db-mode auto \ + --fimo-pvalue 1e-4 --fimo-qvalue 0.05 + +fimo --thresh 1e-4 --qv-thresh \ + --oc fimo_out plant_motifs.meme family.promoter.fa +``` + +### A.2.12 种内/种间共线性和染色体定位 + +```bash +python -m apps.cli.run_synteny \ + --input-dir input/ \ + --workspace "$WS" --prefix "$PREFIX" \ + --gff3 target.representative.gff3 \ + --genome target.genome.fa \ + --gene-id-file "$WS/final/$PREFIX.final_members.id" + +python -m apps.cli.run_interspecies_synteny \ + --input-dir input/ \ + --workspace "$WS" --prefix "$PREFIX" \ + --target-protein target.pep.fa \ + --ref-protein reference.pep.fa \ + --blast-engine diamond --plot-mode both + +python -m apps.cli.run_karyotype \ + --input-dir input/ \ + --workspace "$WS" --prefix "$PREFIX" \ + --genome target.genome.fa --gff3 target.gff3 \ + --gene-id-file "$WS/final/$PREFIX.final_members.id" +``` + +```bash +# MCScanX 所需的相似性和坐标输入 +makeblastdb -in all.pep.fa -dbtype prot -out all_pep_db +blastp -query all.pep.fa -db all_pep_db -evalue 1e-5 \ + -outfmt 6 -num_threads "$CPU" > species.blast +MCScanX species +``` + +### A.2.13 Ka/Ks + +```bash +python -m apps.cli.run_kaks_v2 \ + --output "$WS/kaks" --prefix-name "$PREFIX" \ + --cds1 target.cds.fa --cds2 reference.cds.fa \ + --pep1 target.pep.fa --pep2 reference.pep.fa \ + --pairs ortholog_pairs.tsv +``` + +```bash +muscle -align homologs.pep.fa -output homologs.pep.aln.fa +yn00 yn00.ctl +codeml codeml.ctl +``` + +核心逻辑是“同源 CDS → 蛋白/密码子一致对齐 → `yn00` 估算 Ka、Ks → 必要时 `codeml` 比较选择模型”，不能直接对未校正阅读框的 CDS 做普通核酸比对。 + +### A.2.14 PPI、蛋白性质和表达 + +```bash +python -m apps.cli.run_ppi_v2 \ + --input "$WS/final/$PREFIX.final_members.id" \ + --species \ + --output "$WS/ppi" --prefix "$PREFIX" + +python -m apps.cli.run_pep \ + --workspace "$WS" --prefix "$PREFIX" \ + --fasta "$WS/final/$PREFIX.final_members.pep.fa" --cpu "$CPU" + +python -m apps.cli.run_expression \ + --workspace "$WS" --prefix "$PREFIX" \ + --mode legacy \ + --expression-file expression.tpm.tsv \ + --sample-file samples.tsv \ + --gene-id-file "$WS/final/$PREFIX.final_members.id" +``` + +```python +from Bio.SeqUtils.ProtParam import ProteinAnalysis + +seq = str(record.seq).replace("*", "") +pa = ProteinAnalysis(seq) +row = { + "length": len(seq), + "molecular_weight": pa.molecular_weight(), + "theoretical_pI": pa.isoelectric_point(), + "gravy": pa.gravy(), +} +``` + +STRING 查询的简化逻辑： + +```python +mapped = string_map_identifiers(family_ids, species=taxon_id) +edges = string_network(mapped, required_score=400) +write_tsv(edges, "family.ppi_network.tsv") +``` + +### A.2.15 综合图和报告 + +```bash +python -m apps.cli.run_combo \ + --workspace "$WS" --prefix "$PREFIX" \ + --subfamily confirmed_subfamily.tsv \ + --panels motif,structure,domain + +python -m apps.cli.run_gene_family_report \ + --workspace "$WS" --prefix "$PREFIX" +``` + +## A.3 RNA-seq 阶段：软件与项目入口 + +| 阶段 | 主要软件/包 | 当前项目入口 | +|---|---|---| +| `rnaseq_qc` | fastp；显式备选 cutadapt | `python -m apps.cli.run_rnaseq_qc` | +| `rnaseq_align` | HISAT2、samtools | `python -m apps.cli.run_rnaseq_align` | +| `rnaseq_quant` | featureCounts/Subread；可选 Salmon 或 Kallisto + tximport；辅助 StringTie | `python -m apps.cli.run_rnaseq_quant` | +| `rnaseq_sample_qc` | R：DESeq2、ggplot2、pheatmap | `python -m apps.cli.run_rnaseq_sample_qc` | +| `rnaseq_de` | R：DESeq2、ggplot2 | `python -m apps.cli.run_rnaseq_de` | +| `annotation` | InterProScan、KofamScan | `python -m apps.cli.run_annotation` | +| `rnaseq_go_enrich` | R：clusterProfiler 或项目通用富集代码 | `python -m apps.cli.run_rnaseq_enrich` | +| `rnaseq_kegg_enrich` | R：clusterProfiler；gene2ko | `python -m apps.cli.run_rnaseq_enrich` | +| `rnaseq_fgsea` | R：fgsea、data.table、ggplot2 | `python -m apps.cli.run_rnaseq_fgsea` | +| `rnaseq_wgcna` | R：WGCNA、DESeq2、ggplot2、pheatmap | `python -m apps.cli.run_rnaseq_wgcna` | +| `rnaseq_candidate_family` | Python 证据整合与排序；InterPro/RAG 为可选证据 | `python -m apps.cli.run_rnaseq_candidate_family` | +| `rnaseq_report` | Python Manifest 汇总与 Markdown 生成 | `python -m apps.cli.run_rnaseq_report` | + +### A.3.1 FASTQ 质控 + +```bash +python -m apps.cli.run_rnaseq_qc \ + --input raw_fastq --config samples.yaml --output "$WS" \ + --job-id --prefix "$PREFIX" --cpu "$CPU" \ + --qc-tool fastp --detect-adapter-for-pe +``` + +```bash +fastp \ + -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ + -o sample.clean.R1.fastq.gz -O sample.clean.R2.fastq.gz \ + --detect_adapter_for_pe \ + --cut_front --cut_tail --cut_window_size 4 --cut_mean_quality 20 \ + --qualified_quality_phred 20 --unqualified_percent_limit 30 \ + --n_base_limit 5 --length_required 30 \ + --thread 8 --html sample.fastp.html --json sample.fastp.json +``` + +若明确选择 cutadapt： + +```bash +cutadapt -q 20 -m 50 -j 8 \ + -a AGATCGGAAGAGC -A AGATCGGAAGAGC \ + -o sample.clean.R1.fastq.gz -p sample.clean.R2.fastq.gz \ + sample_R1.fastq.gz sample_R2.fastq.gz +``` + +### A.3.2 HISAT2 比对 + +```bash +python -m apps.cli.run_rnaseq_align \ + --input "$WS/rnaseq_qc/analysis" \ + --hisat2-index genome_index \ + --output "$WS" --job-id --prefix "$PREFIX" \ + --sample-config samples.yaml --cpu "$CPU" +``` + +```bash +# 建索引通常只需对同一参考做一次 +hisat2-build -p "$CPU" genome.fa genome_index + +# 非链特异文库示例；链特异文库由项目按样本配置增加 --rna-strandness +hisat2 -x genome_index \ + -1 sample.clean.R1.fastq.gz -2 sample.clean.R2.fastq.gz \ + -p "$CPU" 2> sample.hisat2.log \ +| samtools sort -@ "$CPU" -o sample.sorted.bam +samtools index -@ "$CPU" sample.sorted.bam +``` + +### A.3.3 featureCounts 与可选转录本定量 + +```bash +python -m apps.cli.run_rnaseq_quant \ + --input "$WS/rnaseq_align/analysis" \ + --gtf annotation.gtf --output "$WS" \ + --job-id --prefix "$PREFIX" --cpu "$CPU" \ + --sample-config samples.yaml --quant-mode dual --abundance-engine auto +``` + +```bash +# -s: 0=non_stranded, 1=forward, 2=reverse；混合链特异性样本分批运行后再合并 +featureCounts -T "$CPU" -p -s 0 \ + -a annotation.gtf -o gene.featurecounts.txt *.sorted.bam +``` + +Salmon + tximport 可选轨道： + +```bash +salmon index -t transcripts.fa -i salmon_index -k 31 +salmon quant -i salmon_index -l A \ + -1 sample.clean.R1.fastq.gz -2 sample.clean.R2.fastq.gz \ + -p "$CPU" --validateMappings -o salmon/sample +Rscript scripts/tximport_abundance.R \ + --engine salmon --sample-files sample_quant_files.tsv \ + --tx2gene tx2gene.tsv --output salmon_gene_level +``` + +由 raw count 计算 FPKM/TPM 的核心公式： + +```python +rpk = counts / (gene_length_bp / 1000.0) +tpm = rpk / rpk.sum() * 1_000_000 +fpkm = counts / ((gene_length_bp / 1000.0) * (library_size / 1_000_000.0)) +``` + +### A.3.4 样本 PCA、相关性与聚类 + +```bash +python -m apps.cli.run_rnaseq_sample_qc \ + --workspace "$WS" --prefix "$PREFIX" \ + --count-matrix "$WS/final/$PREFIX.count_matrix.txt" \ + --config samples.yaml +``` + +```r +library(DESeq2) +library(pheatmap) + +dds <- DESeqDataSetFromMatrix(countData=round(counts), + colData=data.frame(row.names=colnames(counts)), + design=~1) +vsd <- vst(dds, blind=TRUE) +pca <- prcomp(t(assay(vsd))) +pheatmap(cor(assay(vsd)), annotation_col=sample_info) +``` + +### A.3.5 DESeq2 差异表达 + +```bash +python -m apps.cli.run_rnaseq_de \ + --input "$WS/final" --output "$WS" \ + --prefix "$PREFIX" --job-id \ + --config samples.yaml --log2fc 1 --padj 0.05 +``` + +```r +library(DESeq2) + +dds <- DESeqDataSetFromMatrix(countData=round(counts), + colData=sample_info, + design=~condition) +dds <- DESeq(dds) +res <- results(dds, contrast=c("condition", "Treatment", "CK")) +deg <- subset(as.data.frame(res), !is.na(padj) & padj < 0.05 & abs(log2FoldChange) >= 1) +``` + +### A.3.6 InterProScan 与 KofamScan 注释 + +```bash +python -m apps.cli.run_annotation \ + --input proteome.fa --workspace "$WS" --prefix "$PREFIX" \ + --job-id --cpu "$CPU" --annotation-scope both +``` + +```bash +interproscan.sh -i proteome.fa -f TSV,XML,GFF3 \ + -goterms -pa -cpu "$CPU" -o proteome.interproscan.tsv + +exec_annotation -o proteome.kofamscan.detail.tsv \ + --format detail-tsv --cpu "$CPU" proteome.fa +``` + +随后由项目把 InterProScan GO term 转成 `gene2go.txt`，把 KOfam 命中转成 `gene2ko.txt`。 + +### A.3.7 GO/KEGG 过度富集 + +```bash +python -m apps.cli.run_rnaseq_enrich \ + --input "$WS/final" --output "$WS/rnaseq_enrich" \ + --prefix "$PREFIX" --job-id \ + --gene2go gene2go.txt --gene2ko gene2ko.txt +``` + +```r +library(clusterProfiler) + +go <- enricher(gene=deg_ids, + universe=background_ids, + TERM2GENE=gene2go, + pAdjustMethod="BH") +kegg <- enricher(gene=deg_ids, + universe=background_ids, + TERM2GENE=pathway2gene, + pAdjustMethod="BH") +``` + +### A.3.8 FGSEA + +```bash +python -m apps.cli.run_rnaseq_fgsea \ + --workspace "$WS" --prefix "$PREFIX" \ + --comparison Treatment_vs_CK \ + --gene2ko gene2ko.txt \ + --min-size 15 --max-size 500 --nperm 10000 --seed 1 +``` + +```r +library(fgsea) +ranks <- setNames(de_table$log2FoldChange, de_table$gene_id) +ranks <- sort(ranks[is.finite(ranks)], decreasing=TRUE) +fg <- fgsea(pathways=pathway_gene_sets, stats=ranks, + minSize=15, maxSize=500, nperm=10000) +``` + +### A.3.9 WGCNA + +```bash +python -m apps.cli.run_rnaseq_wgcna \ + --workspace "$WS" --prefix "$PREFIX" \ + --expression-matrix "$WS/final/$PREFIX.rnaseq_sample_qc.vst_matrix.tsv" \ + --config samples.yaml \ + --min-mean 1 --max-genes 8000 \ + --min-module-size 30 --merge-cut-height 0.25 +``` + +```r +library(WGCNA) +datExpr <- t(vst_matrix) +powers <- 1:20 +sft <- pickSoftThreshold(datExpr, powerVector=powers, verbose=0) +net <- blockwiseModules(datExpr, + power=sft$powerEstimate, + minModuleSize=30, + mergeCutHeight=0.25, + numericLabels=TRUE) +module_trait_cor <- cor(net$MEs, traits, use="p") +``` + +### A.3.10 候选家族与报告 + +```bash +python -m apps.cli.run_rnaseq_candidate_family \ + --workspace "$WS" --prefix "$PREFIX" \ + --de-files treatment_vs_ck.DE_genes.txt \ + --protein-fasta proteome.fa \ + --interproscan-tsv proteome.interproscan.tsv + +python -m apps.cli.run_rnaseq_report \ + --workspace "$WS" --prefix "$PREFIX" +``` + +```python +# 简化的多证据排序示意，正式权重以项目参数为准 +score = ( + de_score(gene) + + pathway_score(gene) + + wgcna_score(gene) + + domain_score(gene) +) +ranked_genes = sorted(genes, key=lambda g: score[g], reverse=True) +``` + +## A.4 WGS 阶段：软件与项目入口 + +| 阶段 | 主要软件 | 当前项目入口 | +|---|---|---| +| `wgs_qc` | fastp | `python -m apps.cli.run_wgs_qc` | +| `wgs_index` | samtools、BWA-MEM2、GATK | `python -m apps.cli.run_wgs_index` | +| `wgs_align` | BWA-MEM2、samtools | `python -m apps.cli.run_wgs_align` | +| `wgs_markdup` | GATK MarkDuplicates | `python -m apps.cli.run_wgs_markdup` | +| `wgs_bam_qc` | samtools flagstat/depth | `python -m apps.cli.run_wgs_bam_qc` | +| `wgs_haplotype_calling` | GATK HaplotypeCaller | `python -m apps.cli.run_wgs_haplotype_calling` | +| `wgs_gvcf_combine` | GATK CombineGVCFs | `python -m apps.cli.run_wgs_gvcf_combine` | +| `wgs_joint_genotyping` | GATK GenotypeGVCFs | `python -m apps.cli.run_wgs_joint_genotyping` | +| `wgs_variant_select` | GATK SelectVariants；用户确认后的 bcftools 备选 | `python -m apps.cli.run_wgs_variant_select` | +| `wgs_variant_filter` | GATK VariantFiltration/SelectVariants；确认后的 bcftools 备选 | `python -m apps.cli.run_wgs_variant_filter` | +| `wgs_impute` | Beagle | `python -m apps.cli.run_wgs_impute` | +| `wgs_maf_missing_filter` | bcftools | `python -m apps.cli.run_wgs_maf_missing_filter` | + +### A.4.1 WGS clean reads + +```bash +python -m apps.cli.run_wgs_qc \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --wgs-samples wgs_samples.yaml --fastp fastp \ + --threads 8 --max-jobs 2 +``` + +```bash +fastp -i S01_R1.fastq.gz -I S01_R2.fastq.gz \ + -o S01.clean.R1.fastq.gz -O S01.clean.R2.fastq.gz \ + --qualified_quality_phred 5 \ + --unqualified_percent_limit 50 \ + --n_base_limit 10 --length_required 50 \ + --html S01.fastp.html --json S01.fastp.json +``` + +### A.4.2 参考索引 + +```bash +python -m apps.cli.run_wgs_index \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --reference-fasta reference.fa \ + --samtools samtools --bwa-mem2 bwa-mem2 --gatk gatk --threads "$CPU" +``` + +```bash +samtools faidx reference.fa +gatk CreateSequenceDictionary -R reference.fa -O reference.dict +bwa-mem2 index reference.fa +``` + +### A.4.3 比对、排序和索引 + +```bash +python -m apps.cli.run_wgs_align \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --qc-manifest "$WS/final/$PREFIX.wgs_qc.manifest.json" \ + --index-manifest "$WS/final/$PREFIX.wgs_index.manifest.json" \ + --bwa-mem2 bwa-mem2 --samtools samtools \ + --threads 8 --max-jobs 2 +``` + +```bash +bwa-mem2 mem -t "$CPU" \ + -R '@RG\tID:S01\tSM:S01\tPL:ILLUMINA' \ + reference.fa S01.clean.R1.fastq.gz S01.clean.R2.fastq.gz \ +| samtools sort -@ "$CPU" -o S01.sorted.bam +samtools index -@ "$CPU" S01.sorted.bam +``` + +### A.4.4 去重复与 BAM 质控 + +```bash +python -m apps.cli.run_wgs_markdup \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --align-manifest "$WS/final/$PREFIX.wgs_align.manifest.json" \ + --gatk gatk --samtools samtools \ + --max-jobs 2 --remove-duplicates true + +python -m apps.cli.run_wgs_bam_qc \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --markdup-manifest "$WS/final/$PREFIX.wgs_markdup.manifest.json" \ + --index-manifest "$WS/final/$PREFIX.wgs_index.manifest.json" \ + --samtools samtools --coverage-threshold 1 \ + --threads 2 --max-jobs 2 +``` + +```bash +gatk MarkDuplicates \ + -I S01.sorted.bam -O S01.markdup.bam \ + -M S01.markdup.metrics.txt \ + --REMOVE_DUPLICATES true --CREATE_INDEX true \ + --VALIDATION_STRINGENCY LENIENT + +samtools flagstat S01.markdup.bam > S01.flagstat.txt +samtools depth -aa -d 0 S01.markdup.bam > S01.depth.tsv +``` + +`samtools depth -aa` 保留零覆盖位点，才能以整个参考基因组为分母计算覆盖率。 + +### A.4.5 单样本 gVCF 与队列合并 + +```bash +python -m apps.cli.run_wgs_haplotype_calling \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --markdup-manifest "$WS/final/$PREFIX.wgs_markdup.manifest.json" \ + --bam-qc-manifest "$WS/final/$PREFIX.wgs_bam_qc.manifest.json" \ + --index-manifest "$WS/final/$PREFIX.wgs_index.manifest.json" \ + --gatk gatk --native-pair-hmm-threads 2 --max-jobs 2 + +python -m apps.cli.run_wgs_gvcf_combine \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --haplotype-manifest "$WS/final/$PREFIX.wgs_haplotype_calling.manifest.json" \ + --index-manifest "$WS/final/$PREFIX.wgs_index.manifest.json" \ + --gatk gatk +``` + +```bash +gatk HaplotypeCaller \ + -R reference.fa -I S01.markdup.bam \ + -O S01.g.vcf.gz -ERC GVCF \ + --native-pair-hmm-threads 2 + +gatk CombineGVCFs -R reference.fa \ + --variant S01.g.vcf.gz --variant S02.g.vcf.gz \ + -O cohort.combined.g.vcf.gz +``` + +### A.4.6 联合分型 + +```bash +python -m apps.cli.run_wgs_joint_genotyping \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --combine-manifest "$WS/final/$PREFIX.wgs_gvcf_combine.manifest.json" \ + --index-manifest "$WS/final/$PREFIX.wgs_index.manifest.json" \ + --gatk gatk --max-alternate-alleles 6 --calling-threshold 10 +``` + +```bash +gatk GenotypeGVCFs \ + -R reference.fa -V cohort.combined.g.vcf.gz \ + -O cohort.joint.vcf.gz \ + --max-alternate-alleles 6 \ + --standard-min-confidence-threshold-for-calling 10 +``` + +### A.4.7 SNP/INDEL 分离 + +```bash +python -m apps.cli.run_wgs_variant_select \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --joint-manifest "$WS/final/$PREFIX.wgs_joint_genotyping.manifest.json" \ + --method gatk --gatk gatk +``` + +```bash +gatk SelectVariants -R reference.fa -V cohort.joint.vcf.gz \ + --select-type-to-include SNP -O cohort.snp.vcf.gz +gatk SelectVariants -R reference.fa -V cohort.joint.vcf.gz \ + --select-type-to-include INDEL -O cohort.indel.vcf.gz +``` + +### A.4.8 SNP/INDEL 硬过滤 + +```bash +python -m apps.cli.run_wgs_variant_filter \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --select-manifest "$WS/final/$PREFIX.wgs_variant_select.manifest.json" \ + --method gatk --gatk gatk +``` + +```bash +gatk VariantFiltration -R reference.fa -V cohort.snp.vcf.gz \ + --filter-name QD2 --filter-expression 'QD < 2.0' \ + --filter-name QUAL30 --filter-expression 'QUAL < 30.0' \ + --filter-name SOR3 --filter-expression 'SOR > 3.0' \ + --filter-name FS60 --filter-expression 'FS > 60.0' \ + --filter-name MQ40 --filter-expression 'MQ < 40.0' \ + --filter-name MQRankSum-12.5 --filter-expression 'MQRankSum < -12.5' \ + --filter-name ReadPosRankSum-8 --filter-expression 'ReadPosRankSum < -8.0' \ + -O cohort.snp.filtered.vcf.gz + +gatk SelectVariants -R reference.fa -V cohort.snp.filtered.vcf.gz \ + --exclude-filtered -O cohort.snp.pass.vcf.gz + +gatk VariantFiltration -R reference.fa -V cohort.indel.vcf.gz \ + --filter-name QD2 --filter-expression 'QD < 2.0' \ + --filter-name QUAL30 --filter-expression 'QUAL < 30.0' \ + --filter-name FS200 --filter-expression 'FS > 200.0' \ + --filter-name ReadPosRankSum-20 --filter-expression 'ReadPosRankSum < -20.0' \ + -O cohort.indel.filtered.vcf.gz +``` + +GATK 失败时项目不会自动切换到 bcftools；必须由用户确认并重新提交 `--method bcftools`。 + +### A.4.9 Beagle 填补 + +```bash +python -m apps.cli.run_wgs_impute \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --filter-manifest "$WS/final/$PREFIX.wgs_variant_filter.manifest.json" \ + --beagle-jar beagle.jar \ + --memory 8g --window 40 --overlap 2 --ne 100000 \ + --iterations 12 --phase-states 280 --imp-states 1600 --cluster 0.005 +``` + +```bash +java -Xmx8g -jar beagle.jar \ + gt=cohort.snp.pass.vcf.gz out=cohort.imputed \ + window=40 overlap=2 ne=100000 iterations=12 \ + phase-states=280 imp-states=1600 cluster=0.005 +``` + +### A.4.10 MAF 与缺失率终筛 + +```bash +python -m apps.cli.run_wgs_maf_missing_filter \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --impute-manifest "$WS/final/$PREFIX.wgs_impute.manifest.json" \ + --bcftools bcftools --maf 0.05 --max-missing 0.5 +``` + +```bash +# F_MISSING <= 0.5，且次要等位基因频率 MAF >= 0.05 +bcftools view -m2 -M2 -v snps cohort.imputed.vcf.gz \ +| bcftools +fill-tags -Ou -- -t MAF,F_MISSING \ +| bcftools view -i 'MAF>=0.05 && F_MISSING<=0.5' \ + -Oz -o cohort.final.snp.vcf.gz +bcftools index -t cohort.final.snp.vcf.gz +``` + +## A.5 群体遗传与 GWAS：软件和项目入口 + +| 阶段 | 主要软件/包 | 当前项目入口 | +|---|---|---| +| `pop_pca` | bcftools、PLINK、R/ggplot2 | `python -m apps.cli.run_pop_pca` | +| `pop_phylogeny` | TASSEL、FastTree、R/ape | `python -m apps.cli.run_pop_phylogeny` | +| `pop_structure` | bcftools、PLINK、ADMIXTURE、R | `python -m apps.cli.run_pop_structure` | +| `pop_ld_decay` | PopLDdecay、R/ggplot2 | `python -m apps.cli.run_pop_ld_decay` | +| `pop_snp_density` | bcftools、R/CMplot | `python -m apps.cli.run_pop_snp_density` | +| `gwas_fmt` | TASSEL | `python -m apps.cli.run_gwas_tassel` | +| `gwas_kinship` | TASSEL centered IBS | `python -m apps.cli.run_gwas_kinship` | +| `gwas_trait` | Python 表型校验与标准化 | `python -m apps.cli.run_gwas_trait` | +| `gwas_cov` | Python Q/FAM 对齐 | `python -m apps.cli.run_gwas_covariate` | +| `gwas_mlm` | TASSEL MLM/P3D/CMLM | `python -m apps.cli.run_gwas_tassel_mlm` | +| `gwas_plot` | Python/R Manhattan、QQ | `python -m apps.cli.run_gwas_plot` | +| `gwas_emmax` | PLINK/内部格式准备、EMMAX | `python -m apps.cli.run_gwas_emmax` | +| `gwas_emplot` | Python/R Manhattan、QQ | `python -m apps.cli.run_gwas_emmax_plot` | +| `gwas_glm` | TASSEL GLM | `python -m apps.cli.run_gwas_tassel_glm` | + +### A.5.1 PCA + +```bash +python -m apps.cli.run_pop_pca \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf cohort.final.snp.vcf.gz +``` + +```bash +bcftools view -m2 -M2 -v snps cohort.final.snp.vcf.gz \ + -Oz -o population.biallelic.snp.vcf.gz +plink --vcf population.biallelic.snp.vcf.gz \ + --indep-pairwise 50 10 0.2 --out population +plink --vcf population.biallelic.snp.vcf.gz \ + --extract population.prune.in --pca 20 --out population +``` + +### A.5.2 群体样本树 + +```bash +python -m apps.cli.run_pop_phylogeny \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf cohort.final.snp.vcf.gz +``` + +```bash +run_pipeline.pl -Xmx16g \ + -importGuess population.biallelic.snp.vcf.gz \ + -export population -exportType Phylip_Inter +FastTree -nt -gtr population.phy > population.tree.nwk +``` + +```r +library(ape) +tree <- read.tree("population.tree.nwk") +pdf("population.tree.circular.pdf") +plot(tree, type="fan", cex=0.6) +dev.off() +``` + +### A.5.3 ADMIXTURE 群体结构 + +```bash +python -m apps.cli.run_pop_structure \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf cohort.final.snp.vcf.gz --k-min 2 --k-max 10 +``` + +```bash +plink --vcf population.biallelic.snp.vcf.gz \ + --indep-pairwise 50 10 0.2 --make-bed --out population + +for K in 2 3 4 5 6 7 8 9 10; do + admixture --cv=5 -s 1 population.bed "$K" \ + > "population.K${K}.log" 2>&1 +done +``` + +### A.5.4 LD 衰减 + +```bash +python -m apps.cli.run_pop_ld_decay \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf cohort.final.snp.vcf.gz \ + --maf 0.05 --missing 0.2 --max-dist-kb 500 +``` + +```bash +PopLDdecay -InVCF cohort.final.snp.vcf.gz \ + -OutStat population.ld.stat.gz \ + -MAF 0.05 -Miss 0.2 -MaxDist 500 +``` + +### A.5.5 SNP 密度 + +```bash +python -m apps.cli.run_pop_snp_density \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf cohort.final.snp.vcf.gz --chromosome-policy all_vcf_contigs +``` + +```bash +# 生成绘图所需的 Chr、Pos 表 +bcftools query -f '%CHROM\t%POS\n' cohort.final.snp.vcf.gz \ + > population.snp.coordinates.tsv +``` + +```r +library(CMplot) +CMplot(density_table, plot.type="d", bin.size=1e6, + file="png", file.name="snp_density_1Mb") +``` + +### A.5.6 GWAS 格式、表型、Q 协变量与亲缘矩阵 + +```bash +python -m apps.cli.run_gwas_tassel \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf cohort.final.snp.vcf.gz --tassel-pipeline run_pipeline.pl + +python -m apps.cli.run_gwas_trait \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --phenotype phenotype.tsv \ + --sample-column sample_id --trait-column height_cm --trait-name height_cm + +python -m apps.cli.run_gwas_covariate \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --q population.3.Q --fam population.fam --k 3 + +python -m apps.cli.run_gwas_kinship \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf "$WS/final/$PREFIX.gwas_fmt.sorted.vcf" +``` + +TASSEL 核心操作简版： + +```bash +# VCF 排序/导出 HapMap +run_pipeline.pl -Xmx16g \ + -importGuess cohort.final.snp.vcf.gz \ + -SortGenotypeFilePlugin -endPlugin \ + -export gwas.sorted -exportType VCF + +run_pipeline.pl -Xmx16g \ + -importGuess gwas.sorted.vcf \ + -export gwas -exportType Hapmap + +# Centered IBS kinship +run_pipeline.pl -Xmx16g \ + -h gwas.hmp.txt \ + -KinshipPlugin -method Centered_IBS -endPlugin \ + -export gwas.kinship +``` + +协变量生成的核心逻辑： + +```python +sample_ids = read_fam_ids("population.fam") +Q = read_numeric_matrix("population.3.Q") +assert len(sample_ids) == Q.shape[0] and Q.shape[1] >= 2 +covariates = Q[:, :-1] # 建模时省略最后一个祖先成分，避免完全共线 +write_covariates(sample_ids, covariates, "gwas.covariates.txt") +``` + +### A.5.7 TASSEL MLM/CMLM + +```bash +python -m apps.cli.run_gwas_tassel_mlm \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf "$WS/final/$PREFIX.gwas_fmt.sorted.vcf" \ + --kinship "$WS/final/$PREFIX.kinship.txt" \ + --trait "$WS/final/$PREFIX.trait.txt" \ + --covariate "$WS/final/$PREFIX.covariates.txt" \ + --genotype-scope-id \ + --phenotype-scope-id \ + --covariate-scope-id +``` + +```bash +# MLM：P3D，compression=None +run_pipeline.pl -Xmx32g \ + -h gwas.hmp.txt -t gwas.trait.txt -q gwas.covariates.txt \ + -k gwas.kinship.txt \ + -mlm -mlmVarCompEst P3D -mlmCompressionLevel None \ + -export gwas.mlm + +# CMLM：P3D，compression=Optimum +run_pipeline.pl -Xmx32g \ + -h gwas.hmp.txt -t gwas.trait.txt -q gwas.covariates.txt \ + -k gwas.kinship.txt \ + -mlm -mlmVarCompEst P3D -mlmCompressionLevel Optimum \ + -export gwas.cmlm +``` + +### A.5.8 EMMAX + +```bash +python -m apps.cli.run_gwas_emmax \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf "$WS/final/$PREFIX.gwas_fmt.sorted.vcf" \ + --trait "$WS/final/$PREFIX.trait.txt" \ + --vcf-manifest "$WS/final/$PREFIX.gwas_fmt.manifest.json" \ + --trait-manifest "$WS/final/$PREFIX.gwas_trait.manifest.json" \ + --genotype-scope-id \ + --phenotype-scope-id \ + --plink plink --emmax emmax-intel64 --rscript Rscript +``` + +EMMAX 底层逻辑简版： + +```bash +# 由同一 genotype scope 生成 EMMAX 输入和关系矩阵 +emmax-kin-intel64 -v -d 10 gwas_genotype + +# -t 基因型、-p 表型、-k 关系矩阵、-c PCA/固定协变量 +emmax-intel64 \ + -t gwas_genotype \ + -p gwas.phenotype.txt \ + -k gwas_genotype.aIBS.kinf \ + -c gwas.pca.covariates.txt \ + -o gwas.emmax +``` + +### A.5.9 TASSEL GLM/Q-GLM + +```bash +python -m apps.cli.run_gwas_tassel_glm \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --vcf "$WS/final/$PREFIX.gwas_fmt.sorted.vcf" \ + --trait "$WS/final/$PREFIX.trait.txt" \ + --covariate "$WS/final/$PREFIX.covariates.txt" \ + --vcf-manifest "$WS/final/$PREFIX.gwas_fmt.manifest.json" \ + --trait-manifest "$WS/final/$PREFIX.gwas_trait.manifest.json" \ + --covariate-manifest "$WS/final/$PREFIX.gwas_cov.manifest.json" \ + --genotype-scope-id \ + --phenotype-scope-id \ + --covariate-scope-id +``` + +```bash +# 不加 Q 的普通 GLM +run_pipeline.pl -Xmx16g \ + -h gwas.hmp.txt -t gwas.trait.txt \ + -FixedEffectLMPlugin -endPlugin \ + -export gwas.glm + +# 加 Q 的结构校正 GLM +run_pipeline.pl -Xmx16g \ + -h gwas.hmp.txt -t gwas.trait.txt -q gwas.covariates.txt \ + -FixedEffectLMPlugin -endPlugin \ + -export gwas.q_glm +``` + +### A.5.10 Manhattan 与 QQ 图 + +```bash +python -m apps.cli.run_gwas_plot \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --mlm-pvalues mlm.pvalues.tsv --cmlm-pvalues cmlm.pvalues.tsv \ + --genotype-scope-id \ + --phenotype-scope-id \ + --covariate-scope-id \ + --source-model-scope-id \ + --mlm-tested-hypotheses \ + --cmlm-tested-hypotheses \ + --rscript Rscript --log-file "$WS/logs/gwas_plot.log" + +python -m apps.cli.run_gwas_emmax_plot \ + --workspace "$WS" --job-id --prefix "$PREFIX" \ + --emmax-pvalues emmax.pvalues.tsv \ + --emmax-manifest "$WS/final/$PREFIX.gwas_emmax.manifest.json" \ + --genotype-scope-id \ + --phenotype-scope-id \ + --source-model-scope-id \ + --emmax-tested-hypotheses \ + --rscript Rscript --log-file "$WS/logs/gwas_emplot.log" +``` + +```python +import numpy as np + +p = result["p"].clip(lower=np.finfo(float).tiny, upper=1.0) +result["minus_log10_p"] = -np.log10(p) +bonferroni = -np.log10(0.05 / tested_marker_count) + +# QQ 理论值与观测值 +observed = -np.log10(np.sort(p)) +expected = -np.log10((np.arange(1, len(p)+1) - 0.5) / len(p)) +``` + +## A.6 正式运行时不应省略的工程代码 + +上面的简版命令只展示科学计算。当前项目阶段入口还必须完成以下工作： + +```python +# 伪代码：正式阶段入口的共同结构 +upstream = resolve_successful_manifest(required_stage) +verify_path_size_sha256(upstream["artifacts"]) +verify_sample_reference_and_scope_identity(upstream) + +attempt_dir = create_new_attempt_directory(stage, analysis_id) +run_scientific_tool(inputs, parameters, cwd=attempt_dir) +metrics = validate_scientific_outputs(attempt_dir) + +manifest = build_manifest_v1( + stage=stage, + status="success", + inputs=hashed_inputs, + parameters=parameters, + metrics=metrics, + artifacts=validated_artifacts, +) +atomic_publish(attempt_dir, workspace_final_dir) +``` + +不能为了命令简短而省略的关键原则： + +- 不直接用目录中“看起来像”的文件替代权威 Manifest 输入。 +- 不混用不同样本集合、参考基因组、prefix 或 genotype/model scope。 +- 不在工具失败后静默换算法；允许备选软件的阶段必须由用户明确确认。 +- 不把旧任务成功产物复制成新任务结果。 +- 不把空表、损坏 VCF、缺样本 BAM 或未闭合的外部服务任务发布为成功。 +- 每个正式结果都应记录文件大小和 SHA256，保留命令、软件环境和来源信息。 + +--- + +**文档结束。** diff --git a/docs/team/brief_bioflow_dev1_p0.md b/docs/team/brief_bioflow_dev1_p0.md index ac7c79e..05216b6 100644 --- a/docs/team/brief_bioflow_dev1_p0.md +++ b/docs/team/brief_bioflow_dev1_p0.md @@ -31,6 +31,7 @@ G. 阶段 → 镜像映射表（基因家族链约 18 阶段 → bioflow-gf 等 - 适配改动保持最小：除 A/B 清单外一律照抄 dev，发现第 3 处以上需适配的硬编码 → 记入报告交 TL 裁决 ## 【自测要求】 +- 镜像构建：开发/自测阶段的 docker build 由你执行（属自测动作，非发布）；发布级镜像构建与 tagging 由 TL 执行。若 docker build 遇权限/磁盘问题 → 记入报告等 TL - 金标 diff 命令输出为空（24 ID 全等）；manifest 产出合规；镜像内 python -c "import gfa_core.xxx" 加载 .so 成功 - 冒烟脚本沉淀为 bioflow_kernel/tests/smoke_golden.sh - 自测结果写入 docs/team/report_bioflow_dev1.md（完成清单、适配改动清单（逐条）、金标结果、冒烟结果、已知限制） diff --git a/docs/team/brief_bioflow_dev4_p6.md b/docs/team/brief_bioflow_dev4_p6.md index fb84488..863cb55 100644 --- a/docs/team/brief_bioflow_dev4_p6.md +++ b/docs/team/brief_bioflow_dev4_p6.md @@ -7,7 +7,10 @@ ## 职责范围（只动 vstore 侧新增集合与入库脚本；QA 树（8091）零改动；前端引用块渲染与 A3 协作） A. 语料切块入库（vstore 新增独立集合 flow_docs，不混入既有公文/资料集合）： - - 语料源（按优先级）：1) Obsidian《PlantAI_现有分析全流程_基因家族_RNAseq_WGS_GWAS》（79 章，TL 提供导出副本）2) dev 的 gfa_core/schemas/stages/*_contract_v1.md（59 份阶段契约）3) 技能说明（dev docs 内分析方法说明） + - 语料源（按优先级）： + 1) Obsidian《PlantAI_现有分析全流程_基因家族_RNAseq_WGS_GWAS》——TL 已导出副本到仓库 docs/bioflow_corpus/PlantAI_现有分析全流程_基因家族_RNAseq_WGS_GWAS.md（79 章，只读使用；入库脚本以该副本为源，不依赖 Obsidian） + 2) dev 的 gfa_core/schemas/stages/*_contract_v1.md（59 份阶段契约，复制源 /home/xzq/5090/teamWorkspace/zh/PlantAI-dev/，只读） + 3) 技能说明（dev docs 内分析方法说明，只读） - 切块策略：按章/节切块，块含标题路径 + 原文 + 源文件标记；入库脚本 vstore/scripts/ingest_flow_docs.py（幂等可重跑） B. 生信对话内流程咨询：检索 flow_docs 集合 → 答案含引用原文（引用块带源文件与章节定位，UI 可跳转——跳转渲染与 A3 协作，接口约定写入报告） C. 评估集：30-50 条（覆盖各链典型咨询：阶段输入输出/参数/前置条件/软件对照），评测脚本（检索命中 + 答案含引用断言）