Creating implementation plan for integration of FastCall3 and mpileup steps

README.md

@@ -76,11 +76,20 @@ cropability snp -i samples.fa --min-af 0.05 --min-depth 10
 # 计算 LD 矩阵
 cropability ld --n-samples 500 --n-snps 1000
 
+# 运行原生 mpileup（真实 BAM/CRAM 输入）
+cropability pileup -r ref.fa -b sample1.bam sample2.bam -o cohort.mpileup
+
+# 运行变异检测流程（默认 hybrid: 原生 mpileup + 原生 FastCall3 逻辑）
+cropability call-variants -r ref.fa -b sample1.bam sample2.bam -o cohort.vcf --mode hybrid
+
 # 导出 TorchScript 模型（供 Java 调用）
 cropability export --model add --output model.pt
 cropability export --model embedding --output embed.pt
 ```
 
+> NGS 流程默认在 CropAbility 内部执行，不依赖外部 `samtools/FastCall3` 命令。
+> 运行时需要 `pysam`（`pip install "cropability[io]"`）。可选 Rust 后端可用于性能加速。
+
 ### Python API
 
 ```python

cropability/cli/main.py

@@ -9,6 +9,8 @@
   ld        — 计算连锁不平衡矩阵
   gwas      — 全基因组关联分析
   align     — 批量序列比对
