Merge: Fixes issues when ECTyper doesn't find all alleles for Phylotyper

app/modules/phylotyper/ontology.py

@@ -21,7 +21,11 @@
 typing_ontology_version = '<https://www.github.com/superphy/typing/1.0.0>'
 __location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__)))
 
-
+LOCI = {
+    'stx1': [':stx1A', ':stx1B'],
+    'stx2': [':stx2A', ':stx2B'],
+    'eae': [':eae']
+}
 
 @submit
 @prefix
@@ -158,9 +162,9 @@ def generate_graph(uri, loci, values):
     subtype = uri.split(':')[1]
 
     # Check for existance of schema Marker components
-    for l in loci:
-        if not find_object(l, ':Marker'):
-            raise DatabaseError(uri, l)
+    # for l in loci:
+    #     if not find_object(l, ':Marker'):
+    #         raise DatabaseError(uri, l)
 
     # Proceed with creating subtype schema
     graph = Graph()
@@ -214,7 +218,7 @@ def stx1_graph():
 
     """
 
-    return generate_graph('subt:stx1', [':stx1A',':stx1B'], ['a','c','d','untypeable'])
+    return generate_graph('subt:stx1', LOCI['stx1'], ['a','c','d','untypeable'])
 
 
 
@@ -224,7 +228,7 @@ def stx2_graph():
 
     """
 
-    return generate_graph('subt:stx2', [':stx2A',':stx2B'], ['a','b','c','d','e','f','g','untypeable'])
+    return generate_graph('subt:stx2', LOCI['stx2'], ['a','b','c','d','e','f','g','untypeable'])
 
 
 def eae_graph():
@@ -233,7 +237,7 @@ def eae_graph():
 
     """
 
-    return generate_graph('subt:eae', [':eae'],
+    return generate_graph('subt:eae', LOCI['eae'],
         ["alpha-1","alpha-2","beta-1","beta-2","epsilon-1","epsilon-2","eta-1","eta-2",
         "gamma-1","iota-1","iota-2","kappa-1","lambda-1","mu-1","nu-1","omicron-1","pi-1",
         "rho-1","sigma-1","theta-2","xi-1","zeta-1","untypeable"])

app/modules/phylotyper/phylotyper.py

@@ -34,10 +34,14 @@
 redis_url = config.REDIS_URL
 redis_conn = redis.from_url(redis_url)
 
-def _check_tsv(pt_file):
-    pt_results = pd.read_table(pt_file)
-    if pt_results.empty:
-        raise Exception('_check_tsv() failed as pt_results.empty == true for pt_file: {0} with df content: {1}'.format(pt_file, str(pt_results)))
+def _noloci(output_file):
+    """Writes out no loci found.
+    """
+    with open(output_file, 'w') as fh:
+            fh.write('\t'.join(['genome', 'tree_label', 'subtype',
+                                'probability', 'phylotyper_assignment', 'loci']))
+            fh.write('\t'.join(['lcl|query|', 'not applicable', 'not applicable',
+                                'not applicable', 'Subtype loci not found in genome', 'not applicable']))
 
 def phylotyper(uriIsolate, subtype, result_file, id_file=None, job_id=None, job_turtle=None, job_ectyper_datastruct_vf=None):
     """ Wrapper for Phylotyper
