Merge pull request #414 from Ecogenomics/staging

pchaumeil · web-flow · commit 3656d719d170 · 2022-07-06T15:03:56.000+10:00
Staging to Master for 2.1.1
diff --git a/docs/requirements.txt b/docs/requirements.txt
@@ -4,3 +4,6 @@ sphinx-rtd-theme ~= 0.5.0
 recommonmark ~= 0.7.0
 sphinx-sitemap ~= 2.2.0
 nbsphinx ~= 0.8.0
+matplotlib ~= 3.5.2
+linuxdoc == 20211220
+jupyter ~= 1.0.0
diff --git a/docs/src/announcements.rst b/docs/src/announcements.rst
@@ -3,7 +3,7 @@ Announcements
 
 
 GTDB-Tk 2.1.0 available
--------------------
+-----------------------
 
 *May 11, 2022*
 
diff --git a/docs/src/commands/convert_to_itol.rst b/docs/src/commands/convert_to_itol.rst
@@ -0,0 +1,41 @@
+.. _commands/convert_to_itol:
+
+convert_to_itol
+===============
+
+The `convert_to_itol` command will remove internal labels from Newick tree, making it suitable for visualization in `iTOL <http://itol.embl.de/>`_.  
+
+Arguments
+---------
+
+.. argparse::
+   :module: gtdbtk.cli
+   :func: get_main_parser
+   :prog: gtdbtk
+   :path: convert_to_itol
+   :nodefaultconst:
+
+Example
+-------
+
+Input
+^^^^^
+
+
+.. code-block:: bash
+
+    gtdbtk convert_to_itol --input some_tree.tree --output itol.tree
+
+
+Output
+^^^^^^
+
+
+.. code-block:: text
+
+    [2022-06-30 18:44:54] INFO: GTDB-Tk v2.1.0
+    [2022-06-30 18:44:54] INFO: gtdbtk convert_to_itol --input /tmp/decorated.tree --output new.tree
+    [2022-06-30 18:44:54] INFO: Using GTDB-Tk reference data version r207: /gtdbtk-data
+    [2022-06-30 18:44:54] INFO: Convert GTDB-Tk tree to iTOL format
+    [2022-06-30 18:44:54] INFO: Done.
+
diff --git a/docs/src/commands/index.rst b/docs/src/commands/index.rst
@@ -14,6 +14,7 @@ Below is a list of all GTDB-Tk command line options:
    check_install
    classify
    classify_wf
+   convert_to_itol
    de_novo_wf
    decorate
    export_msa
diff --git a/docs/src/conf.py b/docs/src/conf.py
@@ -28,8 +28,8 @@
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
-extensions = ['sphinxarg.ext', 'sphinx.ext.napoleon', 'sphinx.ext.autodoc',
-              'recommonmark', 'sphinx_sitemap', 'nbsphinx']
+extensions = ['sphinxarg.ext', 'sphinx.ext.napoleon', 'sphinx.ext.autodoc', 'linuxdoc.rstFlatTable',
+              'recommonmark', 'sphinx_sitemap', 'nbsphinx','matplotlib.sphinxext.plot_directive']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
diff --git a/docs/src/files/summary.tsv.rst b/docs/src/files/summary.tsv.rst
@@ -12,11 +12,11 @@ Classifications provided by the GTDB-Tk are in the files \<prefix>.bac120.summar
 * fastani_reference_radius: indicates the species-specific ANI circumscription radius of the reference genomes used to determine if a query genome should be classified to the same species as the reference.
 * fastani_taxonomy: indicates the GTDB taxonomy of the above reference genome.
 * fastani_ani: indicates the ANI between the query and above reference genome.
-* fastani_af: indicates the AF between the query and above reference genome.
+* fastani_af: indicates the alignment fraction (AF) between the query and above reference genome.
 * closest_placement_reference: indicates the accession number of the reference genome when a genome is placed on a terminal branch. 
 * closest_placement_taxonomy: indicates the GTDB taxonomy of the above reference genome.
 * closest_placement_ani: indicates the ANI between the query and above reference genome.
-* closest_placement_af: indicates the AF between the query and above reference genome.
+* closest_placement_af: indicates the alignment fraction (AF) between the query and above reference genome.
 * pplacer_taxonomy: indicates the pplacer taxonomy of the query genome.
 * classification_method:	indicates the rule used to classify the genome. This field will be one of: i) ANI, indicating a species assignement was based solely on the calculated ANI and AF with a reference genome; ii) ANI/Placement, indicating a species assignment was made based on both ANI and the placement of the genome in the reference tree; iii) taxonomic classification fully defined by topology, indicating that the classification could be determine based solely on the genome's position in the reference tree; or iv) taxonomic novelty determined using RED, indicating that the relative evolutionary divergence (RED) and placement of the genome in the reference tree were used to determine the classification.
 * note: provides additional information regarding the classification of the genome. Currently this field is only filled out when a species determination is made and indicates if the placement of the genome in the reference tree and closest reference according to ANI/AF are the same (congruent) or different (incongruent). 
