Madhava/libs

.cargo/config.toml

@@ -0,0 +1,2 @@
+[net]
+git-fetch-with-cli = true

.gitignore

@@ -81,3 +81,8 @@ examples/**/.nextflow*
 # Read-only reference material (htslib/rust-htslib source trees used during
 # the CRAM/noodles migration; not part of the build).
 context/
+/repos/
+
+# Local-only cargo env overrides (macOS Xcode CC/CXX/SDKROOT etc.).
+# Belongs per-developer; recreate locally if your toolchain needs it.
+rust/.cargo/config.toml

.gitmodules

@@ -5,4 +5,32 @@
 [submodule "noodles"]
 	path = noodles
 	url = git@github.com:madhavajay/noodles.git
-	branch = madhava/streaming-slice-records
+	branch = madhava/bioscript
+[submodule "vendor/python/pysam"]
+	path = vendor/python/pysam
+	url = https://github.com/pysam-developers/pysam.git
+[submodule "vendor/python/pyfaidx"]
+	path = vendor/python/pyfaidx
+	url = https://github.com/mdshw5/pyfaidx.git
+[submodule "ports/vntyper/vntyper"]
+	path = ports/vntyper/vntyper
+	url = https://github.com/madhavajay/VNtyper.git
+[submodule "ports/vntyper/kestrel"]
+	path = ports/vntyper/kestrel
+	url = https://github.com/paudano/kestrel.git
+[submodule "ports/vntyper/kescases"]
+	path = ports/vntyper/kescases
+	url = https://github.com/paudano/kescases.git
+[submodule "vendor/rust/kestrel-rs"]
+	path = vendor/rust/kestrel-rs
+	url = git@github.com:madhavajay/kestrel-rs.git
+[submodule "vendor/rust/htslib-rs"]
+	path = vendor/rust/htslib-rs
+	url = git@github.com:madhavajay/htslib-rs.git
+	branch = main
+[submodule "vendor/rust/bcftools-rs"]
+	path = vendor/rust/bcftools-rs
+	url = git@github.com:madhavajay/bcftools-rs.git
+[submodule "vendor/rust/samtools-rs"]
+	path = vendor/rust/samtools-rs
+	url = git@github.com:madhavajay/samtools-rs.git

.repoverse.yaml

@@ -0,0 +1,34 @@
+version: 1
+defaults:
+  remote: github
+  revision: main
+  scheme: ssh
+remotes:
+  github:
+    host: github.com
+projects:
+- name: madhavajay/htslib-rs
+  path: repos/htslib-rs
+  revision: main
+- name: madhavajay/noodles
+  path: repos/noodles
+  revision: madhava/bioscript
+provides:
+- madhavajay/htslib-rs
+- madhavajay/noodles
+links:
+- repo: madhavajay/htslib-rs
+  at: vendor/rust/bcftools-rs/repos/htslib-rs
+  branch: main
+- repo: madhavajay/htslib-rs
+  at: vendor/rust/htslib-rs
+  branch: main
+- repo: madhavajay/htslib-rs
+  at: vendor/rust/samtools-rs/repos/htslib-rs
+  branch: main
+- repo: madhavajay/noodles
+  at: noodles
+  branch: madhava/bioscript
+- repo: madhavajay/noodles
+  at: repos/htslib-rs/repos/noodles
+  branch: madhava/bioscript

bioscripts/apol1-new.py

@@ -0,0 +1,97 @@
+from bioscript import pysam
+
+
+G1_SITE_1 = bioscript.variant(
+    rsid="rs73885319",
+    grch37="22:36661906-36661906",
+    grch38="22:36265860-36265860",
+    ref="A",
+    alt="G",
+    kind="snp",
+)
+
+G1_SITE_2 = bioscript.variant(
+    rsid="rs60910145",
+    grch37="22:36662034-36662034",
+    grch38="22:36265988-36265988",
+    ref="T",
+    alt="G",
+    kind="snp",
+)
+
+G2_SITE = bioscript.variant(
+    rsid=["rs71785313", "rs1317778148", "rs143830837"],
+    grch37="22:36662046-36662051",
+    grch38="22:36266000-36266005",
+    ref="I",
+    alt="D",
+    kind="deletion",
+    deletion_length=6,
+    motifs=["TTATAA", "ATAATT"],
+)
+
+
+def count_char(text, needle):
+    if text is None:
+        return 0
+    total = 0
+    for ch in text:
+        if ch == needle:
+            total = total + 1
+    return total
+
+
+def count_non_ref(text, ref):
+    if text is None:
+        return 0
+    total = 0
+    for ch in text:
+        if ch != ref and ch != "-":
+            total = total + 1
+    return total
+
+
+def classify_apol1(genotypes):
+    site1 = genotypes.lookup_variant(G1_SITE_1)
+    site2 = genotypes.lookup_variant(G1_SITE_2)
+    g2 = genotypes.lookup_variant(G2_SITE)
+
+    if site1 is None and site2 is None and g2 is None:
+        return "G-/G-"
+
+    d_count = count_char(g2, "D")
+    site1_variants = count_non_ref(site1, "A")
+    site2_variants = count_non_ref(site2, "T")
+
+    has_g1 = site1_variants > 0 and site2_variants > 0
+    if has_g1:
+        g1_total = site1_variants + site2_variants
+    else:
+        g1_total = 0
+
+    if d_count == 2:
+        return "G2/G2"
+    if d_count == 1:
+        if g1_total >= 2:
+            return "G2/G1"
+        return "G2/G0"
+    if g1_total == 4:
+        return "G1/G1"
+    if g1_total >= 2:
+        return "G1/G0"
+    return "G0/G0"
+
+
+def main():
+    genotypes = bioscript.load_genotypes(input_file)
+    status = classify_apol1(genotypes)
+    rows = [{
+        "participant_id": participant_id,
+        "apol1_status": status,
+    }]
+    bioscript.write_tsv(output_file, rows)
+    print(status)
+
+
+if __name__ == "__main__":
+    main()

