def do_blast(query_fpath, db_fpath, program, out_fpath, params=None,
             remote=False):
    if remote:
        _do_blast_remote(query_fpath, db_fpath, program, out_fpath,
                         params=params)
    else:
        _do_blast_local(query_fpath, db_fpath, program, out_fpath,
                        params=params) class BlasterForFewSubjects(object):
    '''It matches the given SeqRecords against the reads in the file.

    This class uses Blast to do the matching and it is optimized for having
    few SeqRecords to match against a medium sized sequence file.

    In this class the seqs_fpath, that in the API acts as the query, it's
    internally indexed into a blast database and the seqrecords are blasted
    against it. This is done because in the case of having lots of queries to
    blast against few subjects it is more efficient to do the blast backwards.
    This class behaves as if the seqs_fpath were the query and the seqrecords
    the subject although internally the blast is done with the query and
    subject changed.
    '''
    def __init__(self, seqs_fpath, seqs, program, params=None,
                 filters=None, elongate_for_global=False, seqs_type=None):
        '''It inits the class.'''
        self.program = program
        if params is None:
            params = {}
        params['max_target_seqs'] = str(get_setting('PACKET_SIZE'))
        self.params = params
        if filters is None:
            filters = []
        self.filters = filters
        self.elongate_for_global = elongate_for_global
        self._match_parts = self._look_for_blast_matches(seqs_fpath,
                                                         seqs,
                                                         seqs_type)

    def _look_for_blast_matches(self, seq_fpath, oligos, seqs_type):
        'It looks for the oligos in the given sequence files'
        # we need to keep the blast_fhands, because they're temp files and
        # otherwise they might be removed
        temp_dir = TemporaryDir()
        dbpath = os.path.join(temp_dir.name, os.path.basename(seq_fpath))
        seqio([open(seq_fpath)], open(dbpath, 'w'), out_format='fasta',
              copy_if_same_format=False)

        blasts, blast_fhand = _do_blast_2(dbpath, oligos, params=self.params,
                                          program=self.program,
                                          dbtype=seqs_type)
        if self.filters is not None:
            blasts = filter_alignments(blasts, config=self.filters)

        # Which are the regions covered in each sequence?
        indexed_match_parts = {}
        for blast in blasts:
            oligo = blast['query']
            for match in blast['matches']:
                read = match['subject']
                if self.elongate_for_global:
                    elongate_match_parts_till_global(match['match_parts'],
                                                 query_length=oligo['length'],
                                                 subject_length=read['length'],
                                                 align_completely=QUERY)

                # match_parts = [m['match_parts'] for m in blast['matches']]
                match_parts = match['match_parts']
                try:
                    indexed_match_parts[read['name']].extend(match_parts)
                except KeyError:
                    indexed_match_parts[read['name']] = match_parts

        temp_dir.close()
        blast_fhand.close()
        return indexed_match_parts

    def get_matched_segments_for_read(self, read_name):
        'It returns the matched segments for any oligo'
        setting_key = 'DEFAULT_IGNORE_ELONGATION_SHORTER'
        ignore_elongation_shorter = get_setting(setting_key)

        try:
            match_parts = self._match_parts[read_name]
        except KeyError:
            # There was no match in the blast
            return None

        # Any of the match_parts has been elongated?
        elongated_match = False
        for m_p in match_parts:
            if ELONGATED in m_p and m_p[ELONGATED] > ignore_elongation_shorter:
                elongated_match = True
        segments = covered_segments_from_match_parts(match_parts,
                                                     in_query=False)
        return segments, elongated_match def get_or_create_blastdb(blastdb_or_path, dbtype=None, directory=None):
    '''it returns a blast database.

    If it does not exists it creates and if you give it a directory it will
    create it in that directory if it does not exist yet
    '''
    seq_fpath = _get_abs_blastdb_path(blastdb_or_path, dbtype)
    if directory:
        dbname = os.path.basename(seq_fpath)
        dbpath = os.path.join(directory, dbname)
    else:
        dbpath = seq_fpath

    if not _blastdb_exists(dbpath, dbtype=dbtype):
        if not os.path.exists(seq_fpath):
            msg = 'An input sequence is required to create a blastdb'
            raise RuntimeError(msg)
        if seq_fpath != dbpath:
            seqio([open(seq_fpath)], open(dbpath, 'w'), out_format='fasta',
                  copy_if_same_format=False)
        if dbtype is None:
            dbtype = guess_seq_type(open(dbpath))
        _makeblastdb_plus(dbpath, dbtype)
    return dbpath def _blastdb_exists(dbpath, dbtype=None):
    'It checks if a a blast db exists giving its path'
    assert dbtype in (NUCL, PROT, None)
    if dbtype == NUCL:
        exts = ['.nin', '.00.nin']
    elif dbtype == PROT:
        exts = ['.pin', '.00.pin']
    else:
        exts = ['.nin', '.pin', '.00.nin', '.00.pin']
    return any([os.path.exists(dbpath + ext) for ext in exts]) class Blaster(object):
    '''It matches the given SeqRecords against a blast database.
    It needs iterable with seqrecords and a blast dabatase
    '''

    def __init__(self, seqrecords, blastdb, program, dbtype=None, params=None,
                 filters=None, remote=False):
        self.program = program
        if params is None:
            params = {}
        self.params = params
        if filters is None:
            filters = []
        self.filters = filters
        self._remote = remote
        if dbtype not in (NUCL, PROT, None):
            raise ValueError('dbtype must be NUCL, PROT or None (we guess)')
        self._blasts = self._look_for_blast_matches(seqrecords, blastdb,
                                                    dbtype)

    def _look_for_blast_matches(self, seqrecords, blastdb, dbtype):
        'it makes the blast and filters the results'
        blasts, blast_fhand = _do_blast_2(blastdb, seqrecords, self.program,
                                          params=self.params, dbtype=dbtype,
                                          remote=self._remote)
        # print open(blast_fhand.name).read()
        if self.filters is not None:
            blasts = filter_alignments(blasts, config=self.filters)

        blasts = {blast['query']['name']: blast for blast in blasts}
        blast_fhand.close()
        return blasts

    def get_matched_segments(self, seqrecord_name):
        'It returns the matched segments for the given query (SeqRecord) name.'
        try:
            blast = self._blasts[seqrecord_name]
        except KeyError:
            # There was no match in the blast
            return None

        match_parts = [mp for m in blast['matches'] for mp in m['match_parts']]
        segments = covered_segments_from_match_parts(match_parts,
                                                     in_query=False)
        return segments

    @property
    def blasts(self):
        'It returns the blasts alignment results.'
        # we trust the clients not to mess with this dict or its nested
        # structure, to deepcopy would be a waste of resources in a
        # well-behaved world
        return self._blasts class TemporaryDir(object):
    '''This class creates temporary directories '''
    def __init__(self, suffix='', prefix='', directory=None):
        '''It initiates the class.'''
        self._name = tempfile.mkdtemp(suffix=suffix, prefix=prefix,
                                      dir=directory)

    @property
    def name(self):
        'Returns path to the dict'
        return self._name

    def close(self):
        '''It removes the temp dir'''
        if os.path.exists(self._name):
            shutil.rmtree(self._name) def get_setting(key):
    'It returns the value for one setting'
    return _settings[key] TEST_DATA_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), '..',
                                             '..', 'test', 'test_data')) NUCL = 'nucl' SEQITEM = 'seqitem' SEQRECORD = 'seqrecord' class SeqItem(_SeqItem):
    def __new__(cls, name, lines, annotations=None):
