diff --git a/.github/workflows/build_deploy_main.yaml b/.github/workflows/build_deploy_main.yaml index be8a4fb9..6f1ab59a 100644 --- a/.github/workflows/build_deploy_main.yaml +++ b/.github/workflows/build_deploy_main.yaml @@ -3,6 +3,7 @@ on: push: branches: - main + - dev permissions: id-token: write @@ -18,7 +19,7 @@ jobs: runs-on: ubuntu-latest environment: production outputs: - matrix: ${{ steps.get_next_version.outputs.next_version }} + next_version: ${{ steps.get_next_version.outputs.next_version }} steps: - name: 'Checkout repo' uses: actions/checkout@v4 @@ -55,13 +56,12 @@ jobs: environment: production needs: - get_next_version - if: ${{ needs.get_next_version.outputs.matrix != '{"include":[]}' && needs.get_next_version.outputs.matrix != '' }} - strategy: - matrix: ${{ fromJson(needs.get_next_version.outputs.matrix) }} + if: ${{ needs.get_next_version.outputs.next_version != '' && github.event.pull_request.draft == false }} env: DOCKER_BUILDKIT: 1 BUILDKIT_PROGRESS: plain CLOUDSDK_CORE_DISABLE_PROMPTS: 1 + TAG: ${{ needs.get_next_version.outputs.next_version }}${{ github.ref == 'refs/heads/dev' && '_dev' || '' }} steps: - uses: actions/checkout@v4 with: @@ -84,9 +84,8 @@ jobs: - name: 'build image' run: | docker build . -f src/talos/cpg_internal_scripts/CPG_Dockerfile \ - --build-arg VERSION=${{ matrix.tag }} \ - --tag ${{ env.IMAGES_PREFIX }}/${{ matrix.name }}:${{ matrix.tag }} + --tag ${{ env.IMAGES_PREFIX }}/talos:${{ env.TAG }} - name: 'push image' run: | - docker push "${{ env.IMAGES_PREFIX }}/${{ matrix.name }}:${{ matrix.tag }}" + docker push "${{ env.IMAGES_PREFIX }}/talos:${{ env.TAG }}" diff --git a/.github/workflows/dev_auto_build.yaml b/.github/workflows/dev_auto_build.yaml index a8ff6619..d5e5f1ea 100644 --- a/.github/workflows/dev_auto_build.yaml +++ b/.github/workflows/dev_auto_build.yaml @@ -24,7 +24,7 @@ jobs: runs-on: ubuntu-latest environment: production outputs: - matrix: ${{ steps.get_next_version.outputs.next_version }} + next_version: ${{ steps.get_next_version.outputs.next_version }} steps: - name: 'Checkout repo' uses: actions/checkout@v4 @@ -61,9 +61,7 @@ jobs: environment: production needs: - get_next_version - if: ${{ needs.get_next_version.outputs.matrix != '{"include":[]}' && needs.get_next_version.outputs.matrix != '' && github.event.pull_request.draft == false }} - strategy: - matrix: ${{ fromJson(needs.get_next_version.outputs.matrix) }} + if: ${{ needs.get_next_version.outputs.next_version != '' && github.event.pull_request.draft == false }} steps: - name: 'checkout repo' uses: actions/checkout@v4 @@ -88,12 +86,11 @@ jobs: - name: 'build image with get_next_version' run: | docker build . -f src/talos/cpg_internal_scripts/CPG_Dockerfile \ - --build-arg VERSION=${{ matrix.tag }} \ - --tag ${{ env.IMAGES_PREFIX }}/${{ matrix.name }}:${{ matrix.tag }} + --tag ${{ env.IMAGES_PREFIX }}/talos:${{ needs.get_next_version.outputs.next_version }} - name: 'Push to dev artifactory' run: | - docker push "${{ env.IMAGES_PREFIX }}/${{ matrix.name }}:${{ matrix.tag }}" | tee push.log + docker push "${{ env.IMAGES_PREFIX }}/talos:${{ needs.get_next_version.outputs.next_version }}" | tee push.log digest=$(grep 'digest:' push.log | grep -o 'sha256:[0-9a-f]*') echo ":package: Link to image:" >> $GITHUB_STEP_SUMMARY - echo "https://console.cloud.google.com/artifacts/docker/cpg-common/australia-southeast1/images-dev/${{ matrix.name }}/$digest?project=cpg-common" >> $GITHUB_STEP_SUMMARY + echo "https://console.cloud.google.com/artifacts/docker/cpg-common/australia-southeast1/images-dev/talos/$digest?project=cpg-common" >> $GITHUB_STEP_SUMMARY diff --git a/.github/workflows/get_version.py b/.github/workflows/get_version.py index b46b5341..e880bbd3 100644 --- a/.github/workflows/get_version.py +++ b/.github/workflows/get_version.py @@ -1,5 +1,4 @@ import json -import logging import os import re import subprocess @@ -9,11 +8,11 @@ def extract_version_from_file(file_path: str) -> str | None: """ Extract the version from a Dockerfile by searching for a line like: - ARG VERSION=1.0.0 + ENV VERSION=1.0.0 """ with open(file_path) as f: content = f.read() - pattern = re.compile(r'^\s*ARG\s+VERSION\s*=\s*([^\s]+)', re.MULTILINE) + pattern = re.compile(r'^\s*ENV\s+VERSION\s*=\s*([^\s]+)', re.MULTILINE) match = pattern.search(content) return match.group(1) if match else None @@ -27,48 +26,31 @@ def get_next_version_tag(folder: str, version: str) -> str: 'GCP_BASE_IMAGE', 'australia-southeast1-docker.pkg.dev/cpg-common/images', ) - base_image_path_archive = os.environ.get( - 'GCP_BASE_ARCHIVE_IMAGE', - 'australia-southeast1-docker.pkg.dev/cpg-common/images-archive', - ) - full_image_name_prod = f'{base_image_path_prod}/{folder}' - full_image_name_archive = f'{base_image_path_archive}/{folder}' - - tags_list = [] - - logging.basicConfig(level=logging.ERROR) - - for full_image_name in [full_image_name_prod, full_image_name_archive]: - cmd = [ - 'gcloud', - 'container', - 'images', - 'list-tags', - full_image_name, - '--format=json', - ] - try: - result = subprocess.run(cmd, capture_output=True, text=True, check=True) # noqa: S603 - except subprocess.CalledProcessError: - logging.error(f'Failed to list tags for {full_image_name}') - continue - - # If existing tags are found, proceed to determine the next version. - tags_list += json.loads(result.stdout) + cmd = [ + 'gcloud', + 'container', + 'images', + 'list-tags', + f'{base_image_path_prod}/{folder}', + '--format=json', + ] max_suffix = 0 pattern = re.compile(rf'^{re.escape(version)}-(\d+)$') - # If no tags are found, it returns the next version as -1. - for entry in tags_list: - tags = entry.get('tags', []) - for tag in tags: - match = pattern.match(tag) + try: + result = subprocess.run(cmd, capture_output=True, text=True, check=True) # noqa: S603 + for block in json.loads(result.stdout): + if not block.get('tags', []): + continue + match = pattern.match(block['tags'][0]) if match: num = int(match.group(1)) max_suffix = max(max_suffix, num) - new_suffix = max_suffix + 1 - return f'{version}-{new_suffix}' + except subprocess.CalledProcessError as err: + raise RuntimeError('Failed to list tags for the given image') from err + + return f'{version}-{max_suffix + 1}' def main(): @@ -77,17 +59,13 @@ def main(): current_version = extract_version_from_file(dockerfile_name) if current_version is None: # Throw an error here - raise NotImplementedError('The Dockerfile needs to contain a version string in the format "ARG VERSION=x.x.x"') + raise NotImplementedError('The Dockerfile needs to contain a version string in the format "ENV VERSION=x.x.x"') # Determine the next available tag based on current_version. new_tag = get_next_version_tag(container_name, current_version) - include_entries = [{'name': container_name, 'tag': new_tag}] - - # Build the final matrix structure. - matrix = {'include': include_entries} - print(json.dumps(matrix, separators=(',', ':')), file=sys.stderr) - print(json.dumps(matrix, separators=(',', ':')), end='') + print(new_tag, file=sys.stderr) + print(new_tag, end='') main() diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 9dfc4ddd..f5c3b2af 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -28,7 +28,7 @@ repos: # Static type analysis (as much as it's possible in python using type hints) - repo: https://github.com/pre-commit/mirrors-mypy - rev: v1.20.2 + rev: v2.1.0 hooks: - id: mypy args: [--pretty, --show-error-codes, --install-types, --non-interactive] diff --git a/CHANGELOG.md b/CHANGELOG.md index 73c78679..77821964 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 +[11.1.0] - 2026-07 + +### Removed + + * All references to the private build of BCFtools. the `--greedy` behaviour is now released as standard in BCFtools 1.24 + +### Added + + * STR functionality! This is currently specific to STRipy data, with a script to translate STRipy JSON data into a pseudo-joint called VCF + * If there is appetite for STR functionality, this could be expanded to other STR callers. + + + [11.0.0] - 2026-05-15 ### Added @@ -48,8 +61,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 > NOTE! For existing Talos users, obtaining this extra annotation will require the re-run of the annotation workflow. One example gene which is known to be obscured by this default behaviour was RNU2-2, but there may be others. - - [10.0.1] - 2026-03-23 ### Fixed diff --git a/Dockerfile b/Dockerfile index e46962c3..51641f03 100644 --- a/Dockerfile +++ b/Dockerfile @@ -19,49 +19,35 @@ RUN apt-get update && apt-get install -y --no-install-recommends \ FROM base AS bcftools_compiler -ARG BCFTOOLS_VERSION=1.23.1 - -# AS OF 11.0.0, Talos is building BCFtools from a private fork. This fork contains a single change - csq applies annotations -# to both coding and non-coding genes in the event of overlapping genes. By default BCFtools skips non-coding gene annotation -# if a coding transcript consequence was detected, but in practice this is masking clinically relevant non-coding gene variation -# in cases where the non-coding gene overlaps with a non-clinically relevant coding gene. -# The change made is to always search for non-coding variation, but to mask non-coding consequences in otherwise coding transcripts -# if a coding change was already detected. +ARG BCFTOOLS_VERSION=1.24 + RUN apt-get update && apt-get install --no-install-recommends -y \ - autoconf \ gcc \ - git \ libbz2-dev \ libcurl4-openssl-dev \ liblzma-dev \ libssl-dev \ make \ zlib1g-dev && \ - wget https://github.com/samtools/htslib/releases/download/${BCFTOOLS_VERSION}/htslib-${BCFTOOLS_VERSION}.tar.bz2 && \ - tar -xf htslib-${BCFTOOLS_VERSION}.tar.bz2 && \ - cd htslib-${BCFTOOLS_VERSION} && \ - ./configure --enable-libcurl && \ - make && \ - make DESTDIR=/bcftools_install install && \ - cd .. && \ - git clone https://github.com/populationgenomics/bcftools.git && \ - cd bcftools && \ - autoheader && \ - autoconf && \ + rm -rf /var/lib/apt/lists/* && \ + wget https://github.com/samtools/bcftools/releases/download/${BCFTOOLS_VERSION}/bcftools-${BCFTOOLS_VERSION}.tar.bz2 && \ + tar -xf bcftools-${BCFTOOLS_VERSION}.tar.bz2 && \ + cd bcftools-${BCFTOOLS_VERSION} && \ ./configure --enable-libcurl --enable-s3 --enable-gcs && \ make && \ strip bcftools plugins/*.so && \ + make DESTDIR=/bcftools_install install && \ + cd htslib-${BCFTOOLS_VERSION} && \ + make && \ make DESTDIR=/bcftools_install install FROM base AS talos COPY --from=bcftools_compiler /bcftools_install/usr/local/bin/* /usr/local/bin/ COPY --from=bcftools_compiler /bcftools_install/usr/local/libexec/bcftools/* /usr/local/libexec/bcftools/ -COPY --from=bcftools_compiler /bcftools_install/usr/local/lib/ /usr/local/lib/ -RUN ldconfig ARG ECHTVAR_VERSION=v0.2.2 -ARG VERSION=11.0.3 +ENV VERSION=11.1.0 RUN wget -q -O /bin/echtvar "https://github.com/brentp/echtvar/releases/download/${ECHTVAR_VERSION}/echtvar" && \ chmod +x /bin/echtvar diff --git a/README.md b/README.md index 50076a0d..d25b6568 100644 --- a/README.md +++ b/README.md @@ -69,7 +69,7 @@ There are two primary workflows: To build the Docker image: ``` -docker build -t talos:11.0.3 . +docker build -t talos:11.1.0 . ``` ### **2. Download Annotation Resources** diff --git a/docs/getting-started.md b/docs/getting-started.md index 7e29b012..9a57da10 100644 --- a/docs/getting-started.md +++ b/docs/getting-started.md @@ -25,7 +25,7 @@ You will need: Build the Talos Docker image locally: ```bash -docker build -t talos:11.0.3 . +docker build -t talos:11.1.0 . ``` --- @@ -81,7 +81,10 @@ The optional columns (history, ext_ids, seqr_map, mito) can be omitted completel | `seqr_map` | optional | ID mapping for Seqr hyperlinks in the HTML report. | | `mito` | optional | Path to a joint-called mitochondrial VCF. | -### Input types +STR can be handled by Talos, but only data called by STRipy, and aggregated into a joint VCF format by the script `talos/scripts/stripy_json_to_vcf.py`. +This is not yet exposed in the nextflow implementation, but may be in future. + +### Small-variant VCF input types diff --git a/nextflow.config b/nextflow.config index 434e4a32..26faa58e 100644 --- a/nextflow.config +++ b/nextflow.config @@ -16,8 +16,8 @@ params { large_files = "large_files" ref_genome = "${params.large_files}/ref.fa" - // Docker container - "docker build -t talos:11.0.3 ." - container = 'talos:11.0.3' + // Docker container - "docker build -t talos:11.1.0 ." + container = 'talos:11.1.0' // MANE transcript resource mane = "${params.large_files}/MANE.GRCh38.v1.5.summary.txt.gz" @@ -31,6 +31,7 @@ params { // Result files from the CreateRoiFromGff3 nextflow module ensembl_bed = "${params.processed_annotations}/GRCh38.bed" ensembl_merged_bed = "${params.processed_annotations}/GRCh38_merged.bed" + ensembl_symbol_lookup = "${params.processed_annotations}/GRCh38_symbol_to_ensg.json" // AlphaMissense raw data, and reformatted as an echtvar-compatible zip alphamissense_tsv = "${params.large_files}/AlphaMissense_hg38.tsv.gz" diff --git a/nextflow/inputs/config.toml b/nextflow/inputs/config.toml index f01bb6c4..e2554b4a 100644 --- a/nextflow/inputs/config.toml +++ b/nextflow/inputs/config.toml @@ -10,7 +10,7 @@ forced_panels = [ 144, 239] # 3 = Green only # 2 = Amber or Green only # 1 = Red, Amber, or Green -confidence_level = 1 +confidence_level = 3 [[GeneratePanelData.manual_overrides]] # this section permits the manual addition of genes to the panel data diff --git a/nextflow/modules/annotation/AnnotateCsqWithBcftools/main.nf b/nextflow/modules/annotation/AnnotateCsqWithBcftools/main.nf index 25712ed2..3c725b5a 100644 --- a/nextflow/modules/annotation/AnnotateCsqWithBcftools/main.nf +++ b/nextflow/modules/annotation/AnnotateCsqWithBcftools/main.nf @@ -14,6 +14,7 @@ process AnnotateCsqWithBcftools { set -euo pipefail bcftools csq --force -f "${reference}" \ + --greedy 1 \ --local-csq \ -g ${gff3} \ --unify-chr-names 'chr,-,chr' \ diff --git a/nextflow/modules/prep/CreateRoiFromGff3/main.nf b/nextflow/modules/prep/CreateRoiFromGff3/main.nf index 3465b69b..78970be2 100644 --- a/nextflow/modules/prep/CreateRoiFromGff3/main.nf +++ b/nextflow/modules/prep/CreateRoiFromGff3/main.nf @@ -8,9 +8,11 @@ process CreateRoiFromGff3 { // generates a BED file one-gene-per-line, for per-gene annotations // generates an BED with overlapping intervals merged for filtering // sorts both + // also generates a JSON file mapping symbols to ENSG IDs output: path "GRCh38.bed", emit: bed path "GRCh38_merged.bed", emit: merged_bed + path "GRCh38_symbol_to_ensg.json", emit: json_lookup script: """ @@ -19,7 +21,8 @@ process CreateRoiFromGff3 { python -m talos.annotation_scripts.create_roi_from_gff3 \ --gff3 ${gff} \ --unmerged_output unsorted_GRCh38.bed \ - --merged_output unsorted_merged_GRCh38.bed + --merged_output unsorted_merged_GRCh38.bed \ + --json_output GRCh38_symbol_to_ensg.json sort -k1,1V -k2,2n unsorted_GRCh38.bed > GRCh38.bed sort -k1,1V -k2,2n unsorted_merged_GRCh38.bed > GRCh38_merged.bed diff --git a/nextflow_schema.json b/nextflow_schema.json index 2cf64099..098d47c4 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -33,7 +33,7 @@ }, "container": { "type": "string", - "default": "talos:11.0.3", + "default": "talos:11.1.0", "description": "Docker container used by processes." }, "mane": { diff --git a/preparation.nf b/preparation.nf index e4a67977..4f6f78a3 100644 --- a/preparation.nf +++ b/preparation.nf @@ -121,13 +121,15 @@ workflow { // generate the Region-of-interest BED file from Ensembl GFF3 // generates a per-gene BED file with ID annotations // and a overlap-merged version of the same for more efficient region filtering - if (!file(params.ensembl_bed).exists() || !file(params.ensembl_merged_bed).exists()) { + if (!file(params.ensembl_bed).exists() || !file(params.ensembl_merged_bed).exists() || !file(params.ensembl_symbol_lookup).exists()) { CreateRoiFromGff3(ch_gff) ch_bed = CreateRoiFromGff3.out.bed ch_merged_bed = CreateRoiFromGff3.out.merged_bed + ch_symbol_lookup = CreateRoiFromGff3.out.json_lookup } else { ch_merged_bed = channel.fromPath(params.ensembl_merged_bed, checkIfExists: true) ch_bed = channel.fromPath(params.ensembl_bed, checkIfExists: true) + ch_symbol_lookup = channel.fromPath(params.ensembl_symbol_lookup, checkIfExists: true) } // pull and parse the MANE data into a Hail Table @@ -156,6 +158,7 @@ workflow { alphamissense = ch_alphamissense_zip bed = ch_bed merged_bed = ch_merged_bed + symbol_lookup = ch_symbol_lookup clinvar_all = ch_clinvar_all clinvar_pm5 = ch_clinvar_pm5 mitimpact = ch_mitimpact_zip @@ -172,6 +175,8 @@ output { } merged_bed { } + symbol_lookup { + } clinvar_all { } clinvar_pm5 { diff --git a/pyproject.toml b/pyproject.toml index a469abea..b48172d2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = "hatchling.build" name='talos' description='Centre for Population Genomics Variant Prioritisation' readme = "README.md" -version="11.0.3" +version="11.1.0" requires-python = ">=3.10,<3.12" license-files = ["LICENSE"] classifiers=[ @@ -23,7 +23,7 @@ dependencies=[ 'hail~=0.2.137', 'clinvarbitration~=2.2.7', 'cloudpathlib[all]>=0.16.0', - 'cpg-utils>=5.4.1', + 'cpg-utils~=5.7', 'cyvcf2>=0.30.18', 'grpcio-status>=1.48,<1.50', 'httpx>=0.27.0', @@ -57,6 +57,9 @@ docs = [ 'mkdocs-material>=9.5', 'mkdocs-include-markdown-plugin>=6.2', ] +stripy = [ + 'pysam' +] [project.urls] Repository = "https://github.com/populationgenomics/talos" @@ -102,7 +105,8 @@ ignore = [ "PLR0913", # Too many arguments in function (> 5) "C901", # method is too complex (> 10 conditions) "S311", # Standard pseudo-random generators are not suitable for cryptographic purposes - "PLR2004" # Magic value used in comparison, consider replacing `X` with a constant variable. I don't rate this + "PLR2004", # Magic value used in comparison, consider replacing `X` with a constant variable. I don't rate this + "PLW0108", # Lambda may be unnecessary; consider inlining inner function ] [tool.ruff.format] @@ -127,7 +131,7 @@ hail = ["hail", "hailtop"] "test/test_cpg_flow_utils.py" = ['ANN202', 'ARG001'] [tool.bumpversion] -current_version = "11.0.3" +current_version = "11.1.0" parse = "(?P\\d+)\\.(?P\\d+)\\.(?P\\d+)" serialize = ["{major}.{minor}.{patch}"] commit = true diff --git a/src/talos/annotation_scripts/annotated_vcf_into_matrixtable.py b/src/talos/annotation_scripts/annotated_vcf_into_matrixtable.py index 037a9ef5..fca46bfa 100755 --- a/src/talos/annotation_scripts/annotated_vcf_into_matrixtable.py +++ b/src/talos/annotation_scripts/annotated_vcf_into_matrixtable.py @@ -191,7 +191,11 @@ def annotate_all_transcript_consequences( mane[x.transcript]['mane_id'], MISSING_STRING, ), - gene_id=ensgs[mt.locus.contig].get(x.gene, x.gene), + gene_id=hl.if_else( + ensgs.contains(mt.locus.contig), + ensgs[mt.locus.contig].get(x.gene, x.gene), + x.gene, + ), ), mt.transcript_consequences, ), diff --git a/src/talos/annotation_scripts/create_roi_from_gff3.py b/src/talos/annotation_scripts/create_roi_from_gff3.py index 69a6d6ea..84115786 100755 --- a/src/talos/annotation_scripts/create_roi_from_gff3.py +++ b/src/talos/annotation_scripts/create_roi_from_gff3.py @@ -9,6 +9,7 @@ """ import gzip +import json import re from argparse import ArgumentParser @@ -31,7 +32,13 @@ CANONICAL_CONTIGS = [f'chr{x}' for x in list(range(1, 23))] + ['chrX', 'chrY', 'chrM'] -def main(gff3_file: str, unmerged_output: str, merged_output: str, flanking: int = 2000): +def main( + gff3_file: str, + unmerged_output: str, + merged_output: str, + flanking: int = 2000, + json_output: str | None = None, +): """ Read the GFF3 file, and generate a BED file of gene regions, plus padding Args: @@ -39,21 +46,29 @@ def main(gff3_file: str, unmerged_output: str, merged_output: str, flanking: int unmerged_output (str): path to the intended BED output file merged_output (str): path to generate a BED file with merged overlapping rows flanking (int): number of bases to add before and after each gene + json_output (str | None): optional, path to write a Symbol: ENSG dictionary """ - unmerged_lines = generate_bed_lines(gff3_file, unmerged_output, flanking) + unmerged_lines, as_dict = generate_bed_lines(gff3_file, unmerged_output, flanking) merge_output(unmerged_lines, merged_output) + # if a json version was requested, write that out + if json_output is not None: + with open(json_output, 'w') as write_handle: + json.dump(as_dict, write_handle, indent=4) + def generate_bed_lines( gff3_file: str, output: str, flanking: int = FLANKING_REGION, -) -> list[tuple[str, int, int]]: +) -> tuple[list[tuple[str, int, int]], dict[str, str]]: """ Generate the new BED file, and return the lines as a list of lists for merging. """ output_lines: list[tuple[str, int, int]] = [] + output_as_dict: dict[str, str] = {} + # open and iterate over the GFF3 file with gzip.open(gff3_file, 'rt') as handle, open(output, 'w') as write_handle: for line in handle: @@ -85,6 +100,8 @@ def generate_bed_lines( print(f'Failed to extract gene name from {line_as_list[DETAILS_INDEX]}') continue + output_as_dict[gene_name] = gene_id + # write the line to the output output_list = [ f'chr{line_as_list[CHROM_INDEX]}', @@ -100,7 +117,7 @@ def generate_bed_lines( int(line_as_list[END_INDEX]) + flanking, ), ) - return output_lines + return output_lines, output_as_dict def merge_output( @@ -159,19 +176,25 @@ def cli_main(): parser.add_argument( '--merged_output', help='Path to output file, regions merged', - required=False, + required=True, ) parser.add_argument( '--flanking', help='Regions to add to each gene', default=FLANKING_REGION, ) + parser.add_argument( + '--json_output', + help='Write a {Gene Symbol: ENSG, } dictionary in JSON format.', + required=False, + ) args = parser.parse_args() main( gff3_file=args.gff3, unmerged_output=args.unmerged_output, merged_output=args.merged_output, flanking=args.flanking, + json_output=args.json_output, ) diff --git a/src/talos/cpg_internal_scripts/CPG_Dockerfile b/src/talos/cpg_internal_scripts/CPG_Dockerfile index e01bd3ff..093a7e37 100644 --- a/src/talos/cpg_internal_scripts/CPG_Dockerfile +++ b/src/talos/cpg_internal_scripts/CPG_Dockerfile @@ -14,46 +14,32 @@ RUN apt-get update && apt-get install -y --no-install-recommends \ FROM base AS bcftools_compiler -# BCFtools 1.22+ doesn't do any inference on chr-prefixes during CSQ, so the mismatch between fasta/vcf, and the gff3 (no chr-prefix) needs to be handled explicitly -# BCFtools >=1.22 is required internally as it's the first version with a built-in Mitochondrial lookup table. -# AS OF 11.0.3, Talos is building BCFtools from a private fork. This fork contains a single change - csq applies annotations -# to both coding and non-coding genes in the event of overlapping genes. By default BCFtools skips non-coding gene annotation -# if a coding transcript consequence was detected, but in practice this is masking clinically relevant non-coding gene variation -# in cases where the non-coding gene overlaps with a non-clinically relevant coding gene. -# The change made is to always search for non-coding variation, but to mask non-coding consequences in otherwise coding transcripts -# if a coding change was already detected. -ARG BCFTOOLS_VERSION=1.23.1 +ARG BCFTOOLS_VERSION=1.24 RUN apt-get update && apt-get install --no-install-recommends -y \ - autoconf \ gcc \ - git \ libbz2-dev \ libcurl4-openssl-dev \ liblzma-dev \ libssl-dev \ make \ zlib1g-dev && \ - wget https://github.com/samtools/htslib/releases/download/${BCFTOOLS_VERSION}/htslib-${BCFTOOLS_VERSION}.tar.bz2 && \ - tar -xf htslib-${BCFTOOLS_VERSION}.tar.bz2 && \ - cd htslib-${BCFTOOLS_VERSION} && \ - ./configure --enable-libcurl && \ - make && \ - make DESTDIR=/bcftools_install install && \ - cd .. && \ - git clone https://github.com/populationgenomics/bcftools.git && \ - cd bcftools && \ - autoheader && \ - autoconf && \ + rm -rf /var/lib/apt/lists/* && \ + wget https://github.com/samtools/bcftools/releases/download/${BCFTOOLS_VERSION}/bcftools-${BCFTOOLS_VERSION}.tar.bz2 && \ + tar -xf bcftools-${BCFTOOLS_VERSION}.tar.bz2 && \ + cd bcftools-${BCFTOOLS_VERSION} && \ ./configure --enable-libcurl --enable-s3 --enable-gcs && \ make && \ strip bcftools plugins/*.so && \ + make DESTDIR=/bcftools_install install && \ + cd htslib-${BCFTOOLS_VERSION} && \ + make && \ make DESTDIR=/bcftools_install install FROM base AS talos ARG ECHTVAR_VERSION=v0.2.2 -ARG VERSION=11.0.3 +ENV VERSION=11.1.0 COPY --from=bcftools_compiler /bcftools_install/usr/local/bin/* /usr/local/bin/ COPY --from=bcftools_compiler /bcftools_install/usr/local/libexec/bcftools/* /usr/local/libexec/bcftools/ diff --git a/src/talos/cpg_internal_scripts/cpgflow_jobs/annotate_csq_with_bcftools.py b/src/talos/cpg_internal_scripts/cpgflow_jobs/annotate_csq_with_bcftools.py index 9c518ebd..1c042b34 100644 --- a/src/talos/cpg_internal_scripts/cpgflow_jobs/annotate_csq_with_bcftools.py +++ b/src/talos/cpg_internal_scripts/cpgflow_jobs/annotate_csq_with_bcftools.py @@ -44,6 +44,7 @@ def make_bcftools_anno_jobs( f""" bcftools index -t {local_vcf} bcftools csq --force -f {fasta} \\ + --greedy 1 \\ --local-csq \\ --threads 4 \\ -g {gff3_file} \\ diff --git a/src/talos/cpg_internal_scripts/extract_fragmented_vcf_from_mt.py b/src/talos/cpg_internal_scripts/extract_fragmented_vcf_from_mt.py index b39a2a01..a72d34fd 100755 --- a/src/talos/cpg_internal_scripts/extract_fragmented_vcf_from_mt.py +++ b/src/talos/cpg_internal_scripts/extract_fragmented_vcf_from_mt.py @@ -54,6 +54,8 @@ }, } +CANONICAL_CONTIGS = {f'chr{x}' for x in list(range(1, 23))} | {'chrX', 'chrY'} + def filter_mt_to_sgids( mt: hl.MatrixTable, @@ -84,6 +86,12 @@ def filter_mt_to_sgids( return mt.filter_rows(hl.agg.any(mt.GT.is_non_ref())) +def filter_mt_to_contigs(mt: hl.MatrixTable) -> hl.MatrixTable: + """Strict filter to just the non-Mito contigs.""" + contig_literal = hl.literal(CANONICAL_CONTIGS) + return mt.filter_rows(contig_literal.contains(mt.locus.contig)) + + def main( mt_path: str, sg_id_file: str, @@ -109,6 +117,9 @@ def main( # filter the MT to a specific set of SG IDs mt = filter_mt_to_sgids(mt, sg_id_file) + # filter out chrM in this context + mt = filter_mt_to_contigs(mt) + # replace the existing INFO block to just have AC/AN/AF - no other carry-over. Allow for this to be missing. if 'AF' not in mt.info: mt = hl.variant_qc(mt) diff --git a/src/talos/cpg_internal_scripts/talos_stages.py b/src/talos/cpg_internal_scripts/talos_stages.py index 544aeb7b..60b7cfab 100644 --- a/src/talos/cpg_internal_scripts/talos_stages.py +++ b/src/talos/cpg_internal_scripts/talos_stages.py @@ -628,6 +628,26 @@ def queue_jobs(self, cohort: targets.Cohort, inputs: stage.StageInput) -> stage. )['vcf.bgz'] mito_vcf_arg = f'--labelled_mito {mito_vcf} ' + # is this run going to include STR data? + str_vcf_arg = '' + # screen out cohorts/datasets which we don't want to run STRs for + str_cohort = cohort.dataset.name in config.config_retrieve(['workflow', 'str_cohorts'], []) + if str_cohort and ( + str_vcf := query_for_latest_analysis( + dataset=cohort.dataset.name, + analysis_type='vcf', + long_read=config.config_retrieve(['workflow', 'long_read'], False), + stage_name='MakeStripyJointCall', + ) + ): + stripy_vcf = hail_batch.get_batch().read_input_group( + **{ + 'vcf.bgz': str_vcf, + 'vcf.bgz.tbi': f'{str_vcf}.tbi', + }, + )['vcf.bgz'] + str_vcf_arg = f'--str {stripy_vcf} ' + labelled_vcf = hail_batch.get_batch().read_input_group( **{ 'vcf.bgz': hail_inputs, @@ -657,7 +677,7 @@ def queue_jobs(self, cohort: targets.Cohort, inputs: stage.StageInput) -> stage. --labelled_vcf {labelled_vcf} \\ --output {job.output} \\ --panelapp {panelapp_data} \\ - --pedigree {pedigree} {sv_vcf_arg} {history_string} {mito_vcf_arg} + --pedigree {pedigree} {sv_vcf_arg} {history_string} {mito_vcf_arg} {str_vcf_arg} """, ) expected_out = self.expected_outputs(cohort) diff --git a/src/talos/create_talos_html.py b/src/talos/create_talos_html.py index 51eb4cad..4a630bcc 100644 --- a/src/talos/create_talos_html.py +++ b/src/talos/create_talos_html.py @@ -23,7 +23,16 @@ from loguru import logger from talos.config import config_retrieve -from talos.models import PanelApp, PanelDetail, PanelShort, ReportVariant, ResultData, SmallVariant, StructuralVariant +from talos.models import ( + PanelApp, + PanelDetail, + PanelShort, + ReportVariant, + ResultData, + ShortTandemRepeat, + SmallVariant, + StructuralVariant, +) from talos.utils import read_json_from_path JINJA_TEMPLATE_DIR = Path(__file__).absolute().parent / 'templates' @@ -157,7 +166,7 @@ def variant_in_forbidden_gene(variant_obj: ReportVariant, forbidden_genes): if gene_id in forbidden_genes: return True - if isinstance(variant_obj.var_data, StructuralVariant): + if isinstance(variant_obj.var_data, (ShortTandemRepeat, StructuralVariant)): return False # Allow for exclusion by Symbol too @@ -568,6 +577,11 @@ def get_var_change(self) -> str: if isinstance(self.var_data, StructuralVariant): return f'{self.var_data.info["svtype"]} {self.var_data.info["svlen"]}bp' + if isinstance(self.var_data, ShortTandemRepeat): + locus = self.var_data.locus or self.var_data.info.get('locus') + repeats = ', '.join(map(str, self.var_data.info['sample_repeats'])) + return f'STR {locus}, Repeats: {repeats}' + raise ValueError(f'Unknown variant type: {self.var_data.__class__.__name__}') def __init__(self, report_variant: ReportVariant, sample: Sample, ext_labels: list, html_builder: HTMLBuilder): # noqa: PLR0915 diff --git a/src/talos/download_panelapp.py b/src/talos/download_panelapp.py index 939936ae..f356cbca 100644 --- a/src/talos/download_panelapp.py +++ b/src/talos/download_panelapp.py @@ -25,6 +25,8 @@ - attempt to find alternative gene symbols for the ENSG ID - attempt to find alternative Ensembl IDs for the gene symbol - record all variations + +STR-specific feature - parses the PanelApp Repeat Disorders panel """ import asyncio @@ -66,6 +68,7 @@ # if this is a massive result, it returns over a number of pages PANEL_TEMPLATE_URL = f'{PANELS_ENDPOINT}/{{id}}' ACTIVITY_TEMPLATE = f'{PANELS_ENDPOINT}/{{id}}/activities' +STR_TEMPLATE = f'{PANELS_ENDPOINT}/3597/strs?confidence_level=3' MITO_BAD = 'MT' MITO_GOOD = 'M' @@ -311,6 +314,24 @@ def reorganise_mane_data(mane_path: str) -> tuple[dict[str, str], dict[str, str] return ensg_as_primary, symbol_as_primary +def parse_repeat_disorders() -> tuple[set[str], set[str]]: + """Parse panel 3597 - find all genes with a green association to a STR disorder.""" + str_genes: set[str] = set() + str_symbols: set[str] = set() + for each_result in get_json_response(STR_TEMPLATE)['results']: + str_association = each_result['gene_data'] + str_symbols.add(str_association['gene_symbol']) + + for build, content in str_association['ensembl_genes'].items(): + if build.lower() == 'grch38': + # the ensembl version may alter over time, but will be singular + ensembl_data = content[next(iter(content.keys()))] + ensg = ensembl_data['ensembl_id'] + str_genes.add(ensg) + + return str_genes, str_symbols + + def cli_main(): logger.info('Starting PanelApp parsing') parser = ArgumentParser() @@ -408,6 +429,13 @@ async def _fetch_all() -> tuple[dict, dict]: logger.info(f'Removing panel {panel_id} from hpo matching - no green genes') del collected_panel_data.hpos[panel_id] + # query panelapp for the repeat disorders panel + str_genes, str_symbols = parse_repeat_disorders() + + # populate the panelapp object + collected_panel_data.str_genes = str_genes + collected_panel_data.str_symbols = str_symbols + with open(output, 'w') as output_file: output_file.write(collected_panel_data.model_dump_json(indent=4)) diff --git a/src/talos/example_config.toml b/src/talos/example_config.toml index 38f21fcf..679683dd 100644 --- a/src/talos/example_config.toml +++ b/src/talos/example_config.toml @@ -112,6 +112,9 @@ support_categories = ['alphamissense', 'clinvar0star'] #super_logging = true #super_logging_path = "exclusions.jsonl.gz" +# any STR loci with unacceptably low signal-to-noise ratios - this is the 'locus' name from STRipy, not the gene ID +noisy_strs = ['ARX_1', 'ARX_2', 'HOXA13_1', 'HOXA13_2', 'HOXA13_3', 'TCF4', 'RUNX2'] + [RunHailFiltering] # variables affecting how the VCF variants are parsed, and AnalysisVariant objects are populated csq_string = [ "consequence", "gene_id", "gene", "transcript", "mane_id", "mane", "biotype", "dna_change", "amino_acid_change", "codon", "ensp", "am_class", "am_pathogenicity",] diff --git a/src/talos/liftover/lift_1_1_0_to_1_2_0.py b/src/talos/liftover/lift_1_1_0_to_1_2_0.py index 51e44e62..e37da122 100644 --- a/src/talos/liftover/lift_1_1_0_to_1_2_0.py +++ b/src/talos/liftover/lift_1_1_0_to_1_2_0.py @@ -2,8 +2,6 @@ code for lifting over models from 1.1.0 to 1.2.0 """ -from talos.config import config_retrieve - def resultdata(data_dict: dict) -> dict: """ @@ -22,9 +20,9 @@ def resultdata(data_dict: dict) -> dict: _ = variant['var_data'].pop('sample_support') # the list of categories which are being treated as support for this run - variant['var_data']['support_categories'] = config_retrieve(['ValidateMOI', 'support_categories'], []) + variant['var_data']['support_categories'] = [] # the list of categories being ignored for this run - variant['var_data']['ignored_categories'] = config_retrieve(['ValidateMOI', 'ignore_categories'], []) + variant['var_data']['ignored_categories'] = [] data_dict['version'] = '1.2.0' return data_dict diff --git a/src/talos/liftover/lift_2_2_0_to_2_3_0.py b/src/talos/liftover/lift_2_2_0_to_2_3_0.py index a2362a0e..3bfd602c 100644 --- a/src/talos/liftover/lift_2_2_0_to_2_3_0.py +++ b/src/talos/liftover/lift_2_2_0_to_2_3_0.py @@ -3,7 +3,14 @@ """ +def panelapp(data_dict: dict) -> dict: + data_dict |= {'str_genes': set(), 'str_symbols': set()} + data_dict['version'] = '2.3.0' + return data_dict + + def dl_panelapp(data_dict: dict) -> dict: + data_dict |= {'str_genes': set(), 'str_symbols': set()} for _, gene_data in data_dict['genes'].items(): gene_data['panel_confidences'] = {} for _, panel_data in gene_data['panels'].items(): diff --git a/src/talos/models.py b/src/talos/models.py index bfdca7e5..24062209 100644 --- a/src/talos/models.py +++ b/src/talos/models.py @@ -21,6 +21,7 @@ from talos.liftover.lift_2_1_0_to_2_2_0 import dl_panelapp as dl_pa_210_to_220 from talos.liftover.lift_2_1_0_to_2_2_0 import resultdata as rd_210_to_220 from talos.liftover.lift_2_2_0_to_2_3_0 import dl_panelapp as dl_pa_220_to_230 +from talos.liftover.lift_2_2_0_to_2_3_0 import panelapp as pa_220_to_230 from talos.liftover.lift_none_to_1_0_0 import resultdata as rd_none_to_1_0_0 from talos.static_values import get_granular_date @@ -330,13 +331,20 @@ def min_alt_ratio(self, sample: str, threshold: float = 0.2) -> bool: class StructuralVariant(VariantCommon): - """ - placeholder for any methods/data specific to Structural Variants - """ + """Placeholder for any methods/data specific to Structural Variants.""" + + +class ShortTandemRepeat(VariantCommon): + """Sub-class designed to parse STR VCFs into a suitable common format.""" + + locus: str + sample_filter: dict[str, str] = Field(default_factory=dict, exclude=True) + sample_repeats: dict[str, tuple[int, int] | tuple[int]] = Field(default_factory=dict, exclude=True) + sample_repeat_details: dict[str, dict] = Field(default_factory=dict, exclude=True) # register all interchangeable models here -VARIANT_MODELS = SmallVariant | StructuralVariant +VARIANT_MODELS = SmallVariant | StructuralVariant | ShortTandemRepeat class ReportPanel(BaseModel): @@ -438,6 +446,8 @@ class PanelApp(BaseModel): participants: dict[str, ParticipantHPOPanels] = Field(default_factory=dict) version: str = CURRENT_VERSION creation_date: str = Field(default=get_granular_date()) + str_genes: set[str] = Field(default_factory=set) + str_symbols: set[str] = Field(default_factory=set) class DownloadedPanelAppGenePanelDetail(BaseModel): @@ -468,6 +478,8 @@ class DownloadedPanelApp(BaseModel): hpos: dict[int, list[HpoTerm]] = Field(default_factory=dict) version: str = CURRENT_VERSION date: str = Field(default=get_granular_date()) + str_genes: set[str] = Field(default_factory=set) + str_symbols: set[str] = Field(default_factory=set) class ResultMeta(BaseModel): @@ -567,6 +579,7 @@ def __str__(self): PanelApp: { '1.2.0_2.0.0': pa_120_to_200, '2.0.0_2.1.0': pa_200_to_210, + '2.2.0_2.3.0': pa_220_to_230, }, ResultData: { 'None_1.0.0': rd_none_to_1_0_0, diff --git a/src/talos/moi_tests.py b/src/talos/moi_tests.py index e488a3e7..78acd4bf 100644 --- a/src/talos/moi_tests.py +++ b/src/talos/moi_tests.py @@ -14,7 +14,7 @@ from talos.config import config_retrieve from talos.exclusion_log import get_exclusion_logger -from talos.models import VARIANT_MODELS, ReportVariant, SmallVariant, StructuralVariant +from talos.models import VARIANT_MODELS, ReportVariant, ShortTandemRepeat, SmallVariant, StructuralVariant from talos.static_values import get_granular_date from talos.utils import X_CHROMOSOME, CompHetDict @@ -23,6 +23,19 @@ SV_HOMS = {'male_n_homalt', 'female_n_homalt'} +def get_str_var_data(var: VARIANT_MODELS, sample_id: str) -> VARIANT_MODELS: + """For STRs, return a shallow copy with a sample-scoped info dict. Other variant types pass through unchanged.""" + if isinstance(var, ShortTandemRepeat): + copy = var.model_copy(deep=False) + copy.info = { + **var.info, + 'sample_repeats': var.sample_repeats[sample_id], + 'sample_repeat_details': var.sample_repeat_details[sample_id], + } + return copy + return var + + @dataclass class GlobalFilter: """ @@ -46,16 +59,19 @@ class GlobalFilter: small_gnomad_hemi: ClassVar[int] = config_retrieve(['ValidateMOI', 'gnomad_max_hemizygotes']) # filters specific to SVs + # todo the field could have a different name sv_dict: ClassVar[dict[str, float]] = { 'gnomad_v2.1_sv_AF': config_retrieve(['ValidateMOI', 'gnomad_sv_max_af']), } - def too_common(self, variant: SmallVariant | StructuralVariant, applied_moi: str | None = None) -> bool: # noqa: PLR0911 + def too_common( # noqa: PLR0911 + self, variant: SmallVariant | ShortTandemRepeat | StructuralVariant, applied_moi: str | None = None + ) -> bool: """ Check if a variant is too common in the population Args: - variant (SmallVariant | StructuralVariant): the variant to check + variant (SmallVariant | ShortTandemRepeat | StructuralVariant): the variant to check applied_moi (str | None): MOI under which this filter is being applied (for exclusion logging only) Returns: @@ -149,6 +165,9 @@ def too_common(self, variant: SmallVariant | StructuralVariant, applied_moi: str ) return True + elif isinstance(variant, ShortTandemRepeat): + return False + else: raise ValueError('Variant type not recognised') @@ -180,12 +199,14 @@ class DominantFilter: 'gnomad_v2.1_sv_AF': config_retrieve(['ValidateMOI', 'dominant_gnomad_sv_max_af']), } - def too_common(self, variant: SmallVariant | StructuralVariant, applied_moi: str | None = None) -> bool: + def too_common( + self, variant: SmallVariant | ShortTandemRepeat | StructuralVariant, applied_moi: str | None = None + ) -> bool: """ Check if a variant is too common in the population Args: - variant (SmallVariant | StructuralVariant): the variant to check + variant (SmallVariant | ShortTandemRepeat | StructuralVariant): the variant to check applied_moi (str | None): MOI under which this filter is being applied (for exclusion logging only) Returns: @@ -261,6 +282,9 @@ def too_common(self, variant: SmallVariant | StructuralVariant, applied_moi: str ) return True + elif isinstance(variant, ShortTandemRepeat): + return False + else: raise ValueError('Variant type not recognised') @@ -559,6 +583,15 @@ def sample_passes_basic_checks( Returns True if the sample passes; otherwise logs the first failing reason and returns False. Caller is responsible for any additional MOI-specific gates (sex, zygosity). """ + if sample_id not in self.pedigree.participants: + self._log_sample_exclusion( + principal, + sample_id, + 'participant_not_in_pedigree', + details={'allow_support': allow_support}, + ) + return False + if self.pedigree.participants[sample_id].is_not_affected: # not logging, self-evident # self._log_sample_exclusion(principal, sample_id, 'sample_not_affected') @@ -802,7 +835,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, reasons=self.applied_moi, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), @@ -877,7 +910,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, reasons=self.applied_moi, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), @@ -952,7 +985,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), reasons=self.applied_moi, @@ -1017,7 +1050,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, reasons=self.applied_moi, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), @@ -1105,7 +1138,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, reasons=self.applied_moi, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), @@ -1172,7 +1205,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), reasons=self.applied_moi, @@ -1238,7 +1271,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), reasons=self.applied_moi, @@ -1315,7 +1348,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, reasons=self.applied_moi, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), @@ -1387,7 +1420,7 @@ def run( sample=sample_id, family=self.pedigree.participants[sample_id].family_id, gene=principal.info.get('gene_id'), - var_data=principal, + var_data=get_str_var_data(principal, sample_id), categories={key: get_granular_date() for key in principal.category_values(sample_id)}, reasons=self.applied_moi, genotypes=self.get_family_genotypes(variant=principal, sample_id=sample_id), diff --git a/src/talos/run_hail_filtering.py b/src/talos/run_hail_filtering.py index 04d731cb..3018d94a 100644 --- a/src/talos/run_hail_filtering.py +++ b/src/talos/run_hail_filtering.py @@ -778,7 +778,7 @@ def get_csq_from_struct(element: hl.expr.StructExpression) -> hl.expr.StringExpr return hl.delimit([hl.or_else(hl.str(fields.get(f, '')), '') for f in csq_fields], '|') csq = hl.empty_array(hl.tstr) - csq = csq.extend(hl.or_else(tx_expr.map(lambda x: get_csq_from_struct(x)), hl.empty_array(hl.tstr))) # noqa: PLW0108 + csq = csq.extend(hl.or_else(tx_expr.map(lambda x: get_csq_from_struct(x)), hl.empty_array(hl.tstr))) # previous consequence filters may make this caution unnecessary return hl.or_missing(hl.len(csq) > 0, csq) diff --git a/src/talos/scripts/stripy_json_to_vcf.py b/src/talos/scripts/stripy_json_to_vcf.py new file mode 100755 index 00000000..c4ffecff --- /dev/null +++ b/src/talos/scripts/stripy_json_to_vcf.py @@ -0,0 +1,386 @@ +#!/usr/bin/env python3 + +""" +Borrowing _heavily_ from https://github.com/broadinstitute/gatk-sv/blob/7eb2af1feea9f81a390caff211f6832c6e7f7ac5/src/stripy/stripy_to_vcf.py +Has been extended to take multiple single-sample STRipy reports as input, and to characterise the GT field as: +- '0/0' 2 called alleles, no pathogenic +- '1/0' a1 pathogenic, a2 non-pathogenic +- '1/.' a1 pathogenic, a2 missing + +Uses fixed contigs (GRCh38) for the header, instead of relying on a reference input + +This script creates VCFs in the format expected by Talos. Specific fields used: + - INFO/GENE - used to detect the appropriate MOI for this STR + - INFO/SYMBOL - used in variant display + - FORMAT/REPCN - counts detected at each allele + - FORMAT/STR_FILTER - rule out quality-failed STR calls + - FORMAT/DISEASE_DETAILS - dense formatted String, presenting each disease and corresponding thresholds + +This implementation of the STRipy VCF generator requires a symbol: ensg lookup dictionary +This is now provided as an additional input from the prep script/module CreateRoiFromGff3 +""" + +import argparse +import json +import logging +import re +from pathlib import Path + +import pysam + +from cpg_utils import config, to_path + +CONTIG_ORDER = [f'chr{x}' for x in list(range(1, 23))] + ['chrX', 'chrY', 'chrM'] +HEADER_LINES = [ + ('chr1', 248956422), + ('chr2', 242193529), + ('chr3', 198295559), + ('chr4', 190214555), + ('chr5', 181538259), + ('chr6', 170805979), + ('chr7', 159345973), + ('chr8', 145138636), + ('chr9', 138394717), + ('chr10', 133797422), + ('chr11', 135086622), + ('chr12', 133275309), + ('chr13', 114364328), + ('chr14', 107043718), + ('chr15', 101991189), + ('chr16', 90338345), + ('chr17', 83257441), + ('chr18', 80373285), + ('chr19', 58617616), + ('chr20', 64444167), + ('chr21', 46709983), + ('chr22', 50818468), + ('chrX', 156040895), + ('chrY', 57227415), + ('chrM', 16569), +] + +VCF_HEADER = { + 'INFO': [ + {'ID': 'END', 'Number': '1', 'Type': 'Integer', 'Description': 'Stop position of the interval'}, + {'ID': 'SVTYPE', 'Number': '1', 'Type': 'String', 'Description': 'Type of structural variant'}, + {'ID': 'RU', 'Number': '1', 'Type': 'String', 'Description': 'Repeat unit in the reference orientation'}, + {'ID': 'PERIOD', 'Number': '1', 'Type': 'Integer', 'Description': 'Length of the repeat unit'}, + { + 'ID': 'DISEASES', + 'Number': '.', + 'Type': 'String', + 'Description': 'Associated disease symbols for this locus (| separated)', + }, + {'ID': 'LOCUS', 'Number': '1', 'Type': 'String', 'Description': 'Gene/locus identifier from STRipy'}, + {'ID': 'SYMBOL', 'Number': '1', 'Type': 'String', 'Description': 'Gene inferred from STRipy locus identifier'}, + {'ID': 'GENE', 'Number': '1', 'Type': 'String', 'Description': 'ENSG ID inferred from Gene Symbol'}, + ], + 'FORMAT': [ + {'ID': 'GT', 'Number': '1', 'Type': 'String', 'Description': 'Unphased genotype'}, + {'ID': 'REPCN', 'Number': '2', 'Type': 'Float', 'Description': 'Number of repeat units spanned by each allele'}, + { + 'ID': 'DISEASE_DETAILS', + 'Number': '1', + 'Type': 'String', + 'Description': '|-delimited details for each disease, in the form diseaseSymbol__normal__intermediate__pathogenic', # noqa: E501 + }, + { + 'ID': 'REPCI1', + 'Number': '2', + 'Type': 'Integer', + 'Description': '95% CI min,max on repeat counts of first allele', + }, + { + 'ID': 'REPCI2', + 'Number': '2', + 'Type': 'Integer', + 'Description': '95% CI min,max on repeat counts of second allele', + }, + { + 'ID': 'OUTLIER', + 'Number': '2', + 'Type': 'Integer', + 'Description': 'Allelic population outlier flags (0/1) assigned by STRipy', + }, + { + 'ID': 'ZSCORE', + 'Number': '2', + 'Type': 'Float', + 'Description': 'Allelic population Z-scores assigned by STRipy', + }, + {'ID': 'DP', 'Number': '1', 'Type': 'Integer', 'Description': 'Total Depth'}, + {'ID': 'STR_FILTER', 'Number': '.', 'Type': 'String', 'Description': 'Filter status assigned by STRipy'}, + ], +} + + +def get_limited_locus_list() -> list[str] | None: + """Try to get the finite list of loci from config, if available.""" + try: + return config.config_retrieve(['stripy', 'loci_lists', 'default_with_exclusions']) + except config.ConfigError: + print('No CPG config detected, using all loci') + return None + + +def parse_coords(coord: str) -> tuple[str, int, int] | None: + m = re.match(r'(chr)?(?P[^:]+):(?P\d+)-(?P\d+)', str(coord)) + if not m: + return None + chrom = m.group('chrom') + if not chrom.startswith('chr'): + chrom = 'chr' + chrom + return chrom, int(m.group('start')), int(m.group('end')) + + +def _to_float(x: str | float) -> float: + try: + return float(x) + except (TypeError, ValueError): + return float('nan') + + +def _to_outlier(x: str | int) -> int | None: + try: + return int(x) + except (TypeError, ValueError): + return None + + +def _ci_tuple(allele: dict[str, dict[str, int]]) -> tuple[int, int]: + return allele['CI'].get('Min', 0), allele['CI'].get('Max', 0) + + +def parse_disease_ranges(content: dict[str, str | dict[str, int]]) -> str: + """ + Read the CorrespondingDisease content of the locus' JSON dictionary. + This creates a hybrid String of Gene, MOI, normal and intermediate ranges, and pathogenic threshold. + These occur per-disease, and there can be multiple per allele. + + example results: "HD__AD__min9max26__min27max35__36" or "DMD__XLR__min11max33__.__59" + + Args: + content: the dictionary from STRipy's CorrespondingDisease data block + + Returns: + A single aggregated String for all the details + """ + # get the symbol for this disease + disease = content['DiseaseSymbol'] + + # the inheritance pattern + inheritance = content['Inheritance'] + + # parse the normal range, or '.' + normal_range: dict[str, int] | str = content['NormalRange'] + normal = '.' + if isinstance(normal_range, dict): + normal = f'min{normal_range["Min"]}max{normal_range["Max"]}' + + # parse the intermediate range, or '.' + inter_range: dict[str, int] | str = content['IntermediateRange'] + intermediate = '.' + if isinstance(inter_range, dict): + intermediate = f'min{inter_range["Min"]}max{inter_range["Max"]}' + + # add the pathogenic cutoff used for this disease, pull it all together + return f'{disease}__{inheritance}__{normal}__{intermediate}__{content["PathogenicCutoff"]}' + + +def load_sample(json_path: str) -> tuple[str, dict]: + """Load a STRipy JSON and return (sample_name, loci_dict keyed by locus_id).""" + with to_path(json_path).open() as f: + data = json.load(f) + + # hard-wired to pull the SampleID from the input filename - consider reverting + input_file = data.get('JobDetails', {}).get('InputFile', '') + if input_file: + stem = Path(input_file).stem + sample_name = stem.split('__')[0] + else: + sample_name = Path(json_path).stem.split('.')[0] + + loci = {} + for entry in data.get('GenotypingResults', []): + for locus_name, locus in entry.items(): + tl = locus['TargetedLocus'] + + # create a |-delimited string for all ranges for all relevant disease for this locus + disease_deets = '|'.join(parse_disease_ranges(corr) for corr in tl['CorrespondingDisease'].values()) + + # add flexibility if the Alleles aren't populated at all + if 'Alleles' not in locus: + explain_string = f'Parsing {sample_name}, locus {tl["LocusID"]}, no Alleles present - skipping. ' + if 'Filter' in locus: + explain_string += f'Filter: {locus["Filter"]} - ' + logging.debug(explain_string) + continue + + parsed = parse_coords(tl['Coordinates']) + if not parsed: + continue + chrom, start, end = parsed + alleles = locus['Alleles'] + a1 = alleles[0] + a2 = alleles[1] if len(alleles) > 1 else None + motif = tl['Motif'] + loci[locus_name] = { + 'chrom': chrom, + 'pos': start, + 'end': end, + 'id': locus_name, + 'motif': motif, + 'period': len(motif), + 'a1_rep': _to_float(a1['Repeats']), + 'a2_rep': _to_float(a2['Repeats']) if a2 is not None else None, + 'a1_ci': _ci_tuple(a1), + 'a2_ci': _ci_tuple(a2) if a2 is not None else (None, None), + 'a1_range': a1['Range'], + 'a2_range': a2['Range'] if a2 else None, + 'a1_out': _to_outlier(a1['IsPopulationOutlier']), + 'a2_out': _to_outlier(a2['IsPopulationOutlier']) if a2 is not None else None, + 'a1_z': _to_float(a1['PopulationZscore']), + 'a2_z': _to_float(a2['PopulationZscore']) if a2 is not None else None, + 'coverage': locus['Metadata']['Coverage'], + 'filter': locus['Filter'], + 'diseases': '|'.join(sorted({meta['DiseaseSymbol'] for meta in tl['CorrespondingDisease'].values()})), + 'disease_details': disease_deets, + } + + return sample_name, loci + + +def get_header(sample_names: list[str]) -> pysam.VariantHeader: + """Fill the header.""" + header = pysam.VariantHeader() + header.add_meta('source', value='STRipy2VCF') + header.add_meta('ALT', items=[('ID', 'STR'), ('Description', 'Short tandem repeat')]) + for info in VCF_HEADER['INFO']: + header.add_meta('INFO', items=list(info.items())) + for fmt in VCF_HEADER['FORMAT']: + header.add_meta('FORMAT', items=list(fmt.items())) + + # generate the contig header lines + for contig, length in HEADER_LINES: + header.contigs.add(contig, length=length) + + for sample_name in sample_names: + header.add_sample(sample_name) + + return header + + +def convert_range_to_gt(range_val: str | None) -> int | None: + """Converts the per-allele range to an integer, for embedding in GT.""" + if range_val is None: + return None + if range_val.lower() == 'pathogenic': + return 1 + return 0 + + +def get_gene_lookup(map_path: str | None = None) -> dict[str, str]: + """Pull out the dictionary reading logic.""" + lookup: dict[str, str] = {} + if map_path is not None: + with to_path(map_path).open() as handle: + lookup = json.load(handle) + return lookup + + +def write_multisample_vcf( # noqa: PLR0915 + samples: list[tuple[str, dict[str, dict]]], + out_path: str, + gene_map: str | None = None, +) -> None: + """ + Integrate the per-sample data into a union VCF + + Args: + samples: list of (sample_name, loci_dict) tuples + out_path: vcf path to write output to. PySam will write compressed if ending is gz/bgz + gene_map: optional, file path to a JSON dict, mapping gene symbols to gene IDs + """ + + gene_lookup = get_gene_lookup(gene_map) + + header = get_header(sample_names=[sam_bit[0] for sam_bit in samples]) + + only_loci = get_limited_locus_list() + + canonical: dict[str, dict] = {} + + for _, loci in samples: + for locus_id, loc in loci.items(): + # if we got a list of specific loci, only use that finite list - ignore everything else + if (only_loci is not None) and (locus_id not in only_loci): + continue + canonical.setdefault(locus_id, loc) + + # double sort contigs by both chrom and position + sorted_loci = sorted(canonical.values(), key=lambda x: (CONTIG_ORDER.index(x['chrom']), x['pos'])) + + with pysam.VariantFile(out_path, mode='w', header=header) as vf: + for loc in sorted_loci: + rec = header.new_record() + rec.contig = loc['chrom'] + rec.start = loc['pos'] - 1 + rec.stop = loc['end'] + rec.id = str(loc['id']) + rec.ref = 'N' + rec.alts = ('',) + rec.filter.add('PASS') + rec.info['SVTYPE'] = 'STR' + if loc['motif']: + rec.info['RU'] = loc['motif'] + if loc['period']: + rec.info['PERIOD'] = int(loc['period']) + rec.info['DISEASES'] = loc['diseases'] + rec.info['LOCUS'] = loc['id'] + + # extract the gene symbol from STRipy annotations + symbol = loc['id'].split('_')[0] + rec.info['SYMBOL'] = symbol + rec.info['GENE'] = gene_lookup.get(symbol, symbol) + + for sample_name, loci in samples: + s_loc = loci.get(loc['id']) + s = rec.samples[sample_name] + if s_loc is None: + s['GT'] = (None, None) + s['REPCN'] = (None, None) + s['REPCI1'] = (None, None) + s['REPCI2'] = (None, None) + s['OUTLIER'] = (None, None) + s['ZSCORE'] = (None, None) + s['DP'] = None + s['STR_FILTER'] = ['.'] + s['DISEASE_DETAILS'] = '.' + else: + s['GT'] = (convert_range_to_gt(s_loc['a1_range']), convert_range_to_gt(s_loc['a2_range'])) + s['REPCN'] = (s_loc['a1_rep'], s_loc['a2_rep']) + s['REPCI1'] = s_loc['a1_ci'] + s['REPCI2'] = s_loc['a2_ci'] + s['OUTLIER'] = (s_loc['a1_out'], s_loc['a2_out']) + s['ZSCORE'] = (s_loc['a1_z'], s_loc['a2_z']) + s['DP'] = int(s_loc['coverage']) + s['STR_FILTER'] = [str(s_loc['filter'])] if s_loc['filter'] else ['PASS'] + s['DISEASE_DETAILS'] = s_loc['disease_details'] + + vf.write(rec) + + +def main() -> None: + ap = argparse.ArgumentParser(description='Convert STRipy JSON output(s) into a multi-sample VCF.') + ap.add_argument('--json', required=True, nargs='+', help='STRipy JSON report(s); one per sample') + ap.add_argument('--output', required=True, help='Output VCF') + ap.add_argument('--mapping', default=None, help='A JSON dict of Gene Symbol:Gene ID to update annotation') + args = ap.parse_args() + + samples = [load_sample(p) for p in args.json] + + write_multisample_vcf(samples=samples, out_path=args.output, gene_map=args.mapping) + + +if __name__ == '__main__': + main() diff --git a/src/talos/templates/index.html.jinja b/src/talos/templates/index.html.jinja index c6b7400d..62a4e5f1 100644 --- a/src/talos/templates/index.html.jinja +++ b/src/talos/templates/index.html.jinja @@ -453,7 +453,7 @@ var html = '
'; - html += '
'; + html += '
'; html += '
Variant Details
'; html += '
Coordinates
'; @@ -608,6 +608,65 @@ } } html += '
'; + } else if (rowData.variantType === 'ShortTandemRepeat' && variantData.info) { + var strRepeats = variantData.info.sample_repeats; + var strDetails = variantData.info.sample_repeat_details; + + if (strRepeats && strRepeats.length > 0) { + html += '
Repeat Counts
'; + strRepeats.forEach(function(count, i) { + html += 'Allele ' + (i + 1) + ': ' + escapeHtml(String(count)) + '
'; + }); + html += '
'; + } + + if (strDetails && Object.keys(strDetails).length > 0) { + var maxRepeats = strRepeats && strRepeats.length > 0 ? Math.max.apply(null, strRepeats) : null; + html += '
Disease Classification
'; + html += '
'; + html += ''; + html += ''; + + Object.keys(strDetails).forEach(function(disease) { + var d = strDetails[disease]; + var norm = d.norm ? d.norm[0] + '–' + d.norm[1] : '—'; + var inter = d.inter ? d.inter[0] + '–' + d.inter[1] : '—'; + var pathThresh = d.pathogenic && d.pathogenic !== '.' ? '≥' + d.pathogenic : '—'; + + var badge = 'UNCLASSIFIED'; + if (maxRepeats !== null) { + var pathNum = d.pathogenic && d.pathogenic !== '.' ? parseInt(d.pathogenic, 10) : null; + if (pathNum !== null && maxRepeats >= pathNum) { + badge = 'PATHOGENIC'; + } else if (d.inter && maxRepeats >= d.inter[0] && maxRepeats <= d.inter[1]) { + badge = 'INTERMEDIATE'; + } else if (d.norm && maxRepeats >= d.norm[0] && maxRepeats <= d.norm[1]) { + badge = 'NORMAL'; + } + } + + html += ''; + html += ''; + html += ''; + html += ''; + html += ''; + html += ''; + html += ''; + html += ''; + }); + + html += '
DiseaseMOINormalIntermediatePathogenicClassification
' + escapeHtml(disease) + '' + escapeHtml(d.moi || '—') + '' + escapeHtml(norm) + '' + escapeHtml(inter) + '' + escapeHtml(pathThresh) + '' + badge + '
'; + } + + if (rowData.genes) { + var strGenes = rowData.genes.split(','); + var strFirst = strGenes[0] ? strGenes[0].split(':') : ['', '']; + if (strFirst[1]) { + html += '
External Links
'; + html += 'OMIM'; + html += '
'; + } + } } html += ''; diff --git a/src/talos/unified_panelapp_parser.py b/src/talos/unified_panelapp_parser.py index cfbd8c87..cad57359 100644 --- a/src/talos/unified_panelapp_parser.py +++ b/src/talos/unified_panelapp_parser.py @@ -430,6 +430,10 @@ def main(panel_data: str, output_file: str, pedigree_path: str, hpo_file: str | panelapp_data.genes = {key: value for key, value in panelapp_data.genes.items() if key not in genes_to_remove} + # add in the STR content + panelapp_data.str_genes = cached_panelapp.str_genes + panelapp_data.str_symbols = cached_panelapp.str_symbols + # validate and write using pydantic valid_cohort_details = PanelApp.model_validate(panelapp_data) with open(output_file, 'w', encoding='utf-8') as handle: diff --git a/src/talos/utils.py b/src/talos/utils.py index 99f12f85..f353946e 100644 --- a/src/talos/utils.py +++ b/src/talos/utils.py @@ -12,6 +12,7 @@ import zoneinfo from collections import defaultdict from datetime import datetime +from functools import cache from itertools import chain, combinations, combinations_with_replacement, islice from random import choices from typing import TYPE_CHECKING, Any @@ -21,13 +22,16 @@ from loguru import logger from mendelbrot.bcftools_interpreter import TYPES_RE, classify_change from mendelbrot.pedigree_parser import PedigreeParser +from numpy import isnan from tenacity import retry, retry_if_exception_type, stop_after_attempt, wait_exponential_jitter from talos.config import config_retrieve from talos.models import ( VARIANT_MODELS, Coordinates, + PanelApp, ResultData, + ShortTandemRepeat, SmallVariant, StructuralVariant, lift_up_model_version, @@ -75,6 +79,8 @@ # not in use yet, if we want to reduce the number of TranscriptConsequence entries we store, this is a good start BORING_CONSEQUENCES = ['downstream_gene_variant', 'intron_variant', 'upstream_gene_variant'] +STR_RANGE = re.compile(r'min(?P[0-9]+)max(?P[0-9]+)$') + def parse_mane_json_to_dict(mane_json: str) -> dict: """ @@ -134,6 +140,14 @@ def generator_chunks(generator, size): yield list(chain([first], islice(iterator, size - 1))) +@cache +def get_categories_translated(conf_string: str = 'support_categories') -> set[str]: + """One central method to parse categories from config.""" + support_categories = {word.lower() for word in config_retrieve(['ValidateMOI', conf_string], [])} + support_categories.update({translate_category(x) for x in support_categories}) + return support_categories + + @retry( stop=stop_after_attempt(3), wait=wait_exponential_jitter(initial=1, max=5, exp_base=2), @@ -399,6 +413,51 @@ def organise_svdb_doi(info_dict: dict[str, Any]): info_dict['svdb_doi'] = doi_urls +def parse_str_disease_details(detail_string: str) -> dict[str, dict[str, str | tuple[int, int] | None]]: + """ + Parse the pipe-delimited disease details for this sample & locus. + Each block is in the format "Disease__MOI__NormalRange__IntermediateRange__PathogenicThreshold". + Ranges can be missing: "." + Other values can be missing: "NA" + + Args: + detail_string: the single pipe-delimited String + + Returns: + dict of each disease name and its related thresholds and ranges + """ + + results: dict[str, dict] = {} + + # missing value = no details, return an empty dict + if detail_string == '.': + return results + + # parse each disease's data separately + for disease_block in detail_string.split('|'): + # split up the formatted string + gene, moi, norm, inter, pathogenic = disease_block.split('__') + + # try and detect the normal and intermediate ranges + normal_range: None | tuple[int, int] = None + if matchy := re.match(STR_RANGE, norm): + normal_range = (int(matchy.group('min')), int(matchy.group('max'))) + + inter_range: None | tuple[int, int] = None + if matchy := re.match(STR_RANGE, inter): + inter_range = (int(matchy.group('min')), int(matchy.group('max'))) + + # build the + results[gene] = { + 'moi': moi, + 'norm': normal_range, + 'inter': inter_range, + 'pathogenic': pathogenic, + } + + return results + + def organise_de_novo(info_dict: dict[str, Any], alt_depths: dict[str, int], ab_ratios: dict[str, float]) -> None: """ apply some late checking on de novo attributes @@ -424,10 +483,65 @@ def organise_de_novo(info_dict: dict[str, Any], alt_depths: dict[str, int], ab_r # if we detected any failing samples against these rules, fish them out if to_pop: - logger.info(f'Removing de novo status for {len(to_pop)} samples: {", ".join(to_pop)}') + logger.debug(f'Removing de novo status for {len(to_pop)} samples: {", ".join(to_pop)}') info_dict['categorysampledenovo'] = [sam for sam in info_dict['categorysampledenovo'] if sam not in to_pop] +def create_str_variant( + var: 'cyvcf2.Variant', + samples: list[str], +) -> ShortTandemRepeat | None: + """Takes a variant row from an STR 'joint-call', and parses into a Variant object.""" + + ignored_categories = get_categories_translated('ignore_categories') + if 'categorybooleanstr' in ignored_categories: + return None + + coordinates = Coordinates(chrom=var.CHROM.replace('chr', ''), pos=var.POS, ref='STR', alt=var.ID) + + # parse the info dict, and shove in a recognisable category label + boolean_categories = ['categorybooleanstr'] + info: dict[str, Any] = {x.lower(): y for x, y in var.INFO} | {'categorybooleanstr': 1} + + if info.get('locus') in config_retrieve(['ValidateMOI', 'noisy_strs'], []): + return None + + # shuffle the gene ID (ENSG) to match the location Hail classifiers put it + info['gene_id'] = info.pop('gene') + + het_samples, hom_samples = get_non_ref_samples(variant=var, samples=samples) + + # no samples with an STR call, reject the whole variant row + if not (het_samples or hom_samples): + return None + + # get repeat counts for each sample with a variant call + vcf_repcn = dict(zip(samples, var.format('REPCN'), strict=True)) + + disease_details = dict(zip(samples, var.format('DISEASE_DETAILS'), strict=True)) + + # dictionary to record repeat counts for called variants - tuple of length 1 | 2 + repeat_counts: dict[str, tuple[int, int] | tuple[int]] = {} + sample_disease_details: dict = {} + for sample_id in het_samples | hom_samples: + # strip out numpy.nan(float64) values, cast the real values as floats + repeat_counts[sample_id] = tuple(int(value) for value in vcf_repcn[sample_id] if ~isnan(value)) # type: ignore[assignment] + sample_disease_details[sample_id] = parse_str_disease_details(disease_details[sample_id]) + + return ShortTandemRepeat( + coordinates=coordinates, + locus=var.ID, + info=info, + het_samples=het_samples, + hom_samples=hom_samples, + boolean_categories=boolean_categories, + ignored_categories=ignored_categories, + support_categories=get_categories_translated(), + sample_repeats=repeat_counts, + sample_repeat_details=sample_disease_details, + ) + + def create_small_variant( var: 'cyvcf2.Variant', samples: list[str], @@ -454,10 +568,8 @@ def create_small_variant( ), ) - # optionally - ignore some categories from this analysis - # keep it super flexible, lower case and a translated version - ignored_categories = {word.lower() for word in config_retrieve(['ValidateMOI', 'ignore_categories'], [])} - ignored_categories.update({translate_category(x) for x in ignored_categories}) + # optionally - ignore some categories from this analysis. Super flexible, lower case and a translated version. + ignored_categories = get_categories_translated('ignore_categories') # set the class attributes - skipping over categories we've chosen to ignore boolean_categories = [ @@ -471,10 +583,6 @@ def create_small_variant( if key.startswith('categorysample') and key.replace('categorysample', '') not in ignored_categories ] - # the categories to be treated as support-only for this runtime - make it a set - support_categories = {word.lower() for word in config_retrieve(['ValidateMOI', 'support_categories'], [])} - support_categories.update({translate_category(x) for x in support_categories}) - # overwrite with true booleans for cat in boolean_categories: info[cat] = info.get(cat, 0) == 1 @@ -522,7 +630,7 @@ def create_small_variant( boolean_categories=boolean_categories, sample_categories=sample_categories, ignored_categories=ignored_categories, - support_categories=support_categories, + support_categories=get_categories_translated(), phased=phased, alt_depths=alt_depths, depths=depths, @@ -595,8 +703,8 @@ def canonical_contigs_from_vcf(reader) -> set[str]: def gather_gene_dict_from_contig( contig: str, - variant_source: 'cyvcf2.VCFReader', - sv_source: 'cyvcf2.VCFReader | None' = None, + variant_sources: dict[str, 'cyvcf2.VCFReader'], + panelapp: PanelApp, ) -> GeneDict: """ takes a cyvcf2.VCFReader instance, and a specified chromosome @@ -606,8 +714,8 @@ def gather_gene_dict_from_contig( Args: contig (): contig name from VCF header - variant_source (): the VCF reader instance - sv_source (): an optional list of SV VCFs + variant_sources (): dict mapping each variant type to a VCF reader instance + known types: small, sv, mito (chrM only), str Returns: A lookup in the form @@ -631,36 +739,55 @@ def gather_gene_dict_from_contig( contig_dict = defaultdict(list) # iterate over all variants on this contig and store by unique key - # if contig has no variants, prints an error and returns [] - for variant in variant_source(contig): - if (small_variant := create_small_variant(var=variant, samples=variant_source.samples)) is None: + small_samples = variant_sources['small'].samples + for variant in variant_sources['small'](contig): + if (small_variant := create_small_variant(var=variant, samples=small_samples)) is None: continue if small_variant.coordinates.string_format in blacklist: logger.info(f'Skipping blacklisted variant: {small_variant.coordinates.string_format}') continue - - # update the variant count contig_variants += 1 - - # update the gene index dictionary contig_dict[small_variant.info['gene_id']].append(small_variant) # parse the SV VCF if provided, but not a necessary part of processing - if sv_source: + if variant_sources.get('sv'): + sv_samples = variant_sources['sv'].samples structural_variants = 0 - for variant in sv_source(contig): - # create an abstract SV variant - structural_variant = create_structural_variant(var=variant, samples=sv_source.samples) - - # update the variant count + for variant in variant_sources['sv'](contig): + structural_variant = create_structural_variant(var=variant, samples=sv_samples) structural_variants += 1 - - # update the gene index dictionary contig_dict[structural_variant.info['gene_id']].append(structural_variant) - logger.info(f'Contig {contig} contained {structural_variants} SVs') + # parse STR VCF if provided + if variant_sources.get('str'): + # limit STRs to PanelApp green-evidence accepted repeat disorder genes + repeat_disorder_genes = panelapp.str_genes + str_samples = variant_sources['str'].samples + str_variants = 0 + for variant in variant_sources['str'](contig): + if str_var := create_str_variant(var=variant, samples=str_samples): + # skip any which aren't accepted in PanelApp's repeat disorders panel + if str_var.info['gene_id'] not in repeat_disorder_genes: + continue + + str_variants += 1 + contig_dict[str_var.info['gene_id']].append(str_var) + + logger.info(f'Contig {contig} contained {str_variants} STRs') + + # special segment for mito data - if Mito is provided and chrM is being searched, provide variants + if contig == 'chrM' and variant_sources.get('mito'): + mito_samples = variant_sources['mito'].samples + mito_variants = 0 + for variant in variant_sources['mito'](contig): + if mito_var := create_small_variant(var=variant, samples=mito_samples): + mito_variants += 1 + contig_dict[mito_var.info['gene_id']].append(mito_var) + + logger.info(f'Mito VCF contained {mito_variants} variants') + logger.info(f'Contig {contig} contained {contig_variants} variants, in {len(contig_dict)} genes') return contig_dict diff --git a/src/talos/validate_moi.py b/src/talos/validate_moi.py index 87a68792..c62a7ef0 100644 --- a/src/talos/validate_moi.py +++ b/src/talos/validate_moi.py @@ -279,8 +279,7 @@ def count_families(pedigree: PedigreeParser) -> dict: def prepare_results_shell( results_meta: ResultMeta, - small_samples: set[str], - sv_samples: set[str], + source_samples: dict[str, set[str]], pedigree: PedigreeParser, panelapp: PanelApp, ) -> ResultData: @@ -289,8 +288,7 @@ def prepare_results_shell( Args: results_meta (): metadata for the results - small_samples (): samples in the Small VCF - sv_samples (): samples in the SV VCFs + source_samples (): samples indexed by variant type pedigree (): the Pedigree object, already reduced to samples in the callset panelapp (): dictionary of gene data @@ -328,8 +326,8 @@ def prepare_results_shell( if panel_id in panelapp.metadata }, solved=bool(participant.sample_id in solved_cases or participant.family_id in solved_cases), - present_in_small=participant.sample_id in small_samples, - present_in_sv=participant.sample_id in sv_samples, + present_in_small=participant.sample_id in source_samples['small'], + present_in_sv=participant.sample_id in source_samples.get('sv', set()), ), ) @@ -339,8 +337,9 @@ def prepare_results_shell( def cli_main(): parser = ArgumentParser(description='Startup commands for the MOI testing phase of Talos') parser.add_argument('--labelled_vcf', help='Category-labelled VCF') - parser.add_argument('--labelled_sv', help='Category-labelled SV VCF', default=None) - parser.add_argument('--labelled_mito', help='Category-labelled Mito VCF', default=None) + parser.add_argument('--labelled_sv', help='Optional, Category-labelled SV VCF', default=None) + parser.add_argument('--labelled_mito', help='Optional, Category-labelled Mito VCF', default=None) + parser.add_argument('--str', help='Optional, VCF containing STR data', default=None) parser.add_argument('--output', help='Prefix to write JSON results to', required=True) parser.add_argument('--panelapp', help='QueryPanelApp JSON', required=True) parser.add_argument('--pedigree', help='Path to PED file', required=True) @@ -354,6 +353,7 @@ def cli_main(): pedigree=args.pedigree, labelled_sv=args.labelled_sv, labelled_mito=args.labelled_mito, + str_vcf=args.str, previous=args.previous, ) @@ -365,6 +365,7 @@ def main( pedigree: str, labelled_sv: str | None = None, labelled_mito: str | None = None, + str_vcf: str | None = None, previous: str | None = None, ): """ @@ -378,6 +379,7 @@ def main( labelled_vcf (str): VCF output from Hail Labelling stage labelled_sv (str | None): optional second VCF (SV) labelled_mito (str | None): optional Mitochondrial VCF + str_vcf (str | None): optional STR 'joint-call' VCF output (str): location to write output file panelapp_path (str): location of PanelApp data JSON pedigree (str): location of PED file @@ -398,6 +400,7 @@ def main( # initialise the optional MOI-stage exclusion logger; no-op when disabled in config exclusion_logger = get_exclusion_logger() + # shove everything in a try-except to make sure the global logger is closed down try: panelapp: PanelApp = read_json_from_path( panelapp_path, @@ -413,26 +416,31 @@ def main( result_list: list[ReportVariant] = [] - # collect all sample IDs from each VCF type - small_vcf_samples: set[str] = set() - sv_vcf_samples: set[str] = set() - - # open the small variant VCF using a cyvcf2 reader + # open the small variant VCF using a cyvcf2 reader - mandatory data source vcf_opened = VCFReader(labelled_vcf) - small_vcf_samples.update(set(vcf_opened.samples)) + source_samples: dict[str, set[str]] = {'small': set(vcf_opened.samples)} + source_vcfs: dict[str, VCFReader] = {'small': vcf_opened} # optional SV behaviour - sv_opened = None if labelled_sv: sv_opened = VCFReader(labelled_sv) - sv_vcf_samples = set(sv_opened.samples) - - all_samples = small_vcf_samples.union(sv_vcf_samples) + source_samples['sv'] = set(sv_opened.samples) + source_vcfs['sv'] = sv_opened - mito_opened = None + # optional mito behaviour if labelled_mito: mito_opened = VCFReader(labelled_mito) - all_samples |= set(mito_opened.samples) + source_vcfs['mito'] = mito_opened + source_samples['mito'] = set(mito_opened.samples) + + # optional str behaviour + if str_vcf: + str_opened = VCFReader(str_vcf) + source_vcfs['str'] = str_opened + source_samples['str'] = set(str_opened.samples) + + # collect all unique samples across the various sources + all_samples = set.union(*source_samples.values()) # parse the pedigree from the file ped = PedigreeParser(pedigree) @@ -454,21 +462,10 @@ def main( # assemble {gene: [var1, var2, ..]} contig_dict = gather_gene_dict_from_contig( contig=contig, - variant_source=vcf_opened, - sv_source=sv_opened, - ) - - result_list.extend( - apply_moi_to_variants( - variant_dict=contig_dict, - moi_lookup=moi_lookup, - panelapp_data=panelapp.genes, - pedigree=ped, - ), + variant_sources=source_vcfs, + panelapp=panelapp, ) - if mito_opened: - contig_dict = gather_gene_dict_from_contig('chrM', variant_source=mito_opened) result_list.extend( apply_moi_to_variants( variant_dict=contig_dict, @@ -488,8 +485,7 @@ def main( # create a shell to store results in, adds participant metadata results_model = prepare_results_shell( results_meta=results_meta, - small_samples=small_vcf_samples, - sv_samples=sv_vcf_samples, + source_samples=source_samples, pedigree=ped, panelapp=panelapp, ) diff --git a/src/talos/version.py b/src/talos/version.py index d1bb7002..77f1fe77 100644 --- a/src/talos/version.py +++ b/src/talos/version.py @@ -1 +1 @@ -__version__ = '11.0.3' +__version__ = '11.1.0' diff --git a/test/test_unified_panelapp_parser.py b/test/test_unified_panelapp_parser.py index a6c7cd61..d2a3e4de 100644 --- a/test/test_unified_panelapp_parser.py +++ b/test/test_unified_panelapp_parser.py @@ -101,7 +101,6 @@ def test_match_participants_to_panels(): 'genes': {}, 'participants': { 'sam1': ParticipantHPOPanels( - external_id='', family_id='', hpo_terms=[HpoTerm(id='HP:1', label='1'), HpoTerm(id='HP:2', label='2')], panels={99, 1, 3, 5}, @@ -109,7 +108,6 @@ def test_match_participants_to_panels(): matched_phenotypes=set(), ), 'sam2': ParticipantHPOPanels( - external_id='', family_id='', hpo_terms=[HpoTerm(id='HP:1', label='1')], panels={99, 1, 3}, @@ -118,6 +116,8 @@ def test_match_participants_to_panels(): ), }, 'version': CURRENT_VERSION, + 'str_genes': set(), + 'str_symbols': set(), } diff --git a/test/test_utils.py b/test/test_utils.py index 74ab45a4..3caccb53 100644 --- a/test/test_utils.py +++ b/test/test_utils.py @@ -26,6 +26,11 @@ FIVE_EXPECTED = 5 +class FakePanelApp: + def __init__(self): + self.str_genes: set[str] = set() + + def test_coord_sorting(): """ check that coord sorting methods work @@ -132,8 +137,7 @@ def test_gene_dict(two_trio_variants_vcf): gene = ENSG00000075043 """ reader = VCFReader(two_trio_variants_vcf) - contig = 'chr20' - var_dict = gather_gene_dict_from_contig(contig=contig, variant_source=reader) + var_dict = gather_gene_dict_from_contig(contig='chr20', variant_sources={'small': reader}, panelapp=FakePanelApp()) assert len(var_dict) == 1 assert 'ENSG00000075043' in var_dict assert len(var_dict['ENSG00000075043']) == TWO_EXPECTED @@ -163,7 +167,7 @@ def test_phased_dict(phased_vcf_path): gene = ENSG00000075043 """ reader = VCFReader(phased_vcf_path) - var_dict = gather_gene_dict_from_contig(contig='chr20', variant_source=reader) + var_dict = gather_gene_dict_from_contig(contig='chr20', variant_sources={'small': reader}, panelapp=FakePanelApp()) assert len(var_dict) == ONE_EXPECTED assert 'ENSG00000075043' in var_dict assert len(var_dict['ENSG00000075043']) == TWO_EXPECTED diff --git a/test/test_validate_methods.py b/test/test_validate_methods.py index c13cb224..d4887ebc 100644 --- a/test/test_validate_methods.py +++ b/test/test_validate_methods.py @@ -49,10 +49,10 @@ def test_results_shell(pedigree_path: str): 'mother_2': {'external_id': 'mother_2', 'hpo_terms': []}, }, ) + source_samples = {'small': {'male'}, 'sv': {'female'}} shell = prepare_results_shell( results_meta=ResultMeta(), - small_samples={'male'}, - sv_samples={'female'}, + source_samples=source_samples, pedigree=PedigreeParser(pedigree_path=pedigree_path), panelapp=panelapp, ) diff --git a/uv.lock b/uv.lock index 81e1331e..f335b60e 100644 --- a/uv.lock +++ b/uv.lock @@ -591,7 +591,7 @@ wheels = [ [[package]] name = "cpg-utils" -version = "5.4.2" +version = "5.7.0" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "boto3" }, @@ -600,13 +600,15 @@ dependencies = [ { name = "deprecated" }, { name = "frozendict" }, { name = "google-auth" }, + { name = "google-auth-oauthlib" }, { name = "google-cloud-artifact-registry" }, { name = "google-cloud-secret-manager" }, + { name = "metamist" }, { name = "requests" }, { name = "tabulate" }, { name = "toml" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/f1/55/e069ef106f576e3e7916f4801d56dc2efd867387ff4fc0b6593815a54f47/cpg_utils-5.4.2.tar.gz", hash = "sha256:eb6f0808655e836171db279b44db3bff66b4e7a9b3a8629caf630483e7cc5bc6", size = 43646, upload-time = "2025-08-13T01:56:47.044Z" } +sdist = { url = "https://files.pythonhosted.org/packages/01/9e/145494610bd318645f834a0d8f7e7b0f819f77337b334254bb9ccfa89bd3/cpg_utils-5.7.0.tar.gz", hash = "sha256:aba177361a73003836de4a8c310dadccf6e7bff16760bf51bbb263a133d2d421", size = 51003, upload-time = "2026-06-30T21:32:18.12Z" } [[package]] name = "cryptography" @@ -2320,6 +2322,26 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/5a/dc/491b7661614ab97483abf2056be1deee4dc2490ecbf7bff9ab5cdbac86e1/pyreadline3-3.5.4-py3-none-any.whl", hash = "sha256:eaf8e6cc3c49bcccf145fc6067ba8643d1df34d604a1ec0eccbf7a18e6d3fae6", size = 83178, upload-time = "2024-09-19T02:40:08.598Z" }, ] +[[package]] +name = "pysam" +version = "0.24.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/09/9b/e856ad525ed0847810cca936ee551b8e0a709525fea34e0e0e723173f62d/pysam-0.24.0.tar.gz", hash = "sha256:db0f86c15532ef5dad263748324f45d9a639668e3497d8cabce54ef47a1a78d9", size = 5201690, upload-time = "2026-04-27T12:26:50.879Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/25/ed/5b642a89de9a4cf5ce21f70b4c2c2cd10720ef12cb6ef60b005ecf9dcb08/pysam-0.24.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:d06a94e2886cb8affba77bc4fe337513f29387113d708d14609765aa8bb143d3", size = 6217780, upload-time = "2026-04-27T12:25:00.598Z" }, + { url = "https://files.pythonhosted.org/packages/c1/39/2c9aea082383c8d41d0d8b8f9907dc6273f8cee47a8c83d0c8bad02488b4/pysam-0.24.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:4d3260795e066f659ea94e6b7cac07727ef2877b874f69029b6a2e48551d72bd", size = 5946970, upload-time = "2026-04-27T12:25:02.783Z" }, + { url = "https://files.pythonhosted.org/packages/12/0b/d9c1155257f2c2d0d8a33d15011964bf0c9ad927f85fa57f00965b0c1cfc/pysam-0.24.0-cp310-cp310-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:9b299f6128601356e25ca94040eef9dadb1e196aa2f36736c93948aaa821c4bc", size = 22114024, upload-time = "2026-04-27T12:25:04.779Z" }, + { url = "https://files.pythonhosted.org/packages/21/48/92dddc8df65b576c9d30752650c89301b5222d4ac10187724796cedfd723/pysam-0.24.0-cp310-cp310-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:98a40440b16bbf67b696914a10acc241a0cd0fa02dcb14aee8e614eb2c5d8b37", size = 22675883, upload-time = "2026-04-27T12:25:07.582Z" }, + { url = "https://files.pythonhosted.org/packages/07/ff/eee9caac42c255150cc7eae19afa71bd2bb8347fd9895196f7daae2ccd62/pysam-0.24.0-cp310-cp310-musllinux_1_2_aarch64.whl", hash = "sha256:8e4f5f1ebf3bdcad4c3b8f3a656ec2afdc44f0efc96ce8523858556e95d8ca9a", size = 22056750, upload-time = "2026-04-27T12:25:10.127Z" }, + { url = "https://files.pythonhosted.org/packages/d7/7f/da483a98b11690f3c4b24f4aa12b7e0d0d1e19da8b8ab45785c129240eae/pysam-0.24.0-cp310-cp310-musllinux_1_2_x86_64.whl", hash = "sha256:fcce1edc48a9d8610958296c0aaf3e28723d143113bcd9b7c13b5cd821300d2a", size = 22364034, upload-time = "2026-04-27T12:25:12.709Z" }, + { url = "https://files.pythonhosted.org/packages/68/b3/1d586b4db6ff7220ac1df4b9588fae7c2622680ca51d53aa7cf32ee73604/pysam-0.24.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:48524b2132158ffeca7ca21c52490e92dca6687446faad95522db2354e376641", size = 8711987, upload-time = "2026-04-27T12:25:15.415Z" }, + { url = "https://files.pythonhosted.org/packages/6d/29/5fe6f53416cfc17cbe1440f5f55b655e5dcb64b67e18cae9b0a38366db0b/pysam-0.24.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:58e347406cec254bd68ac0f50da1e772273c3c891f65eab4b54641f2a89e8644", size = 8441919, upload-time = "2026-04-27T12:25:17.984Z" }, + { url = "https://files.pythonhosted.org/packages/1a/75/69379edfcd26c88f749e29b5d3aee9044fbd201c1a678371e7ccd9010937/pysam-0.24.0-cp311-cp311-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:ed1f202fe27f2bd677b3d210616890a8cbf043043ea89404ada3c121a5e9a20e", size = 25122944, upload-time = "2026-04-27T12:25:20.982Z" }, + { url = "https://files.pythonhosted.org/packages/13/6d/2e2a80cda66b152ec143ceb1784e56d8abc6abeeeb65f8ef6aa52c731038/pysam-0.24.0-cp311-cp311-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:8de0b19109da0c4aa8e44c6fef2fae97f0e616f3cdf5be0f1b74f7946b8de262", size = 25697041, upload-time = "2026-04-27T12:25:24.566Z" }, + { url = "https://files.pythonhosted.org/packages/91/70/852f6917e9d7592482db2c43ccd91fa57158cf20ba07c31e9a477af84b9d/pysam-0.24.0-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:d528fc09a44c85b7d4018ba8be543e1cc2e7a4d06eeb779a28ac34947e182f5a", size = 25069351, upload-time = "2026-04-27T12:25:27.617Z" }, + { url = "https://files.pythonhosted.org/packages/b0/6b/0596765a59253958488fc7a8816df347d815ac3850a7a9cab785c1b1e184/pysam-0.24.0-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:1e9f20388c80dd0e4c09dfd7fe38adc9fdf771ec665e071e1b7350dc8da2d692", size = 25378507, upload-time = "2026-04-27T12:25:30.504Z" }, +] + [[package]] name = "pyspark" version = "3.5.3" @@ -2728,7 +2750,7 @@ wheels = [ [[package]] name = "talos" -version = "11.0.0" +version = "11.1.0" source = { editable = "." } dependencies = [ { name = "clinvarbitration" }, @@ -2762,6 +2784,9 @@ docs = [ { name = "mkdocs-include-markdown-plugin" }, { name = "mkdocs-material" }, ] +stripy = [ + { name = "pysam" }, +] test = [ { name = "bump-my-version" }, { name = "pre-commit" }, @@ -2776,7 +2801,7 @@ requires-dist = [ { name = "clinvarbitration", specifier = "~=2.2.7" }, { name = "cloudpathlib", extras = ["all"], specifier = ">=0.16.0" }, { name = "cpg-flow", marker = "extra == 'cpg'", specifier = "~=1.3" }, - { name = "cpg-utils", specifier = ">=5.4.1" }, + { name = "cpg-utils", specifier = "~=5.7" }, { name = "cyvcf2", specifier = ">=0.30.18" }, { name = "google-cloud-storage", marker = "extra == 'cpg'", specifier = "==2.14.0" }, { name = "grpcio-status", specifier = ">=1.48,<1.50" }, @@ -2792,6 +2817,7 @@ requires-dist = [ { name = "pendulum", specifier = ">=3.1.0" }, { name = "pre-commit", marker = "extra == 'test'" }, { name = "pydantic", specifier = ">=2.5.2" }, + { name = "pysam", marker = "extra == 'stripy'" }, { name = "pytest", marker = "extra == 'test'" }, { name = "pytest-httpx", marker = "extra == 'test'" }, { name = "pytest-xdist", marker = "extra == 'test'", specifier = ">=3.6.0" }, @@ -2801,7 +2827,7 @@ requires-dist = [ { name = "tenacity", specifier = ">=9.0.0" }, { name = "toml", specifier = "==0.10.2" }, ] -provides-extras = ["cpg", "test", "docs"] +provides-extras = ["cpg", "test", "docs", "stripy"] [[package]] name = "tenacity"