added check for invalid min-/maxDist; works with fasta36.3.8g

Author: trvinh

fdog/data/conda_requirements.yml

@@ -1,6 +1,6 @@
 blast
 hmmer
-fasta3
+fasta3=36.3.8i
 clustalw
 mafft
-muscle
+muscle=5.1

fdog/libs/alignment.py

@@ -140,7 +140,10 @@ def calc_aln_score(fa1, fa2, aln_strategy = 'local', debugCore = False):
     Return dictionary {gene_id:aln_score}
     """
     fdog_path = os.path.realpath(__file__).replace('/libs/alignment.py','')
-    fasta36_options = '%s %s -s BP62 -m 9 -d 0 -z -1 -E 100' % (fa1, fa2)
+    os.symlink(fa1, 'fasta36_1.fa')
+    os.symlink(fa2, 'fasta36_2.fa')
+    # fasta36_options = '%s %s -s BP62 -m 9 -d 0 -z -1 -E 100' % (fa1, fa2)
+    fasta36_options = 'fasta36_1.fa fasta36_2.fa -s BP62 -m 9 -d 0 -z -1 -E 100'
     fdog_path = os.path.realpath(__file__).replace('/libs/alignment.py','')
     fasta36_bin = check_fasta36_executable(fdog_path)
     if aln_strategy == 'global':
@@ -173,4 +176,6 @@ def calc_aln_score(fa1, fa2, aln_strategy = 'local', debugCore = False):
                 if re.search('\(\s+\d+\)', l):
                     l = re.sub(r'\(\s+','(', l)
                 aln_score[gene_id] = aln_score[gene_id] + int(l.split()[2])
+    os.remove('fasta36_1.fa')
+    os.remove('fasta36_2.fa')
     return(aln_score)

fdog/libs/preparation.py

@@ -20,11 +20,13 @@
 from pathlib import Path
 from Bio import SeqIO
 from Bio.Blast.Applications import NcbiblastpCommandline
+from ete3 import NCBITaxa
 
 import fdog.libs.zzz as general_fn
 import fdog.libs.fasta as fasta_fn
 import fdog.libs.blast as blast_fn
 import fdog.libs.output as output_fn
+import fdog.libs.tree as tree_fn
 
 
 ##### FUNCTIONS FOR DATA/INPUT PREPARATION #####
@@ -117,6 +119,42 @@ def check_blast_version(corepath, refspec):
             'ERROR: Error running blast (probably conflict with BLAST DBs versions)\n%s'
             % (NcbiblastpCommandline(query = query, db = blast_db)))
 
+def check_ranks_core_taxa(corepath, minDist, maxDist):
+    """ Check if all core taxa have a valid minDist and maxDist tax ID
+    Return 2 dictionaries of taxa for invalid minDist and maxDist, where
+    keys is taxon name and value is the next valid rank
+    """
+    invalid_minDist = []
+    invalid_maxDist = []
+    ncbi = NCBITaxa()
+    rank_list = ['species', 'genus', 'family', 'order', 'class', 'phylum', 'kingdom', 'superkingdom']
+    suggest_minIndex = rank_list.index(minDist)
+    suggest_maxIndex = rank_list.index(maxDist)
+    for f in os.listdir(corepath):
+        if os.path.isdir(f'{corepath}/{f}'):
+            id = f.split('@')[1]
+            lineage = ncbi.get_lineage(id)
+            ranks = ncbi.get_rank(lineage)
+            if len(general_fn.matching_elements(ranks, minDist)) < 1:
+                invalid_minDist.append(f)
+                index_minDist = rank_list.index(minDist) + 1
+                while index_minDist < len(rank_list):
+                    if len(general_fn.matching_elements(ranks, rank_list[index_minDist])) > 0:
+                        if index_minDist > suggest_minIndex:
+                            suggest_minIndex = index_minDist
+                        break
+                    index_minDist += 1
+            if len(general_fn.matching_elements(ranks, maxDist)) < 1:
+                invalid_maxDist.append(f)
+                index_maxDist = rank_list.index(maxDist) + 1
+                while index_maxDist < len(rank_list):
+                    if len(general_fn.matching_elements(ranks, rank_list[index_maxDist])) > 0:
+                        if index_maxDist > suggest_maxIndex:
+                            suggest_maxIndex = index_maxDist
+                        break
+                    index_maxDist += 1
+    return(invalid_minDist, invalid_maxDist, rank_list[suggest_minIndex], rank_list[suggest_maxIndex])
+
 
 def get_seed_id_from_fa(core_fa, refspec):
     """ Get seed ID from core ortholog fasta file