def get_title(seq):
def get_description(seq):
def get_name(seq):
def get_file_format(seq):
def _break():
def _is_fastq_plus_line(line, seq_name):
def _get_seqitem_quals(seq):
def get_str_seq(seq):
def get_length(seq):
def _get_seqitem_qualities(seqwrap):
def get_int_qualities(seq):
def _int_quals_to_str_quals(int_quals, out_format):
def get_str_qualities(seq, out_format=None):
def get_annotations(seq):
def _copy_seqrecord(seqrec, seq=None, name=None, id_=None):
def _copy_seqitem(seqwrapper, seq=None, name=None):
def copy_seq(seqwrapper, seq=None, name=None):
def _slice_seqitem(seqwrap, start, stop):
def slice_seq(seq, start=None, stop=None):
def assing_kind_to_seqs(kind, seqs, file_format):
SANGER_QUALS = {chr(i): i - 33 for i in range(33, 127)}
ILLUMINA_QUALS = {chr(i): i - 64 for i in range(64, 127)}
SANGER_STRS = {i - 33: chr(i)  for i in range(33, 127)}
ILLUMINA_STRS = {i - 64: chr(i) for i in range(64, 127)}
import unittest
import os.path
from tempfile import NamedTemporaryFile
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq
from crumbs.blast import (do_blast, BlasterForFewSubjects,
                          get_or_create_blastdb, _blastdb_exists, Blaster)
from crumbs.utils.file_utils import TemporaryDir
from crumbs.settings import get_setting
from crumbs.utils.test_utils import TEST_DATA_DIR
from crumbs.utils.tags import NUCL, SEQITEM, SEQRECORD
from crumbs.seq.seq import SeqWrapper, SeqItem, assing_kind_to_seqs

# You should have received a copy of the GNU General Public License
# along with ngs_crumbs. If not, see <http://www.gnu.org/licenses/>.




TITANIUM_LINKER = get_setting('TITANIUM_LINKER')
FLX_LINKER = get_setting('FLX_LINKER')

# pylint: disable=R0201
# pylint: disable=R0904
# pylint: disable=C0111