@@ -79,41 +83,38 @@ def phylotyper(uriIsolate, subtype, result_file, id_file=None, job_id=None, job_
     loci_results = ontology.schema_query(uri)
     loci = [ gu(l['locus']) for l in sorted(loci_results, key=lambda k: k['i'])]
 
-    # Get alleles for this genome
-    markerseqs = MarkerSequences(loci, job_id, job_turtle, job_ectyper_datastruct_vf, redis_conn)
-    fasta = markerseqs.fasta(uriIsolate)
-
-    temp_dir = mkdtemp(prefix='pt'+subtype, dir=config.DATASTORE)
+    # Temp files.
+    temp_dir = mkdtemp(prefix='pt' + subtype, dir=config.DATASTORE)
     query_file = os.path.join(temp_dir, 'query.fasta')
     output_file = os.path.join(temp_dir, 'subtype_predictions.tsv')
 
-    if fasta:
-        # Run phylotyper
-        with open(query_file, 'w') as fh:
-            fh.write(fasta)
-
-        subprocess.check_call(['phylotyper', 'genome', '--noplots',
-                         subtype,
-                         temp_dir,
-                         query_file])
-
+    # Get alleles for this genome
+    markerseqs = MarkerSequences(loci, job_id, job_turtle, job_ectyper_datastruct_vf, redis_conn)
+    # Validation all the alleles required are in this genome.
+    if not markerseqs.validate(subtype):
+        _noloci(output_file)
     else:
-        # No loci
-        # raise Exception('phylotyper.phylotyper() could not retrieve reference sequences for loci: {0}, uriIsolate: {1}, subtype: {2}'.format(
-        #     str(loci),
-        #     str(uriIsolate),
-        #     subtype
-        # ))
-        # Report no loci status in output
-        with open(output_file, 'w') as fh:
-            fh.write('\t'.join(['genome','tree_label','subtype','probability','phylotyper_assignment','loci']))
-            fh.write('\t'.join(['lcl|query|','not applicable','not applicable','not applicable','Subtype loci not found in genome','not applicable']))
+        fasta = markerseqs.fasta(uriIsolate)
+        if fasta:
+            # Run phylotyper
+            with open(query_file, 'w') as fh:
+                fh.write(fasta)
+            subprocess.check_call(
+                ['phylotyper',
+                'genome',
+                '--noplots',
+                subtype,
+                temp_dir,
+                query_file],
+                stderr=subprocess.STDOUT)
+
+        else:
+            # No loci
+            _noloci(output_file)
 
     shutil.move(output_file, result_file)
     shutil.rmtree(temp_dir)
 
-    # _check_tsv(result_file)
-
     return result_file
 
 

app/modules/phylotyper/sequences.py

@@ -10,6 +10,8 @@
 from middleware.decorators import submit, prefix, tojson
 from middleware.graphers import turtle_utils
 from routes.job_utils import fetch_job
+from modules.phylotyper.ontology import stx1_graph, stx2_graph, eae_graph, LOCI
+from middleware.graphers.turtle_utils import generate_uri as gu
 
 __location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__)))
 
@@ -136,8 +138,10 @@ def __init__(self, markers=[':stx2A',':stx2B'], job_id=None, job_turtle=None, jo
         g = Graph()
         ontology_turtle_file = os.path.join(__location__, 'superphy_subtyping.ttl')
         g.parse(ontology_turtle_file, format="turtle")
+        # Add Phylotyper ontology graphs.
+        g = g + stx1_graph() + stx2_graph() + eae_graph()
         # Retrieve and merge graphs from pre-req. jobs.
-        self.graph = g + fetch_job(job_id, redis_conn).result + fetch_job(job_turtle, redis_conn).result + fetch_job(job_ectyper_datastruct_vf, redis_conn).result
+        self.graph = g + fetch_job(job_id, redis_conn).result +  fetch_job(job_turtle, redis_conn).result + fetch_job(job_ectyper_datastruct_vf, redis_conn).result
 
     def sequences(self, genome_uri):
         """Retrieve sequences for object alleles
@@ -197,6 +201,54 @@ def fasta(self, genome_uri):
 
         return fasta_string
 
+    @prefix
+    def _subtype_query(self):
+        """
+        Queries for a specific URI of given type
+
+        Returns:
+            dictionary
+
+        """
+        query = '''
+            SELECT ?subtype
+            WHERE {{
+                ?region a faldo:Region ; :hasPart ?subtype .
+                ?subtype a :VirulenceFactor .
+            }}
+            '''
+
+        return query
+
+    def _find_object(self, uri):
+        """
+        Returns true if URI is already in database
+
+        Args:
+            uri(str): URI with prefix defined in config.py
+            rdftype(str): the URI linked by a rdf:type relationship to URI
+
+        """
+
+        query = self._subtype_query()
+
+        query_result = self.graph.query(query)
+
+        l = [tup[0].toPython() for tup in query_result]
+        full_uri = str(gu(uri))
+
+        return full_uri in l
+
+    def validate(self, subtype):
+        """Checks that the MakerSequence.graph has all the alleles required
+        for phylotyper analysis. Returns False if not (& Phylotyper should
+        not be run).
+        """
+        # Check for existance of schema Marker components
+        for l in LOCI[subtype]:
+            if not self._find_object(l):
+                return False
+        return True
 
 if __name__=='__main__':
     import argparse