@@ -147,11 +185,18 @@ def identify_seed_id(seqFile, refspec, corepath, debug, silentOff):
     # otherwise, perform blast search
     blast_xml = blast_fn.do_blastsearch(seqFile, refspec_db, evalBlast = 0.001)
     blast_out = blast_fn.parse_blast_xml(blast_xml)
+    if len(blast_out['hits']) < 1:
+        print(f'ERROR: Cannot find seed sequence {blast_out["query"]} in genome of reference species!')
+        print(f'You can check it by running:\nblastp -query {seqFile} -db {corepath}/{refspec}/{refspec} -evalue 0.001 -outfmt 7')
+        sys.exit()
     for hit in blast_out['hits']:
         if blast_out['hits'][hit]['align_len'] == blast_out['query_len']:
+            print("BEST BLAST HIT")
             return(hit)
         elif abs(int(blast_out['hits'][hit]['align_len']) - int(blast_out['query_len'])) < 10:
             output_fn.print_stdout(silentOff, 'WARNING: Found seed sequence shorter/longer than input!')
             return(hit)
         else:
-            sys.exit('ERROR: Cannot find seed sequence in genome of reference species for %s!' % blast_out['query'])
+            print(f'ERROR: Cannot find seed sequence {blast_out["query"]} in genome of reference species!')
+            print(f'You can check it by running:\nblastp -query {seqFile} -db {corepath}/{refspec}/{refspec} -evalue 0.001 -outfmt 7')
+            sys.exit()

fdog/runMulti.py

@@ -320,12 +320,25 @@ def main():
 
     begin = time.time()
     ##### Check and group parameters
+    print('Preparing & Checking...')
     (inFol, hmmpath, corepath, searchpath, annopath) = prepare_fn.check_input(
                     [inFol, refspec, outpath, hmmpath,
                     corepath, searchpath, annopath, pathFile])
     pathArgs = [outpath, hmmpath, corepath, searchpath, annopath]
     prepare_fn.check_blast_version(corepath, refspec)
 
+    (invalid_minDist, invalid_maxDist, suggested_minRank, suggested_maxRank) = prepare_fn.check_ranks_core_taxa(corepath, minDist, maxDist)
+    if len(invalid_minDist) > 0 or len(invalid_maxDist) > 0:
+        print(f'Invalid {minDist} (--minDist) for {len(invalid_minDist)} species:\n{invalid_minDist}')
+        print(f'Invalid {maxDist} (--maxDist) for {len(invalid_maxDist)} species:\n{invalid_maxDist}')
+        if not minDist == "genus" and not maxDist == "kingdom":
+            print(f'Please consider setting --minDist and --maxDist with these valid ranks:\n--minDist {suggested_minRank} --maxDist {suggested_maxRank}')
+            sys.exit()
+        else:
+            print(f'WARNING: --minDist and --maxDist will be automatically changed to {suggested_minRank} and {suggested_maxRank}')
+            minDist = suggested_minRank
+            maxDist = suggested_maxRank
+
     if not fasOff:
         check_fas = fas_fn.check_fas_executable()
         if check_fas == 0:

fdog/runSingle.py

@@ -189,6 +189,7 @@ def main():
 
     begin = time.time()
     ##### Check and group parameters
+    print('Preparing & Checking...')
     if seqFile == 'infile.fa':
         fdogPath = os.path.realpath(__file__).replace('/runSingle.py','')
         seqFile = '%s/data/infile.fa' % fdogPath
@@ -200,6 +201,18 @@ def main():
 
     prepare_fn.check_blast_version(corepath, refspec)
 
+    (invalid_minDist, invalid_maxDist, suggested_minRank, suggested_maxRank) = prepare_fn.check_ranks_core_taxa(corepath, minDist, maxDist)
+    if len(invalid_minDist) > 0 or len(invalid_maxDist) > 0:
+        print(f'Invalid {minDist} (--minDist) for {len(invalid_minDist)} species:\n{invalid_minDist}')
+        print(f'Invalid {maxDist} (--maxDist) for {len(invalid_maxDist)} species:\n{invalid_maxDist}')
+        if not minDist == "genus" and not maxDist == "kingdom":
+            print(f'Please consider setting --minDist and --maxDist with these valid ranks:\n--minDist {suggested_minRank} --maxDist {suggested_maxRank}')
+            sys.exit()
+        else:
+            print(f'WARNING: --minDist and --maxDist will be automatically changed to {suggested_minRank} and {suggested_maxRank}')
+            minDist = suggested_minRank
+            maxDist = suggested_maxRank
+
     if not fasOff:
         check_fas = fas_fn.check_fas_executable()
         if check_fas == 0:

setup.py

@@ -26,7 +26,7 @@
 
 setup(
     name="fdog",
-    version="0.1.17",
+    version="0.1.18",
     python_requires='>=3.7.0',
     description="Feature-aware Directed OrtholoG search tool",
     long_description=long_description,