+  pileup    — 运行原生 mpileup（CropAbility 内置）
+  call-variants — 运行 NGS 变异检测流程
   export    — 导出 TorchScript 模型
   benchmark — GPU 性能基准测试
 """
@@ -199,6 +201,60 @@ def cmd_ld(args: argparse.Namespace) -> int:
     return 0
 
 
+def cmd_pileup(args: argparse.Namespace) -> int:
+    """运行原生 mpileup 并输出位点汇总文本。"""
+    from cropability.genomics.pipeline import QCThresholds, VariantPipeline
+
+    qc = QCThresholds(
+        min_depth=args.min_depth,
+        min_base_quality=args.min_baseq,
+        min_mapping_quality=args.min_mapq,
+        min_alt_freq=args.min_af,
+    )
+    pipeline = VariantPipeline()
+    report = pipeline.run(
+        mode="mpileup",
+        reference=args.reference,
+        bam_files=args.bam,
+        output=args.output,
+        qc=qc,
+        regions=args.region,
+        dry_run=args.dry_run,
+    )
+    print("native pileup completed")
+    print(f"  output: {report['mpileup']['output']}")
+    return 0
+
+
+def cmd_call_variants(args: argparse.Namespace) -> int:
+    """运行原生变异检测流程（mpileup / fastcall3 / hybrid）。"""
+    from cropability.genomics.pipeline import QCThresholds, VariantPipeline
+
+    qc = QCThresholds(
+        min_depth=args.min_depth,
+        min_base_quality=args.min_baseq,
+        min_mapping_quality=args.min_mapq,
+        min_alt_freq=args.min_af,
+    )
+    pipeline = VariantPipeline()
+    report = pipeline.run(
+        mode=args.mode,
+        reference=args.reference,
+        bam_files=args.bam,
+        output=args.output,
+        qc=qc,
+        regions=args.region,
+        mpileup_output=args.mpileup_output,
+        dry_run=args.dry_run,
+    )
+    print(f"pipeline mode={report['mode']} completed")
+    if "mpileup" in report:
+        print(f"  mpileup: {report['mpileup']['output']}")
+    if "fastcall3" in report:
+        print(f"  vcf: {report['fastcall3']['output']}")
+    return 0
+
+
 # ---------------------------------------------------------------------------
 # 主解析器
 # ---------------------------------------------------------------------------
@@ -243,6 +299,37 @@ def build_parser() -> argparse.ArgumentParser:
     ld.add_argument("--n-samples", type=int, default=200)
     ld.add_argument("--n-snps", type=int, default=500)
 
+    # pileup
+    pileup = sub.add_parser("pileup", help="运行原生 mpileup")
+    pileup.add_argument("-r", "--reference", required=True, help="参考基因组 FASTA")
+    pileup.add_argument("-b", "--bam", required=True, nargs="+", help="输入 BAM/CRAM 文件列表")
+    pileup.add_argument("-o", "--output", required=True, help="输出位点汇总路径（TSV）")
+    pileup.add_argument("--region", default=None, help="区域过滤，如 chr1:1-100000")
+    pileup.add_argument("--min-depth", type=int, default=10, help="最小深度阈值")
+    pileup.add_argument("--min-baseq", type=int, default=20, help="最小碱基质量")
+    pileup.add_argument("--min-mapq", type=int, default=20, help="最小比对质量")
+    pileup.add_argument("--min-af", type=float, default=0.05, help="最小替代等位频率")
+    pileup.add_argument("--dry-run", action="store_true", help="仅打印命令，不实际执行")
+
+    # call-variants
+    call = sub.add_parser("call-variants", help="运行变异检测流程")
+    call.add_argument("-r", "--reference", required=True, help="参考基因组 FASTA")
+    call.add_argument("-b", "--bam", required=True, nargs="+", help="输入 BAM/CRAM 文件列表")
+    call.add_argument("-o", "--output", required=True, help="输出 VCF（或 mpileup）路径")
+    call.add_argument(
+        "--mode",
+        choices=["mpileup", "fastcall3", "hybrid"],
+        default="hybrid",
+        help="流程模式",
+    )
+    call.add_argument("--mpileup-output", default=None, help="hybrid 模式的中间 mpileup 输出路径")
+    call.add_argument("--region", default=None, help="区域过滤，如 chr1:1-100000")
+    call.add_argument("--min-depth", type=int, default=10, help="最小深度阈值")
+    call.add_argument("--min-baseq", type=int, default=20, help="最小碱基质量")
+    call.add_argument("--min-mapq", type=int, default=20, help="最小比对质量")
+    call.add_argument("--min-af", type=float, default=0.05, help="最小替代等位频率")
+    call.add_argument("--dry-run", action="store_true", help="仅打印命令，不实际执行")
+
     return parser
 
 
@@ -252,6 +339,8 @@ def build_parser() -> argparse.ArgumentParser:
     "snp": cmd_snp,
     "export": cmd_export,
     "ld": cmd_ld,
+    "pileup": cmd_pileup,
+    "call-variants": cmd_call_variants,
 }
 
 

cropability/genomics/__init__.py

@@ -9,10 +9,13 @@
 - alignment: 序列比对评分
 """
 
-from cropability.genomics.variant import VariantCaller, SNPResult
-from cropability.genomics.ld import LDCalculator, LDResult
-from cropability.genomics.gwas import GWASEngine, GWASResult
 from cropability.genomics.alignment import SmithWatermanGPU
+from cropability.genomics.fastcall3 import FastCall3Config, FastCall3Runner, FastCall3RunResult
+from cropability.genomics.gwas import GWASEngine, GWASResult
+from cropability.genomics.ld import LDCalculator, LDResult
+from cropability.genomics.pileup import MpileupParser, PileupRecord, PileupSample, PileupSiteSummary
+from cropability.genomics.pipeline import PipelineConfig, QCThresholds, VariantPipeline
+from cropability.genomics.variant import SNPResult, VariantCaller
 
 __all__ = [
     "VariantCaller",
@@ -22,4 +25,14 @@
     "GWASEngine",
     "GWASResult",
     "SmithWatermanGPU",
+    "FastCall3Runner",
+    "FastCall3Config",
+    "FastCall3RunResult",
+    "MpileupParser",
+    "PileupRecord",
+    "PileupSample",
+    "PileupSiteSummary",
+    "VariantPipeline",
+    "PipelineConfig",
+    "QCThresholds",
 ]