crest4.classify

  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3
  4"""
  5Written by Lucas Sinclair.
  6GNUv3 Licensed.
  7Contact at www.sinclair.bio
  8"""
  9
 10# Built-in modules #
 11import os, multiprocessing
 12
 13# Internal modules #
 14import crest4
 15import crest4.databases
 16from crest4.query import Query
 17from crest4.databases import CrestDatabase
 18
 19# First party modules #
 20from functools import cached_property
 21from autopaths.file_path import FilePath
 22from autopaths.dir_path  import DirectoryPath
 23
 24# Constants #
 25all_db_choices = ('midori253darn', 'silvamod138pr2', 'mitofish', 'ssuome')
 26
 27###############################################################################
 28class Classify:
 29    """
 30    This is the main object offered by the `crest4` package.
 31    It enables you to automatically assign taxonomic names to DNA sequences
 32    obtained from environmental sequencing.
 33
 34    After creating a new instance of a `Classify` object, you can simply call
 35    it to get your input data processed and the assignments output file
 36    generated. Examples are included in the `README.md` file of this package,
 37    which can directly be seen on the GitHub page at:
 38
 39    https://github.com/xapple/crest4/
 40    """
 41
 42    def __init__(self,
 43                 fasta,
 44                 search_algo = 'blast',
 45                 num_threads = 1,
 46                 search_db   = 'ssuome',
 47                 output_dir  = None,
 48                 search_hits = None,
 49                 min_score   = None,
 50                 score_drop  = 2.0,
 51                 min_smlrty  = True,
 52                 otu_table   = None,
 53                 ):
 54        """
 55        Args:
 56
 57            fasta: The path to a single FASTA file as a string.
 58                   These are the sequences that will be taxonomically
 59                   classified.
 60
 61            search_algo: The algorithm used for the sequence similarity search
 62                         that will be run to match the sequences against the
 63                         database chosen. Either `blast` or `vsearch`. No
 64                         other values are currently supported. By default,
 65                         `blast`.
 66
 67            num_threads: The number of processors to use for the sequence
 68                         similarity search. By default, parallelism is turned
 69                         off, and this value is 1. If you pass the value
 70                         `True,` we will run as many processes as there are
 71                         CPUs but no more than 32.
 72
 73            search_db: The database used for the sequence similarity search.
 74                       Either `ssuome`, `silvamod138pr2`, 'mitofish',
 75                       or `midori253darn`. By default, `ssuome`.
 76                       Optionally, the user can provide a custom database by
 77                       specifying the full path to a directory containing all
 78                       required files under `search_db`.
 79                       See the README for more information.
 80
 81            output_dir: The directory into which all the classification
 82                        results will be written to. This defaults to a
 83                        directory with the same name as the original FASTA
 84                        file and a `.crest4` suffix appended.
 85
 86            search_hits: The path where the search results will be stored.
 87                         This defaults to the output directory. However,
 88                         if the search operation has already been completed
 89                         beforehand, specify the path here to skip the
 90                         sequence similarity search step and go directly to
 91                         the taxonomy step. If a hits file exists in the output
 92                         directory and this option is not specified, it is
 93                         deleted and regenerated.
 94
 95            min_score: The minimum bit-score for a search hit to be considered
 96                       when using BLAST as the search algorithm. All hits below
 97                       this score are ignored. When using VSEARCH, this value
 98                       instead indicates the minimum identity between two
 99                       sequences for the hit to be considered.
100                       The default is `155` for BLAST and `0.75` for VSEARCH.
101
102            score_drop: Determines the range of hits to retain and the range
103                        to discard based on a drop in percentage from the score
104                        of the best hit. Any hit below the following value:
105                        "(100 - score_drop)/100 * best_hit_score" is ignored.
106                        By default, `2.0`.
107
108            min_smlrty: Determines if the minimum similarity filter is turned
109                        on or off. Pass any value like `False` to turn it off.
110                        The minimum similarity filter prevents classification
111                        to higher ranks when a minimum rank-identity is not
112                        met. The default is `True`.
113
114            otu_table: Optionally, one can specify the path to an existing OTU
115                       table in CSV or TSV format when running `crest4`.
116                       The sequence names in the OTU table must be rows and
117                       have to match the names in the FASTA file. The column,
118                       on the other hand, provide your samples names.
119                       When this option is used, then two extra output files
120                       are generated. Firstly, a table summarizing the
121                       assignment counts per taxa. Secondly, a table
122                       propagating the sequence counts upwards
123                       in a cumulative fashion.
124                       """
125        # Save attributes #
126        self.fasta       = fasta
127        self.search_algo = search_algo
128        self.num_threads = num_threads
129        self.search_db   = search_db
130        self.output_dir  = output_dir
131        self.search_hits = search_hits
132        self.min_score   = min_score
133        self.score_drop  = score_drop
134        self.min_smlrty  = min_smlrty
135        self.otu_table   = otu_table
136        # Assign default values and change others #
137        self.transform()
138        # Validate attributes #
139        self.validate()
140
141    def transform(self):
142        """
143        This method will replace empty attributes with defaults when this is
144        needed and will convert others to proper types.
145        """
146        # The fasta should be a FASTA object #
147        if self.fasta is not None:
148            from fasta import FASTA
149            self.fasta = FASTA(self.fasta)
150        # Default for the number of threads #
151        if not isinstance(self.num_threads, int):
152            if self.num_threads is True:
153                self.num_threads = min(multiprocessing.cpu_count(), 32)
154            elif self.num_threads is False:
155                self.num_threads = 1
156            elif self.num_threads.lower() == 'true':
157                self.num_threads = min(multiprocessing.cpu_count(), 32)
158        # The database is always a string #
159        self.search_db = str(self.search_db)
160        # Default for the output directory #
161        if self.output_dir is None:
162            self.output_dir = self.fasta + '.crest4/'
163        self.output_dir = DirectoryPath(self.output_dir)
164        # The search hits is a file somewhere if passed #
165        if self.search_hits is not None:
166            self.search_hits = FilePath(self.search_hits)
167        # Default for the search hits file if not passed #
168        if self.search_hits is None:
169            self.search_hits = FilePath(self.output_dir + 'search.hits')
170            self.search_hits.remove()
171        # Default for the minimum score #
172        if self.min_score is None:
173            if self.search_algo == 'blast':
174                self.min_score = 155
175            if self.search_algo == 'vsearch':
176                self.min_score = 0.75
177        # The minimum score has to be a number, not a string #
178        try:
179            self.min_score = float(self.min_score)
180        except (ValueError, TypeError):
181            msg = "The minimum score value must be numerical (not '%s')."
182            raise ValueError(msg % self.min_score)
183        # The score drop has to be a float, not a string #
184        try:
185            self.score_drop = float(self.score_drop)
186        except (ValueError, TypeError):
187            msg = "The score drop value must be numerical (not '%s')."
188            raise ValueError(msg % self.score_drop)
189        # Turn off the minimum similarity filter if the user passed any value #
190        if self.min_smlrty is not True: self.min_smlrty = False
191        # The OTU table is a file somewhere if passed #
192        if self.otu_table is not None:
193            self.otu_table = FilePath(self.otu_table)
194            self.otu_table.must_exist()
195
196    def validate(self):
197        """
198        This method will raise an Exception if any of the arguments passed by
199        the user are illegal.
200        """
201        # The fasta should exist if passed #
202        if self.fasta is not None:
203            self.fasta.must_exist()
204        # Either the FASTA file or the hits file has to contain something #
205        if not self.fasta and not self.search_hits:
206            msg = "Neither the FASTA file at '%s' nor the search hits file" \
207                  " at '%s' contain any data. Cannot proceed."
208            raise Exception(msg % (self.fasta, self.search_hits))
209        # Check the search algorithm #
210        if self.search_algo not in ('blast', 'vsearch'):
211            msg = "The search algorithm '%s' is not supported."
212            raise ValueError(msg % self.search_algo)
213        # The search database is a known entry or exists on the filesystem #
214        if self.search_db not in all_db_choices:
215            if not os.path.exists(self.search_db):
216                msg = "The search database '%s' is not supported."
217                raise ValueError(msg % self.search_db)
218        # Check the minimum score value is above zero #
219        if self.min_score < 0.0:
220            msg = "The minimum score cannot be smaller than zero ('%s')."
221            raise ValueError(msg % self.min_score)
222        # Check the minimum score value is below one #
223        if self.min_score > 1.0:
224            if self.search_algo == 'vsearch':
225                msg = "The minimum score cannot be more than 1.0 when" \
226                      " using VSEARCH ('%s') because it represents the" \
227                      " the minimum identity between two sequences."
228                raise ValueError(msg % self.min_score)
229        # Check the score drop value #
230        if self.score_drop < 0.0:
231            msg = "The score drop value cannot be smaller than zero ('%s')."
232            raise ValueError(msg % self.min_score)
233        if self.score_drop > 100.0:
234            msg = "The score drop value cannot be over 100 ('%s')."
235            raise ValueError(msg % self.min_score)
236
237    def __repr__(self):
238        """A simple representation of this object to avoid memory addresses."""
239        return "<%s object on '%s'>" % (self.__class__.__name__, self.fasta)
240
241    @cached_property
242    def database(self):
243        """
244        Retrieve the database object that the user has selected.
245        This can be either a standard database or a custom-specified path.
246        """
247        if self.search_db not in all_db_choices:
248            # Take the absolute path #
249            user_path = os.path.abspath(self.search_db)
250            # Make the object #
251            return CrestDatabase(custom_path=user_path)
252        else:
253            return getattr(crest4.databases, self.search_db)
254
255    #------------------------------ Searching --------------------------------#
256    @cached_property
257    def seqsearch(self):
258        """
259        An object representing the sequence similarity search.
260        Makes use of the `seqsearch` module. For reference:
261
262        * Setting `-outfmt` to 5 means XML output.
263        * Setting `-outfmt` to 6 means tabular output.
264        * Setting `-outfmt` to 7 means tabular output with comments.
265        """
266        # Initialize #
267        params = {}
268        # If the user chose BLAST then we have to specify tabular output #
269        if self.search_algo == 'blast':
270            params = {'-outfmt': '7 qseqid sseqid bitscore length nident'}
271        # In case the user chose VSEARCH we specify the minimum identity
272        # and the minimum sequence match length
273        if self.search_algo == 'vsearch':
274            params = {'--id':      self.min_score,
275                      '--mincols': 25}
276        # Build and return the object #
277        from seqsearch.search import SeqSearch
278        return SeqSearch(input_fasta = self.fasta,
279                         database    = self.database,
280                         seq_type    = 'nucl',
281                         algorithm   = self.search_algo,
282                         filtering   = {'max_targets': 100},
283                         num_threads = self.num_threads,
284                         out_path    = self.search_hits,
285                         params      = params)
286
287    def search(self):
288        """A method to launch the sequence similarity search."""
289        # Launch the search algorithm #
290        return self.seqsearch.run()
291
292    #----------------------------- Assigning ---------------------------------#
293    @cached_property
294    def score_frac(self):
295        """
296        Using the parameter `self.score_drop`, which is a percentage (e.g., 2)
297        indicating a drop, we compute the minimum remaining amount of score
298        allowed, as a fraction (e.g., 0.98).
299        """
300        return 1 - (self.score_drop / 100)
301
302    @cached_property
303    def queries(self):
304        """
305        Parses the output of the sequence search program used which returns
306        a list of `QueryResult` objects from the `SearchIO.parse` module.
307        The list contains one Query object per sequence that was originally
308        inputted. Use these objects to access the taxonomic assignments.
309        """
310        # Check if the search has been done already #
311        if not self.search_hits: self.search()
312        # Iterate on the sequence search results #
313        result = [Query(self, query) for query in self.seqsearch.results]
314        # VSEARCH entirely forgets about sequences that had no hits.
315        # Instead of still listing them in the output like BLAST.
316        # So we have to add them back to the list in this awkward manner
317        if self.search_algo == 'vsearch':
318            reported_names = set(query.name for query in result)
319            for seq in self.fasta:
320                if seq.id not in reported_names:
321                    q = type('FakeQuery', (), {'hits': [], 'id': seq.id})
322                    result.append(Query(self, q))
323        # Return #
324        return result
325
326    @cached_property
327    def queries_by_id(self):
328        """
329        References the same Query objects as the `queries` property above,
330        except that this time they are in a dictionary with the query ids
331        (i.e., the original fasta ids) as keys instead of in a list.
332        """
333        return {query.name: query for query in self.queries}
334
335    #------------------------------- Outputs ---------------------------------#
336    @cached_property
337    def out_file(self):
338        """
339        The path to the file that will contain the taxonomic assignments
340        for every sequence.
341        """
342        # Make sure that the output directory exists #
343        self.output_dir.create_if_not_exists()
344        # Return #
345        return self.output_dir + "assignments.txt"
346
347    @cached_property
348    def otu_info(self):
349        """An object giving access to the OTU table information and methods."""
350        from crest4.otu_tables import InfoFromTableOTUs
351        return InfoFromTableOTUs(self, self.otu_table)
352
353    def __call__(self):
354        """Generate outputs."""
355        # Intro message #
356        print('Running crest4 version ' + crest4.__version__)
357        # Iterate #
358        self.out_file.writelines(query.tax_string for query in self.queries)
359        # Special case where an OTU table was passed #
360        if self.otu_table:
361            path_by_rank    = self.output_dir + 'otus_by_rank.tsv'
362            path_cumulative = self.output_dir + 'otus_cumulative.tsv'
363            otus_by_rank    = self.otu_info.otus_by_rank
364            otus_cumulative = self.otu_info.otus_cumulative
365            otus_by_rank.to_csv(path_by_rank, index=False, sep='\t')
366            otus_cumulative.to_csv(path_cumulative, index=False, sep='\t')
367        # Print a success message #
368        msg = "Classification ran successfully. Results are placed in '%s'."
369        print(msg % self.out_file)
370        # Return #
371        return self.out_file