class BlastTest(unittest.TestCase):
    'It tests the blast infrastructure'

    @staticmethod
    def test_blastdb():
        'It creates a blast database.'
        db_name = 'arabidopsis_genes'
        seq_fpath = os.path.join(TEST_DATA_DIR, db_name)
        db_dir = TemporaryDir(prefix='blast_dbs_')
        try:
            db_path1 = get_or_create_blastdb(seq_fpath, directory=db_dir.name,
                                    dbtype='nucl')
            db_path = os.path.join(db_dir.name, db_name)
            assert 'CATAGGGTCACCAATGGC' in open(db_path1).read(100)
            assert db_path1 == db_path
            assert os.path.exists(db_path)
            index_fpath = os.path.join(db_dir.name, db_name + '.nsq')
            assert os.path.exists(index_fpath)
        finally:
            db_dir.close()

    def test_blast_search(self):
        'It does a blast search'
        db_name = 'arabidopsis_genes'
        seq_fpath = os.path.join(TEST_DATA_DIR, db_name)
        db_dir = TemporaryDir(prefix='blast_dbs_')
        try:
            db_fpath = get_or_create_blastdb(seq_fpath, directory=db_dir.name,
                                              dbtype='nucl')
            query_fhand = NamedTemporaryFile()
            query_fhand.write(open(seq_fpath).read(200))
            query_fhand.flush()
            out_fhand = NamedTemporaryFile()
            do_blast(seq_fpath, db_fpath, program='blastn',
                     out_fpath=out_fhand.name)
            assert '</BlastOutput>' in open(out_fhand.name).read()
        finally:
            db_dir.close()

    def xtest_remote_blast(self):
        'It does a remote blast search'
        seq_fhand = NamedTemporaryFile()
        fasta = '>seq1\nCTAGTCTAGTCGTAGTCATGGCTGTAGTCTAGTCTACGATTCGTATCAGTTGT'
        fasta += 'GTGACATCGATCATGTTGTATTGTGTACTATACACACACGTAGGTCGACTATCGTAGC\n'
        seq_fhand.write(fasta)
        seq_fhand.flush()
        out_fhand = NamedTemporaryFile()
        do_blast(seq_fhand.name, 'nt', program='blastn',
                 out_fpath=out_fhand.name, remote=True)
        assert '</BlastOutput>' in open(out_fhand.name).read()

        # fail if outfmt is not xml
        try:
            params = {'outfmt': 'txt'}
            do_blast(seq_fhand.name, 'nt', program='blastn',
                out_fpath=out_fhand.name, remote=True, params=params)
            self.fail()
        except RuntimeError:
            pass

    @staticmethod
    def test_get_or_create_blastdb():
        blastdb = os.path.join(TEST_DATA_DIR, 'arabidopsis_genes')

        directory = TemporaryDir()
        assert not _blastdb_exists(blastdb, NUCL)
        get_or_create_blastdb(blastdb, NUCL, directory.name)
        new_blast_path = os.path.join(directory.name,
                                      os.path.basename(blastdb))
        assert _blastdb_exists(new_blast_path, NUCL)
        get_or_create_blastdb(blastdb, NUCL, directory.name)
        assert _blastdb_exists(new_blast_path, NUCL)
        directory.close()

        # already exists
        blastdb = os.path.join(TEST_DATA_DIR, 'blastdbs', 'arabidopsis_genes')
        assert _blastdb_exists(blastdb, NUCL)
        get_or_create_blastdb(blastdb, NUCL)
        assert _blastdb_exists(blastdb, NUCL)


def create_a_matepair_file():
    'It creates a matepair fasta file'

    seq_5 = 'CTAGTCTAGTCGTAGTCATGGCTGTAGTCTAGTCTACGATTCGTATCAGTTGTGTGAC'
    seq_3 = 'ATCGATCATGTTGTATTGTGTACTATACACACACGTAGGTCGACTATCGTAGCTAGT'
    mate_seq = seq_5 + TITANIUM_LINKER + seq_3
    mate_fhand = NamedTemporaryFile(suffix='.fasta')
    mate_fhand.write('>seq1\n' + mate_seq + '\n')
    mate_fhand.flush()
    return mate_fhand


class BlastMaterTest(unittest.TestCase):
    'It tests the splitting of mate pairs'
    def test_matching_segments(self):
        'It tests the detection of oligos in sequence files'
        seq_5 = 'CTAGTCTAGTCGTAGTCATGGCTGTAGTCTAGTCTACGATTCGTATCAGTTGTGTGAC'
        mate_fhand = create_a_matepair_file()

        linkers = [SeqItem('titan', ['>titan\n', TITANIUM_LINKER + '\n']),
                   SeqItem('flx', ['>flx\n', FLX_LINKER + '\n'])]
        linkers = assing_kind_to_seqs(SEQITEM, linkers, 'fasta')

        expected_region = (len(seq_5), len(seq_5 + TITANIUM_LINKER) - 1)

        matcher = BlasterForFewSubjects(mate_fhand.name, linkers,