Simplified code to run MAFFT and verify results.

Author: rfm-targa

CHEWBBACA/AlleleCallEvaluator/evaluate_calls.py

@@ -419,7 +419,7 @@ def main(input_files, schema_directory, output_directory, annotations,
 
             protein_files = [r[1] for r in results]
 
-            # Align sequences with MAFFT
+            # Run MAFFT to compute MSA
             print('\nDetermining the MSA for each locus...')
             mafft_outdir = fo.join_paths(temp_directory, ['alignment_files'])
             fo.create_directory(mafft_outdir)
@@ -436,15 +436,23 @@ def main(input_files, schema_directory, output_directory, annotations,
                                                common_args,
                                                mw.call_mafft)
 
-            results = mo.map_async_parallelizer(inputs,
-                                                mo.function_helper,
-                                                cpu_cores,
-                                                show_progress=True)
+            mafft_results = mo.map_async_parallelizer(inputs,
+                                                      mo.function_helper,
+                                                      cpu_cores,
+                                                      show_progress=True)
+
+            # Identify cases where MAFFT failed
+            mafft_failed = [r[0] for r in mafft_results if r[1] is False]
+            if len(mafft_failed) > 0:
+                print(f'\nCould not determine MSA for {len(mafft_failed)} loci.')
+
+            # Get files that were created by MAFFT
+            mafft_successful = [r[0] for r in mafft_results if r[1] is True]
 
             print('\nCreating file with the full cgMLST alignment...', end='')
             # Concatenate all alignment files and index with BioPython
             concat_aln = fo.join_paths(mafft_outdir, ['cgMLST_concat.fasta'])
-            fo.concatenate_files(mafft_outfiles, concat_aln)
+            fo.concatenate_files(mafft_successful, concat_aln)
             # Index file
             concat_index = fao.index_fasta(concat_aln)
             sample_alignment_files = []
@@ -495,25 +503,15 @@ def main(input_files, schema_directory, output_directory, annotations,
     report_html_file = fo.join_paths(output_directory, [ct.ALLELECALL_REPORT_BASENAME])
     fo.write_to_file(report_html, report_html_file, 'w', '\n')
 
-    # Copy JS bundle to output directory
+    # Copy the JS bundle to the output directory
+    # When chewBBACA is installed
     script_path = os.path.dirname(os.path.abspath(__file__))
-    parent_dir = os.path.dirname(script_path)
-    try:
-        fo.copy_file(fo.join_paths(script_path,
-                                   ['report_template_components',
-                                    'src',
-                                    'bundles',
-                                    'AlleleCallEvaluator',
-                                    'report_bundle.js']),
-                     output_directory)
-    except Exception as e:
-        fo.copy_file(fo.join_paths(parent_dir,
-                                   ['report_template_components',
-                                    'src',
-                                    'bundles',
-                                    'AlleleCallEvaluator',
-                                    'report_bundle.js']),
-                     output_directory)
+    # For development
+    if script_path.endswith('CHEWBBACA') is False:
+        script_path = os.path.dirname(script_path)
+
+    fo.copy_file(fo.join_paths(script_path, [ct.ALLELECALL_EVALUATOR_BUNDLE]),
+                 output_directory)
 
     # Delete all temporary files
     fo.delete_directory(temp_directory)

CHEWBBACA/SchemaEvaluator/evaluate_schema.py

@@ -502,23 +502,27 @@ def locus_report(locus_file, locus_data, annotation_columns,
             mafft_outfile = os.path.basename(mafft_infile)
             mafft_outfile = mafft_outfile.replace('.fasta', '_aligned.fasta')
             mafft_outfile = fo.join_paths(output_directory, [mafft_outfile])
-            mafft_outfile_exists = mw.call_mafft(mafft_infile, mafft_outfile)
-            # Get MSA data
-            alignment_text = fo.read_file(mafft_outfile)
-            msa_data['sequences'] = alignment_text
-
-            # Get Tree data
-            # Get the phylocanvas data
-            tree_file = mafft_outfile.replace('_aligned.fasta', '.fasta.tree')
-            phylo_data = fo.read_file(tree_file)
-
-            # Start by substituting greatest value to avoid substituting
-            # smaller values contained in greater values
-            for i in range(locus_data[1], 0, -1):
-                phylo_data = phylo_data.replace(f'{i}_', '')
-
-            phylo_data = phylo_data.replace('\n', '')
-            phylo_data = {"phylo_data": phylo_data}
+            # Run MAFFT to compute MSA
+            mafft_results = mw.call_mafft(mafft_infile, mafft_outfile)
+            if mafft_results[1] is True:
+                # Get MSA data
+                alignment_text = fo.read_file(mafft_outfile)
+                msa_data['sequences'] = alignment_text
+
+                # Get Tree data
+                # Get the phylocanvas data
+                tree_file = mafft_outfile.replace('_aligned.fasta', '.fasta.tree')
+                phylo_text = fo.read_file(tree_file)
+
+                # Start by substituting greatest value to avoid substituting
+                # smaller values contained in greater values
+                for i in range(locus_data[1], 0, -1):
+                    phylo_text = phylo_text.replace(f'{i}_', '')
+
+                phylo_text = phylo_text.replace('\n', '')
+                phylo_data['phylo_data'] = phylo_text
+            else:
+                print(f'Could not compute MSA and phylogenetic tree for locus {locus}')
 
     locus_columns = ct.LOCUS_COLUMNS.format(minimum_length,
                                             locus_data[14],
@@ -772,50 +776,20 @@ def main(schema_directory, output_directory, genes_list, annotations,
 
     # Copy the JS bundle files to the respective directories
     # JS bundle used by schema report
-    script_path = os.path.dirname(os.path.abspath(__file__))
-    parent_dir = os.path.dirname(script_path)
     # When chewBBACA is installed
-    try:
-        shutil.copy(fo.join_paths(script_path,
-                                  ['report_template_components',
-                                   'src',
-                                   'bundles',
-                                   'SchemaEvaluator',
-                                   'schema_report',
-                                   'report_bundle.js']),
-                    output_directory)
+    script_path = os.path.dirname(os.path.abspath(__file__))
     # For development
-    except Exception as e:
-        shutil.copy(fo.join_paths(parent_dir,
-                                  ['report_template_components',
-                                   'src',
-                                   'bundles',
-                                   'SchemaEvaluator',
-                                   'schema_report',
-                                   'report_bundle.js']),
-                    output_directory)
+    if script_path.endswith('CHEWBBACA') is False:
+        script_path = os.path.dirname(script_path)
+
+    shutil.copy(fo.join_paths(script_path, [ct.SCHEMA_EVALUATOR_SCHEMA_BUNDLE]),
+                output_directory)
+
     if loci_reports is True:
-        # JS bundle used by loci reports
-        try:
-            shutil.copy(fo.join_paths(script_path,
-                                      ['report_template_components',
-                                       'src',
-                                       'bundles',
-                                       'SchemaEvaluator',
-                                       'loci_reports',
-                                       'report_bundle.js']),
-                        html_dir)
-        except Exception as e:
-            shutil.copy(fo.join_paths(parent_dir,
-                                      ['report_template_components',
-                                       'src',
-                                       'bundles',
-                                       'SchemaEvaluator',
-                                       'loci_reports',
-                                       'report_bundle.js']),
-                        html_dir)
+        shutil.copy(fo.join_paths(script_path, [ct.SCHEMA_EVALUATOR_LOCI_BUNDLE]),
+                    html_dir)
 
     # Delete all temporary files
     fo.delete_directory(temp_directory)
 
-    print(f'\nResults available in {output_directory}.')
+    print(f'\nResults available in {os.path.abspath(output_directory)}.')

CHEWBBACA/utils/blast_wrapper.py

@@ -48,15 +48,17 @@ def make_blast_db(makeblastdb_path, input_fasta, output_path, db_type,
         A list with the warnings and errors raised by BLAST.
     """
     # use '-parse-seqids' to be able to retrieve/align sequences by identifier
-    blastdb_cmd = [makeblastdb_path, '-in', input_fasta,
+    makedb_cmd = [makeblastdb_path, '-in', input_fasta,
                    '-out', output_path, '-parse_seqids',
                    '-dbtype', db_type]
 
-    makedb_cmd = subprocess.Popen(blastdb_cmd,
-                                  stdout=subprocess.PIPE,
-                                  stderr=subprocess.PIPE)
+    makedb_process = subprocess.Popen(makedb_cmd,
+                                      stdout=subprocess.PIPE,
+                                      stderr=subprocess.PIPE)
+
+    #stdout, stderr = makedb_process.communicate()
 
-    stderr = makedb_cmd.stderr.readlines()
+    stderr = makedb_process.stderr.readlines()
 
     # ignore errors/warnings provided to `ignore`
     if len(stderr) > 0:

CHEWBBACA/utils/constants.py

@@ -13,6 +13,7 @@
 
 
 import sys
+import shutil
 import platform
 
 
@@ -124,12 +125,9 @@
 BLAST_MAJOR = 2
 BLAST_MINOR = 9
 
-# BLASTp and makeblastdb aliases for Linux and Windows
-BLASTP_ALIAS = 'blastp.exe' if platform.system() == 'Windows' else 'blastp'
-MAKEBLASTDB_ALIAS = 'makeblastdb.exe' if platform.system() == 'Windows' else 'makeblastdb'
-
-# Prodigal alias (must be in PATH)
-PRODIGAL_PATH = 'prodigal'
+# Paths to BLASTp and makeblastdb executables in Linux and Windows
+BLASTP_ALIAS = 'blastp.exe' if platform.system() == 'Windows' else shutil.which('blastp')
+MAKEBLASTDB_ALIAS = 'makeblastdb.exe' if platform.system() == 'Windows' else shutil.which('makeblastdb')
 
 # BLAST warnings to be ignored
 IGNORE_RAISED = ['Warning: [blastp] To obtain better run time '
@@ -138,6 +136,9 @@
                  '<OUT_FILE_NAME> and use <OUT_FILE_NAME> as the '
                  'argument to -seqidlist']
 
+# Path to MAFFT executable
+MAFFT_ALIAS = shutil.which('mafft')
+
 # Replacements for genome and loci identifiers
 CHAR_REPLACEMENTS = [("|", "_"), ("_", "-"), ("(", ""),
                      (")", ""), ("'", ""), ("\"", ""),
@@ -306,6 +307,15 @@
 PRESENCE_ABSENCE_BASENAME = 'presence_absence.tsv'
 MISSING_LOCI_BASENAME = 'missing_loci_stats.tsv'
 
+# Relative path to the JS bundles used by SchemaEvaluator
+# Main page
+SCHEMA_EVALUATOR_SCHEMA_BUNDLE = 'report_template_components/src/bundles/SchemaEvaluator/schema_report/report_bundle.js'
+# Loci pages
+SCHEMA_EVALUATOR_LOCI_BUNDLE = 'report_template_components/src/bundles/SchemaEvaluator/loci_reports/report_bundle.js'
+
+# Relative path to the JS bundle used by AlleleCallEvaluator
+ALLELECALL_EVALUATOR_BUNDLE = 'report_template_components/src/bundles/AlleleCallEvaluator/report_bundle.js'
+
 # Do not use list of strings as constants if the strings include formatting
 # placeholders. Multiple references to the list of strings will have the same
 # id and altering the strings with format will not change the list id. In

CHEWBBACA/utils/mafft_wrapper.py

@@ -13,35 +13,50 @@
 
 
 import os
-import shutil
 import subprocess
 
 try:
-    from utils import iterables_manipulation as im
+    from utils import (constants as ct,
+                       iterables_manipulation as im)
 except ModuleNotFoundError:
-    from CHEWBBACA.utils import iterables_manipulation as im
+    from CHEWBBACA.utils import (constants as ct,
+                                 iterables_manipulation as im)
 
 
 def call_mafft(input_file, output_file):
-    """Call MAFFT to generate an alignment.
+    """Call MAFFT to compute a MSA.
 
     Parameters
     ----------
     input_file : str
         Path to a FASTA file with the sequences to align.
+    output_file : str
+        Path to the output file created by MAFFT.
 
     Returns
     -------
-    Path to the file with the computed MSA if successful, False otherwise.
+    output_file : str
+        Path to the output file.
+    outfile_exists : bool
+        True if the output file was created, False otherwise.
+    stdout : bytes
+        MAFFT stdout.
+    stderr : bytes
+        MAFFT stderr.
+
+    True if the output file was created, False otherwise.
     """
-    mafft_cmd = [shutil.which('mafft'), '--thread', '1', '--treeout', '--retree', '1',
+    mafft_cmd = [ct.MAFFT_ALIAS, '--thread', '1', '--treeout', '--retree', '1',
                  '--maxiterate', '0', input_file, '>', output_file]
     mafft_cmd = ' '.join(mafft_cmd)
 
+    # Must include subprocess.PIPE to get stdout and stderr
     mafft_cmd = subprocess.Popen(mafft_cmd,
-                                 shell=True,
                                  stdout=subprocess.PIPE,
-                                 stderr=subprocess.PIPE)
-    mafft_cmd.wait()
+                                 stderr=subprocess.PIPE,
+                                 shell=True)
+    stdout, stderr = mafft_cmd.communicate()
 
-    return os.path.exists(output_file)
+    outfile_exists = os.path.exists(output_file)
+
+    return [output_file, outfile_exists, stdout, stderr]