aboutsummaryrefslogtreecommitdiff
path: root/gn3
diff options
context:
space:
mode:
Diffstat (limited to 'gn3')
-rw-r--r--gn3/api/search.py246
1 files changed, 207 insertions, 39 deletions
diff --git a/gn3/api/search.py b/gn3/api/search.py
index 83951a1..650208b 100644
--- a/gn3/api/search.py
+++ b/gn3/api/search.py
@@ -1,13 +1,16 @@
"""Search using Xapian index."""
from collections import namedtuple
+import gzip
import json
from functools import partial, reduce
+from pathlib import Path
from typing import Callable
import urllib.parse
from flask import abort, Blueprint, current_app, jsonify, request
from pymonad.maybe import Just, Maybe, Nothing
+from pymonad.tools import curry
import xapian
from gn3.monads import MonadicDict
@@ -17,6 +20,7 @@ search = Blueprint("search", __name__)
ChromosomalPosition = namedtuple("ChromosomalPosition", "chromosome position")
ChromosomalInterval = namedtuple("ChromosomalInterval", "chromosome start end")
+IntervalLiftoverFunction = Callable[[ChromosomalInterval], Maybe[ChromosomalInterval]]
FieldProcessorFunction = Callable[[str], xapian.Query]
@@ -35,6 +39,28 @@ def combine_queries(operator: int, *queries: xapian.Query) -> xapian.Query:
return reduce(partial(xapian.Query, operator), queries)
+def query_subqueries(query: xapian.Query) -> list[xapian.Query]:
+ """Return list of child queries in query."""
+ return [query.get_subquery(i) for i in range(query.get_num_subqueries())]
+
+
+def query_terms(query: xapian.Query) -> list[str]:
+ """Return list of terms in query."""
+ # Unfortunately, the TermIterator from python xapian bindings seems
+ # buggy. So, we resort to traversing the query tree.
+ # python xapian bindings do not expose xapian.Query.LEAF_TERM. This is
+ # most likely a bug.
+ leaf_type = 100
+ if query.get_type() == leaf_type:
+ # We have no choice but to access the protected _get_terms_begin method.
+ # pylint: disable=protected-access
+ return [query._get_terms_begin().get_term().decode("utf-8")]
+ else:
+ return reduce(lambda result, subquery: result + query_terms(subquery),
+ query_subqueries(query),
+ [])
+
+
class FieldProcessor(xapian.FieldProcessor):
"""
Field processor for use in a xapian query parser.
@@ -52,6 +78,68 @@ class FieldProcessor(xapian.FieldProcessor):
return self.proc(query)
+def field_processor_or(*field_processors: FieldProcessorFunction) -> FieldProcessorFunction:
+ """Combine field processors using the OR operator."""
+ return (lambda query:
+ combine_queries(xapian.Query.OP_OR,
+ *[field_processor(query)
+ for field_processor in field_processors]))
+
+
+def liftover(chain_file: Path, position: ChromosomalPosition) -> Maybe[ChromosomalPosition]:
+ """Liftover chromosomal position using chain file."""
+ # The chain file format is described at
+ # https://genome.ucsc.edu/goldenPath/help/chain.html
+ def split_chain_line(line: str) -> tuple[ChromosomalInterval, ChromosomalInterval, int, str]:
+ (_, _, target_chromosome, _, _, target_start, target_end,
+ query_chromosome, query_size, query_strand, query_start, query_end, _) = line.split()
+ return (ChromosomalInterval(target_chromosome.removeprefix("chr"),
+ int(target_start), int(target_end)),
+ ChromosomalInterval(query_chromosome.removeprefix("chr"),
+ int(query_start), int(query_end)),
+ int(query_size), query_strand)
+
+ with gzip.open(chain_file, "rt") as file:
+ for line in file:
+ if line.startswith('chain'):
+ (target, query, query_size, query_strand) = split_chain_line(line)
+ if (target.chromosome == position.chromosome
+ and target.start <= position.position < target.end):
+ transformed_position = query.start + position.position - target.start
+ if query_strand == '-':
+ transformed_position = query_size - 1 - transformed_position
+ return Just(ChromosomalPosition(query.chromosome, transformed_position))
+ return Nothing
+
+
+def liftover_interval(chain_file: str, interval: ChromosomalInterval) -> ChromosomalInterval:
+ """
+ Liftover interval using chain file.
+
+ This function lifts over the interval by merely lifting the start and end
+ points. This is simplistic and not strictly correct, but will do for the
+ purposes of synteny in search.
+ """
+ point_liftover = partial(liftover, chain_file)
+ return (Maybe.apply(curry(2,
+ lambda start, end:
+ ChromosomalInterval(start.chromosome,
+ Just(start.position),
+ Just(end.position))))
+ .to_arguments(interval_start(interval).bind(point_liftover),
+ interval_end(interval).bind(point_liftover))
+ # In '-' strands, the order of the interval may be reversed. Work
+ # around.
+ # KLUDGE: Inspecting the value contained in a monad is usually bad
+ # practice. But, in this case, there doesn't seem to be a way
+ # around.
+ .map(lambda interval: ChromosomalInterval(interval.chromosome,
+ Just(min(interval.start.value, interval.end.value)),
+ Just(max(interval.start.value, interval.end.value)))
+ if interval.start.is_just() and interval.end.is_just()
+ else interval))
+
+
def parse_range(range_string: str) -> tuple[Maybe[str], Maybe[str]]:
"""Parse xapian range strings such as start..end."""
start, end = range_string.split("..")
@@ -65,8 +153,9 @@ def apply_si_suffix(location: str) -> int:
return int(float(location[:-1])*10**suffixes.get(location[-1].lower(), 0))
-def parse_location_field(species: str, species_prefix: str,
+def parse_location_field(lifted_species: str, species_prefix: str,
chromosome_prefix: str, location_slot: int,
+ liftover_function: IntervalLiftoverFunction,
query: bytes) -> xapian.Query:
"""Parse location shorthands and return a xapian query.
@@ -82,51 +171,130 @@ def parse_location_field(species: str, species_prefix: str,
*[location.map(apply_si_suffix)
for location in parse_range(location)])
+ def make_query(species: str, interval: ChromosomalInterval) -> xapian.Query:
+ # TODO: Convert the xapian index to use bases instead of megabases.
+ to_megabases = lambda x: str(float(x)/1e6)
+ return combine_queries(xapian.Query.OP_AND,
+ xapian.Query(species_prefix + species),
+ xapian.Query(chromosome_prefix + interval.chromosome),
+ xapian.NumberRangeProcessor(location_slot)
+ (interval.start.maybe("", to_megabases),
+ interval.end.maybe("", to_megabases)))
try:
interval = split_query(query.decode("utf-8"))
except ValueError:
return xapian.Query(xapian.Query.OP_INVALID)
- return combine_queries(xapian.Query.OP_AND,
- xapian.Query(species_prefix + species),
- xapian.Query(chromosome_prefix + interval.chromosome),
- xapian.NumberRangeProcessor(location_slot)
- (interval.start.maybe("", str),
- interval.end.maybe("", str)))
+ return (liftover_function(interval)
+ .maybe(xapian.Query.MatchNothing,
+ partial(make_query, lifted_species)))
+
+
+def parse_synteny_field(synteny_prefix: str, query: bytes) -> xapian.Query:
+ """Parse synteny field and return a xapian query."""
+ if query.decode("utf-8") in ["on", "off"]:
+ return xapian.Query(synteny_prefix + query.decode("utf-8"))
+ else:
+ return xapian.Query(xapian.Query.OP_INVALID)
+
+def is_synteny_on(synteny_prefix: str, query: xapian.Query) -> bool:
+ """Check if synteny search is requested in query."""
+ return synteny_prefix + "on" in query_terms(query)
-def parse_query(query: str):
+
+def remove_synteny_field(synteny_prefix: str, query: xapian.Query,
+ parent_operator: int = xapian.Query.OP_AND) -> xapian.Query:
+ """Return a new query with the synteny field removed."""
+ # Note that this function only supports queries that exclusively use the
+ # AND, OR, FILTER, RANGE and INVALID operators.
+ # python xapian bindings do not expose xapian.Query.LEAF_TERM. This is
+ # most likely a bug.
+ leaf_type = 100
+ # Handle leaf node.
+ if query.get_type() == leaf_type:
+ if not any(term.startswith(synteny_prefix) for term in query_terms(query)):
+ return query
+ elif parent_operator in [xapian.Query.OP_AND, xapian.Query.OP_FILTER]:
+ return xapian.Query.MatchAll
+ elif parent_operator == xapian.Query.OP_OR:
+ return xapian.Query.MatchNothing
+ else:
+ raise ValueError("Unexpected operator in query", query.get_type())
+ # Recurse on non-leaf nodes with the AND, OR or FILTER operators as root.
+ elif query.get_type() in (xapian.Query.OP_AND, xapian.Query.OP_OR, xapian.Query.OP_FILTER):
+ return combine_queries(query.get_type(),
+ *[remove_synteny_field(synteny_prefix, subquery, query.get_type())
+ for subquery in query_subqueries(query)])
+ # Return other supported non-leaf nodes verbatim.
+ elif query.get_type() in [xapian.Query.OP_VALUE_RANGE, xapian.Query.OP_INVALID]:
+ return query
+ # Raise an exception on unsupported non-leaf nodes.
+ else:
+ raise ValueError("Unexpected operator in query", query.get_type())
+
+
+def parse_query(synteny_files_directory: Path, query: str):
"""Parse search query using GeneNetwork specific field processors."""
- queryparser = xapian.QueryParser()
- queryparser.set_stemmer(xapian.Stem("en"))
- queryparser.set_stemming_strategy(queryparser.STEM_SOME)
- species_prefix = "XS"
- chromosome_prefix = "XC"
- queryparser.add_boolean_prefix("author", "A")
- queryparser.add_boolean_prefix("species", species_prefix)
- queryparser.add_boolean_prefix("group", "XG")
- queryparser.add_boolean_prefix("tissue", "XI")
- queryparser.add_boolean_prefix("dataset", "XDS")
- queryparser.add_boolean_prefix("symbol", "XY")
- queryparser.add_boolean_prefix("chr", chromosome_prefix)
- queryparser.add_boolean_prefix("peakchr", "XPC")
- queryparser.add_prefix("description", "XD")
- range_prefixes = ["mean", "peak", "mb", "peakmb", "additive", "year"]
- for i, prefix in enumerate(range_prefixes):
- queryparser.add_rangeprocessor(xapian.NumberRangeProcessor(i, prefix + ":"))
-
- # Add field processors for location shorthands.
- species_shorthands = {"Hs": "human",
- "Mm": "mouse",
- "Rn": "rat"}
- for shorthand, species in species_shorthands.items():
- queryparser.add_boolean_prefix(
- shorthand, FieldProcessor(partial(parse_location_field,
- species,
- species_prefix,
- chromosome_prefix,
- range_prefixes.index("mb"))))
- return queryparser.parse_query(query)
+ synteny_prefix = "XSYN"
+
+ def make_query_parser(synteny: bool) -> xapian.QueryParser:
+ queryparser = xapian.QueryParser()
+ queryparser.set_stemmer(xapian.Stem("en"))
+ queryparser.set_stemming_strategy(queryparser.STEM_SOME)
+ species_prefix = "XS"
+ chromosome_prefix = "XC"
+ queryparser.add_boolean_prefix("author", "A")
+ queryparser.add_boolean_prefix("species", species_prefix)
+ queryparser.add_boolean_prefix("group", "XG")
+ queryparser.add_boolean_prefix("tissue", "XI")
+ queryparser.add_boolean_prefix("dataset", "XDS")
+ queryparser.add_boolean_prefix("symbol", "XY")
+ queryparser.add_boolean_prefix("chr", chromosome_prefix)
+ queryparser.add_boolean_prefix("peakchr", "XPC")
+ queryparser.add_prefix("description", "XD")
+ queryparser.add_prefix("synteny", FieldProcessor(partial(parse_synteny_field,
+ synteny_prefix)))
+ range_prefixes = ["mean", "peak", "mb", "peakmb", "additive", "year"]
+ for i, prefix in enumerate(range_prefixes):
+ queryparser.add_rangeprocessor(xapian.NumberRangeProcessor(i, prefix + ":"))
+
+ # Add field processors for location shorthands.
+ species_shorthands = {"Hs": "human",
+ "Mm": "mouse"}
+ for shorthand, species in species_shorthands.items():
+ field_processors = [partial(parse_location_field,
+ species,
+ species_prefix,
+ chromosome_prefix,
+ range_prefixes.index("mb"),
+ Just)]
+ # If human and synteny is requested, add liftover.
+ # With synteny search, we search for the same gene sequences
+ # across different species. But, the same gene sequences may be
+ # present in very different chromosomal positions in different
+ # species. So, we first liftover.
+ if shorthand == "Hs" and synteny:
+ chain_files = {"mouse": "hg19ToMm10-chains.over.chain.gz"}
+ for lifted_species, chain_file in chain_files.items():
+ field_processors.append(
+ partial(parse_location_field,
+ lifted_species,
+ species_prefix,
+ chromosome_prefix,
+ range_prefixes.index("mb"),
+ partial(liftover_interval,
+ synteny_files_directory / chain_file)))
+ queryparser.add_boolean_prefix(
+ shorthand,
+ FieldProcessor(field_processor_or(*field_processors)))
+ return queryparser
+
+ return remove_synteny_field(
+ synteny_prefix,
+ make_query_parser(is_synteny_on(synteny_prefix,
+ make_query_parser(False).parse_query(query)))
+ .parse_query(query))
@search.route("/")
@@ -143,7 +311,7 @@ def search_results():
if results_per_page > maximum_results_per_page:
abort(400, description="Requested too many search results")
- query = parse_query(querystring)
+ query = parse_query(Path(current_app.config["DATA_DIR"]) / "synteny", querystring)
traits = []
# pylint: disable=invalid-name
with xapian_database(current_app.config["XAPIAN_DB_PATH"]) as db: