Export and re-import managed variants with bcftools liftover

CHANGELOG.md

@@ -27,6 +27,7 @@ About changelog [here](https://keepachangelog.com/en/1.0.0/)
 - An icon to the Gens-to-cases link to distinguish it from Gens-to-individuals links on the SV variant page (#6410)
 - Display both GnomAD and GnomAD non-UK-Biobank links on variants pages for build 38 SNVs (#6413)
 - Variant dismiss tags shown on pages can be customised in the institute settings page by admin users (#6416)
+- A `--liftover-from` option to the `export managed` command line, so managed variants can be exported and imported again in another genome build (#6225)
 ### Changed
 - Replaced Ensembl rest liftover service with liftover API from the Broad Institute (#6293)
 - Avoid fetching genes and panels multiple times when loading variants (#6350)

scout/commands/export/variant.py

@@ -8,10 +8,11 @@
 from flask.cli import with_appcontext
 from xlsxwriter import Workbook
 
-from scout.constants import CALLERS, DATE_DAY_FORMATTER
+from scout.constants import BUILDS, CALLERS, DATE_DAY_FORMATTER
 from scout.constants.managed_variant import MANAGED_CATEGORIES
 from scout.constants.variants_export import VERIFIED_VARIANTS_HEADER
 from scout.export.variant import (
+    export_lift_over_managed_variants,
     export_managed_variants,
     export_verified_variants,
 )
@@ -124,8 +125,15 @@ def verified(collaborator, test, outpath=None):
 @collaborator_option
 @build_option
 @json_option
+@click.option(
+    "--liftover-from",
+    type=click.Choice(BUILDS),
+    help="Perform liftover on coordinates and export as managed variants infile.",
+)
 @with_appcontext
-def managed(collaborator: str, category: Tuple[str], build: str, json: bool):
+def managed(
+    collaborator: str, category: Tuple[str], build: str, json: bool, liftover_from: Optional[str]
+):
     """Export managed variants for a collaborator in VCF or JSON format"""
     LOG.info("Running scout export managed variants")
     adapter = store
@@ -141,6 +149,9 @@ def managed(collaborator: str, category: Tuple[str], build: str, json: bool):
     if json:
         click.echo(json_lib.dumps([var for var in variants], default=bson_handler))
         return
+    if liftover_from:
+        export_lift_over_managed_variants(managed_variants=variants, liftover_from=liftover_from)
+        return
 
     print_vcf(variants=variants, build=build, export_category="MANAGED")
 

scout/export/variant.py

@@ -1,13 +1,16 @@
 # -*- coding: utf-8 -*-
 import logging
 import urllib.parse
-from typing import List
+from typing import Iterable, List
+
+import click
 
 from scout.adapter.mongo.base import MongoAdapter
 from scout.constants import CHROMOSOME_INTEGERS
-from scout.constants.managed_variant import MANAGED_CATEGORIES
+from scout.constants.managed_variant import MANAGED_CATEGORIES, MANAGED_VARIANTS_INFILE_HEADER
 from scout.constants.query_terms import GT_NO_ALT_CALL
 from scout.models.managed_variant import ManagedVariant
+from scout.utils.broad_liftover_client import BroadLiftoverApiClient
 
 LOG = logging.getLogger(__name__)
 
@@ -24,6 +27,74 @@
     return sorted(variants, key=sort_key)
 
 
+def export_lift_over_managed_variants(managed_variants: Iterable, liftover_from: str):
+    """Perform liftover over a list of managed variants and print a list of lines formatted as a managed variants upload infile."""
+
+    export_lines = [MANAGED_VARIANTS_INFILE_HEADER]
+    client = BroadLiftoverApiClient()
+
+    lifted_build = "38" if liftover_from == "37" else "37"
+
+    nfailed = 0
+    nprocessed = 0
+
+    for nprocessed, variant_obj in enumerate(managed_variants, 1):
+        if nprocessed % 50 == 0:
+            LOG.info(f"Processed {nprocessed} variants")
+
+        category = variant_obj.get("category", "snv")
+        if category not in ["snv", "cancer_snv"]:
+            continue
+
+        build = variant_obj.get("build")
+
+        if build == lifted_build:
+            chrom = variant_obj["chromosome"]
+            pos = variant_obj["position"]
+            end = variant_obj.get("end", pos)
+            ref = variant_obj.get("reference")
+            alt = variant_obj.get("alternative")
+        else:
+            result = client.liftover(
+                build_from=liftover_from,
+                chrom=variant_obj.get("chromosome"),
+                start=variant_obj.get("position"),
+                end=variant_obj.get("end"),
+                ref=variant_obj.get("reference", ""),
+                alt=variant_obj.get("alternative", ""),
+            )
+
+            if "error" in result:
+                nfailed += 1
+                LOG.error(result)
+                continue
+
+            chrom = result["output_chrom"].replace("chr", "")
+            pos = result["output_pos"]
+            end = result.get("output_end") or result.get("output_pos")
+            ref = result["output_ref"]
+            alt = result["output_alt"]
+
+        sub_category = variant_obj.get("sub_category", "snv")
+
+        desc = variant_obj.get("description")
+        if "(causatives" not in (desc or ""):
+            description = f"{desc} (managed, build{liftover_from})"
+        else:
+            description = desc
+
+        institutes = ",".join(variant_obj.get("institute") or [])
+
+        export_lines.append(
+            f"{chrom};{pos};{end};{ref};{alt};"
+            f"{category};{sub_category};{lifted_build};{description};;{institutes}"
+        )
+
+    LOG.info(f"Done. Total processed: {nprocessed} - total failed: {nfailed}")
+    for line in export_lines:
+        click.echo(line)
+
+
 def export_managed_variants(
     adapter: MongoAdapter,
     institute: str = None,

tests/export/test_export_variants.py

@@ -1,7 +1,9 @@
-# -*- coding: utf-8 -*-
+import responses
 
+from scout.constants.managed_variant import MANAGED_VARIANTS_INFILE_HEADER
 from scout.constants.variants_export import MT_EXPORT_HEADER
-from scout.export.variant import export_mt_variants
+from scout.export.variant import export_lift_over_managed_variants, export_mt_variants
+from scout.utils.broad_liftover_client import LIFTOVER_URL
 
 
 def test_export_mt_variants(case_obj, real_populated_database):
@@ -26,7 +28,60 @@ def test_export_mt_variants(case_obj, real_populated_database):
     for sample in samples:
         sample_lines = export_mt_variants(variants=mt_variants, sample_id=sample["individual_id"])
 
-        # check that rows to write to excel corespond to number of variants
+        # check that rows to write to excel correspond to number of variants
         assert len(sample_lines) == len(mt_variants)
-        # check that cols to write to excel corespond to fields of excel header
+        # check that cols to write to excel correspond to fields of Excel header
         assert len(sample_lines[0]) == len(MT_EXPORT_HEADER)
+
+
+@responses.activate
+def test_export_lift_over_managed_variants(broad_bcftools_liftover_response, capsys):
+    """Test the function lifts over and formats managed variants into a managed variants infile."""
+
+    # GIVEN an Iterable with managed variants:
+    managed_variant = {
+        "chromosome": "8",
+        "position": 141310715,
+        "end": 141310715,
+        "reference": "T",
+        "alternative": "G",
+        "category": "snv",
+        "build": "37",
+    }
+
+    managed_variants = [managed_variant]
+
+    url = (
+        f"{LIFTOVER_URL}/"
+        "?hg=hg19-to-hg38"
+        "&format=variant"
+        f"&chrom={managed_variant['chromosome']}"
+        f"&pos={managed_variant['position']}"
+        f"&end={managed_variant['end']}"
+        f"&ref={managed_variant['reference']}"
+        f"&alt={managed_variant['alternative']}"
+    )
+
+    # GIVEN a mocked call to the liftover service (BCFTools via Broad institute's API)
+    resp = broad_bcftools_liftover_response
+    responses.add(
+        responses.GET,
+        url,
+        json=resp,
+        status=200,
+    )
+
+    # WHEN exporting the managed variants
+    export_lift_over_managed_variants(
+        managed_variants=managed_variants,
+        liftover_from="37",
+    )
+
+    # THEN a header and a lifted variant line should be printed
+    out = capsys.readouterr().out.splitlines()
+
+    assert out[0] == MANAGED_VARIANTS_INFILE_HEADER
+    assert (
+        out[1]
+        == f"{resp['chrom'].replace('chr','')};{resp['output_pos']};{resp['output_pos']};{resp['output_ref']};{resp['output_alt']};snv;snv;38;None (managed, build37);;"
+    )

tests/server/blueprints/managed_variants/test_managed_variants_views.py

@@ -66,7 +66,7 @@ def test_add_and_remove_managed_variants(app, mocker, mock_redirect):
         resp = client.post(url_for("managed_variants.add_managed_variant"), data=add_form_data)
         # THEN the status code should still be redirect
         assert resp.status_code == 302
-        # THEN the database should still contatain only one variant
+        # THEN the database should still contain only one variant
         assert sum(1 for i in store.managed_variant_collection.find()) == 1
 
         # WHEN requesting to remove the selected variant