diff --git a/docs/src/index.rst b/docs/src/index.rst
@@ -46,6 +46,7 @@ We encourage you to cite GTDB-Tk and the third-party dependencies as described i
    :caption: Running GTDB-Tk
    :maxdepth: 1
 
+   performance/index
    commands/index
    files/index
    examples/classify_wf
diff --git a/docs/src/installing/index.rst b/docs/src/installing/index.rst
@@ -38,7 +38,7 @@ Hardware requirements
      - ~65 GB
      - ~1 hour / 1,000 genomes @ 64 CPUs
    * - Bacteria
-     - ~320 GB ( 55GB for divide-and-conquer)
+     - ~55GB (320 GB when using --full_tree)
      - ~65 GB
      - ~1 hour / 1,000 genomes @ 64 CPUs
 
diff --git a/docs/src/performance/accuracy.rst b/docs/src/performance/accuracy.rst
@@ -0,0 +1,58 @@
+.. _performance/Accuracy:
+
+
+Accuracy
+========
+
+The similarity of GTDB-Tk v1 and v2 classifications were first assessed using 16,710 bacterial genomes from the GEMs dataset (Nayfach et al., 2021) that represent novel taxa relative to GTDB R07-RS207.
+| Only 12 genomes (0.07%) did not have identical classifications between GTDB-Tk v1 and the divide-and-conquer approach used in GTDB-Tk v2.
+| The majority of incongruence was due to genomes being over- (6 genomes) or under-classified (4 genomes) by a single taxonomic rank. Only 2 genomes had conflicting taxonomic assignments, and these were both relatively poor-quality genomes assigned as new classes in alternative phyla.
+
+.. flat-table:: Table 1. Novelty of GEM genomes relative to GTDB R07-RS207 based on GTDB-Tk v1 classifications.
+   :header-rows: 2
+
+   * -
+     -
+     - :cspan:`4` GTDB-Tk v2 classifications relative to GTDB-Tk v1 classifications
+   * - Toxon Novelty
+     - No genomes
+     - Congruent
+     - Conflict
+     - Underclassified
+     - Overclassified
+   * - Novel phylum
+     - 3
+     - 2
+     - 0
+     - 0
+     - 1
+   * - Novel class
+     - 42
+     - 35
+     - 2
+     - 2
+     - 2
+   * - Novel order
+     - 144
+     - 143
+     - 0
+     - 0
+     - 1
+   * - Novel family
+     - 543
+     - 540
+     - 0
+     - 1
+     - 2
+   * - Novel genus
+     - 3,222
+     - 3,219
+     - 0
+     - 1
+     - 0
+   * - Novel species
+     - 12,756
+     - 12,756
+     - 0
+     - 0
+     - 0
diff --git a/docs/src/performance/index.rst b/docs/src/performance/index.rst
@@ -0,0 +1,12 @@
+.. _performance:
+
+############################
+  Performance and Accuracy
+############################
+
+
+.. toctree::
+   :maxdepth: 1
+
+   performance
+   accuracy
diff --git a/docs/src/performance/performance.rst b/docs/src/performance/performance.rst
@@ -0,0 +1,77 @@
+.. _performance/Performance:
+
+Performance
+===========
+
+
+| GTDB-Tk v2 runs 22% to 35% faster when processing 1000 genomes with 1 to 64 CPUs (**Fig.1**) and is >40% faster when processing 5,000 genomes using 32 CPUs (**Fig.2**).
+| The tests below were run on a machine with 4 AMD EPYC 7402 24-Core Processor and 512 GB of RAM.
+
+.. plot::
+
+      import matplotlib.pyplot as plt
+      from matplotlib.ticker import ScalarFormatter
+      import numpy as np
+
+      color_setup = ['#5a6855','red','color']
+      setup = color_setup
+      cpus_list = [1, 8, 16, 32, 64]
+      split_list = [x/60 for x in [2571, 457, 309, 230, 216]]
+      nosplit_list= [x/60 for x in [3980, 638, 398, 295, 279]]
+      values = ['1', '8', '16', '32','64']
+
+      fig, ax = plt.subplots()
+
+      plt.ylim(3,70)
+      plt.xticks(cpus_list, values)
+
+      ax.scatter(cpus_list, split_list, label="GTDB-Tk v2", marker="s", s=30, color=setup[0])
+      ax.plot(cpus_list, split_list,linestyle='dashed', color=setup[0])
+      ax.scatter(cpus_list, nosplit_list, label="GTDB-Tk v1",color=setup[1])
+      ax.plot(cpus_list, nosplit_list,linestyle=':', color=setup[1])
+
+      ax.set_yscale('log')
+      ax.set_yticks([3,5,7,10,20,30,50,70])
+      ax.yaxis.set_major_formatter(ScalarFormatter())
+
+      plt.ylabel('Runtime (h)')
+      plt.xlabel('Number of CPUs')
+      plt.title('Fig.1: GTDB-Tk runtime for 1000 genomes (log scale)')
+
+      # show a legend on the plot
+      plt.legend(loc=1, prop={'size': 12},frameon=False)
+      plt.show()
+
+**Fig. 1**: Processing time (log scale) for 1,000 randomly selected GEM MAGs for increasing numbers of CPUs.
+
+.. plot::
+
+      color_setup = ['#5a6855','red','color']
+      setup = color_setup
+      pool_size_list=[10, 50, 100, 200, 500, 1000, 2000, 5000]
+      split_list=[x/60 for x in [88, 150, 160, 169, 195, 235, 312, 558]]
+      nosplit_list=[x/60 for x in [137, 151, 163, 180, 219, 280, 416, 934]]
+      values = [10,500,1000,1500,2000,2500,3000,3500,4000,4500,5000]
+
+      plt.scatter(pool_size_list, split_list, label="GTDB-Tk v2" , marker="s", s=30, color=setup[0])
+      plt.plot(pool_size_list, split_list,linestyle='dashed', color=setup[0])
+      plt.scatter(pool_size_list, nosplit_list, label="GTDB-Tk v1",color=setup[1])
+      plt.plot(pool_size_list, nosplit_list,linestyle=':', color=setup[1])
+
+      # naming the x axis
+
+      plt.xticks(values)
+      plt.yticks([1,3,5,7,9,11,13,15])
+
+      # naming the axis
+      plt.ylabel('Runtime (h)')
+      plt.xlabel('Number of genomes')
+      # giving a title to my graph
+      plt.title('Fig.2: GTDB-Tk runtime with 32CPUs')
+
+      # show a legend on the plot
+      plt.legend(loc=2, prop={'size': 12},frameon=False)
+
+      # function to show the plot
+      plt.show()
+**Fig. 2**: Processing time with 32 CPUs on increasing numbers of randomly selected GEM MAGs.
diff --git a/gtdbtk/__init__.py b/gtdbtk/__init__.py
@@ -29,4 +29,4 @@
 __status__ = 'Production'
 __title__ = 'GTDB-Tk'
 __url__ = 'https://github.com/Ecogenomics/GTDBTk'