bioscripts/apol1-pysam-proof.py

@@ -0,0 +1,64 @@
+from bioscript import pysam
+
+
+APOL1_SITES = [
+    {
+        "key": "G1_SITE_1",
+        "chrom": "22",
+        "start": 36265859,
+        "stop": 36265860,
+        "ref": "A",
+        "alt": "G",
+    },
+    {
+        "key": "G1_SITE_2",
+        "chrom": "22",
+        "start": 36265987,
+        "stop": 36265988,
+        "ref": "T",
+        "alt": "G",
+    },
+    {
+        "key": "G2_SITE",
+        "chrom": "22",
+        "start": 36265999,
+        "stop": 36266005,
+        "ref": "TTATAA",
+        "alt": "<DEL:6>",
+    },
+]
+
+
+def count_region_reads(bam, site):
+    total = 0
+    for read in bam.fetch(site["chrom"], site["start"], site["stop"]):
+        if not read.is_unmapped:
+            total = total + 1
+    return total
+
+
+def main():
+    bam = pysam.AlignmentFile(
+        input_file,
+        "rc",
+        reference_filename=reference_file,
+        index_filename=input_index,
+    )
+    rows = []
+    for site in APOL1_SITES:
+        rows.append(
+            {
+                "participant_id": participant_id,
+                "variant_key": site["key"],
+                "chrom": site["chrom"],
+                "start": str(site["start"]),
+                "stop": str(site["stop"]),
+                "depth": str(count_region_reads(bam, site)),
+                "proof_status": "region_fetch_only",
+            }
+        )
+    bioscript.write_tsv(output_file, rows)
+
+
+if __name__ == "__main__":
+    main()

bioscripts/examples/vntyper/assay.yaml

@@ -0,0 +1,84 @@
+schema: bioscript:assay:1.0
+version: "1.0"
+name: vntyper_muc1
+label: MUC1 VNTR (VNtyper)
+summary: >
+  Advanced assay that genotypes the MUC1 VNTR frameshift (ADTKD-MUC1) from an
+  aligned genome. It slices the MUC1 region, extracts reads, runs mapping-free
+  Kestrel genotyping, and applies VNtyper post-processing to call the
+  pathogenic cytosine insertion. Requires aligned input (BAM/CRAM) — it cannot
+  run from SNP-chip or VCF genotypes.
+tags:
+  - type:risk
+  - gene:MUC1
+  - kind:vntr
+support:
+  input:
+    - bam
+    - cram
+members:
+  - kind: variant
+    path: muc1-vntr.yaml
+    version: "1.0"
+analyses:
+  - id: vntyper_muc1
+    kind: bioscript
+    path: vntyper.py
+    output_format: tsv
+    label: MUC1 VNTR genotype
+    derived_from:
+      - muc1-vntr.yaml
+    assets:
+      - id: muc1_reference
+        path: assets/muc1_motifs.fa
+    emits:
+      - key: vntyper_outcome
+        label: MUC1 VNTR outcome
+        value_type: string
+        format: badge
+      - key: vntyper_status
+        label: MUC1 VNTR status
+        value_type: string
+      - key: vntyper_confidence
+        label: Kestrel confidence
+        value_type: string
+        format: badge
+      - key: vntyper_variant
+        label: Called variant
+        value_type: string
+      - key: vntyper_alt_depth
+        label: Alternate-variant depth
+        value_type: string
+    logic:
+      source:
+        name: VNtyper / Kestrel
+        url: https://github.com/hassansaei/VNtyper
+      description: >
+        The MUC1 region is sliced from the aligned input, converted to FASTQ,
+        and genotyped with mapping-free Kestrel against the MUC1 motif
+        reference. VNtyper motif/frameshift/confidence post-processing selects
+        the called variant; vntyper_status is positive when a valid
+        High_Precision* / High_Precision frameshift passes the filters.
+findings:
+  - schema: bioscript:pgx-label:1.0
+    id: muc1_vntr_positive_finding
+    label: MUC1 VNTR pathogenic frameshift detected
+    authority_type: clinical_annotation
+    binding:
+      source: analysis
+      analysis_id: vntyper_muc1
+      key: vntyper_status
+      operator: equals
+      value: positive
+    regulatory_sources:
+      - "ClinVar"
+    pgx_action_level: "Informative"
+    evidence:
+      source: VNtyper
+      kind: method_annotation
+      id: ADTKD-MUC1
+      url: https://github.com/hassansaei/VNtyper
+    notes: >
+      A positive call indicates a MUC1 VNTR frameshift consistent with
+      ADTKD-MUC1; confirm with an orthogonal method (e.g. SNaPshot or
+      long-read sequencing).