(run_dir: Path, args: argparse.Namespace, rows: list[dict[str, str]])
| 643 | |
| 644 | |
| 645 | def execute(run_dir: Path, args: argparse.Namespace, rows: list[dict[str, str]]) -> dict[str, Any]: |
| 646 | reference = args.reference_fasta.expanduser().resolve() |
| 647 | annotation_vcf = args.annotation_vcf.expanduser().resolve() if args.annotation_vcf else None |
| 648 | results: dict[str, Any] = {"ok": True, "steps": []} |
| 649 | (run_dir / "qc").mkdir(parents=True, exist_ok=True) |
| 650 | (run_dir / "variants").mkdir(parents=True, exist_ok=True) |
| 651 | if not (Path(str(reference) + ".fai")).exists(): |
| 652 | faidx = run_cmd(["samtools", "faidx", str(reference)], run_dir, timeout=600) |
| 653 | write_json(run_dir / "logs" / "samtools_faidx.json", faidx) |
| 654 | results["steps"].append({"name": "samtools_faidx", "ok": faidx.get("ok")}) |
| 655 | results["ok"] = bool(results["ok"] and faidx.get("ok")) |
| 656 | |
| 657 | for row in rows: |
| 658 | sample = row["sample"] |
| 659 | alignment = Path(row["alignment"]) |
| 660 | quickcheck = run_cmd(["samtools", "quickcheck", "-v", str(alignment)], run_dir, timeout=300) |
| 661 | write_json(run_dir / "logs" / f"{sample}.quickcheck.json", quickcheck) |
| 662 | flagstat = run_cmd(["samtools", "flagstat", str(alignment)], run_dir, timeout=600) |
| 663 | write_json(run_dir / "logs" / f"{sample}.flagstat.json", flagstat) |
| 664 | write_text(run_dir / "qc" / f"{sample}.flagstat.txt", flagstat.get("stdout_tail", "")) |
| 665 | idxstats = run_cmd(["samtools", "idxstats", str(alignment)], run_dir, timeout=600) |
| 666 | write_json(run_dir / "logs" / f"{sample}.idxstats.json", idxstats) |
| 667 | write_text(run_dir / "qc" / f"{sample}.idxstats.tsv", idxstats.get("stdout_tail", "")) |
| 668 | region_qc = run_region_qc(run_dir, sample, alignment, args.region, args.callable_min_depth) |
| 669 | write_json(run_dir / "logs" / f"{sample}.coverage.json", region_qc["coverage"]) |
| 670 | write_json(run_dir / "logs" / f"{sample}.depth.json", region_qc["depth"]) |
| 671 | write_json(run_dir / "qc" / f"{sample}.callability.json", region_qc["callability"]) |
| 672 | vcf = run_dir / "variants" / f"{sample}.vcf.gz" |
| 673 | call = run_cmd( |
| 674 | ["bash", "-c", bcftools_call_command(reference, alignment, vcf, args.region)], |
| 675 | run_dir, |
| 676 | timeout=3600, |
| 677 | ) |
| 678 | write_json(run_dir / "logs" / f"{sample}.bcftools_call.json", call) |
| 679 | write_text(run_dir / "logs" / f"{sample}.bcftools_call.log", call.get("stdout_tail", "")) |
| 680 | mq_header_fix = ( |
| 681 | normalize_mq_header(run_dir, sample, vcf) |
| 682 | if call.get("ok") |
| 683 | else {"ok": False, "skipped": True} |
| 684 | ) |
| 685 | write_json(run_dir / "logs" / f"{sample}.mq_header_fix.json", mq_header_fix) |
| 686 | pre_annotation_index = ( |
| 687 | run_cmd(["bcftools", "index", "-t", str(vcf)], run_dir, timeout=600) |
| 688 | if call.get("ok") and mq_header_fix.get("ok") |
| 689 | else {"ok": False, "skipped": True} |
| 690 | ) |
| 691 | write_json(run_dir / "logs" / f"{sample}.pre_annotation_index.json", pre_annotation_index) |
| 692 | final_vcf, annotation_result = ( |
| 693 | annotate_vcf(run_dir, sample, vcf, annotation_vcf, args.annotation_columns) |
| 694 | if pre_annotation_index.get("ok") |
| 695 | else (vcf, {"ok": False, "skipped": True}) |
| 696 | ) |
| 697 | write_json(run_dir / "logs" / f"{sample}.annotation.json", annotation_result) |
| 698 | filtered_vcf, filter_result = ( |
| 699 | filter_vcf(run_dir, final_vcf, args.filter_min_qual, args.filter_min_site_dp) |
| 700 | if call.get("ok") and mq_header_fix.get("ok") and annotation_result.get("ok") |
| 701 | else (final_vcf, {"ok": False, "skipped": True}) |
| 702 | ) |
no test coverage detected