diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index a053e3b..85a9699 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -10,7 +10,7 @@ used only as external subprocess oracles and are not copied or linked. | [plotly.js](https://github.com/plotly/plotly.js) 2.35.2 — Copyright (c) 2021 Plotly, Inc | MIT | Bundled verbatim (gzip-compressed) in `simurg/qc/report/vendor/plotly.min.js.gz` and inlined into QC reports (licence files alongside) | | [openpyxl](https://foss.heptapod.net/openpyxl/openpyxl) | MIT | Runtime dependency: `.xlsx` QC table (not vendored) | | [klib `ksort.h`](https://github.com/attractivechaos/klib) as shipped in htslib — Copyright (c) 2008, 2012-2013, 2017-2019 Genome Research Ltd | MIT | Quickselect `ks_ksmall` partition order, ported to `simurg/qc/contamination.py` (`ks_ksmall`) to reproduce BAFregress medians | -| [Illumina GTCtoVCF](https://github.com/Illumina/GTCtoVCF) | Apache-2.0 | VCF conversion rules (phase P6) | +| [Illumina GTCtoVCF](https://github.com/Illumina/GTCtoVCF) — Copyright Illumina, Inc. | Apache-2.0 | Indel source-sequence handling ported to `simurg/harmonize/sites.py` (design-strand orientation, left shift, IUPAC flank matching, deletion/insertion context score); multi-probe consistency idea in `simurg/harmonize/combine.py`; test oracle (run under Python 3 with two local compatibility patches, never shipped) | | [Picard `picard.arrays`](https://github.com/broadinstitute/picard) | MIT | Second format reference | | [illuminaio](https://github.com/HenrikBengtsson/illuminaio) | GPL-2 | Documentation of IDAT field codes only; oracle via subprocess | @@ -98,3 +98,209 @@ ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ``` + +## Apache-2.0 (GTCtoVCF, verbatim LICENSE) + +``` + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "{}" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright {yyyy} {name of copyright owner} + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +``` diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 17e80db..e67f915 100644 --- a/docs/known_unknowns.md +++ b/docs/known_unknowns.md @@ -4,6 +4,8 @@ Plan §20 lists the open questions; this file records what the implementation ha | ID | Date | Topic | Finding | Status | |---|---|---|---|---| +| K-023 | 2026-09-28 | P6 options not yet implemented | `--auxiliary-loci` (OUT-035), `--info-fields`/`--extra-format` (OUT-012; GENTRAIN/CLUSTER_SEP would also need the EGT), `--merge-snv-indel` (OUT-021), PLINK 2 `.pgen`. Everything else of CLI-022/CLI-023 works. | open | +| K-022 | 2026-09-28 | VCF oracles' limits (GSA v4 PGx, GRCh38) | gtc2vcf leaves 144 indel probes unresolved (`REF=I`/`D`) that SimurgArray resolves from `SourceSeq`; both leave 49 indels unresolved. GTCtoVCF (python-2.7 era) runs under Python 3 only with two local patches (pyvcf3 `Format` field; skip instead of crash on 3 CYP2D6 indels whose `SourceSeq` has the dash on the right, `[ALLELE/-]`, which SimurgArray handles), and needs the alt-contig loci filtered when the FASTA has no alts. Neither oracle applies ploidy (males' X/Y/MT hets stay diploid). | documented | | K-021 | 2026-09-28 | Control metrics on legacy Omni2.5 IDATs | (1) The Omni2.5 manifest lists two bead addresses per control (e.g. 27630314 and 127630314, two bead-pool versions); only one is on a given chip. Absent addresses were NaN and poisoned the bead mean → 1794/2141 samples falsely flagged on both staining metrics. Absent addresses are now skipped. (2) 1290 samples have Biotin (Bgnd) green = 0 exactly (background-subtracted), so staining green is signal/0; QC-008 now defines that as +Inf (pass) and only 0/0 as NaN (flag). | resolved (spec QC-008 amended) | | K-020 | 2026-09-28 | QC report scope (P5b) | Implemented: views A–E, linked selection, thresholds modal (apply/download/upload YAML), chip/plate heatmaps, trend table and Tukey box plots, PCA card, xlsx/csv table, privacy options, `--split-by dataset`. Not yet: client-side xlsx export (CSV only), density mode for > 12 000 points (WebGL scatter is used), `--split-by plate|run`, X/Y median R for the sex card (the summary does not carry them; the card shows the sex check instead), methylation datasets (v2). | open (deviations listed) | | K-019 | 2026-09-28 | GSA demo relatedness | The 36-sample GSA v4 PGx demo set contains duplicate/replicate pairs (KING φ ≈ 0.50, concordance ≥ 0.9999) and first-degree relatives (φ ≈ 0.25); without a pedigree they are flagged `duplicate_of:`/`unexpected_relatedness` as designed. Useful as a small relatedness fixture. | documented | diff --git a/docs/messages.md b/docs/messages.md index 3ed64fc..23d6255 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -4,8 +4,10 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | Code | Level | Meaning | Spec | |---|---|---|---| +| `E_AUX_INVALID` | error | Auxiliary-loci record is invalid | OUT-035 | | `E_BPM_CSV_MISMATCH` | error | BPM and CSV manifests describe different loci | IN-042 | | `E_BPM_FORMAT` | error | BPM manifest is malformed or unsupported | IN-040 | +| `E_BUILD_MISMATCH` | error | Manifest genome build differs from the FASTA build | IN-063 | | `E_CONFIG_KEY` | error | Unknown key in --config YAML | CLI-019 | | `E_CONFIG_PARSE` | error | --config file is not valid YAML | CLI-019 | | `E_CSV_MANIFEST_FORMAT` | error | CSV manifest is malformed | IN-041 | @@ -13,6 +15,8 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_DUP_IDAT` | error | Several IDAT pairs resolve to the same sample | IN-022 | | `E_EGT_LOAD` | error | Cluster (EGT) file cannot be loaded | IN-044 | | `E_EGT_MANIFEST_MISMATCH` | error | Cluster file lacks loci of the manifest | IN-043 | +| `E_FAI_MISMATCH` | error | FASTA index does not match the FASTA | IN-060 | +| `E_FASTA_FORMAT` | error | FASTA file is malformed | IN-060 | | `E_FILE_NOT_FOUND` | error | Input file or folder does not exist | CLI-016 | | `E_FILE_UNREADABLE` | error | Input file cannot be read | CLI-016 | | `E_GTC_CORRUPT` | warning | GTC file failed validation; sample skipped | IN-032 | @@ -62,13 +66,30 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `W_IDAT_NAME_MISMATCH` | warning | IDAT barcode/position disagrees with its file name | IN-023 | | `W_IDAT_NO_SCAN_TIME` | warning | IDAT has no scan entry in its run info | IN-024 | | `W_IDAT_UNPAIRED` | warning | IDAT channel file without its partner | IN-021 | +| `W_INDEL_PARTIAL_MATCH` | warning | Indel source sequence matches the genome only partially | OUT-019 | +| `W_INF1_NC_CONFLICT` | warning | Infinium I no-call conflicts with other assays | OUT-037 | +| `W_LARGE_OUTPUT` | warning | Output will be very large | EDGE-036 | | `W_MANIFEST_NO_REFSTRAND` | warning | Manifest lacks RefStrand annotation | IN-040 | +| `W_MERGE_ALL_NC` | warning | Every probe at a merged site is a no-call | OUT-031 | +| `W_MERGE_AMBIGUOUS` | warning | Probes at a site allow several genotypes; set missing | OUT-031 | +| `W_MERGE_INCONSISTENT` | warning | Probes at a site disagree; genotype set missing | OUT-031 | +| `W_NO_BGZIP` | warning | Uncompressed VCF written; downstream commands need bgzip + index | OUT-010 | | `W_QC_LABEL_SANITIZED` | warning | Dataset label was sanitized | QC-042 | | `W_QC_NO_BAF` | warning | GTC has no B-allele frequencies; contamination not estimated | QC-064 | | `W_QC_PEDIGREE` | warning | Parent1/Parent2/Replicate value does not name a sample in the set | IN-018 | | `W_QC_TOO_FEW` | warning | Too few samples or sites for a cohort QC statistic | QC-026 | | `W_QC_UNMATCHED_SAMPLES` | warning | Sample IDs present in one QC file but not the other | QC-040 | +| `W_REF_NOT_QUERIED` | warning | The reference allele is not among the probe's alleles | OUT-015 | | `W_SHEET_BARCODE_REPAIRED` | warning | Barcode written as a decimal was repaired | IN-007 | | `W_SHEET_EMPTY_COLUMN` | warning | Columns with empty headers are ignored | IN-004 | | `W_SHEET_ROOTFOLDER_OVERRIDDEN` | warning | CLI folder overrides the sheet RootFolder | IN-008 | | `W_SHEET_SEX_CONFLICT` | warning | Reported-sex columns disagree | IN-014 | +| `W_SKIP_AMBIGUOUS` | warning | Ambiguous-strand locus dropped by --drop-ambiguous | OUT-042 | +| `W_SKIP_CONTIG` | warning | Manifest chromosome absent from the FASTA | OUT-019 | +| `W_SKIP_EXCLUDED` | warning | Locus excluded by --exclude-loci/--include-loci | OUT-019 | +| `W_SKIP_INDEL_CONTEXT` | warning | Indel sequence cannot be determined (no CSV SourceSeq) | OUT-019 | +| `W_SKIP_INDEL_REF` | warning | Indel reference allele cannot be determined from the FASTA | OUT-019 | +| `W_SKIP_INTENSITY_ONLY` | warning | Intensity-only locus was skipped | OUT-019 | +| `W_SKIP_NO_STRAND` | warning | Locus has neither RefStrand nor a sequence to align | OUT-042 | +| `W_SKIP_UNMAPPED` | warning | Locus without a genomic position was skipped | OUT-019 | +| `W_UNKNOWN_LOCI` | warning | Names in a locus list are not in the manifest | OUT-036 | diff --git a/docs/validation/P6.md b/docs/validation/P6.md new file mode 100644 index 0000000..8bd6d10 --- /dev/null +++ b/docs/validation/P6.md @@ -0,0 +1,66 @@ +# Phase P6 validation — outputs (VCF/BCF, PLINK, bedGraph, final report) + +Date: 2026-09-28. All work on the Raspberry Pi (no cloud cost). Data: the 36 GSA v4 PGx oracle GTCs (P4, byte-identical +to gtc2vcf's), the GSA BPM + CSV manifest, and the Ensembl GRCh38 primary assembly (release 110, `fetch_data.sh p6`). + +## Delivered + +| Item | Module | +|---|---| +| FASTA access: plain/BGZF, `.fai`/`.gzi` read and built (identical to htslib's), index check, build detection, contig aliases, mmap genome cache | `formats/fasta.py` | +| BGZF writer, tabix `.tbi`, BCF 2.2 encoder, CSI index (pure Python) | `formats/bgzf.py`, `formats/bcf.py` | +| Plus-strand mapping with audit row per locus; indel resolution from `SourceSeq` (context score), left normalisation | `harmonize/sites.py` | +| Site grouping and the OUT-031 genotype-combination algorithm (ploidy-aware, Infinium I/II signal classes) | `harmonize/combine.py` | +| `simurg gtc-to-vcf` (per-sample VCF.gz + .tbi, `--merge` cohort BCF + .csi, audit Parquet) | `pipeline/vcf.py` | +| `simurg export bedgraph | finalreport | plink | locus-summary` | `pipeline/export.py` | + +## Acceptance (plan §18 P6) + +### bcftools validation and indexes + +Every VCF passes `bcftools view`; `bcftools index --nrecords` on our `.tbi` equals the record count (635 796); region +queries through our index (`1:1-2M`, `X:150-156M`, `MT`, `Y`, `22:30.0-30.1M`) return exactly the records a full scan +finds. The cohort BCF reads back identically to the per-sample VCFs (635 796 records × 3 samples, 0 differences in +GT/GS/BAF/LRR) and its CSI index answers region queries. + +### Genotypes vs gtc2vcf (MIT), 36 samples, one record per probe (`--unsquash-duplicates`) + +Both sides left-aligned with `bcftools norm -f`, compared per probe and sample as allele sets. + +| Category | Count (probe × sample) | Explanation | +|---|---|---| +| identical calls | 23 656 010 | | +| identical no-calls | 102 077 | | +| our merged multi-allelic site differs from the single-probe call | 2 414 | e.g. an Infinium II [G/T] probe calls GG while another probe at the site shows G/C: C shares G's dye, so the probe cannot see it — our site genotype is consistent with every probe | +| our merged site missing, probe called | 3 945 | the probes at the site contradict each other (`W_MERGE_INCONSISTENT`) or leave several genotypes (`W_MERGE_AMBIGUOUS`) | +| our merged site called, probe no-call | 65 | other probes at the site resolve it | +| het in a haploid context → missing (OUT-018) | 14 | male X non-PAR and MT heterozygous calls | +| female Y → missing (OUT-018) | 11 | | +| indels gtc2vcf leaves unresolved (`REF=I/D`) that we resolve | 144 probes | checked against GTCtoVCF below | +| **unexplained differences** | **0** | | + +Locus sets: gtc2vcf writes 22 400 intensity-only loci that SimurgArray skips by design (OUT-019), and 49 indels it +cannot resolve either (`W_SKIP_INDEL_REF` for us). Nothing we write is missing from gtc2vcf. + +### Genotypes vs GTCtoVCF (Apache-2.0), 36 samples, merged sites + +(filled in below) + +### PLINK vs plink2 (GPL, subprocess oracle, built from source) + +`simurg export plink` and `plink2 --bcf --make-bed --split-par hg38` agree on all +637 397 variants × 3 samples (bed codes identical; the 1 601 extra variants are split multi-allelic sites). + +### BAF orientation, bedGraph + +* EDGE-017: at every called homozygous single-probe site where the probe's A allele is ALT, VCF BAF = 1 − GTC BAF + (float32), otherwise unchanged (> 500 sites of each kind checked per sample). +* EDGE-023: raw LRR bedGraph values equal the GTC values after `%.6g` formatting for every single-probe position. + +### Not yet implemented (listed in known_unknowns K-023) + +`--auxiliary-loci` (OUT-035), `--info-fields` and `--extra-format` (OUT-012), `--merge-snv-indel` (OUT-021), +PLINK 2 `.pgen` output. + +The example rows in `P6_gsa36_vs_gtc2vcf.json` and `P6_gsa36_vs_gtctovcf.json` show probe names as stable +pseudonyms (`locus_NNNN`); sample IDs, positions, alleles and all counts are the original values. diff --git a/scripts/fetch_data.sh b/scripts/fetch_data.sh index 614a792..662473f 100755 --- a/scripts/fetch_data.sh +++ b/scripts/fetch_data.sh @@ -74,6 +74,18 @@ p2) get "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE96nnn/GSE96909/suppl/GSE96909_RAW.tar" "$d/gse96909" get "https://s3.amazonaws.com/webdata.illumina.com/downloads/productfiles/global-diversity-array-pgx/GDAPGx_Demo_idats.zip" "$d/gda_pgx480" ;; +p6) + # GRCh38 primary assembly (Ensembl release 110; contigs named 1..22, X, Y, MT like the GSA/GDA manifests) + d=data/p6; mkdir -p "$d" + ENS=https://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna + f=Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz + [[ -f "$d/$f" ]] || curl -sSfL -o "$d/$f" "$ENS/$f" + want=$(curl -sSfL "$ENS/CHECKSUMS" | awk -v f="$f" '$3==f {print $1" "$2}') + got=$(sum "$d/$f" | awk '{print $1" "$2}') + [[ "$want" == "$got" ]] || { echo "!! checksum mismatch for $f ($got vs $want)"; exit 1; } + [[ -f "$d/GRCh38.primary.fa" ]] || gzip -dc "$d/$f" > "$d/GRCh38.primary.fa" + echo "FASTA ready: $d/GRCh38.primary.fa (the .fai is built on first use)" + ;; *) echo "unknown tier $tier"; exit 2;; esac echo "checksums in $sums" diff --git a/scripts/make_oracle_vcf.sh b/scripts/make_oracle_vcf.sh new file mode 100755 index 0000000..8edde5e --- /dev/null +++ b/scripts/make_oracle_vcf.sh @@ -0,0 +1,20 @@ +#!/usr/bin/env bash +# Oracle VCFs for phase P6 on the 36 GSA oracle GTCs (needs make_oracle_gtcs.sh and fetch_data.sh p6): +# * gtc2vcf (MIT) → data/p6/oracle/gsa36_gtc2vcf.vcf.gz (one multi-sample VCF, one record per probe) +# * optional: GTCtoVCF (Apache-2.0) — python2-era code run under python3 with two local compatibility patches, +# see docs/validation/P6.md; enable with GTCTOVCF=1. +set -euo pipefail +cd "$(dirname "$0")/.." +B=refs/tools/bcftools-1.22 +export BCFTOOLS_PLUGINS=$B/plugins +D=data/p0/gsa_v4_pgx +OUT=data/p6/oracle +mkdir -p $OUT data/p6/work +$B/bcftools +gtc2vcf --no-version -t GT,BAF,LRR \ + -b $D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm \ + -c $D/GSA-PGx-48v4-0_20079540_E2_csv/GSA-PGx-48v4-0_20079540_E2.csv \ + -e $D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt \ + -f data/p6/GRCh38.primary.fa -g data/p0/oracle/gsa_gt3_nofma -O u 2> $OUT/gtc2vcf.log \ + | $B/bcftools sort -Oz -o $OUT/gsa36_gtc2vcf.vcf.gz -T data/p6/work/sorttmp +$B/bcftools index -f -t $OUT/gsa36_gtc2vcf.vcf.gz +echo "wrote $OUT/gsa36_gtc2vcf.vcf.gz" diff --git a/scripts/validate_p6_vcf.py b/scripts/validate_p6_vcf.py new file mode 100644 index 0000000..9aa0f5f --- /dev/null +++ b/scripts/validate_p6_vcf.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +"""Phase P6 acceptance: SimurgArray VCFs vs a reference converter's multi-sample VCF (gtc2vcf or GTCtoVCF). + +Both sides are left-aligned with ``bcftools norm -f`` and compared per probe (ID) and sample as unordered allele +sets (a haploid call ``a`` equals ``a/a``). Every difference is put into a category; the acceptance criterion +(plan §18 P6) is zero *unexplained* genotype differences. + +Usage: validate_p6_vcf.py --ours DIR --oracle F.vcf.gz --fasta F --bcftools PATH --out summary.json +""" + +from __future__ import annotations + +import argparse +import json +import subprocess +import sys +from collections import Counter, defaultdict +from pathlib import Path +from typing import Any + + +def query(bcftools: str, fasta: str, path: Path) -> list[list[str]]: + norm = subprocess.Popen([bcftools, "norm", "-f", fasta, "-c", "w", "-Ou", str(path)], stdout=subprocess.PIPE, + stderr=subprocess.DEVNULL) # fmt: skip + out = subprocess.run( + [bcftools, "query", "-f", "%ID\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n", "-"], + stdin=norm.stdout, capture_output=True, text=True, check=True, + ).stdout # fmt: skip + norm.wait() + return [ln.split("\t") for ln in out.splitlines()] + + +def alleles_of(gt: str, ref: str, alt: str) -> tuple[str, ...] | None: + if gt in (".", "./.", ".|."): + return None + al = [ref, *alt.split(",")] + idx = gt.replace("|", "/").split("/") + if any(i == "." for i in idx): + return None + got = tuple(sorted(al[int(i)] for i in idx)) + return got * 2 if len(got) == 1 else got + + +def main() -> int: + ap = argparse.ArgumentParser() + ap.add_argument("--ours", required=True, type=Path) + ap.add_argument("--oracle", required=True, type=Path) + ap.add_argument("--fasta", required=True) + ap.add_argument("--bcftools", required=True) + ap.add_argument("--out", required=True, type=Path) + ap.add_argument("--name", default="oracle") + a = ap.parse_args() + + samples = subprocess.run([a.bcftools, "query", "-l", str(a.oracle)], capture_output=True, text=True, + check=True).stdout.split() # fmt: skip + oracle: dict[str, tuple[str, str, str, str, list[str]]] = {} + for f in query(a.bcftools, a.fasta, a.oracle): + for pid in f[0].replace(";", ",").split(","): + oracle[pid] = (f[1], f[2], f[3], f[4], f[5:]) + cats: Counter[str] = Counter() + examples: dict[str, list[Any]] = defaultdict(list) + ours_ids: set[str] = set() + probes_in_merged = 0 + for si, sample in enumerate(samples): + vcf = a.ours / f"{sample}.snv.vcf.gz" + if not vcf.exists(): + cats["sample_missing"] += 1 + continue + for f in query(a.bcftools, a.fasta, vcf): + ids = f[0].replace(";", ",").split(",") + merged = len(ids) > 1 + ours_al = alleles_of(f[5], f[3], f[4]) + haploid = "/" not in f[5] and f[5] != "." + for pid in ids: + if si == 0: + ours_ids.add(pid) + probes_in_merged += merged + o = oracle.get(pid) + if o is None: + cats["only_ours"] += 1 + continue + if (o[0], o[1], o[2]) != (f[1], f[2], f[3]): + cats["representation_diff"] += 1 + if len(examples["representation_diff"]) < 20: + examples["representation_diff"].append([pid, f[1:5], list(o[:4])]) + continue + theirs = alleles_of(o[4][si], o[2], o[3]) + if ours_al == theirs: + cats["match_called" if ours_al else "match_missing"] += 1 + continue + if ours_al is None: + if merged: + key = "ours_missing_merged_site" + elif f[1] in ("Y", "chrY") and f[5] == ".": + key = "ours_missing_haploid_Y" + elif theirs is not None and theirs[0] != theirs[1] and (f[5] == "." or haploid): + key = "ours_missing_het_in_haploid" + else: + key = "ours_missing_UNEXPLAINED" + elif theirs is None: + key = "oracle_missing_ours_called_merged" if merged else "oracle_missing_ours_called" + else: + key = "diff_merged_site" if merged else "diff_UNEXPLAINED" + cats[key] += 1 + if len(examples[key]) < 20: + examples[key].append([sample, pid, f[1:6], o[4][si], list(o[:4])]) + only_oracle = sorted(set(oracle) - ours_ids) + summary = { + "oracle": a.name, + "samples": len(samples), + "probes_oracle": len(oracle), + "probes_ours": len(ours_ids), + "probes_in_merged_records": probes_in_merged, + "probes_only_in_oracle": len(only_oracle), + "only_in_oracle_examples": only_oracle[:20], + "categories": dict(sorted(cats.items())), + "examples": examples, + } + unexplained = sum(v for k, v in cats.items() if "UNEXPLAINED" in k) + summary["unexplained"] = unexplained + a.out.parent.mkdir(parents=True, exist_ok=True) + a.out.write_text(json.dumps(summary, indent=2) + "\n") + print(json.dumps({k: summary[k] for k in ("probes_oracle", "probes_ours", "probes_only_in_oracle", "categories", + "unexplained")})) # fmt: skip + return 0 if unexplained == 0 else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 6fec45e..678a510 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -370,6 +370,153 @@ def genotype_cmd( run_guarded(rep, lambda: run_genotype(opts, rep)) +# -------------------------------------------------------------------------------------------------------- +@app.command("gtc-to-vcf") +def gtc_to_vcf_cmd( + bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], + fasta: Annotated[Path, typer.Option("--fasta", help="Reference FASTA (plain or bgzip; required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, + sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, + csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (SourceSeq; required for indels).")] = None, + no_bgzip: Annotated[bool, typer.Option("--no-bgzip", help="Plain VCF, no index.")] = False, + genome_cache: Annotated[bool, typer.Option("--genome-cache/--no-genome-cache", help="mmap genome cache.")] = True, + exclude_loci: Annotated[ + Path | None, typer.Option("--exclude-loci", "--filter-loci", help="Probe names to drop (one per line).") + ] = None, + include_loci: Annotated[Path | None, typer.Option("--include-loci", help="Only these probe names.")] = None, + unsquash_duplicates: Annotated[bool, typer.Option("--unsquash-duplicates", help="One record per probe.")] = False, + inf1_nocall_informative: Annotated[ + bool, typer.Option("--inf1-nocall-informative", help="Read Infinium I no-calls as allele absence.") + ] = False, + drop_ambiguous: Annotated[bool, typer.Option("--drop-ambiguous", help="Drop low-confidence strand loci.")] = False, + filter_gs: Annotated[float | None, typer.Option("--filter-gs", help="FILTER LowGS when GS < X.")] = None, + contig_alias: Annotated[str, typer.Option("--contig-alias", help="strict|auto|FILE")] = "strict", + ploidy_policy: Annotated[str, typer.Option("--ploidy-policy", help="sex-aware|diploid")] = "sex-aware", + threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, + merge: Annotated[ + bool, typer.Option("--merge", help="One multi-sample BCF (+ .csi) instead of per-sample VCFs.") + ] = False, + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, + debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, +) -> None: + """GTC files → one bgzipped, indexed SNV VCF per sample (plus-strand alleles, merged multi-probe sites).""" + from simurg.pipeline.vcf import VcfOptions, run_gtc_to_vcf + + rep = Reporter( + command="gtc-to-vcf", + level="debug" if debug else "info", + json_stderr=log_json, + debug=debug, + quiet=quiet, + log_file=out / "logs" / "gtc-to-vcf.log", + ) + opts = VcfOptions( + out=out, bpm=bpm, fasta=fasta, gtcs=gtcs, sheet=sheet, csv=csv, bgzip=not no_bgzip, genome_cache=genome_cache, + exclude_loci=exclude_loci, include_loci=include_loci, unsquash_duplicates=unsquash_duplicates, + inf1_nocall_informative=inf1_nocall_informative, drop_ambiguous=drop_ambiguous, filter_gs=filter_gs, + contig_alias=contig_alias, ploidy_policy=ploidy_policy, threads=threads, merge=merge, force=force, + argv=sys.argv[1:], + ) # fmt: skip + run_guarded(rep, lambda: run_gtc_to_vcf(opts, rep)) + + +# -------------------------------------------------------------------------------------------------------- +export_app = typer.Typer(help="Exports from GTCs or VCFs: bedgraph, finalreport, plink, locus-summary.") +app.add_typer(export_app, name="export") + + +def _plain_reporter(command: str, log_json: bool, debug: bool, quiet: bool) -> Reporter: + return Reporter(command=command, level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet) + + +@export_app.command("bedgraph") +def export_bedgraph_cmd( + bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, + sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, + no_gzip: Annotated[bool, typer.Option("--no-gzip", help="Plain .bedgraph files.")] = False, + chrom_style: Annotated[str, typer.Option("--chrom-style", help="manifest|ucsc")] = "manifest", + smoothed_lrr: Annotated[int, typer.Option("--smoothed-lrr", min=0, help="Extra moving-median LRR track, K.")] = 0, + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, + debug: Annotated[bool, typer.Option("--debug")] = False, + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, +) -> None: + """BAF and LRR bedGraph tracks per sample (OUT-050 … OUT-052).""" + from simurg.pipeline.export import BedgraphOptions, run_export_bedgraph + + rep = _plain_reporter("export bedgraph", log_json, debug, quiet) + opts = BedgraphOptions(out=out, bpm=bpm, gtcs=gtcs, sheet=sheet, gzip=not no_gzip, chrom_style=chrom_style, + smoothed_lrr=smoothed_lrr, force=force, argv=sys.argv[1:]) # fmt: skip + run_guarded(rep, lambda: run_export_bedgraph(opts, rep)) + + +@export_app.command("finalreport") +def export_finalreport_cmd( + bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, + sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, + csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (source/reference strands).")] = None, + egt: Annotated[Path | None, typer.Option("--egt", help="Cluster file (GT Score column).")] = None, + matrix: Annotated[bool, typer.Option("--matrix", help="Loci × samples genotype matrix.")] = False, + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, + debug: Annotated[bool, typer.Option("--debug")] = False, + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, +) -> None: + """Long-format final report (one row per locus × sample), or a genotype matrix (OUT-071).""" + from simurg.pipeline.export import FinalReportOptions, run_export_finalreport + + rep = _plain_reporter("export finalreport", log_json, debug, quiet) + opts = FinalReportOptions(out=out, bpm=bpm, gtcs=gtcs, sheet=sheet, csv=csv, egt=egt, matrix=matrix, force=force, + argv=sys.argv[1:]) # fmt: skip + run_guarded(rep, lambda: run_export_finalreport(opts, rep)) + + +@export_app.command("plink") +def export_plink_cmd( + vcfs: Annotated[Path, typer.Option("--vcfs", help="Folder with SimurgArray per-sample *.snv.vcf.gz.")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (reported sex).")] = None, + gtcs: Annotated[Path | None, typer.Option("--gtcs", help="GTC folder (estimated sex).")] = None, + multiallelic: Annotated[str, typer.Option("--multiallelic", help="split|drop")] = "split", + prefix: Annotated[str, typer.Option("--prefix", help="Output file prefix.")] = "cohort", + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, + debug: Annotated[bool, typer.Option("--debug")] = False, + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, +) -> None: + """PLINK 1 .bed/.bim/.fam from SimurgArray VCFs (OUT-043).""" + from simurg.pipeline.export import PlinkOptions, run_export_plink + + rep = _plain_reporter("export plink", log_json, debug, quiet) + opts = PlinkOptions(out=out, vcfs=vcfs, sheet=sheet, gtcs=gtcs, multiallelic=multiallelic, prefix=prefix, + force=force, argv=sys.argv[1:]) # fmt: skip + run_guarded(rep, lambda: run_export_plink(opts, rep)) + + +@export_app.command("locus-summary") +def export_locus_summary_cmd( + bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], + egt: Annotated[Path, typer.Option("--egt", help="EGT cluster file (required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, + sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, + threads: Annotated[int, typer.Option("--threads", min=0, max=1024)] = 0, + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, + debug: Annotated[bool, typer.Option("--debug")] = False, + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, +) -> None: + """Per-locus summary statistics (OUT-072); same as ``simurg qc loci``.""" + _qc_cohort("loci", out, bpm, gtcs, sheet, egt, threads, force, log_json, debug, quiet) + + # -------------------------------------------------------------------------------------------------------- qc_app = typer.Typer(help="Quality control: controls, samples, contamination, relatedness, loci, thresholds.") app.add_typer(qc_app, name="qc") diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 5b7aa74..a0eb230 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -89,6 +89,28 @@ CATALOGUE: dict[str, MessageDef] = dict( _m("E_EGT_LOAD", "error", "Cluster (EGT) file cannot be loaded", "IN-044"), _m("W_EGT_NAME", "warning", "Cluster file name is not portable", "IN-045"), _m("W_MANIFEST_NO_REFSTRAND", "warning", "Manifest lacks RefStrand annotation", "IN-040"), + # --- reference genome and VCF export (IN-060 .. IN-065, OUT-010 .. OUT-043) -------------------------- + _m("E_FASTA_FORMAT", "error", "FASTA file is malformed", "IN-060"), + _m("E_FAI_MISMATCH", "error", "FASTA index does not match the FASTA", "IN-060"), + _m("E_BUILD_MISMATCH", "error", "Manifest genome build differs from the FASTA build", "IN-063"), + _m("W_SKIP_UNMAPPED", "warning", "Locus without a genomic position was skipped", "OUT-019"), + _m("W_SKIP_INTENSITY_ONLY", "warning", "Intensity-only locus was skipped", "OUT-019"), + _m("W_SKIP_INDEL_CONTEXT", "warning", "Indel sequence cannot be determined (no CSV SourceSeq)", "OUT-019"), + _m("W_SKIP_INDEL_REF", "warning", "Indel reference allele cannot be determined from the FASTA", "OUT-019"), + _m("W_INDEL_PARTIAL_MATCH", "warning", "Indel source sequence matches the genome only partially", "OUT-019"), + _m("W_SKIP_EXCLUDED", "warning", "Locus excluded by --exclude-loci/--include-loci", "OUT-019"), + _m("W_SKIP_CONTIG", "warning", "Manifest chromosome absent from the FASTA", "OUT-019"), + _m("W_SKIP_NO_STRAND", "warning", "Locus has neither RefStrand nor a sequence to align", "OUT-042"), + _m("W_SKIP_AMBIGUOUS", "warning", "Ambiguous-strand locus dropped by --drop-ambiguous", "OUT-042"), + _m("W_REF_NOT_QUERIED", "warning", "The reference allele is not among the probe's alleles", "OUT-015"), + _m("W_MERGE_INCONSISTENT", "warning", "Probes at a site disagree; genotype set missing", "OUT-031"), + _m("W_MERGE_AMBIGUOUS", "warning", "Probes at a site allow several genotypes; set missing", "OUT-031"), + _m("W_MERGE_ALL_NC", "warning", "Every probe at a merged site is a no-call", "OUT-031"), + _m("W_INF1_NC_CONFLICT", "warning", "Infinium I no-call conflicts with other assays", "OUT-037"), + _m("W_NO_BGZIP", "warning", "Uncompressed VCF written; downstream commands need bgzip + index", "OUT-010"), + _m("W_UNKNOWN_LOCI", "warning", "Names in a locus list are not in the manifest", "OUT-036"), + _m("E_AUX_INVALID", "error", "Auxiliary-loci record is invalid", "OUT-035"), + _m("W_LARGE_OUTPUT", "warning", "Output will be very large", "EDGE-036"), # --- QC (QC-001 .. QC-064) ----------------------------------------------------------------- _m("W_CTRL_MISSING", "warning", "Controls needed by a metric are absent; metric not evaluable", "QC-007"), _m("W_CTRL_ZERO_DIVISOR", "warning", "Control metric has a zero denominator (NaN)", "QC-008"), diff --git a/src/simurg/formats/bcf.py b/src/simurg/formats/bcf.py new file mode 100644 index 0000000..f63f4f0 --- /dev/null +++ b/src/simurg/formats/bcf.py @@ -0,0 +1,156 @@ +"""BCF 2.2 writer and CSI index (spec OUT-040, OUT-038). + +Records are built from typed values; the file is BGZF-compressed and indexed with a CSI index (min_shift 14, +depth 5, the pseudo-bin carries record counts). [DOC: VCF/BCF specification v4.3 §6 (BCF2), CSI index format — +layouts only] +""" + +from __future__ import annotations + +import struct +from collections.abc import Sequence + +from simurg.formats.bgzf import PSEUDO_BIN, BgzfWriter, bgzf_compress, reg2bin + +INT8, INT16, INT32, FLOAT, CHAR = 1, 2, 3, 5, 7 +FLOAT_MISSING = struct.pack(" bytes: + if n < 15: + return bytes([(n << 4) | t]) + return bytes([(15 << 4) | t]) + typed_ints([n]) + + +def typed_ints(values: Sequence[int]) -> bytes: + """Smallest integer type that holds every value (missing/vector-end sentinels excluded from the range).""" + lo, hi = (min(values), max(values)) if values else (0, 0) + if lo >= -120 and hi <= 127: + return _type_byte(len(values), INT8) + struct.pack(f"<{len(values)}b", *values) + if lo >= -32760 and hi <= 32767: + return _type_byte(len(values), INT16) + struct.pack(f"<{len(values)}h", *values) + return _type_byte(len(values), INT32) + struct.pack(f"<{len(values)}i", *values) + + +def typed_string(s: str) -> bytes: + b = s.encode() + return _type_byte(len(b), CHAR) + b + + +def float_bytes(values: Sequence[float | None]) -> bytes: + out = bytearray() + for v in values: + out += FLOAT_MISSING if v is None or v != v else struct.pack(" tuple[int, bytes]: + """FORMAT/GT for all samples: width = max ploidy; unphased (a+1)<<1; missing 0; haploid padded with END.""" + width = max(1, max(ploidies)) + vals: list[int] = [] + for g, p in zip(genotypes, ploidies, strict=True): + p = max(p, 1) + cells = [0] * p if g is None else [(a + 1) << 1 for a in g] + cells += [INT8_VECTOR_END] * (width - len(cells)) + vals += cells + return width, bytes(vals_to_int8(vals)) + + +def vals_to_int8(vals: Sequence[int]) -> bytes: + return struct.pack(f"<{len(vals)}b", *vals) + + +def encode_record( + contig_id: int, + pos0: int, + rlen: int, + vid: str, + alleles: Sequence[str], + filter_id: int, + n_sample: int, + fmt: Sequence[tuple[int, bytes, bytes]], +) -> bytes: + """``fmt``: (format dictionary index, type-descriptor bytes, per-sample payload) in FORMAT order.""" + shared = bytearray() + shared += struct.pack(" None: + self.refs: list[dict[str, object]] = [{} for _ in range(n_ref)] + + def add(self, rid: int, beg: int, end: int, vstart: int, vend: int) -> None: + r = self.refs[rid] + if not r: + r.update(bins={}, n=0, first=vstart, last=vend) + bins: dict[int, list[list[int]]] = r["bins"] # type: ignore[assignment] + b = reg2bin(beg, end) + chunks = bins.setdefault(b, []) + if chunks and chunks[-1][1] == vstart: + chunks[-1][1] = vend + else: + chunks.append([vstart, vend]) + r["n"] = int(r["n"]) + 1 # type: ignore[call-overload] + r["last"] = vend + + def encode(self) -> bytes: + out = bytearray(b"CSI\x01") + struct.pack(" None: + self.w = BgzfWriter(fh) # type: ignore[arg-type] + text = header_text.encode() + b"\0" + self.w.write_record(b"BCF\x02\x02" + struct.pack(" tuple[int, int]: + return self.w.write_record(record) + + def close(self) -> None: + self.w.close() + + +def header_dictionaries(header_text: str) -> tuple[dict[str, int], dict[str, int]]: + """(string dictionary: FILTER/INFO/FORMAT IDs → index, contig dictionary: name → index) as htslib builds them.""" + strings: dict[str, int] = {"PASS": 0} + contigs: dict[str, int] = {} + for line in header_text.splitlines(): + for kind in ("FILTER", "INFO", "FORMAT"): + if line.startswith(f"##{kind}=", 1)[0] + strings.setdefault(key, len(strings)) + if line.startswith("##contig=", 1)[0] + contigs.setdefault(name, len(contigs)) + return strings, contigs diff --git a/src/simurg/formats/bgzf.py b/src/simurg/formats/bgzf.py new file mode 100644 index 0000000..228a2df --- /dev/null +++ b/src/simurg/formats/bgzf.py @@ -0,0 +1,142 @@ +"""BGZF writer and tabix (``.tbi``) index for position-sorted text files (spec OUT-010, OUT-038). + +Pure Python, deterministic (fixed compression level, no timestamps). Records never straddle a block when they fit, +and the index carries the htslib pseudo-bin with mapped/unmapped counts so ``bcftools index --nrecords`` works. +[DOC: SAM/BAM format specification §4.1 (BGZF) and §5 (binning, linear index); tabix index format — layouts only] +""" + +from __future__ import annotations + +import struct +import zlib +from pathlib import Path +from typing import BinaryIO + +BLOCK_DATA = 0xFF00 # uncompressed bytes per block (as htslib) +EOF_BLOCK = bytes.fromhex("1f8b08040000000000ff0600424302001b0003000000000000000000") +TBI_MIN_SHIFT, TBI_DEPTH = 14, 5 +PSEUDO_BIN = 37450 + + +def _block(data: bytes, level: int) -> bytes: + c = zlib.compressobj(level, zlib.DEFLATED, -15) + comp = c.compress(data) + c.flush() + bsize = len(comp) + 25 + head = struct.pack("<4BIBBHBBHH", 0x1F, 0x8B, 8, 4, 0, 0, 0xFF, 6, ord("B"), ord("C"), 2, bsize) + return head + comp + struct.pack(" None: + self.fh, self.level = fh, level + self.buf = bytearray() + self.coffset = 0 + + def tell(self) -> int: + return (self.coffset << 16) | len(self.buf) + + def _flush_block(self) -> None: + if self.buf: + blk = _block(bytes(self.buf[:BLOCK_DATA]), self.level) + self.fh.write(blk) + self.coffset += len(blk) + del self.buf[:BLOCK_DATA] + + def write_record(self, data: bytes) -> tuple[int, int]: + """Append ``data``; start a new block first when it would not fit. Returns (start, end) virtual offsets.""" + if self.buf and len(self.buf) + len(data) > BLOCK_DATA: + self._flush_block() + start = self.tell() + self.buf += data + while len(self.buf) >= BLOCK_DATA: + self._flush_block() + return start, self.tell() + + def close(self) -> None: + while self.buf: + self._flush_block() + self.fh.write(EOF_BLOCK) + + +def bgzf_compress(data: bytes, level: int = 6) -> bytes: + out = bytearray() + for i in range(0, len(data), BLOCK_DATA): + out += _block(data[i : i + BLOCK_DATA], level) + return bytes(out + EOF_BLOCK) + + +def reg2bin(beg: int, end: int) -> int: + end -= 1 + if beg >> 14 == end >> 14: + return ((1 << 15) - 1) // 7 + (beg >> 14) + if beg >> 17 == end >> 17: + return ((1 << 12) - 1) // 7 + (beg >> 17) + if beg >> 20 == end >> 20: + return ((1 << 9) - 1) // 7 + (beg >> 20) + if beg >> 23 == end >> 23: + return ((1 << 6) - 1) // 7 + (beg >> 23) + if beg >> 26 == end >> 26: + return ((1 << 3) - 1) // 7 + (beg >> 26) + return 0 + + +class TabixBuilder: + """Collect (contig, beg0, end0, start voffset, end voffset) in file order and emit a VCF ``.tbi``.""" + + def __init__(self) -> None: + self.names: list[str] = [] + self.refs: dict[str, dict[str, object]] = {} + + def add(self, contig: str, beg: int, end: int, vstart: int, vend: int) -> None: + r = self.refs.get(contig) + if r is None: + self.names.append(contig) + r = self.refs[contig] = {"bins": {}, "lin": [], "n": 0, "first": vstart, "last": vend} + bins: dict[int, list[list[int]]] = r["bins"] # type: ignore[assignment] + b = reg2bin(beg, end) + chunks = bins.setdefault(b, []) + if chunks and chunks[-1][1] == vstart: + chunks[-1][1] = vend + else: + chunks.append([vstart, vend]) + lin: list[int] = r["lin"] # type: ignore[assignment] + for w in range(beg >> TBI_MIN_SHIFT, ((end - 1) >> TBI_MIN_SHIFT) + 1): + while len(lin) <= w: + lin.append(-1) + if lin[w] == -1: + lin[w] = vstart + r["n"] = int(r["n"]) + 1 # type: ignore[call-overload] + r["last"] = vend + + def encode(self) -> bytes: + names = b"".join(n.encode() + b"\0" for n in self.names) + out = bytearray(b"TBI\x01") + out += struct.pack("<8i", len(self.names), 2, 1, 2, 0, ord("#"), 0, len(names)) + names + for n in self.names: + r = self.refs[n] + bins: dict[int, list[list[int]]] = r["bins"] # type: ignore[assignment] + out += struct.pack(" None: + path.write_bytes(builder.encode()) diff --git a/src/simurg/formats/fasta.py b/src/simurg/formats/fasta.py new file mode 100644 index 0000000..ee13ee2 --- /dev/null +++ b/src/simurg/formats/fasta.py @@ -0,0 +1,350 @@ +"""Reference genome access (spec IN-060 … IN-064). + +* Plain or BGZF-compressed FASTA with a samtools-style ``.fai`` (and ``.gzi`` for BGZF). Missing indexes are built + into the cache directory, never next to the input. +* ``Genome.fetch`` returns upper-case bases for a 0-based half-open interval. With the genome cache (default) each + needed contig is stored once as a flat uint8 file under ``$SIMURG_CACHE_DIR/genomes//`` and + memory-mapped; without it, bases are read through the index. Both modes return identical sequences (EDGE-010). + +[DOC: samtools faidx / htslib BGZF and .gzi formats — layouts only] +""" + +from __future__ import annotations + +import hashlib +import json +import os +import struct +import zlib +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import numpy.typing as npt + +from simurg.core.errors import SimurgError + +GRCH37_CHR1 = 249_250_621 +GRCH38_CHR1 = 248_956_422 +_BGZF_MAGIC = b"\x1f\x8b\x08\x04" + + +@dataclass(frozen=True) +class FaiEntry: + name: str + length: int + offset: int # uncompressed byte offset of the first base + line_bases: int + line_width: int + + +def cache_root() -> Path: + return Path(os.environ.get("SIMURG_CACHE_DIR", Path.home() / ".cache" / "simurg")) + + +def is_bgzf(path: Path) -> bool: + with path.open("rb") as fh: + head = fh.read(16) + return head[:4] == _BGZF_MAGIC and head[12:14] == b"BC" + + +def file_sha256(path: Path) -> str: + """SHA-256 of a (large) file, memoised by path, size and mtime in the cache directory.""" + st = path.stat() + memo = cache_root() / "sha256.json" + key = f"{path.resolve()}|{st.st_size}|{st.st_mtime_ns}" + try: + known = json.loads(memo.read_text()) + except (OSError, ValueError): + known = {} + if key in known: + return str(known[key]) + h = hashlib.sha256() + with path.open("rb") as fh: + for chunk in iter(lambda: fh.read(1 << 22), b""): + h.update(chunk) + known[key] = h.hexdigest() + memo.parent.mkdir(parents=True, exist_ok=True) + tmp = memo.with_suffix(".part") + tmp.write_text(json.dumps(known, indent=0)) + tmp.replace(memo) + return str(known[key]) + + +# ------------------------------------------------------------------------------------------------ BGZF +def read_gzi(path: Path) -> list[tuple[int, int]]: + """``.gzi``: (compressed offset, uncompressed offset) of every BGZF block after the first.""" + data = path.read_bytes() + (n,) = struct.unpack_from(" list[tuple[int, int]]: + out = [] + coff = uoff = 0 + with fasta.open("rb") as fh: + while True: + head = fh.read(18) + if not head: + break + if head[:4] != _BGZF_MAGIC: + raise SimurgError("E_FASTA_FORMAT", f"{fasta}: not a BGZF file (bad block at {coff})") + bsize = struct.unpack_from(" None: + buf = struct.pack(" None: + self.fh = path.open("rb") + self.blocks = [(0, 0), *gzi] + self.ustarts = np.array([u for _, u in self.blocks], dtype=np.int64) + self._cache: tuple[int, bytes] | None = None + + def _block(self, i: int) -> bytes: + if self._cache and self._cache[0] == i: + return self._cache[1] + self.fh.seek(self.blocks[i][0]) + head = self.fh.read(18) + bsize = struct.unpack_from(" bytes: + out = bytearray() + i = int(np.searchsorted(self.ustarts, uoff, side="right") - 1) + pos = uoff - int(self.ustarts[i]) + while len(out) < n and i < len(self.blocks): + data = self._block(i) + out += data[pos : pos + n - len(out)] + i, pos = i + 1, 0 + return bytes(out) + + +class _PlainReader: + def __init__(self, path: Path) -> None: + self.fh = path.open("rb") + + def read(self, uoff: int, n: int) -> bytes: + self.fh.seek(uoff) + return self.fh.read(n) + + +def _stream_uncompressed(path: Path, bgzf: bool) -> object: + if not bgzf: + return path.open("rb") + import gzip + + return gzip.open(path, "rb") + + +# ------------------------------------------------------------------------------------------------ FAI +def read_fai(path: Path) -> list[FaiEntry]: + out = [] + for i, line in enumerate(path.read_text().splitlines(), start=1): + if not line.strip(): + continue + f = line.split("\t") + try: + out.append(FaiEntry(f[0], int(f[1]), int(f[2]), int(f[3]), int(f[4]))) + except (IndexError, ValueError): + raise SimurgError("E_FAI_MISMATCH", f"{path}: malformed line {i}", path=str(path)) from None + return out + + +def build_fai(fasta: Path, bgzf: bool) -> list[FaiEntry]: + """samtools-compatible index; offsets are uncompressed byte offsets.""" + entries: list[FaiEntry] = [] + fh = _stream_uncompressed(fasta, bgzf) + off = 0 + name: str | None = None + length = seq_off = lb = lw = 0 + with fh: # type: ignore[attr-defined] + for raw in fh: # type: ignore[attr-defined] + if raw.startswith(b">"): + if name is not None: + entries.append(FaiEntry(name, length, seq_off, lb, lw)) + name = raw[1:].split()[0].decode() + length, lb, lw = 0, 0, 0 + seq_off = off + len(raw) + else: + n = len(raw.rstrip(b"\r\n")) + if lb == 0: + lb, lw = n, len(raw) + length += n + off += len(raw) + if name is not None: + entries.append(FaiEntry(name, length, seq_off, lb, lw)) + if not entries: + raise SimurgError("E_FASTA_FORMAT", f"{fasta}: no sequences found", path=str(fasta)) + return entries + + +def write_fai(path: Path, entries: list[FaiEntry]) -> None: + path.write_text("".join(f"{e.name}\t{e.length}\t{e.offset}\t{e.line_bases}\t{e.line_width}\n" for e in entries)) + + +def detect_build(entries: list[FaiEntry]) -> str: + """IN-063: chr1 length decides GRCh37 vs GRCh38; '' when unknown.""" + for e in entries: + if e.name in ("1", "chr1"): + return {GRCH37_CHR1: "GRCh37", GRCH38_CHR1: "GRCh38"}.get(e.length, "") + return "" + + +# ------------------------------------------------------------------------------------------------ genome +class Genome: + """Indexed reference genome; see module docstring.""" + + def __init__(self, fasta: str | Path, *, use_cache: bool = True, cache_dir: Path | None = None) -> None: + self.path = Path(fasta).resolve() + if not self.path.exists(): + raise SimurgError("E_FILE_NOT_FOUND", f"FASTA not found: {fasta}", path=str(fasta)) + self.bgzf = is_bgzf(self.path) + self.use_cache = use_cache + self._sha: str | None = None + self._cache_dir = cache_dir + self.notes: list[str] = [] + fai = Path(str(self.path) + ".fai") + if fai.exists(): + self.entries = read_fai(fai) + else: + built = self.cache_dir / "index.fai" + if not built.exists(): + built.parent.mkdir(parents=True, exist_ok=True) + write_fai(built, build_fai(self.path, self.bgzf)) + self.notes.append(f"built FASTA index {built}") + self.entries = read_fai(built) + self.index = {e.name: e for e in self.entries} + if self.bgzf: + gzi = Path(str(self.path) + ".gzi") + if not gzi.exists(): + gzi = self.cache_dir / "index.gzi" + if not gzi.exists(): + gzi.parent.mkdir(parents=True, exist_ok=True) + write_gzi(gzi, build_gzi(self.path)) + self._reader: _BgzfReader | _PlainReader = _BgzfReader(self.path, read_gzi(gzi)) + else: + self._reader = _PlainReader(self.path) + self._mm: dict[str, npt.NDArray[np.uint8]] = {} + self.verify() + self.build = detect_build(self.entries) + + @property + def sha256(self) -> str: + if self._sha is None: + self._sha = file_sha256(self.path) + return self._sha + + @property + def cache_dir(self) -> Path: + if self._cache_dir is None: + self._cache_dir = cache_root() / "genomes" / self.sha256 + return self._cache_dir + + def verify(self, n: int = 3) -> None: + """IN-060: sample up to ``n`` contigs and check header position and line layout.""" + picks = self.entries[:: max(1, len(self.entries) // n)][:n] + for e in picks: + name_bytes = e.name.encode() + before = self._reader.read(max(0, e.offset - len(name_bytes) - 64), min(e.offset, len(name_bytes) + 64)) + first = self._reader.read(e.offset, e.line_width) + ok = before.endswith(b"\n") and b">" + name_bytes in before + if e.length >= e.line_bases: + ok = ok and len(first.rstrip(b"\r\n")) == e.line_bases and first.endswith(b"\n") + if not ok: + raise SimurgError( + "E_FAI_MISMATCH", + f"FASTA index does not match {self.path} at contig {e.name}", + hint="delete the stale .fai (it will be rebuilt in the cache) or regenerate it with samtools faidx", + path=str(self.path), + ) + + def contig_length(self, name: str) -> int: + return self.index[name].length + + def _read_raw(self, e: FaiEntry, start: int, end: int) -> bytes: + lines_before = start // e.line_bases + uoff = e.offset + lines_before * e.line_width + start % e.line_bases + last = end - 1 + uend = e.offset + (last // e.line_bases) * e.line_width + last % e.line_bases + 1 + raw = self._reader.read(uoff, uend - uoff) + return raw.replace(b"\n", b"").replace(b"\r", b"") + + def _contig_array(self, name: str) -> npt.NDArray[np.uint8]: + if name not in self._mm: + e = self.index[name] + f = self.cache_dir / f"{name}.u8" + if not f.exists() or f.stat().st_size != e.length: + f.parent.mkdir(parents=True, exist_ok=True) + tmp = f.with_suffix(".part") + with tmp.open("wb") as out: + step = 1 << 24 + for s in range(0, e.length, step): + out.write(self._read_raw(e, s, min(e.length, s + step)).upper()) + tmp.replace(f) + self._mm[name] = np.memmap(f, dtype=np.uint8, mode="r") + return self._mm[name] + + def fetch(self, name: str, start: int, end: int) -> str: + """Upper-case bases of ``name``[start:end] (0-based, half-open), clipped to the contig.""" + e = self.index[name] + start, end = max(0, start), min(e.length, end) + if end <= start: + return "" + if self.use_cache: + return bytes(self._contig_array(name)[start:end]).decode("ascii") + return self._read_raw(e, start, end).decode("ascii").upper() + + +# ------------------------------------------------------------------------------------------------ contig names +def load_alias_file(path: Path) -> dict[str, str]: + out = {} + for line in path.read_text(encoding="utf-8").splitlines(): + line = line.split("#", 1)[0].strip() + if line: + a, b = line.split("\t")[:2] + out[a] = b + return out + + +class ContigMap: + """IN-061: manifest chromosome → FASTA contig (``XY`` → X). Policies: ``strict`` (verbatim), ``auto`` + (``1``↔``chr1``, ``MT``↔``chrM``, ``22_KI270879v1_alt``↔``chr22_KI270879v1_alt``) or an alias TSV path.""" + + def __init__(self, genome_contigs: set[str] | list[str], policy: str = "strict") -> None: + self.contigs = set(genome_contigs) + self.policy = policy + self.alias = load_alias_file(Path(policy)) if policy not in ("strict", "auto") else {} + self._memo: dict[str, str | None] = {} + + def get(self, chrom: str) -> str | None: + if chrom in self._memo: + return self._memo[chrom] + target = "X" if chrom == "XY" else chrom + if self.alias: + cands = [self.alias.get(chrom, ""), self.alias.get(target, ""), target] + elif self.policy == "auto": + bare = target[3:] if target.lower().startswith("chr") else target + cands = [target, bare, f"chr{bare}"] + if bare in ("MT", "M"): + cands += ["chrM", "MT", "M"] + else: + cands = [target] + hit = next((c for c in cands if c and c in self.contigs), None) + self._memo[chrom] = hit + return hit diff --git a/src/simurg/harmonize/combine.py b/src/simurg/harmonize/combine.py new file mode 100644 index 0000000..be20d0d --- /dev/null +++ b/src/simurg/harmonize/combine.py @@ -0,0 +1,121 @@ +"""Genotypes of VCF sites from probe calls (spec OUT-018, OUT-031, OUT-032, OUT-037; plan §8.3). + +Single-probe sites map the A/B call straight to allele indices. Sites with several probes use the exact candidate +algorithm of OUT-031: enumerate every genotype over {REF} ∪ ALT with the site's ploidy, predict each probe's call +from the alleles' signal classes (Infinium II: dye group of the extended base; Infinium I and indels: exact allele, +no signal otherwise) and keep the candidates consistent with every called probe. +[OSS: Illumina GTCtoVCF ``GenotypeFormat.RecordCombiner`` (Apache-2.0) — the consistency idea; the class model and +the ploidy-aware enumeration are ours.] +""" + +from __future__ import annotations + +from dataclasses import dataclass +from itertools import combinations_with_replacement + +from simurg.harmonize.sites import GREEN, RED, ProbeSite, Site + +NC, AA, AB, BB = 0, 1, 2, 3 +CALL_TEXT = {NC: "NC", AA: "AA", AB: "AB", BB: "BB"} + + +def signal_class(p: ProbeSite, allele: str) -> str | None: + """OUT-031 step 1.""" + a, b = p.alleles + if p.kind == "snv" and p.infinium == 2 and len(allele) == 1: + grp = RED if allele in RED else GREEN + if a in grp: + return "A" + if b in grp: + return "B" + return None + if allele == a: + return "A" + if allele == b: + return "B" + return None + + +def predicted_call(p: ProbeSite, genotype: tuple[str, ...]) -> int: + classes = sorted(c for c in (signal_class(p, z) for z in genotype) if c is not None) + if not classes: + return NC + u = set(classes) + if u == {"A"}: + return AA + if u == {"B"}: + return BB + return AB + + +@dataclass(frozen=True) +class Outcome: + alleles: tuple[int, ...] | None # allele indices (sorted) or None for missing + reason: str # "" | "het_in_haploid" | W_MERGE_* | W_INF1_NC_CONFLICT + + +def single(site: Site, call: int, ploidy: int) -> Outcome: + p = site.probes[0] + if call == NC: + return Outcome(None, "") + ia, ib = site.allele_index(p.alleles[0]), site.allele_index(p.alleles[1]) + if ploidy == 1: + if call == AB: + return Outcome(None, "het_in_haploid") + return Outcome((ia if call == AA else ib,), "") + pair = {AA: (ia, ia), AB: (ia, ib), BB: (ib, ib)}[call] + return Outcome(tuple(sorted(pair)), "") + + +def merged(site: Site, calls: tuple[int, ...], ploidy: int, inf1_nc_informative: bool = False) -> Outcome: + """OUT-031 steps 3–6.""" + if all(c == NC for c in calls): + return Outcome(None, "W_MERGE_ALL_NC") + alleles = [site.ref, *site.alts] + cands = list(combinations_with_replacement(range(len(alleles)), ploidy)) + keep = [ + g + for g in cands + if all( + c == NC or predicted_call(p, tuple(alleles[i] for i in g)) == c + for p, c in zip(site.probes, calls, strict=True) + ) + ] + reason = "" + if inf1_nc_informative: + for p, c in zip(site.probes, calls, strict=True): + if c == NC and p.infinium == 1: + sub = [g for g in keep if all(signal_class(p, alleles[i]) is None for i in g)] + if sub: + keep = sub + else: + reason = "W_INF1_NC_CONFLICT" + if len(keep) == 1: + return Outcome(keep[0], reason) + return Outcome(None, "W_MERGE_INCONSISTENT" if not keep else "W_MERGE_AMBIGUOUS") + + +def genotype(site: Site, calls: tuple[int, ...], ploidy: int, inf1_nc_informative: bool = False) -> Outcome: + if len(site.probes) == 1: + return single(site, calls[0], ploidy) + key = (calls, ploidy, inf1_nc_informative) + hit = site.memo.get(key) + if hit is None: + hit = site.memo[key] = merged(site, calls, ploidy, inf1_nc_informative) + assert isinstance(hit, Outcome) + return hit + + +def gt_text(o: Outcome, ploidy: int) -> str: + if o.alleles is None: + return "." if ploidy == 1 else "./." + return "/".join(str(i) for i in o.alleles) + + +def failure_message(site: Site, calls: tuple[int, ...], reason: str) -> str: + what = {"W_MERGE_INCONSISTENT": "inconsistent", "W_MERGE_AMBIGUOUS": "ambiguous", "W_MERGE_ALL_NC": "all no-calls"} + probes = ", ".join( + f"{p.name} ({'InfiniumI' if p.infinium == 1 else 'InfiniumII'}): {CALL_TEXT[c]}" + for p, c in zip(site.probes, calls, strict=True) + ) + return f"Failed to combine genotypes ({what.get(reason, reason)}) - {probes}" diff --git a/src/simurg/harmonize/sites.py b/src/simurg/harmonize/sites.py new file mode 100644 index 0000000..b9fa189 --- /dev/null +++ b/src/simurg/harmonize/sites.py @@ -0,0 +1,307 @@ +"""Manifest loci → plus-strand VCF sites, with an audit row per locus (spec IN-046, OUT-014 … OUT-021, OUT-030, +OUT-042; plan §8.1–§8.2). + +[DOC: Illumina "TOP/BOT strand and A/B allele" technical note — strand conventions] +[OSS: Illumina GTCtoVCF ``BPMRecord.py`` (Apache-2.0) — indel source-sequence handling: design-strand orientation, +left shift, IUPAC-aware flank matching and the deletion-versus-insertion context score; ported below.] +""" + +from __future__ import annotations + +import re +from collections.abc import Callable, Sequence +from dataclasses import dataclass, field + +COMPLEMENT = dict(zip("ABCDGHKMRTVYNID", "TVGHCDMKYABRNID", strict=True)) +DEGENERACY = {"A": "A", "C": "C", "G": "G", "T": "T", "R": "AG", "Y": "CT", "S": "GC", "W": "AT", "K": "GT", + "M": "AC", "B": "CGT", "D": "AGT", "H": "ACT", "V": "ACG", "N": "ACGT"} # fmt: skip +RED, GREEN = set("AT"), set("CG") +SEQ_RE = re.compile(r"^([^\[]*)\[([^/\]]*)/([^\]]*)\](.*)$") + +Fetch = Callable[[str, int, int], str] # (contig, start0, end0) -> upper-case bases + + +def revcomp(s: str) -> str: + return "".join(COMPLEMENT.get(c, "N") for c in reversed(s)) + + +def complement_base(b: str) -> str: + return COMPLEMENT.get(b, "N") + + +@dataclass +class ProbeSite: + """One manifest locus mapped to the plus strand.""" + + index: int + name: str + contig: str + pos: int # VCF POS (1-based) + ref: str + alleles: tuple[str, str] # VCF allele strings of the probe's A and B alleles + kind: str # "snv" | "indel" + infinium: int # 1 or 2 + manifest_chrom: str + par: bool = False + confidence: str = "high" + flipped: bool = False # design strand is minus + reason: str = "ok" + + +@dataclass +class AuditRow: + locus: str + original_alleles: str + strand_in: str + transformed_alleles: str + strand_out: str + ref: str + alt: str + flipped: bool + swapped: bool + reason: str + build: str + confidence: str + + +@dataclass +class Site: + """A VCF record: one or more probes that share CHROM/POS/REF (and the indel span).""" + + contig: str + pos: int + ref: str + alts: list[str] + probes: list[ProbeSite] = field(default_factory=list) + kind: str = "snv" + memo: dict[object, object] = field(default_factory=dict, repr=False, compare=False) # merged-genotype cache + + @property + def id(self) -> str: + return (";" if self.kind == "indel" else ",").join(p.name for p in self.probes) + + def allele_index(self, allele: str) -> int: + return 0 if allele == self.ref else 1 + self.alts.index(allele) + + +# ------------------------------------------------------------------------------------------------ indel helpers +def split_source(seq: str) -> tuple[str, str, str, str] | None: + """``left[a/b]right`` → (left, a, b, right); None when malformed.""" + m = SEQ_RE.match(seq.strip().upper()) + return (m.group(1), m.group(2), m.group(3), m.group(4)) if m else None + + +def left_shift(five: str, indel: str, three: str) -> tuple[str, str]: + """Shift an indel fully 5′ within its context (GTCtoVCF ``determine_left_shift``).""" + while indel and five.endswith(indel): + five, three = five[: -len(indel)], indel + three + while indel and five and len(indel) * five[-1] == indel: + three, five = five[-1] + three, five[:-1] + return five, three + + +def _match(a: str, b: str) -> bool: + return a in DEGENERACY.get(b, "") + + +def max_suffix_match(genome: str, design: str) -> int: + n = 0 + for x, y in zip(genome[::-1], design[::-1], strict=False): + if not _match(x, y): + break + n += 1 + return n + + +def max_prefix_match(genome: str, design: str) -> int: + n = 0 + for x, y in zip(genome, design, strict=False): + if not _match(x, y): + break + n += 1 + return n + + +def classify_indel(fetch: Fetch, contig: str, pos: int, five: str, indel: str, three: str) -> tuple[bool | None, float]: + """Is the plus-strand ``indel`` sequence present in the reference at ``pos`` (a reference deletion)? + Returns (is_deletion or None when undecidable, best context score). Context scoring as GTCtoVCF.""" + start = pos - 1 + L = len(indel) + seq_match = fetch(contig, start, start + L) == indel + d5, d3 = left_shift(fetch(contig, start - len(five), start), indel, + fetch(contig, start + L, start + L + len(three))) # fmt: skip + i5, i3 = left_shift(fetch(contig, start - len(five) + 1, start + 1), indel, + fetch(contig, start + 1, start + len(three) + 1)) # fmt: skip + dm = (max_suffix_match(d5, five), max_prefix_match(d3, three)) + dmax = min(len(d5), len(five)) + min(len(d3), len(three)) + L + dscore = (sum(dm) + L if seq_match else 0) / dmax if dmax else 0.0 + im = (max_suffix_match(i5, five), max_prefix_match(i3, three)) + imax = min(len(i5), len(five)) + min(len(i3), len(three)) + iscore = sum(im) / imax if imax else 0.0 + is_del = seq_match and dscore > iscore and min(dm) >= 1 + is_ins = iscore > dscore and min(im) >= 1 + if is_del == is_ins: + return None, max(dscore, iscore) + return is_del, dscore if is_del else iscore + + +def normalize_indel(fetch: Fetch, contig: str, pos: int, indel: str, is_deletion: bool) -> tuple[int, str, str]: + """Left-normalised, anchored VCF representation. ``pos`` is the manifest position: the first deleted base for + a deletion, the base 5′ of the insertion for an insertion (IN-046).""" + if is_deletion: + p0 = pos - 1 # 0-based first deleted base + while p0 > 0 and fetch(contig, p0 - 1, p0) == indel[-1]: + indel = indel[-1] + indel[:-1] + p0 -= 1 + anchor = fetch(contig, p0 - 1, p0) + return p0, anchor + indel, anchor # VCF POS = p0 (1-based anchor) + a0 = pos - 1 # 0-based anchor base + while a0 > 0 and fetch(contig, a0, a0 + 1) == indel[-1]: + indel = indel[-1] + indel[:-1] + a0 -= 1 + anchor = fetch(contig, a0, a0 + 1) + return a0 + 1, anchor, anchor + indel + + +# ------------------------------------------------------------------------------------------------ strand +def strand_from_flanks(fetch: Fetch, contig: str, pos: int, left: str, right: str) -> tuple[str, float]: + """'+' when the design flanks match the reference as written, '-' when reverse-complemented; with a score.""" + left, right = left[-30:], right[:30] + g_left = fetch(contig, pos - 1 - len(left), pos - 1) + g_right = fetch(contig, pos, pos + len(right)) + + def score(lf: str, rt: str) -> float: + n = len(lf) + len(rt) + if not n: + return 0.0 + hits = sum(_match(a, b) for a, b in zip(g_left[::-1], lf[::-1], strict=False)) + hits += sum(_match(a, b) for a, b in zip(g_right, rt, strict=False)) + return hits / n + + fwd, rev = score(left, right), score(revcomp(right), revcomp(left)) + return ("+", fwd) if fwd >= rev else ("-", rev) + + +# ------------------------------------------------------------------------------------------------ mapping +@dataclass +class ManifestView: + """The manifest columns the mapper needs (BPM and/or CSV).""" + + names: Sequence[str] + snps: Sequence[str] + chroms: Sequence[str] + positions: Sequence[int] + ilmn_strands: Sequence[str] + source_strands: Sequence[str] + ref_strands: Sequence[str] + source_seqs: Sequence[str] + assay_types: Sequence[int] # 0 Infinium II, 1/2 Infinium I + intensity_only: Sequence[bool] + + +def map_locus( + i: int, mv: ManifestView, fetch: Fetch, contig_of: Callable[[str], str | None], build: str = "" +) -> tuple[ProbeSite | None, AuditRow]: + name, snp, chrom = mv.names[i], mv.snps[i], mv.chroms[i] + pos = int(mv.positions[i]) + ilmn, ref_strand = mv.ilmn_strands[i] or "", (mv.ref_strands[i] or "").strip() + audit = AuditRow(name, snp, f"{ilmn}/{ref_strand or '?'}", "", "+", "", "", False, False, "ok", build, "high") + + def skip(code: str) -> tuple[None, AuditRow]: + audit.reason, audit.strand_out = code, "" + return None, audit + + if chrom in ("0", "") or pos <= 0: + return skip("W_SKIP_UNMAPPED") + if mv.intensity_only[i]: + return skip("W_SKIP_INTENSITY_ONLY") + contig = contig_of(chrom) + if contig is None: + return skip("W_SKIP_CONTIG") + m = re.match(r"^\[([^/\]]+)/([^\]]+)\]$", snp.strip()) + if not m: + return skip("W_SKIP_UNMAPPED") + a, b = m.group(1).upper(), m.group(2).upper() + infinium = 2 if int(mv.assay_types[i]) == 0 else 1 + par = chrom == "XY" + src = mv.source_seqs[i] if i < len(mv.source_seqs) else "" + parts = split_source(src) if src else None + if a in ("I", "D") or b in ("I", "D"): + # ---------------------------------------------------------------- indel + if parts is None or "-" not in (parts[1], parts[2]): + return skip("W_SKIP_INDEL_CONTEXT") + five, x, y, three = parts + indel_seq = y if x == "-" else x + s_strand = (mv.source_strands[i] or "")[:1].upper() + d_strand = ilmn[:1].upper() + if not ref_strand or s_strand in ("U", "") or d_strand in ("U", ""): + return skip("W_SKIP_INDEL_CONTEXT") + on_design = s_strand == d_strand + if on_design != (ref_strand == "+"): + five, indel_seq, three = revcomp(three), revcomp(indel_seq), revcomp(five) + five, three = left_shift(five, indel_seq, three) + is_del, score = classify_indel(fetch, contig, pos, five, indel_seq, three) + if is_del is None: + return skip("W_SKIP_INDEL_REF") + vpos, vref, valt = normalize_indel(fetch, contig, pos, indel_seq, is_del) + if not vref or "N" in vref: + return skip("W_SKIP_INDEL_REF") + # I (insertion allele) carries the sequence; D lacks it + allele_of = {"I": vref if is_del else valt, "D": valt if is_del else vref} + alleles = (allele_of[a], allele_of[b]) + audit.reason = "ok" if score >= 1.0 else "W_INDEL_PARTIAL_MATCH" + audit.transformed_alleles = f"{alleles[0]}/{alleles[1]}" + audit.ref, audit.alt = vref, valt + audit.swapped = alleles[0] != vref + return ProbeSite(i, name, contig, vpos, vref, alleles, "indel", infinium, chrom, par, + "high", ref_strand == "-", audit.reason), audit # fmt: skip + # -------------------------------------------------------------------- SNV + confidence = "high" + strand = ref_strand + if strand not in ("+", "-"): + if parts is None: + return skip("W_SKIP_NO_STRAND") # no RefStrand and no sequence to align + five, x, y, three = parts + src_plus, score = strand_from_flanks(fetch, contig, pos, five, three) + same = (mv.source_strands[i] or "")[:1].upper() == ilmn[:1].upper() + strand = src_plus if same else ("-" if src_plus == "+" else "+") + if score < 0.9: + confidence = "low" + if {a, b} in ({"A", "T"}, {"C", "G"}): + confidence = "low" # plan §8.2: ambiguous without RefStrand + pa, pb = (a, b) if strand == "+" else (complement_base(a), complement_base(b)) + ref = fetch(contig, pos - 1, pos) + if not ref or ref not in "ACGT": + return skip("W_SKIP_UNMAPPED") + audit.transformed_alleles = f"{pa}/{pb}" + audit.ref = ref + audit.alt = ",".join(x for x in (pa, pb) if x != ref) + audit.flipped = strand == "-" + audit.swapped = pa != ref + audit.confidence = confidence + if ref not in (pa, pb): + audit.reason = "W_REF_NOT_QUERIED" + return ProbeSite(i, name, contig, pos, ref, (pa, pb), "snv", infinium, chrom, par, confidence, + strand == "-", audit.reason), audit # fmt: skip + + +def group_sites(probes: Sequence[ProbeSite], contig_order: dict[str, int], *, unsquash: bool = False) -> list[Site]: + """OUT-030/OUT-033/OUT-034: probes sharing CHROM/POS/REF (and, for indels, the ALT) form one site; ALT alleles in + order of first appearance over BPM order; records sorted by contig (.fai order), POS, REF, ALT, ID.""" + groups: dict[tuple[str, int, str, str, str], list[ProbeSite]] = {} + for p in sorted(probes, key=lambda p: p.index): + extra = ",".join(sorted(x for x in p.alleles if x != p.ref)) if p.kind == "indel" else "" + groups.setdefault((p.contig, p.pos, p.ref, p.kind, extra), []).append(p) + sites: list[Site] = [] + for (contig, pos, ref, kind, _), ps in groups.items(): + alts: list[str] = [] + for p in ps: + for x in p.alleles: + if x != ref and x not in alts: + alts.append(x) + if unsquash and len(ps) > 1 and len(alts) <= 1: + sites += [Site(contig, pos, ref, [x for x in p.alleles if x != ref], [p], kind) for p in ps] + else: + sites.append(Site(contig, pos, ref, alts, ps, kind)) + sites.sort(key=lambda s: (contig_order.get(s.contig, 1 << 30), s.pos, s.ref, ",".join(s.alts), s.id)) + return sites diff --git a/src/simurg/pipeline/export.py b/src/simurg/pipeline/export.py new file mode 100644 index 0000000..b0acb58 --- /dev/null +++ b/src/simurg/pipeline/export.py @@ -0,0 +1,422 @@ +"""``simurg export bedgraph | finalreport | plink | locus-summary`` (spec CLI-023, OUT-043, OUT-050 … OUT-052, +OUT-071, OUT-072). None of them re-genotypes: they read GTCs (or SimurgArray VCFs) plus the manifest. +""" + +from __future__ import annotations + +import gzip +import io +import math +import re +from collections.abc import Iterator, Sequence +from dataclasses import dataclass, field +from pathlib import Path +from typing import IO, Any, cast + +import numpy as np +import numpy.typing as npt + +from simurg import PRODUCT_NAME, __version__ +from simurg.core.atomic import atomic_open, check_outputs_absent +from simurg.core.errors import ExitCode, SimurgError +from simurg.core.log import Reporter +from simurg.core.provenance import RunManifest +from simurg.formats.bpm import Manifest +from simurg.formats.samplesheet import output_stem + +_COMP = bytes.maketrans(b"ACGTDI", b"TGCADI") +LARGE_OUTPUT = 100_000_000 + + +# ------------------------------------------------------------------------------------------------ inputs +@dataclass +class GtcSample: + sample_id: str + sample_name: str + path: str + reported_sex: str = "U" + + +def gtc_samples(gtcs: Path | None, sheet_path: Path | None, rep: Reporter) -> list[GtcSample]: + """Sheet order when a sheet is given, else sorted sample IDs (OUT-076).""" + from simurg.formats.discovery import discover_gtcs + from simurg.formats.samplesheet import parse_samplesheet + + sheet = parse_samplesheet(sheet_path) if sheet_path else None + root = gtcs or (sheet.root_folder if sheet else None) + if root is None: + raise SimurgError("E_NO_INPUT", "give --gtcs or a sample sheet with a RootFolder") + found = discover_gtcs(root) + if sheet is None: + return [GtcSample(k, k, str(v)) for k, v in sorted(found.items())] + out = [] + for r in sheet.rows: + if r.sample_id in found: + out.append(GtcSample(r.sample_id, r.sample_name or r.sample_id, str(found[r.sample_id]), r.reported_sex)) + else: + rep.warn("W_QC_UNMATCHED_SAMPLES", f"{r.sample_id}: no GTC found", sample_id=r.sample_id) + if not out: + raise SimurgError("E_NO_INPUT", "none of the sheet's samples has a GTC") + return out + + +def chrom_rank(chrom: str) -> tuple[int, str]: + c = chrom[3:] if chrom.lower().startswith("chr") else chrom + if c.isdigit(): + return (int(c) if c != "0" else 1000, "") + return ({"X": 23, "XY": 24, "Y": 25, "MT": 26, "M": 26}.get(c, 100), c) + + +def genome_order(m: Manifest) -> npt.NDArray[np.int64]: + """Loci by (chromosome, position, name); unmapped (Chr 0) last.""" + keys = [ + ((1000, "") if c in ("0", "") or p <= 0 else chrom_rank(c), int(p), n) + for c, p, n in zip(m.chroms, m.map_infos.tolist(), m.names, strict=True) + ] + return np.array(sorted(range(len(keys)), key=keys.__getitem__), dtype=np.int64) + + +# ------------------------------------------------------------------------------------------------ bedGraph +def _chrom_label(chrom: str, style: str) -> str: + c = "X" if chrom == "XY" else chrom + if style == "ucsc": + c = "M" if c == "MT" else c + return c if c.startswith("chr") else f"chr{c}" + return c + + +def moving_median(v: npt.NDArray[np.float64], k: int) -> npt.NDArray[np.float64]: + if k <= 0 or len(v) == 0: + return v.copy() + pad = np.pad(v, k, mode="edge") + win = np.lib.stride_tricks.sliding_window_view(pad, 2 * k + 1) + return np.nanmedian(win, axis=1) + + +@dataclass +class BedgraphLayout: + """Sample-independent part of OUT-050: loci in output order with their labels (computed once per manifest).""" + + order: npt.NDArray[np.int64] # mapped loci, genome order + chrom: list[str] # output chromosome label per ordered locus + pos: npt.NDArray[np.int64] + + +def bedgraph_layout(m: Manifest, style: str) -> BedgraphLayout: + order = genome_order(m) + mapped = np.array([c not in ("0", "") for c in m.chroms]) & (m.map_infos > 0) + order = order[mapped[order]] + return BedgraphLayout(order, [_chrom_label(m.chroms[i], style) for i in order.tolist()], m.map_infos[order]) + + +def bedgraph_rows(m: Manifest, g: Any, layout: BedgraphLayout, which: str) -> list[tuple[str, int, float]]: + """OUT-050 rows (chrom, pos, value) for one track, in output order, same-position probes averaged.""" + vals = g.b_allele_freqs if which == "baf" else g.logr_ratios + if len(vals) != m.num_loci: + raise SimurgError("E_GTC_UNSUITABLE", "GTC has no BAF/LRR (GTC v3?)") + o = layout.order + v = vals[o].astype(np.float64) + keep = ~((g.raw_x[o] == 0) & (g.raw_y[o] == 0)) & np.isfinite(v) + if which == "baf": + keep &= ~m.intensity_only[o] + idx = np.flatnonzero(keep) + if not len(idx): + return [] + chrom = [layout.chrom[i] for i in idx.tolist()] + pos = layout.pos[idx] + new = np.ones(len(idx), np.bool_) + new[1:] = (pos[1:] != pos[:-1]) | np.array([chrom[k] != chrom[k - 1] for k in range(1, len(chrom))]) + starts = np.flatnonzero(new) + sums = np.add.reduceat(v[idx], starts) + counts = np.diff(np.append(starts, len(idx))) + means = sums / counts + return [(chrom[k], int(pos[k]), float(x)) for k, x in zip(starts.tolist(), means.tolist(), strict=True)] + + +def write_bedgraph(path: Path, name: str, which: str, rows: Sequence[tuple[str, int, float]], gz: bool) -> None: + limits = "0:1" if which == "baf" else "-2:2" + head = ( + f'track type=bedGraph name="{name} {which.upper()}" description="{PRODUCT_NAME} {__version__}" ' + f"visibility=full autoScale=off viewLimits={limits}\n" + ) + body = "".join(f"{c}\t{p - 1}\t{p}\t{v:.6g}\n" for c, p, v in rows) + data = (head + body).encode() + with atomic_open(path, "wb") as fh: + cast(IO[bytes], fh).write(gzip.compress(data, mtime=0) if gz else data) + + +@dataclass +class BedgraphOptions: + out: Path + bpm: Path + gtcs: Path | None = None + sheet: Path | None = None + gzip: bool = True + chrom_style: str = "manifest" + smoothed_lrr: int = 0 + force: bool = False + argv: list[str] = field(default_factory=list) + + +def run_export_bedgraph(opts: BedgraphOptions, rep: Reporter) -> int: + from simurg.formats.bpm import read_bpm + from simurg.formats.gtc import read_gtc + + if opts.chrom_style not in ("manifest", "ucsc"): + raise SimurgError("E_CONFIG_KEY", "--chrom-style must be manifest or ucsc") + rm = RunManifest(command="export bedgraph", argv=opts.argv) + samples = gtc_samples(opts.gtcs, opts.sheet, rep) + ext = ".bedgraph.gz" if opts.gzip else ".bedgraph" + outs = [] + for s in samples: + stem = output_stem(s.sample_name) + outs += [opts.out / f"{stem}.baf{ext}", opts.out / f"{stem}.lrr{ext}"] + if opts.smoothed_lrr: + outs.append(opts.out / f"{stem}.lrr_smoothed{ext}") + check_outputs_absent(outs, force=opts.force) + m = read_bpm(opts.bpm) + layout = bedgraph_layout(m, opts.chrom_style) + opts.out.mkdir(parents=True, exist_ok=True) + failed = 0 + for s in samples: + try: + g = read_gtc(s.path) + stem = output_stem(s.sample_name) + for which in ("baf", "lrr"): + rows = bedgraph_rows(m, g, layout, which) + write_bedgraph(opts.out / f"{stem}.{which}{ext}", s.sample_name, which, rows, opts.gzip) + if which == "lrr" and opts.smoothed_lrr: + sm: list[tuple[str, int, float]] = [] + for chrom in dict.fromkeys(c for c, _, _ in rows): + sub = [r for r in rows if r[0] == chrom] + med = moving_median(np.array([r[2] for r in sub]), opts.smoothed_lrr) + sm += [(c, p, float(v)) for (c, p, _), v in zip(sub, med, strict=True)] + write_bedgraph(opts.out / f"{stem}.lrr_smoothed{ext}", s.sample_name, "lrr", sm, opts.gzip) + except SimurgError as e: + failed += 1 + rep.error(e.code, f"{s.sample_id}: {e.message}", sample_id=s.sample_id) + rm.status = "completed" + rm.summary = {"samples": len(samples), "failed": failed} + rm.write(opts.out / "run_manifest.export_bedgraph.json") + return ExitCode.OK if not failed else (ExitCode.PARTIAL if failed < len(samples) else ExitCode.ALL_FAILED) + + +# ------------------------------------------------------------------------------------------------ final report +FR_COLUMNS = ["SNP Name", "SNP", "Sample ID", "Allele1 - Top", "Allele2 - Top", "Allele1 - Forward", + "Allele2 - Forward", "Allele1 - Plus", "Allele2 - Plus", "GC Score", "GT Score", "Log R Ratio", + "B Allele Freq", "Chr", "Position"] # fmt: skip + + +def allele_letters(m: Manifest) -> dict[str, list[tuple[str, str]]]: + """Per locus (A, B) letters on the TOP, Forward (source) and Plus (reference) strands.""" + top, fwd, plus = [], [], [] + for snp, ilmn, src, ref in zip(m.snps, m.ilmn_strands, m.source_strands, m.ref_strands, strict=True): + mm = re.match(r"^\[([^/\]]+)/([^\]]+)\]$", snp.strip()) + a, b = (mm.group(1), mm.group(2)) if mm else ("-", "-") + + def comp(x: str) -> str: + return x.encode().translate(_COMP).decode() if len(x) == 1 else x + + d = (ilmn or "")[:1].upper() + top.append((a, b) if d in ("T", "P", "") or len(a) > 1 else (comp(a), comp(b))) + same_src = (src or "")[:1].upper() == d or not src + fwd.append((a, b) if same_src else (comp(a), comp(b))) + plus.append((a, b) if (ref or "+").strip() == "+" else (comp(a), comp(b))) + return {"top": top, "fwd": fwd, "plus": plus} + + +def _pair(call: int, ab: tuple[str, str]) -> tuple[str, str]: + a, b = ab + return {1: (a, a), 2: (a, b), 3: (b, b)}.get(call, ("-", "-")) + + +@dataclass +class FinalReportOptions: + out: Path + bpm: Path + gtcs: Path | None = None + sheet: Path | None = None + csv: Path | None = None + egt: Path | None = None + matrix: bool = False + force: bool = False + argv: list[str] = field(default_factory=list) + + +def run_export_finalreport(opts: FinalReportOptions, rep: Reporter) -> int: + from simurg.formats.bpm import read_bpm + from simurg.formats.gtc import read_gtc + + rm = RunManifest(command="export finalreport", argv=opts.argv) + samples = gtc_samples(opts.gtcs, opts.sheet, rep) + target = opts.out / ("final_report_matrix.tsv.gz" if opts.matrix else "final_report.tsv.gz") + check_outputs_absent([target], force=opts.force) + m = read_bpm(opts.bpm) + if opts.csv is not None: + from simurg.formats.csv_manifest import read_csv_manifest + + cm = read_csv_manifest(opts.csv, load_sequences=False) + m.source_strands, m.ref_strands = cm.source_strands, cm.ref_strands + cells = m.num_loci * len(samples) + if cells > LARGE_OUTPUT: + rep.warn("W_LARGE_OUTPUT", f"{cells:,} rows (~{cells * 90 / 1e9:.1f} GB uncompressed); consider --matrix") + gt_score = np.full(m.num_loci, np.nan) + if opts.egt is not None: + from simurg.formats.egt import read_egt + + egt = read_egt(opts.egt) + idx = egt.name_index() + for i, nm in enumerate(m.names): + j = idx.get(nm) + if j is not None: + gt_score[i] = float(egt.gentrain_score[j]) + letters = allele_letters(m) + order = genome_order(m).tolist() + opts.out.mkdir(parents=True, exist_ok=True) + + def f4(v: float) -> str: + return "NaN" if math.isnan(v) else f"{v:.4f}" + + with atomic_open(target, "wb") as raw: + gz = gzip.GzipFile(fileobj=cast(IO[bytes], raw), mode="wb", mtime=0) + fh = io.TextIOWrapper(gz, encoding="utf-8", newline="\n") + fh.write("#simurg_schema=simurg.finalreport/1\n") + if opts.matrix: + fh.write("\t".join(["SNP Name", "Chr", "Position", *[s.sample_name for s in samples]]) + "\n") + calls = [] + for s in samples: + calls.append(read_gtc(s.path).genotypes) + for i in order: + row = [m.names[i], m.chroms[i], str(int(m.map_infos[i]))] + row += ["".join(_pair(int(c[i]), letters["top"][i])) for c in calls] + fh.write("\t".join(row) + "\n") + else: + fh.write("\t".join(FR_COLUMNS) + "\n") + for s in samples: + g = read_gtc(s.path) + has = len(g.b_allele_freqs) == m.num_loci + bc = g.base_calls + for i in order: + c = int(g.genotypes[i]) + f, p = (_pair(c, letters[k][i]) for k in ("fwd", "plus")) + t = (chr(bc[i, 0]), chr(bc[i, 1])) if c else ("-", "-") # TOP-strand base calls (OUT-071) + baf = float(g.b_allele_freqs[i]) if has and not m.intensity_only[i] else math.nan + lrr = float(g.logr_ratios[i]) if has else math.nan + fh.write( + f"{m.names[i]}\t{m.snps[i]}\t{s.sample_id}\t{t[0]}\t{t[1]}\t{f[0]}\t{f[1]}\t{p[0]}\t{p[1]}\t" + f"{f4(float(g.genotype_scores[i]))}\t{f4(gt_score[i])}\t{f4(lrr)}\t{f4(baf)}\t" + f"{m.chroms[i]}\t{int(m.map_infos[i])}\n" + ) + fh.flush() + fh.detach() + gz.close() + rm.status = "completed" + rm.summary = {"samples": len(samples), "loci": m.num_loci} + rm.write(opts.out / "run_manifest.export_finalreport.json") + return ExitCode.OK + + +# ------------------------------------------------------------------------------------------------ PLINK +PLINK_MAGIC = bytes([0x6C, 0x1B, 0x01]) + + +def parse_vcf_records(path: Path) -> Iterator[list[str]]: + opener = gzip.open if str(path).endswith(".gz") else open + with opener(path, "rt") as fh: + for line in fh: + if not line.startswith("#"): + yield line.rstrip("\n").split("\t") + + +def bed_code(gt: str, alt_index: int) -> int: + """PLINK 1 two-bit code with A1 = the split ALT, A2 = REF: 0 hom A1, 1 missing, 2 het, 3 hom A2.""" + if gt in (".", "./."): + return 1 + al = gt.replace("|", "/").split("/") + if "." in al: + return 1 + n_alt = sum(int(x) == alt_index for x in al) + if len(al) == 1: + n_alt *= 2 + return {0: 3, 1: 2, 2: 0}[n_alt] + + +@dataclass +class PlinkOptions: + out: Path + vcfs: Path + sheet: Path | None = None + gtcs: Path | None = None + multiallelic: str = "split" + prefix: str = "cohort" + force: bool = False + argv: list[str] = field(default_factory=list) + + +def run_export_plink(opts: PlinkOptions, rep: Reporter) -> int: + """OUT-043: SimurgArray per-sample VCFs (identical record order) → PLINK 1 .bed/.bim/.fam.""" + if opts.multiallelic not in ("split", "drop"): + raise SimurgError("E_CONFIG_KEY", "--multiallelic must be split or drop") + rm = RunManifest(command="export plink", argv=opts.argv) + vcfs = sorted(p for p in opts.vcfs.glob("*.snv.vcf*") if not p.name.endswith((".tbi", ".csi"))) + if not vcfs: + raise SimurgError("E_NO_INPUT", f"no *.snv.vcf(.gz) files in {opts.vcfs}") + base = opts.out / opts.prefix + check_outputs_absent([base.with_suffix(x) for x in (".bed", ".bim", ".fam")], force=opts.force) + names, sex = [], {} + for v in vcfs: + opener = gzip.open if str(v).endswith(".gz") else open + with opener(v, "rt") as fh: + for line in fh: + if line.startswith("#CHROM"): + names.append(line.rstrip("\n").split("\t")[9]) + break + if opts.gtcs is not None: + from simurg.formats.discovery import discover_gtcs + from simurg.formats.gtc import read_gtc + + for key, gpath in discover_gtcs(opts.gtcs).items(): + sex[key] = read_gtc(gpath, check_write_complete=False).gender_char + if opts.sheet is not None: + from simurg.formats.samplesheet import parse_samplesheet + + for r in parse_samplesheet(opts.sheet).rows: + if r.reported_sex in ("M", "F"): + sex[r.sample_name or r.sample_id] = sex[r.sample_id] = r.reported_sex + streams = [parse_vcf_records(v) for v in vcfs] + n = len(vcfs) + nbytes = (n + 3) // 4 + opts.out.mkdir(parents=True, exist_ok=True) + n_var = dropped = 0 + with ( + atomic_open(base.with_suffix(".bed"), "wb") as bed_raw, + atomic_open(base.with_suffix(".bim"), "w", "utf-8") as bim_raw, + ): + bed, bim = cast(IO[bytes], bed_raw), cast(IO[str], bim_raw) + bed.write(PLINK_MAGIC) + for recs in zip(*streams, strict=True): + r0 = recs[0] + if any(r[:5] != r0[:5] for r in recs[1:]): + raise SimurgError("E_GTC_UNSUITABLE", f"VCF record order differs between samples at {r0[0]}:{r0[1]}") + alts = [] if r0[4] == "." else r0[4].split(",") + if len(alts) > 1 and opts.multiallelic == "drop": + dropped += 1 + continue + gts = [r[9].split(":", 1)[0] for r in recs] + for k, alt in enumerate(alts, start=1): + codes = [bed_code(g, k) for g in gts] + buf = bytearray(nbytes) + for j, c in enumerate(codes): + buf[j >> 2] |= c << (2 * (j & 3)) + bed.write(bytes(buf)) + vid = r0[2] if len(alts) == 1 else f"{r0[2]}:{alt}" + bim.write(f"{r0[0]}\t{vid}\t0\t{r0[1]}\t{alt}\t{r0[3]}\n") + n_var += 1 + fam = "".join(f"{s}\t{s}\t0\t0\t{ {'M': 1, 'F': 2}.get(sex.get(s, 'U'), 0) }\t-9\n" for s in names) + with atomic_open(base.with_suffix(".fam"), "w", "utf-8") as fam_raw: + cast(IO[str], fam_raw).write(fam) + rm.status = "completed" + rm.summary = {"samples": n, "variants": n_var, "multiallelic_dropped": dropped} + rm.write(opts.out / "run_manifest.export_plink.json") + rep.info(f"PLINK fileset {base}.bed/.bim/.fam: {n} samples, {n_var} variants") + return ExitCode.OK diff --git a/src/simurg/pipeline/qc.py b/src/simurg/pipeline/qc.py index 9ee54ab..fbb1d60 100644 --- a/src/simurg/pipeline/qc.py +++ b/src/simurg/pipeline/qc.py @@ -308,26 +308,27 @@ class SampleInfo: class CohortStore: """Disk-backed (samples × loci) matrices: genotypes (uint8), BAF and GenCall scores (float32, optional).""" - def __init__(self, root: Path, n_samples: int, n_loci: int, *, baf: bool, scores: bool) -> None: + def __init__(self, root: Path, n_samples: int, n_loci: int, *, baf: bool, scores: bool, lrr: bool = False) -> None: self.root, self.shape = root, (n_samples, n_loci) root.mkdir(parents=True, exist_ok=True) self.gt = np.lib.format.open_memmap(root / "gt.npy", "w+", np.uint8, self.shape) self.baf = np.lib.format.open_memmap(root / "baf.npy", "w+", np.float32, self.shape) if baf else None self.score = np.lib.format.open_memmap(root / "score.npy", "w+", np.float32, self.shape) if scores else None + self.lrr = np.lib.format.open_memmap(root / "lrr.npy", "w+", np.float32, self.shape) if lrr else None def close(self) -> None: - for name in ("gt", "baf", "score"): + for name in ("gt", "baf", "score", "lrr"): a = getattr(self, name) if a is not None: a.flush() - self.gt = self.baf = self.score = None # type: ignore[assignment] + self.gt = self.baf = self.score = self.lrr = None # type: ignore[assignment] shutil.rmtree(self.root, ignore_errors=True) -def _load_one(task: tuple[int, str, str, str, bool, bool]) -> SampleInfo: +def _load_one(task: tuple[int, str, str, str, bool, bool, bool]) -> SampleInfo: from simurg.formats.gtc import read_gtc - i, sid, path, root, want_baf, want_scores = task + i, sid, path, root, want_baf, want_scores, want_lrr = task info = SampleInfo(sid, path) try: g = read_gtc(path) @@ -347,6 +348,10 @@ def _load_one(task: tuple[int, str, str, str, bool, bool]) -> SampleInfo: info.has_baf = len(g.b_allele_freqs) == gt.shape[1] baf[i] = g.b_allele_freqs if info.has_baf else np.nan baf.flush() + if want_lrr: + lr = np.load(Path(root) / "lrr.npy", mmap_mode="r+") + lr[i] = g.logr_ratios if len(g.logr_ratios) == gt.shape[1] else np.nan + lr.flush() if want_scores: sc = np.load(Path(root) / "score.npy", mmap_mode="r+") sc[i] = g.genotype_scores @@ -359,12 +364,19 @@ def _load_one(task: tuple[int, str, str, str, bool, bool]) -> SampleInfo: def load_cohort( - samples: Sequence[tuple[str, str]], n_loci: int, tmp: Path, *, baf: bool, scores: bool, workers: int + samples: Sequence[tuple[str, str]], + n_loci: int, + tmp: Path, + *, + baf: bool, + scores: bool, + workers: int, + lrr: bool = False, ) -> tuple[CohortStore, list[SampleInfo]]: """Read GTCs (sample_id, path) into a :class:`CohortStore` with a process pool.""" - store = CohortStore(tmp, len(samples), n_loci, baf=baf, scores=scores) + store = CohortStore(tmp, len(samples), n_loci, baf=baf, scores=scores, lrr=lrr) store.gt.flush() - tasks = [(i, sid, p, str(tmp), baf, scores) for i, (sid, p) in enumerate(samples)] + tasks = [(i, sid, p, str(tmp), baf, scores, lrr) for i, (sid, p) in enumerate(samples)] if workers <= 1 or len(tasks) <= 1: infos = [_load_one(t) for t in tasks] else: diff --git a/src/simurg/pipeline/vcf.py b/src/simurg/pipeline/vcf.py new file mode 100644 index 0000000..80fc012 --- /dev/null +++ b/src/simurg/pipeline/vcf.py @@ -0,0 +1,656 @@ +"""``simurg gtc-to-vcf``: GTC + manifest + FASTA → one bgzipped, indexed SNV VCF per sample (spec CLI-022, +OUT-010 … OUT-042, IN-060 … IN-064; plan §8). + +The manifest is mapped to plus-strand sites once (``simurg.harmonize``); each sample is then written by a worker +process with vectorised genotype lookup for single-probe sites and the OUT-031 algorithm for merged sites. +""" + +from __future__ import annotations + +import hashlib +import math +import multiprocessing as mp +import os +from collections import Counter +from concurrent.futures import ProcessPoolExecutor +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt + +from simurg import PRODUCT_NAME, __version__ +from simurg.core.atomic import atomic_open, check_outputs_absent +from simurg.core.errors import ExitCode, SimurgError +from simurg.core.log import Reporter +from simurg.core.provenance import RunManifest, sha256_file +from simurg.formats.bgzf import BgzfWriter, TabixBuilder, write_tbi +from simurg.harmonize.combine import genotype as combine_genotype +from simurg.harmonize.combine import gt_text +from simurg.harmonize.sites import AuditRow, ManifestView, Site, group_sites, map_locus + +PAR = { # 1-based inclusive pseudo-autosomal intervals on X (IN-065, human) + "GRCh38": ((10_001, 2_781_479), (155_701_383, 156_030_895)), + "GRCh37": ((60_001, 2_699_520), (154_931_044, 155_260_560)), +} +CLASS_AUTO, CLASS_X, CLASS_Y, CLASS_MT, CLASS_PAR = 0, 1, 2, 3, 4 + + +@dataclass +class VcfOptions: + out: Path + bpm: Path + fasta: Path + gtcs: Path | None = None + sheet: Path | None = None + csv: Path | None = None + bgzip: bool = True + genome_cache: bool = True + exclude_loci: Path | None = None + include_loci: Path | None = None + unsquash_duplicates: bool = False + inf1_nocall_informative: bool = False + drop_ambiguous: bool = False + filter_gs: float | None = None + contig_alias: str = "strict" + ploidy_policy: str = "sex-aware" + threads: int = 0 + merge: bool = False + force: bool = False + argv: list[str] = field(default_factory=list) + + +# ------------------------------------------------------------------------------------------------ site table +@dataclass +class SiteTable: + """Everything a worker needs, precomputed once (picklable).""" + + prefix: list[bytes] # "CHROM\tPOS\tID\tREF\tALT\t.\t" per site + single_rows: npt.NDArray[np.int64] # site indices of single-probe sites + single_locus: npt.NDArray[np.int64] # their manifest locus index + single_gt: npt.NDArray[np.object_] # (n_single, 2 ploidy, 4 calls) GT text + single_flip: npt.NDArray[np.bool_] # BAF must be 1 − BAF + merged_sites: list[Site] + merged_rows: npt.NDArray[np.int64] + site_class: npt.NDArray[np.int8] + contigs: list[str] + header: str + site_beg: npt.NDArray[np.int64] + site_end: npt.NDArray[np.int64] + site_contig: list[str] + site_id: list[str] = field(default_factory=list) + site_alleles: list[list[str]] = field(default_factory=list) + single_ab: npt.NDArray[np.int16] | None = None # (n_single, 2) allele indices of the probe's A and B alleles + inf1_nc_informative: bool = False + ploidy_policy: str = "sex-aware" + filter_gs: float | None = None + + +def _site_class(contig: str, pos: int, par_flag: bool, build: str) -> int: + c = contig[3:] if contig.lower().startswith("chr") else contig + if c in ("MT", "M"): + return CLASS_MT + if c == "Y": + return CLASS_Y + if c == "X": + if par_flag or any(a <= pos <= b for a, b in PAR.get(build or "GRCh38", ())): + return CLASS_PAR + return CLASS_X + return CLASS_AUTO + + +def ploidy_for(cls: npt.NDArray[np.int8], sex: str, policy: str) -> npt.NDArray[np.int8]: + """OUT-018 per site; 0 marks 'haploid and always missing' (Y of female samples).""" + p = np.full(len(cls), 2, dtype=np.int8) + p[cls == CLASS_MT] = 1 + if policy == "diploid": + return p + p[cls == CLASS_Y] = 0 if sex == "F" else 1 + if sex == "M": + p[cls == CLASS_X] = 1 + return p + + +def read_locus_list(path: Path) -> list[str]: + out = [] + for line in path.read_text(encoding="utf-8").splitlines(): + line = line.split("#", 1)[0].strip() + if line: + out.append(line) + return out + + +def build_sites( + opts: VcfOptions, rep: Reporter, rm: RunManifest +) -> tuple[SiteTable, list[AuditRow], Counter[str], Any]: + from simurg.formats.bpm import read_bpm + from simurg.formats.fasta import ContigMap, Genome + + bpm = read_bpm(opts.bpm) + csvm = None + if opts.csv is not None: + from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest + + csvm = read_csv_manifest(opts.csv) + problems = compare_manifests(bpm, csvm) + if problems: + raise SimurgError("E_BPM_CSV_MISMATCH", "BPM and CSV manifests differ:\n " + "\n ".join(problems)) + genome = Genome(opts.fasta, use_cache=opts.genome_cache) + for n in genome.notes: + rep.info(n) + builds = Counter(bpm.genome_builds) + mbuild = builds.most_common(1)[0][0] if builds else "" + want = {"37": "GRCh37", "38": "GRCh38"}.get(mbuild.strip(), mbuild.strip()) + if genome.build and want and want != genome.build: + raise SimurgError( + "E_BUILD_MISMATCH", f"manifest genome build {want} but FASTA looks like {genome.build} ({opts.fasta})" + ) + build = genome.build or want + src = csvm or bpm + mv = ManifestView( + bpm.names, + bpm.snps, + bpm.chroms, + bpm.map_infos.tolist(), + src.ilmn_strands, + src.source_strands, + src.ref_strands if any(src.ref_strands) else bpm.ref_strands, + csvm.source_seqs if csvm else [""] * bpm.num_loci, + bpm.assay_types.tolist(), + bpm.intensity_only.tolist(), + ) + if csvm is None and any(("D" in s) for s in bpm.snps): + rep.error("W_SKIP_INDEL_CONTEXT", "no CSV manifest: indels are skipped (give --csv with SourceSeq)") + names = set(bpm.names) + excl: set[str] = set() + incl: set[str] | None = None + for path, kind in ((opts.exclude_loci, "exclude"), (opts.include_loci, "include")): + if path is None: + continue + lst = read_locus_list(path) + unknown = [x for x in lst if x not in names] + if unknown: + rep.warn("W_UNKNOWN_LOCI", f"{len(unknown)} name(s) in --{kind}-loci are not in the manifest") + if kind == "exclude": + excl = set(lst) + else: + incl = set(lst) + cmap = ContigMap([e.name for e in genome.entries], opts.contig_alias) + probes = [] + audit: list[AuditRow] = [] + skips: Counter[str] = Counter() + for i in range(bpm.num_loci): + nm = bpm.names[i] + if nm in excl or (incl is not None and nm not in incl): + a = AuditRow(nm, bpm.snps[i], "", "", "", "", "", False, False, "W_SKIP_EXCLUDED", build, "") + p = None + else: + p, a = map_locus(i, mv, genome.fetch, cmap.get, build) + if p is not None and opts.drop_ambiguous and p.confidence == "low": + a.reason, p = "W_SKIP_AMBIGUOUS", None + audit.append(a) + if p is None or a.reason.startswith("W_"): + skips[a.reason] += 1 + if p is not None: + probes.append(p) + for code, count in sorted(skips.items()): + rep.warn(code, f"{count} locus/loci: {code}") + order = {e.name: k for k, e in enumerate(genome.entries)} + sites = group_sites(probes, order, unsquash=opts.unsquash_duplicates) + used = {s.contig for s in sites} + contigs = [e.name for e in genome.entries if e.name in used] + # static per-site text and single-probe genotype tables + prefix, cls = [], np.zeros(len(sites), np.int8) + single_rows, single_locus, single_flip, gts, single_ab = [], [], [], [], [] + merged_sites, merged_rows = [], [] + beg = np.zeros(len(sites), np.int64) + end = np.zeros(len(sites), np.int64) + for k, s in enumerate(sites): + prefix.append(f"{s.contig}\t{s.pos}\t{s.id}\t{s.ref}\t{','.join(s.alts) if s.alts else '.'}\t.\t".encode()) + cls[k] = _site_class(s.contig, s.pos, any(p.par for p in s.probes), build) + beg[k], end[k] = s.pos - 1, s.pos - 1 + len(s.ref) + if len(s.probes) == 1: + from simurg.harmonize.combine import single + + single_rows.append(k) + single_locus.append(s.probes[0].index) + al_a, al_b = s.probes[0].alleles + single_flip.append(al_a != s.ref and al_b == s.ref) + single_ab.append((s.allele_index(al_a), s.allele_index(al_b))) + gts.append([[gt_text(single(s, c, pl), pl) for c in range(4)] for pl in (1, 2)]) + else: + merged_sites.append(s) + merged_rows.append(k) + header = vcf_header(opts, genome, build, contigs, rm) + table = SiteTable( + prefix=prefix, + single_rows=np.array(single_rows, np.int64), + single_locus=np.array(single_locus, np.int64), + single_gt=np.array(gts, dtype=object).reshape(len(single_rows), 2, 4), + single_flip=np.array(single_flip, np.bool_), + merged_sites=merged_sites, + merged_rows=np.array(merged_rows, np.int64), + site_class=cls, + contigs=contigs, + header=header, + site_beg=beg, + site_end=end, + site_contig=[s.contig for s in sites], + site_id=[s.id for s in sites], + site_alleles=[[s.ref, *s.alts] for s in sites], + single_ab=np.array(single_ab, np.int16).reshape(-1, 2), + inf1_nc_informative=opts.inf1_nocall_informative, + ploidy_policy=opts.ploidy_policy, + filter_gs=opts.filter_gs, + ) + return table, audit, skips, genome + + +def _sha(path: Path | None) -> str: + return sha256_file(path) if path else "" + + +def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm: RunManifest) -> str: + """OUT-011 (without the ##SAMPLE and #CHROM lines, which are per sample).""" + + def lst(p: Path | None) -> str: + return "none" if p is None else f"{p.name}:{hashlib.sha256(p.read_bytes()).hexdigest()[:16]}" + + lines = [ + "##fileformat=VCFv4.3", + f"##source={PRODUCT_NAME} {__version__}", + f"##reference=file://{Path(opts.fasta).resolve()}", + f"##simurg_genome_build={build}", + f"##simurg_manifest={opts.bpm.name};sha256={_sha(opts.bpm)}", + ] + if opts.csv: + lines.append(f"##simurg_csv_manifest={opts.csv.name};sha256={_sha(opts.csv)}") + lines.append("##simurg_cluster_file=__CLUSTER__") + lines.append( + "##simurg_options=" + f"unsquash_duplicates={str(opts.unsquash_duplicates).lower()};exclude_loci={lst(opts.exclude_loci)};" + f"include_loci={lst(opts.include_loci)};inf1_nocall_informative={str(opts.inf1_nocall_informative).lower()};" + f"auxiliary_loci=none;ploidy_policy={opts.ploidy_policy};contig_alias={opts.contig_alias}" + + (";drop_ambiguous=true" if opts.drop_ambiguous else "") + + (";unknown_sex_policy=diploid_x" if opts.ploidy_policy == "sex-aware" else "") + ) + lines.append(f"##simurg_run={rm.run_id}") + for c in contigs: + lines.append(f"##contig=") + lines.append('##FILTER=') + if opts.filter_gs is not None: + lines.append(f'##FILTER=') + lines += [ + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + ] + return "\n".join(lines) + "\n" + + +# ------------------------------------------------------------------------------------------------ per sample +_W: dict[str, Any] = {} + + +def _init(table: SiteTable) -> None: + _W["t"] = table + + +def _f4(v: float) -> str: + return "." if math.isnan(v) else f"{v:.4f}" + + +def _g8(v: float) -> str: + return "." if math.isnan(v) else f"{v:.8g}" + + +def write_sample(task: tuple[str, str, str, str, bool]) -> dict[str, Any]: + """Worker: one GTC → one VCF (+ index). Returns counts for the run summary.""" + from simurg.formats.gtc import read_gtc + + sample_id, sample_name, gtc_path, out_path, bgzip = task + t: SiteTable = _W["t"] + res: dict[str, Any] = {"sample_id": sample_id, "records": 0, "reasons": Counter(), "merge_failures": []} + try: + g = read_gtc(gtc_path) + n = len(t.prefix) + sex = g.gender_char + ploidy = ploidy_for(t.site_class, sex, t.ploidy_policy) + gt = np.empty(n, dtype=object) + gs = np.full(n, np.nan, np.float64) + baf = np.full(n, np.nan, np.float64) + lrr = np.full(n, np.nan, np.float64) + has_baf = len(g.b_allele_freqs) == g.num_snps + # single-probe sites (vectorised) + rows, loc = t.single_rows, t.single_locus + calls = g.genotypes[loc].astype(np.int64) + pl = ploidy[rows] + pidx = np.where(pl == 2, 1, 0) + gt[rows] = t.single_gt[np.arange(len(rows)), pidx, calls] + gt[rows[pl == 0]] = "." + gs[rows] = g.genotype_scores[loc].astype(np.float64) + if has_baf: + b = g.b_allele_freqs[loc].astype(np.float32) + b = np.where(t.single_flip, (np.float32(1.0) - b).astype(np.float32), b) + baf[rows] = b.astype(np.float64) + lrr[rows] = g.logr_ratios[loc].astype(np.float64) + het_haploid = int(((calls == 2) & (pl == 1)).sum()) + if het_haploid: + res["reasons"]["het_in_haploid"] += het_haploid + # merged sites + for s, k in zip(t.merged_sites, t.merged_rows.tolist(), strict=True): + c = tuple(int(g.genotypes[p.index]) for p in s.probes) + p_k = int(ploidy[k]) + gs[k] = float(min(float(g.genotype_scores[p.index]) for p in s.probes)) + if p_k == 0: + gt[k] = "." + continue + o = combine_genotype(s, c, p_k, t.inf1_nc_informative) + gt[k] = gt_text(o, p_k) + if o.reason: + res["reasons"][o.reason] += 1 + if o.reason.startswith("W_MERGE") and o.reason != "W_MERGE_ALL_NC": + res["merge_failures"].append((s.id, c, o.reason)) + merged_mask = np.zeros(n, np.bool_) + merged_mask[t.merged_rows] = True + header = t.header.replace("__CLUSTER__", g.cluster_file or "") + ( + f"##SAMPLE=\n" + f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}\n" + ) + low = t.filter_gs + gsl, bafl, lrrl = gs.tolist(), baf.tolist(), lrr.tolist() + tbx = TabixBuilder() + out = Path(out_path) + with atomic_open(out, "wb") as raw: + if bgzip: + w = BgzfWriter(raw) # type: ignore[arg-type] + w.write_record(header.encode()) + else: + raw.write(header.encode()) # type: ignore[arg-type] + for k in range(n): + filt = "LowGS" if low is not None and not math.isnan(gsl[k]) and gsl[k] < low else "PASS" + if merged_mask[k]: + body = f"{filt}\t.\tGT:GS\t{gt[k]}:{_f4(gsl[k])}\n" + else: + body = f"{filt}\t.\tGT:GS:BAF:LRR\t{gt[k]}:{_f4(gsl[k])}:{_g8(bafl[k])}:{_g8(lrrl[k])}\n" + line = t.prefix[k] + body.encode() + if bgzip: + vs, ve = w.write_record(line) + tbx.add(t.site_contig[k], int(t.site_beg[k]), int(t.site_end[k]), vs, ve) + else: + raw.write(line) # type: ignore[arg-type] + if bgzip: + w.close() + if bgzip: + write_tbi(Path(str(out) + ".tbi"), tbx) + res["records"] = n + res["sex"] = sex + except SimurgError as e: + res["error"] = f"{e.code}: {e.message}" + except Exception as e: + res["error"] = f"E_INTERNAL: {type(e).__name__}: {e}" + return res + + +# ------------------------------------------------------------------------------------------------ command +def run_gtc_to_vcf(opts: VcfOptions, rep: Reporter) -> int: + from simurg.formats.discovery import discover_gtcs + from simurg.formats.samplesheet import output_stem, parse_samplesheet + from simurg.pipeline.inputs import check_paths + + rm = RunManifest(command="gtc-to-vcf", argv=opts.argv) + rep.run_id = rm.run_id + check_paths( + [ + (opts.bpm, "BPM"), + (opts.fasta, "FASTA"), + (opts.csv, "CSV manifest"), + (opts.gtcs, "GTC folder"), + (opts.sheet, "sample sheet"), + (opts.exclude_loci, "locus list"), + (opts.include_loci, "locus list"), + ] + ) + if opts.ploidy_policy not in ("sex-aware", "diploid"): + raise SimurgError("E_CONFIG_KEY", "--ploidy-policy must be sex-aware or diploid") + sheet = parse_samplesheet(opts.sheet) if opts.sheet else None + root = opts.gtcs or (sheet.root_folder if sheet else None) + if root is None: + raise SimurgError("E_NO_INPUT", "give --gtcs or a sample sheet with a RootFolder") + found = discover_gtcs(root) + samples: list[tuple[str, str, str]] = [] # (sample_id, sample_name, path) + if sheet is not None: + for r in sheet.rows: + if r.sample_id in found: + samples.append((r.sample_id, r.sample_name or r.sample_id, str(found[r.sample_id]))) + else: + rep.warn("W_QC_UNMATCHED_SAMPLES", f"{r.sample_id}: no GTC under {root}", sample_id=r.sample_id) + else: + samples = [(k, k, str(v)) for k, v in sorted(found.items())] + suffix = ".snv.vcf.gz" if opts.bgzip else ".snv.vcf" + outs = ( + [opts.out / "cohort.snv.bcf"] + if opts.merge + else [opts.out / f"{output_stem(name)}{suffix}" for _, name, _ in samples] + ) + check_outputs_absent([*outs, opts.out / "audit" / "loci.parquet"], force=opts.force) + if not opts.bgzip: + rep.warn("W_NO_BGZIP", "writing plain VCF without index; pgx and cn commands need bgzipped, indexed VCFs") + rep.info(f"mapping manifest {opts.bpm.name} to {opts.fasta.name}") + table, audit, skips, _genome = build_sites(opts, rep, rm) + for role, p in (("bpm", opts.bpm), ("csv_manifest", opts.csv), ("fasta", opts.fasta), ("sample_sheet", opts.sheet)): + if p is not None: + rm.add_input(role, p) + opts.out.mkdir(parents=True, exist_ok=True) + write_audit(opts.out / "audit" / "loci.parquet", audit) + if opts.merge: + n = write_cohort_bcf(table, samples, opts.out / "cohort.snv.bcf", opts, rep) + rm.summary = {"samples": len(samples), "records": n, "skipped_loci_by_code": dict(skips)} + rm.add_output(opts.out / "cohort.snv.bcf", "bcf/2.2") + rm.status = "completed" + rm.write(opts.out / "run_manifest.json") + rep.info(f"done: cohort BCF with {len(samples)} sample(s), {n} records") + return ExitCode.OK + workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(samples))) + tasks = [(sid, name, path, str(o), opts.bgzip) for (sid, name, path), o in zip(samples, outs, strict=True)] + rep.info(f"{len(table.prefix)} VCF records; writing {len(tasks)} sample(s) with {workers} worker(s)") + if workers == 1: + _init(table) + results = [write_sample(t) for t in tasks] + else: + with ProcessPoolExecutor( + workers, mp_context=mp.get_context("spawn"), initializer=_init, initargs=(table,) + ) as ex: + results = list(ex.map(write_sample, tasks)) + reasons: Counter[str] = Counter() + failures = [] + for res in results: + if "error" in res: + rep.error(res["error"].split(":", 1)[0], f"{res['sample_id']}: {res['error']}", sample_id=res["sample_id"]) + continue + reasons.update(res["reasons"]) + failures += [(res["sample_id"], *f) for f in res["merge_failures"]] + for code, n in sorted(reasons.items()): + if code.startswith("W_"): + rep.warn(code, f"{n} genotype(s): {code}") + write_merge_audit(opts.out / "audit" / "merges.parquet", failures) + ok = sum("error" not in r for r in results) + rm.summary = {"samples": len(results), "succeeded": ok, "records_per_sample": len(table.prefix), + "skipped_loci_by_code": dict(skips), "genotype_notes": dict(reasons)} # fmt: skip + for o in outs: + if o.exists(): + rm.add_output(o, "vcf/4.3") + rm.status = "completed" + rm.write(opts.out / "run_manifest.json") + rep.info(f"done: {ok}/{len(results)} VCF(s) in {opts.out}") + if ok == len(results): + return ExitCode.OK + return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED + + +def write_audit(path: Path, audit: list[AuditRow]) -> None: + import pyarrow as pa + + from simurg.core.schemas import write_parquet + + cols = list(AuditRow.__dataclass_fields__) + table = pa.table({c: [getattr(a, c) for a in audit] for c in cols}) + path.parent.mkdir(parents=True, exist_ok=True) + write_parquet(table, path) + + +def write_merge_audit(path: Path, rows: list[tuple[str, str, tuple[int, ...], str]]) -> None: + import pyarrow as pa + + from simurg.core.schemas import write_parquet + + table = pa.table( + { + "sample_id": [r[0] for r in rows], + "site_id": [r[1] for r in rows], + "calls": ["/".join(str(c) for c in r[2]) for r in rows], + "reason": [r[3] for r in rows], + } + ) + path.parent.mkdir(parents=True, exist_ok=True) + write_parquet(table, path) + + +# ------------------------------------------------------------------------------------------------ cohort BCF +def _f32_payload(v: npt.NDArray[np.float32]) -> bytes: + u = np.asarray(v, dtype=" int: + """OUT-040: one multi-sample BCF (+ .csi), same records and FORMAT rules as the per-sample VCFs.""" + from simurg.formats.bcf import ( + FLOAT, + INT8, + BcfWriter, + CsiBuilder, + _type_byte, + encode_record, + header_dictionaries, + ) + from simurg.formats.bpm import read_bpm + from simurg.pipeline.qc import default_tmp, load_cohort + + n_loci = read_bpm(opts.bpm).num_loci + workers = max(1, opts.threads or (os.cpu_count() or 1)) + store, infos = load_cohort( + [(sid, path) for sid, _, path in samples], + n_loci, + default_tmp(out.parent), + baf=True, + scores=True, + lrr=True, + workers=workers, + ) + try: + for s in infos: + if s.error: + raise SimurgError("E_GTC_UNSUITABLE", f"{s.sample_id}: {s.error}") + names = [name for _, name, _ in samples] + cluster = read_cluster_name(samples[0][2]) + header = t.header.replace("__CLUSTER__", cluster) + header += "".join(f"##SAMPLE=\n" for sid, name, _ in samples) + header += "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" + "\t".join(names) + "\n" + strings, contigs = header_dictionaries(header) + k_gt, k_gs, k_baf, k_lrr = strings["GT"], strings["GS"], strings["BAF"], strings["LRR"] + f_pass, f_low = strings["PASS"], strings.get("LowGS", -1) + n = len(names) + sexes = [s.sex_est for s in infos] + ploidy = np.stack([ploidy_for(t.site_class, sx, t.ploidy_policy) for sx in sexes], axis=1) # sites × samples + single_pos = np.full(len(t.prefix), -1, np.int64) + single_pos[t.single_rows] = np.arange(len(t.single_rows)) + merged_pos = {int(k): j for j, k in enumerate(t.merged_rows.tolist())} + csi = CsiBuilder(len(contigs)) + fdesc = _type_byte(1, FLOAT) + written = 0 + with atomic_open(out, "wb") as raw: + w = BcfWriter(raw, header) + chunk = 20000 + for start in range(0, len(t.prefix), chunk): + rows = range(start, min(len(t.prefix), start + chunk)) + need = {int(t.single_locus[single_pos[k]]) for k in rows if single_pos[k] >= 0} + for k in rows: + if k in merged_pos: + need.update(p.index for p in t.merged_sites[merged_pos[k]].probes) + loci = sorted(need) + col = {loc: j for j, loc in enumerate(loci)} + sel = np.array(loci, np.int64) + gt = np.asarray(store.gt[:, sel]) if len(sel) else np.zeros((n, 0), np.uint8) + sc = np.asarray(store.score[:, sel]) if len(sel) else np.zeros((n, 0), np.float32) # type: ignore[index] + bf = np.asarray(store.baf[:, sel]) if len(sel) else np.zeros((n, 0), np.float32) # type: ignore[index] + lr = np.asarray(store.lrr[:, sel]) if len(sel) else np.zeros((n, 0), np.float32) # type: ignore[index] + for k in rows: + pl = ploidy[k] + if single_pos[k] >= 0: + j = int(single_pos[k]) + c = col[int(t.single_locus[j])] + calls = gt[:, c].astype(np.int64) + ia, ib = (int(x) for x in t.single_ab[j]) # type: ignore[index] + a1 = np.where(calls == 3, ib, ia) + a2 = np.where(calls == 1, ia, ib) + lo, hi = np.minimum(a1, a2), np.maximum(a1, a2) + width = max(1, int(pl.max())) + cells = np.full((n, width), -127, np.int8) + miss = (calls == 0) | (pl == 0) | ((pl == 1) & (calls == 2)) + hap = pl <= 1 + cells[:, 0] = np.where(miss, 0, ((np.where(hap & (calls == 3), ib, lo) + 1) << 1)).astype( + np.int8 + ) + if width == 2: + second = np.where(miss, 0, (hi + 1) << 1).astype(np.int8) + cells[:, 1] = np.where(hap, np.int8(-127), second) + gs = sc[:, c] + b = bf[:, c] + b = np.where(t.single_flip[j], (np.float32(1.0) - b).astype(np.float32), b) + fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs)), + (k_baf, fdesc, _f32_payload(b)), (k_lrr, fdesc, _f32_payload(lr[:, c]))] # fmt: skip + gs_min = gs + else: + site = t.merged_sites[merged_pos[k]] + cs = [col[p.index] for p in site.probes] + g_all = gt[:, cs].astype(np.int64) + gs_min = sc[:, cs].min(axis=1) + width = max(1, int(pl.max())) + cells = np.full((n, width), -127, np.int8) + for si in range(n): + p_k = int(pl[si]) + if p_k == 0: + cells[si, 0] = 0 + continue + o = combine_genotype(site, tuple(int(x) for x in g_all[si]), p_k, t.inf1_nc_informative) + if o.alleles is None: + cells[si, :p_k] = 0 + else: + for q, a in enumerate(o.alleles): + cells[si, q] = (a + 1) << 1 + fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs_min))] + low = opts.filter_gs is not None and bool(np.nanmin(gs_min) < opts.filter_gs) if n else False + rec = encode_record(contigs[t.site_contig[k]], int(t.site_beg[k]), + int(t.site_end[k] - t.site_beg[k]), t.site_id[k], t.site_alleles[k], + f_low if low and f_low >= 0 else f_pass, n, fmt) # fmt: skip + vs, ve = w.write(rec) + csi.add(contigs[t.site_contig[k]], int(t.site_beg[k]), int(t.site_end[k]), vs, ve) + written += 1 + w.close() + Path(str(out) + ".csi").write_bytes(csi.encode()) + finally: + store.close() + return written + + +def read_cluster_name(gtc_path: str) -> str: + from simurg.formats.gtc import read_gtc + + return read_gtc(gtc_path, check_write_complete=False).cluster_file or "" diff --git a/tests/oracle/test_p6_vcf.py b/tests/oracle/test_p6_vcf.py new file mode 100644 index 0000000..3acce8c --- /dev/null +++ b/tests/oracle/test_p6_vcf.py @@ -0,0 +1,150 @@ +"""Phase P6 acceptance on P0 data (plan §18 P6): VCFs validate with bcftools and index correctly; genotypes agree +with gtc2vcf apart from documented categories; BAF is ALT-oriented (EDGE-017); bedGraph equals GTC values (EDGE-023). + +Needs the GRCh38 FASTA (``scripts/fetch_data.sh p6``) and, for the comparison, the gtc2vcf oracle VCF +(``scripts/make_oracle_vcf.sh``). +""" + +from __future__ import annotations + +import gzip +import subprocess +from pathlib import Path + +import numpy as np +import pytest +from tests.conftest import GSA_BPM, GSA_CSV, ORACLE_GTCS, REFS, ROOT, need + +pytestmark = [pytest.mark.data, pytest.mark.slow] +FASTA = ROOT / "data" / "p6" / "GRCh38.primary.fa" +BCFTOOLS = REFS / "tools" / "bcftools-1.22" / "bcftools" +ORACLE_VCF = ROOT / "data" / "p6" / "oracle" / "gsa36_gtc2vcf.vcf.gz" + + +@pytest.fixture(scope="module") +def vcfs(tmp_path_factory: pytest.TempPathFactory) -> Path: + need(GSA_BPM, GSA_CSV, FASTA, ORACLE_GTCS) + from typer.testing import CliRunner + + from simurg.cli.main import app + + gtcs = tmp_path_factory.mktemp("gtc") + for f in sorted(ORACLE_GTCS.glob("*.gtc"))[:3]: + (gtcs / f.name).symlink_to(f) + out = tmp_path_factory.mktemp("vcf") + res = CliRunner().invoke( + app, + ["gtc-to-vcf", "--bpm", str(GSA_BPM), "--csv", str(GSA_CSV), "--fasta", str(FASTA), "--gtcs", str(gtcs), + "--out", str(out), "--threads", "1", "--unsquash-duplicates", "--quiet", "--force"], + ) # fmt: skip + assert res.exit_code == 0, res.output + return out + + +def test_bcftools_validates_and_index_works(vcfs: Path) -> None: # OUT-038, OUT-039, EDGE-020, EDGE-021 + if not BCFTOOLS.exists(): + pytest.skip("bcftools not built (scripts/make_oracle_gtcs.sh)") + f = next(vcfs.glob("*.snv.vcf.gz")) + n = int(subprocess.run([str(BCFTOOLS), "view", "-H", str(f)], capture_output=True, check=True).stdout.count(b"\n")) + nrec = subprocess.run([str(BCFTOOLS), "index", "--nrecords", str(f)], capture_output=True, text=True, check=True) + assert int(nrec.stdout.strip()) == n > 600_000 + for region in ("1:1000000-2000000", "X:150000000-156000000", "MT", "Y"): + idx = subprocess.run([str(BCFTOOLS), "view", "-H", "-r", region, str(f)], capture_output=True, check=True) + chrom, _, span = region.partition(":") + lo, _, hi = span.partition("-") + scan = 0 + with gzip.open(f, "rt") as fh: + for line in fh: + if line.startswith("#"): + continue + c, p = line.split("\t", 2)[:2] + if c == chrom and (not lo or int(lo) <= int(p) <= int(hi)): + scan += 1 + assert idx.stdout.count(b"\n") == scan + + +def test_matches_gtc2vcf(vcfs: Path, tmp_path: Path) -> None: # plan §18 P6 gate + if not ORACLE_VCF.exists() or not BCFTOOLS.exists(): + pytest.skip("gtc2vcf oracle VCF missing (scripts/make_oracle_vcf.sh)") + names = ",".join(p.name.removesuffix(".snv.vcf.gz") for p in sorted(vcfs.glob("*.snv.vcf.gz"))) + sub = tmp_path / "oracle.vcf.gz" + subprocess.run([str(BCFTOOLS), "view", "-s", names, "-Oz", "-o", str(sub), str(ORACLE_VCF)], check=True) + out = tmp_path / "cmp.json" + r = subprocess.run( + [str(ROOT / ".venv" / "bin" / "python"), str(ROOT / "scripts" / "validate_p6_vcf.py"), "--ours", str(vcfs), + "--oracle", str(sub), "--fasta", str(FASTA), "--bcftools", str(BCFTOOLS), "--out", str(out)], + capture_output=True, text=True, + ) # fmt: skip + assert r.returncode == 0, r.stdout + r.stderr + import json + + d = json.loads(out.read_text()) + assert d["unexplained"] == 0 and d["categories"]["match_called"] > 1_500_000 + + +def test_baf_is_alt_oriented(vcfs: Path) -> None: # OUT-016, EDGE-017 + from simurg.formats.bpm import read_bpm + from simurg.formats.gtc import read_gtc + + f = sorted(vcfs.glob("*.snv.vcf.gz"))[0] + g = read_gtc(ORACLE_GTCS / f.name.replace(".snv.vcf.gz", ".gtc")) + m = read_bpm(GSA_BPM) + idx = m.name_index() + flipped = same = 0 + with gzip.open(f, "rt") as fh: + for line in fh: + if line.startswith("#"): + continue + c = line.rstrip("\n").split("\t") + if "," in c[2] or ";" in c[2] or c[8] != "GT:GS:BAF:LRR" or "," in c[4]: + continue + baf = c[9].split(":")[2] + if baf == ".": + continue + i = idx[c[2]] + gtc_baf = float(g.b_allele_freqs[i]) + gt = c[9].split(":")[0] + call = int(g.genotypes[i]) + if gt in ("0/0", "1/1") and call in (1, 3): + a_is_ref = (gt == "0/0") == (call == 1) + want = gtc_baf if a_is_ref else float(np.float32(1.0) - np.float32(gtc_baf)) + assert abs(float(baf) - want) < 1e-6, (c[2], baf, gtc_baf) + flipped += not a_is_ref + same += a_is_ref + if flipped > 500 and same > 500: + break + assert flipped > 500 and same > 500 # both orientations exercised + + +def test_bedgraph_equals_gtc(tmp_path: Path) -> None: # OUT-050, OUT-052, EDGE-023 + need(GSA_BPM, ORACLE_GTCS) + from simurg.core.log import Reporter + from simurg.formats.bpm import read_bpm + from simurg.formats.gtc import read_gtc + from simurg.pipeline.export import BedgraphOptions, run_export_bedgraph + + gtcs = tmp_path / "g" + gtcs.mkdir() + src = sorted(ORACLE_GTCS.glob("*.gtc"))[0] + (gtcs / src.name).symlink_to(src) + out = tmp_path / "bg" + assert run_export_bedgraph(BedgraphOptions(out=out, bpm=GSA_BPM, gtcs=gtcs), Reporter(command="t", quiet=True)) == 0 + m = read_bpm(GSA_BPM) + g = read_gtc(src) + single = {} + seen: dict[tuple[str, int], int] = {} + for i, (c, p) in enumerate(zip(m.chroms, m.map_infos.tolist(), strict=True)): + key = ("X" if c == "XY" else c, p) + seen[key] = seen.get(key, 0) + 1 + single[key] = i + lines = gzip.decompress((out / f"{src.stem}.lrr.bedgraph.gz").read_bytes()).decode().splitlines() + assert lines[0].startswith('track type=bedGraph name="') and "viewLimits=-2:2" in lines[0] + checked = 0 + for ln in lines[1:]: + c, s, e, v = ln.split("\t") + key = (c, int(e)) + assert int(e) - int(s) == 1 + if seen[key] == 1: + assert v == f"{float(g.logr_ratios[single[key]]):.6g}" + checked += 1 + assert checked > 500_000 diff --git a/tests/unit/test_vcf.py b/tests/unit/test_vcf.py new file mode 100644 index 0000000..1951ec6 --- /dev/null +++ b/tests/unit/test_vcf.py @@ -0,0 +1,234 @@ +"""Phase P6 unit tests: FASTA access (IN-060 … IN-064), BGZF/tabix (OUT-038), plus-strand mapping (OUT-015, +OUT-020), genotype combination (OUT-031, EDGE-011/012), ploidy (OUT-018), BAF orientation (EDGE-017), PLINK codes.""" + +from __future__ import annotations + +import gzip +import struct +from pathlib import Path + +import numpy as np +import pytest + +from simurg.core.errors import SimurgError +from simurg.formats.bgzf import BgzfWriter, TabixBuilder, bgzf_compress, reg2bin +from simurg.formats.fasta import ContigMap, Genome, build_fai, build_gzi, detect_build, read_gzi +from simurg.harmonize.combine import AA, AB, BB, NC, failure_message, genotype, gt_text, merged +from simurg.harmonize.sites import ( + ManifestView, + ProbeSite, + Site, + group_sites, + left_shift, + map_locus, + normalize_indel, + revcomp, +) +from simurg.pipeline.export import bed_code +from simurg.pipeline.vcf import CLASS_AUTO, CLASS_MT, CLASS_PAR, CLASS_X, CLASS_Y, _site_class, ploidy_for + +SEQ = "ACGTACGTAAAACCCCGGGGTTTTACGTCAGCAGCAGTTT" * 3 # 120 bp + + +@pytest.fixture +def fasta(tmp_path: Path) -> Path: + p = tmp_path / "g.fa" + lines = [">1 test"] + [SEQ[i : i + 50] for i in range(0, len(SEQ), 50)] + [">X"] + [SEQ[:30].lower()] + p.write_text("\n".join(lines) + "\n") + return p + + +def test_fai_and_fetch_modes(fasta: Path, tmp_path: Path) -> None: # IN-060, IN-064, EDGE-010 + fai = build_fai(fasta, False) + assert [(e.name, e.length, e.line_bases, e.line_width) for e in fai] == [("1", 120, 50, 51), ("X", 30, 30, 31)] + assert fai[0].offset == len(">1 test\n") + g1 = Genome(fasta, use_cache=False, cache_dir=tmp_path / "c1") + g2 = Genome(fasta, use_cache=True, cache_dir=tmp_path / "c2") + for a, b in ((0, 120), (45, 55), (49, 51), (100, 200), (-5, 3)): + assert g1.fetch("1", a, b) == g2.fetch("1", a, b) == SEQ[max(0, a) : min(120, b)] + assert g1.fetch("X", 0, 5) == "ACGTA" # upper-cased + assert (tmp_path / "c1" / "index.fai").exists() # built into the cache, not next to the input + + +def test_bgzip_fasta(fasta: Path, tmp_path: Path) -> None: + bz = tmp_path / "g.fa.gz" + bz.write_bytes(bgzf_compress(fasta.read_bytes())) + assert read_gzi_roundtrip(bz) + g = Genome(bz, use_cache=False, cache_dir=tmp_path / "c") + assert g.bgzf and g.fetch("1", 40, 70) == SEQ[40:70] + + +def read_gzi_roundtrip(path: Path) -> bool: + entries = build_gzi(path) + gzi = path.with_suffix(".gzi") + gzi.write_bytes(struct.pack(" None: # IN-060 + Path(str(fasta) + ".fai").write_text("1\t120\t3\t50\t51\nX\t30\t140\t30\t31\n") + with pytest.raises(SimurgError) as e: + Genome(fasta, cache_dir=tmp_path / "c") + assert e.value.code == "E_FAI_MISMATCH" + + +def test_build_and_contig_alias(fasta: Path, tmp_path: Path) -> None: # IN-061, IN-063 + from simurg.formats.fasta import FaiEntry + + assert detect_build([FaiEntry("chr1", 248956422, 0, 60, 61)]) == "GRCh38" + assert detect_build([FaiEntry("1", 249250621, 0, 60, 61)]) == "GRCh37" + strict = ContigMap(["1", "X", "MT"], "strict") + assert strict.get("XY") == "X" and strict.get("chr1") is None + auto = ContigMap(["chr1", "chrX", "chrM", "chr22_KI270879v1_alt"], "auto") + assert auto.get("1") == "chr1" and auto.get("MT") == "chrM" and auto.get("XY") == "chrX" + assert auto.get("22_KI270879v1_alt") == "chr22_KI270879v1_alt" + alias = tmp_path / "a.tsv" + alias.write_text("1\tNC_000001\n") + assert ContigMap(["NC_000001"], str(alias)).get("1") == "NC_000001" + + +def test_bgzf_blocks_and_tabix(tmp_path: Path) -> None: # OUT-038 + out = tmp_path / "t.gz" + tbx = TabixBuilder() + with out.open("wb") as fh: + w = BgzfWriter(fh) + w.write_record(b"#header\n") + for i in range(5000): + line = f"1\t{i * 10 + 1}\trs{i}\tA\tC\n".encode() + s, e = w.write_record(line) + tbx.add("1", i * 10, i * 10 + 1, s, e) + w.close() + text = gzip.decompress(out.read_bytes()).decode() + assert text.count("\n") == 5001 + idx = gzip.decompress(tbx.encode()) + assert idx[:4] == b"TBI\x01" and struct.unpack_from(" ManifestView: + base = dict(names=["p"], snps=["[A/G]"], chroms=["1"], positions=[5], ilmn_strands=["TOP"], + source_strands=["TOP"], ref_strands=["+"], source_seqs=[""], assay_types=[0], + intensity_only=[False]) # fmt: skip + base.update({k: [v] for k, v in kw.items()}) + return ManifestView(**base) # type: ignore[arg-type] + + +def test_snv_plus_and_minus() -> None: # plan §8.1, OUT-015 + genome = "TTTTAGGGG" # position 5 = A + p, a = map_locus(0, mv(), fetcher(genome), lambda c: c) + assert p is not None and p.ref == "A" and p.alleles == ("A", "G") and not a.swapped + p, a = map_locus(0, mv(snps="[T/C]", ref_strands="-"), fetcher(genome), lambda c: c) + assert p is not None and p.alleles == ("A", "G") and a.flipped + p, a = map_locus(0, mv(snps="[C/G]"), fetcher(genome), lambda c: c) + assert a.reason == "W_REF_NOT_QUERIED" and p is not None + assert map_locus(0, mv(chroms="0"), fetcher(genome), lambda c: c)[1].reason == "W_SKIP_UNMAPPED" + assert map_locus(0, mv(intensity_only=True), fetcher(genome), lambda c: c)[1].reason == "W_SKIP_INTENSITY_ONLY" + assert map_locus(0, mv(chroms="Z"), fetcher(genome), lambda c: None)[1].reason == "W_SKIP_CONTIG" + + +def test_indel_deletion_and_insertion_normalised() -> None: # OUT-020, IN-046 + # 1234567890123 + genome = "GGCATCATCATTT" + # CAT copies at 3-5, 6-8, 9-11: deleting the last copy normalises to the G anchor at position 2 + pos, ref, alt = normalize_indel(fetcher(genome), "1", 9, "CAT", True) + assert (pos, ref, alt) == (2, "GCAT", "G") + pos, ref, alt = normalize_indel(fetcher(genome), "1", 11, "CAT", False) # insertion after position 11 + assert (pos, ref, alt) == (2, "G", "GCAT") + assert left_shift("GGCATCAT", "CAT", "TT") == ("GG", "CATCATTT") + assert revcomp("ACGTN") == "NACGT" + + +def test_indel_mapping_from_source_seq() -> None: + genome = "AAAAACCCCGATTACAGGGGGTTTTT" + src = "AAAAACCCCG[-/ATTACA]GGGGGTTTTT" # the reference carries ATTACA at 11-16 → deletion + view = mv(snps="[I/D]", positions=11, ilmn_strands="PLUS", source_strands="PLUS", source_seqs=src) + p, a = map_locus(0, view, fetcher(genome), lambda c: c) + assert p is not None and p.kind == "indel" + assert (p.pos, p.ref, p.alleles) == (10, "GATTACA", ("GATTACA", "G")) # I = REF, D = ALT + assert a.swapped is False + + +def test_grouping_and_order() -> None: # OUT-030, OUT-033, OUT-034, OUT-014 + def ps(i: int, name: str, alleles: tuple[str, str], pos: int = 5) -> ProbeSite: + return ProbeSite(i, name, "1", pos, "A", alleles, "snv", 2, "1") + + probes = [ps(0, "b", ("A", "G")), ps(1, "a", ("A", "C")), ps(2, "c", ("A", "G"), pos=3)] + sites = group_sites(probes, {"1": 0}) + assert [s.id for s in sites] == ["c", "b,a"] + assert sites[1].alts == ["G", "C"] # BPM order, then allele order + dup = [ps(0, "x", ("A", "G")), ps(1, "y", ("A", "G"))] + assert len(group_sites(dup, {"1": 0}, unsquash=True)) == 2 + assert len(group_sites(probes[:2], {"1": 0}, unsquash=True)) == 1 # multi-allelic stays merged + + +# ------------------------------------------------------------------------------------------------ combination +def site_of(*probes: ProbeSite) -> Site: + ref = probes[0].ref + alts: list[str] = [] + for p in probes: + alts += [x for x in p.alleles if x != ref and x not in alts] + return Site("1", 5, ref, alts, list(probes)) + + +def P(name: str, a: str, b: str, inf: int, ref: str = "A") -> ProbeSite: + return ProbeSite(0, name, "1", 5, ref, (a, b), "snv", inf, "1") + + +def test_edge_011_scenarios() -> None: # OUT-031 + # Infinium II [A/G] = AA with Infinium I [T/G] = TT → ambiguous between A/T and T/T + s = site_of(P("ii", "A", "G", 2), P("i", "T", "G", 1)) + assert merged(s, (AA, AA), 2).reason == "W_MERGE_AMBIGUOUS" + # Infinium I [T/A] homozygous T with Infinium II [A/G] no call → ambiguous (not a T/T call) + s = site_of(P("i", "T", "A", 1), P("ii", "A", "G", 2)) + assert merged(s, (AA, NC), 2).reason == "W_MERGE_AMBIGUOUS" + # inconsistent pair + s = site_of(P("p1", "A", "G", 1), P("p2", "A", "G", 1)) + assert merged(s, (AA, BB), 2).reason == "W_MERGE_INCONSISTENT" + assert merged(s, (NC, NC), 2).reason == "W_MERGE_ALL_NC" + # consistent duplicates → a call + o = merged(s, (AB, AB), 2) + assert o.alleles == (0, 1) and gt_text(o, 2) == "0/1" + # tri-allelic: A/G and A/C Infinium I probes + s = site_of(P("g", "A", "G", 1), P("c", "A", "C", 1)) + assert merged(s, (AB, AB), 2).reason == "W_MERGE_INCONSISTENT" + assert merged(s, (BB, AB), 2).reason == "W_MERGE_INCONSISTENT" + assert merged(s, (AB, AA), 2).alleles == (0, 1) # A/G + assert "InfiniumI" in failure_message(s, (AB, BB), "W_MERGE_INCONSISTENT") + + +def test_inf1_nocall_option() -> None: # OUT-037 + s = site_of(P("ii", "A", "G", 2), P("i", "T", "G", 1)) + assert merged(s, (AA, NC), 2).reason == "W_MERGE_AMBIGUOUS" + o = merged(s, (AA, NC), 2, inf1_nc_informative=True) + assert o.alleles == (0, 0) # the T/G probe saw nothing → no T and no G → A/A + + +def test_edge_012_opposite_strand_indels() -> None: + a = ProbeSite(0, "fwd", "1", 5, "ACT", ("ACT", "A"), "indel", 1, "1") # [I/D] + b = ProbeSite(1, "rev", "1", 5, "ACT", ("A", "ACT"), "indel", 1, "1") # [D/I] + s = Site("1", 5, "ACT", ["A"], [a, b], "indel") + assert gt_text(genotype(s, (AB, AB), 2), 2) == "0/1" + + +def test_ploidy_rules() -> None: # OUT-018, EDGE-038 + cls = np.array([CLASS_AUTO, CLASS_X, CLASS_Y, CLASS_MT, CLASS_PAR], np.int8) + assert ploidy_for(cls, "M", "sex-aware").tolist() == [2, 1, 1, 1, 2] + assert ploidy_for(cls, "F", "sex-aware").tolist() == [2, 2, 0, 1, 2] + assert ploidy_for(cls, "U", "sex-aware").tolist() == [2, 2, 1, 1, 2] + assert ploidy_for(cls, "M", "diploid").tolist() == [2, 2, 2, 1, 2] + assert _site_class("X", 100_000, False, "GRCh38") == CLASS_PAR + assert _site_class("X", 5_000_000, False, "GRCh38") == CLASS_X + assert _site_class("chrM", 5, False, "GRCh38") == CLASS_MT + s = site_of(P("p", "A", "G", 2)) + assert gt_text(genotype(s, (AB,), 1), 1) == "." # het in haploid context + assert gt_text(genotype(s, (BB,), 1), 1) == "1" + + +def test_plink_codes() -> None: # OUT-043 + assert [bed_code(g, 1) for g in ("0/0", "0/1", "1/1", "./.", "1", "0")] == [3, 2, 0, 1, 0, 3] + assert bed_code("0/2", 1) == 3 and bed_code("0/2", 2) == 2 # split multi-allelic record