seqsearch.search.vsearch

  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3
  4"""
  5Written by Lucas Sinclair.
  6MIT Licensed.
  7Contact at www.sinclair.bio
  8"""
  9
 10# Built-in modules #
 11
 12# Internal modules #
 13from autopaths.file_path import FilePath
 14from seqsearch.search.core import CoreSearch
 15
 16# Third party modules #
 17from seqsearch import sh
 18
 19###############################################################################
 20class VSEARCHquery(CoreSearch):
 21    """
 22    A vsearch job.
 23
 24    Example commands:
 25
 26        vsearch --usearch_global queries.fas --db queries.fas --id 0.6
 27                --blast6out results.blast6 --maxaccepts 0 --maxrejects 0
 28
 29        ./vsearch --usearch_global queries.fsa --db database.fsa --id 0.9
 30                  --alnout alnout.txt
 31
 32    Some interesting options form the manual follow:
 33
 34    --maxaccepts <positive integer>
 35    > Maximum number of hits to accept before stopping the search. The default
 36    value is 1. This option works in pair with --maxrejects. The search process
 37    sorts target sequences by decreasing number of k-mers they have in common
 38    with the query sequence, using that information as a proxy for sequence
 39    similarity. After pairwise alignments, if the first target sequence passes
 40    the acceptation criteria, it is accepted as best hit and the search process
 41    stops for that query. If --maxaccepts is set to a higher value, more hits
 42    are accepted.
 43    """
 44
 45    extension = 'vsearchout'
 46
 47    @property
 48    def command(self):
 49        # Executable #
 50        if self.executable: cmd = [self.executable.path]
 51        else:               cmd = ['vsearch']
 52        # Other parameters
 53        cmd += ['--usearch_global', self.query,
 54                '--db',             self.db,
 55                '--blast6out',      self.out_path,
 56                '--threads',        self.cpus]
 57        # Options #
 58        for k,v in self.params.items(): cmd += [k, v]
 59        # Return #
 60        return list(map(str, cmd))
 61
 62    #-------------------------------- RUNNING --------------------------------#
 63    def run(self, verbose=False):
 64        """Simply run the VSEARCH search locally."""
 65        # Check the executable is available #
 66        if self.executable:
 67            self.executable.must_exist()
 68        else:
 69            from plumbing.check_cmd_found import check_cmd
 70            check_cmd('vsearch')
 71        # Create the output directory if it doesn't exist #
 72        self.out_path.directory.create_if_not_exists()
 73        # Optionally print the command #
 74        if verbose:
 75            print("Running VSEARCH command:\n    %s" % ' '.join(self.command))
 76        # Run it #
 77        cmd    = sh.Command(self.command[0])
 78        result = cmd(self.command[1:], _out=self._out, _err=self._err)
 79        # Return #
 80        return result
 81
 82    #----------------------------- PARSE RESULTS -----------------------------#
 83    @property
 84    def results(self):
 85        """
 86        Parse the results and yield biopython SearchIO entries.
 87
 88        Beware:
 89        Some databases are not unique on the id, and this causes the parser to
 90        complain about duplicate entries and raise exceptions such as:
 91
 92            ValueError: The ID or alternative IDs of Hit 'DQ448783' exists
 93            in this QueryResult.
 94
 95        Summary of the columns:
 96        https://www.metagenomics.wiki/tools/blast/blastn-output-format-6
 97
 98            qseqid sseqid pident length mismatch gapopen qstart qend sstart send
 99            evalue bitscore
100
101        Warning: Unlike BLAST results, if a sequence got no hits it is NOT
102                 reported at all in VSEARCH. The number of entries yielded
103                 will not match the number of sequences at input.
104        """
105        # Import parsing library #
106        from Bio import SearchIO
107        # Avoid the warning #
108        import warnings
109        warnings.filterwarnings("ignore", 'BiopythonDeprecationWarning')
110        # Iterate over lines #
111        with open(self.out_path, 'rt') as handle:
112            for entry in SearchIO.parse(handle, 'blast-tab'):
113                yield entry
114
115###############################################################################
116class VSEARCHdb(FilePath):
117    """A VSEARCH database one can search against."""
118
119    def __repr__(self):
120        return '<%s on "%s">' % (self.__class__.__name__, self.path)
121
122    def __init__(self, path):
123        # Call parent constructor #
124        FilePath.__init__(self, path)
125        # Check if the corresponding FASTA exists #
126        self.fasta_path = self.replace_extension('fasta')
127
128    def __bool__(self):
129        """Does the indexed database actually exist?"""
130        return bool(self.replace_extension('udb'))
131
132    def create_if_not_exists(self, *args, **kwargs):
133        """If the indexed database has not been generated, generate it."""
134        if not self: return self.makedb(*args, **kwargs)
135
136    def makedb(self, output=None, stdout=None, verbose=False):
137        # We need a corresponding fasta that exists #
138        if not self.fasta_path:
139            raise Exception("No fasta file at '%s'." % self.fasta_path)
140        # Message #
141        if verbose:
142            msg = "Calling `makeudb_usearch` from '%s' to '%s'."
143            print(msg % (self.fasta_path, self))
144        # The output #
145        if output is None: output = self
146        # Options #
147        options = ['--makeudb_usearch', self.fasta_path, '--output', output]
148        # Call the program #
149        sh.vsearch(*options, _out=stdout)
150        # Return #
151        return output