-__version__ = '2.1.0'
+__version__ = '2.1.1'
diff --git a/gtdbtk/config/config.py b/gtdbtk/config/config.py
@@ -47,7 +47,7 @@
 
 BACKBONE_PPLACER_REF_PKG = 'gtdbtk_package_backbone.refpkg'
 BACKBONE_RED_FILE = os.path.join(BACKBONE_RED_DIR, 'backbone_red_value.tsv')
-CLASS_LEVEL_PPLACER_REF_PKG = os.path.join(CLASS_LEVEL_PPLACER_DIR, 'gtdbtk.package.{iter}.refpkg')
+CLASS_LEVEL_PPLACER_REF_PKG = 'gtdbtk.package.{iter}.refpkg'
 CLASS_LEVEL_RED_FILE = os.path.join(CLASS_LEVEL_RED_DIR, 'red_value_{iter}.tsv')
 
 RED_DIST_BAC_DICT = ''
@@ -218,7 +218,7 @@
 PPLACER_AR53_REF_PKG = f"gtdb_{VERSION_DATA}_ar53.refpkg"
 PPLACER_RPS23_REF_PKG = f"gtdb_{VERSION_DATA}_rps23.refpkg"
 PPLACER_MIN_RAM_BAC_FULL = 320
-PPLACER_MIN_RAM_BAC_SPLIT = 50
+PPLACER_MIN_RAM_BAC_SPLIT = 55
 PPLACER_MIN_RAM_ARC = 40
 
 # Fastani configuration
diff --git a/gtdbtk/external/pfam_search.py b/gtdbtk/external/pfam_search.py
@@ -85,14 +85,15 @@ def _workerThread(self, queueIn, queueOut, n_skipped):
         """Process each data item in parallel."""
         try:
             while True:
-                gene_file = queueIn.get(block=True, timeout=None)
-                if gene_file is None:
+                queue_next = queueIn.get(block=True, timeout=None)
+                if queue_next is None:
                     break
 
-                genome_dir, filename = os.path.split(gene_file)
-                genome_id = filename.replace(self.protein_file_suffix, '')
-                output_hit_file = os.path.join(self.output_dir, genome_id, filename.replace(self.protein_file_suffix,
-                                                                                            self.pfam_suffix))
+                genome_id, gene_file = queue_next
+
+                output_hit_file = os.path.join(self.output_dir, genome_id, '{}{}'.format(genome_id, self.pfam_suffix))
+                output_tophit_file = os.path.join(self.output_dir, genome_id,
+                                                  '{}{}'.format(genome_id, self.pfam_top_hit_suffix))
 
                 # Check if this has already been processed.
                 out_files = (output_hit_file, TopHitPfamFile.get_path(self.output_dir, genome_id))
diff --git a/gtdbtk/external/tigrfam_search.py b/gtdbtk/external/tigrfam_search.py
@@ -20,6 +20,7 @@
 import os
 import subprocess
 
+from gtdbtk.biolib_lite.common import make_sure_path_exists
 from gtdbtk.exceptions import GTDBTkExit
 from gtdbtk.io.marker.tophit import TopHitTigrFile
 from gtdbtk.tools import sha256, file_has_checksum, tqdm_log
@@ -101,18 +102,25 @@ def _topHit(self, tigrfam_file):
     def _workerThread(self, queueIn, queueOut, n_skipped):
         """Process each data item in parallel."""
         while True:
-            gene_file = queueIn.get(block=True, timeout=None)
-            if gene_file is None:
+            queue_next = queueIn.get(block=True, timeout=None)
+            if queue_next is None:
                 break
 
-            assembly_dir, filename = os.path.split(gene_file)
-            genome_id = filename.replace(self.protein_file_suffix, '')
-            genome_dir = os.path.join(self.output_dir, genome_id)
-            output_hit_file = os.path.join(genome_dir, filename.replace(self.protein_file_suffix,
-                                                                        self.tigrfam_suffix))
+            genome_id, gene_file = queue_next
+            output_hit_file = os.path.join(self.output_dir, genome_id, '{}{}'.format(genome_id, self.tigrfam_suffix))
+            # output_tophit_file = os.path.join(self.output_dir, genome_id,
+            #                                   '{}{}'.format(genome_id, self.tigrfam_top_hit_suffix))
+            hmmsearch_out = os.path.join(self.output_dir, genome_id, '{}_tigrfam.out'.format(genome_id))
 
-            hmmsearch_out = os.path.join(genome_dir, filename.replace(self.protein_file_suffix,
-                                                                      '_tigrfam.out'))
+
+            # assembly_dir, filename = os.path.split(gene_file)
+            # genome_id = filename.replace(self.protein_file_suffix, '')
+            # genome_dir = os.path.join(self.output_dir, genome_id)
+            # output_hit_file = os.path.join(genome_dir, filename.replace(self.protein_file_suffix,
+            #                                                             self.tigrfam_suffix))
+
+            # hmmsearch_out = os.path.join(genome_dir, filename.replace(self.protein_file_suffix,
+            #                                                           '_tigrfam.out'))
 
             # Check if this has already been processed.
             out_files = (output_hit_file, hmmsearch_out, TopHitTigrFile.get_path(self.output_dir, genome_id))
@@ -122,6 +130,11 @@ def _workerThread(self, queueIn, queueOut, n_skipped):
                     n_skipped.value += 1
 
             else:
+
+                genome_dir = os.path.join(self.output_dir, genome_id)
+                hmmsearch_out = os.path.join(genome_dir, '{}_tigrfam.out'.format(genome_id))
+                make_sure_path_exists(genome_dir)
+
                 args = ['hmmsearch', '-o', hmmsearch_out, '--tblout', output_hit_file,
                         '--noali', '--notextw', '--cut_nc', '--cpu',
                         str(self.cpus_per_genome), self.tigrfam_hmms, gene_file]
diff --git a/gtdbtk/io/prodigal/tln_table_summary.py b/gtdbtk/io/prodigal/tln_table_summary.py
@@ -53,4 +53,4 @@ def read(self):
         with open(self.path, 'r') as fh:
             for line in fh.readlines():
                 gid, tbl = line.strip().split('\t')
-                self.genomes[gid] = int(tbl)
+                self.genomes[gid] = str(tbl)
diff --git a/gtdbtk/markers.py b/gtdbtk/markers.py
diff --git a/scripts/restructure_reps_fna_folder.sh b/scripts/restructure_reps_fna_folder.sh
