NCBIWWW = create_fake_class(MSG + BIO) def seqio(in_fhands, out_fhand, out_format, copy_if_same_format=True):
    'It converts sequence files between formats'
    if out_format not in get_setting('SUPPORTED_OUTPUT_FORMATS'):
        raise IncompatibleFormatError("This output format is not supported")

    in_formats = [get_format(fhand) for fhand in in_fhands]

    if len(in_fhands) == 1 and in_formats[0] == out_format:
        if copy_if_same_format:
            copyfileobj(in_fhands[0], out_fhand)
        else:
            rel_symlink(in_fhands[0].name, out_fhand.name)
    else:
        seqs = _read_seqrecords(in_fhands)
        try:
            write_seqrecs(seqs, out_fhand, out_format)
        except ValueError, error:
            if error_quality_disagree(error):
                raise MalformedFile(str(error))
            if 'No suitable quality scores' in str(error):
                msg = 'No qualities available to write output file'
                raise IncompatibleFormatError(msg)
            raise
    flush_fhand(out_fhand) def guess_seq_type(fhand):
    '''It guesses if the file is nucleotide or protein'''
    rna = set(ambiguous_rna_letters)
    dna = set(ambiguous_dna_letters)
    rna_dna = rna.union(dna)

    protein = set(extended_protein_letters)
    only_prot = list(protein.difference(rna_dna))

    chunk_size = 1024
    chunk = peek_chunk_from_file(fhand, chunk_size)
    if not chunk:
        raise UnknownFormatError('The file is empty')
    fhand_ = cStringIO.StringIO(chunk)
    total_letters = 0
    nucleotides = 0
    for seq in read_seqs([fhand_]):
        for letter in get_str_seq(seq):
            total_letters += 1
            if letter in ('gcatnuGCATNU'):
                nucleotides += 1
            if letter in only_prot:
                return 'prot'
    nucl_freq = nucleotides / total_letters
    if nucl_freq > 0.8:
        return 'nucl'

    raise RuntimeError('unable to guess the seq type') def write_seqs(seqs, fhand=None, file_format=None):
    'It writes the given sequences'
    if fhand is None:
        fhand = NamedTemporaryFile(suffix='.' + file_format.replace('-', '_'))

    seqs, seqs2 = tee(seqs)
    try:
        seq = seqs2.next()
    except StopIteration:
        # No sequences to write, so we're done
        return fhand
    del seqs2
    seq_class = seq.kind
    if seq_class == SEQITEM:
        _write_seqitems(seqs, fhand, file_format)
    elif seq_class == SEQRECORD:
        seqs = (seq.object for seq in seqs)
        _write_seqrecords(seqs, fhand, file_format)
    else:
        raise ValueError('Unknown class for seq: ' + seq_class)
    return fhand def check_process_finishes(process, binary, stdout=None, stderr=None):
    'It checks that the given process finishes OK, otherwise raises an Error'

    stdout_msg, stderr_msg = process.communicate()

    returncode = process.returncode
    if returncode == 0:
        return
    elif returncode is None:
        msg = 'The process for {:s} is still running with PID'
        msg = msg.format(binary, process.PID)
        raise RuntimeError(msg)

    if stdout and not stdout_msg:
        stdout.flush()
        stdout_msg = open(stdout.name).read()
    if stderr and not stderr_msg:
        stderr.flush()
        stderr_msg = open(stderr.name).read()
    msg = '{:s} had a problem running\n'.format(binary)
    if stdout_msg:
        msg += 'stdout:\n{:s}\n'.format(stdout_msg)
    if stderr_msg:
        msg += 'stderr:\n{:s}\n'.format(stderr_msg)
    raise ExternalBinaryError(msg) def popen(*args, **kwargs):
    'It spawns a process with subprocess.Popen'
    try:
        return Popen(*args, **kwargs)
    except OSError, error:
        msg = 'The binary "%s" is not in the path.' % args[0][0]
        if 'No such file' in str(error):
            raise MissingBinaryError(msg) def get_binary_path(binary_name):
    '''It return the path to the proper binary. It looks on platform and
    architecture to decide it.

    Fails if there is not binary for that architecture
    '''
    if get_setting('USE_EXTERNAL_BIN_PREFIX'):
        ext_binary_name = get_setting('EXTERNAL_BIN_PREFIX') + binary_name
        if os.path.exists(ext_binary_name):
            return ext_binary_name

    if not get_setting('ADD_PATH_TO_EXT_BIN'):
        # I have to check if the binary is on my current directory.
        # If it is there use it, else assumes that it is on the path
        if os.path.exists(os.path.join(os.getcwd(), ext_binary_name)):
            return os.path.join(os.getcwd(), ext_binary_name)
        # return binary_name

    system = platform.system().lower()
    if system == 'windows':
        binary_name += '.exe'
    arch = platform.architecture()[0]

    join = os.path.join

    third_party_path = get_setting('THIRD_PARTY_DIR')

    binary_path = os.path.abspath(join(third_party_path, 'bin', system, arch,
                                       binary_name))
    if os.path.exists(binary_path):
        return binary_path
    elif arch == '64bit':
        arch = '32bit'
        binary_path = os.path.abspath(join(third_party_path, 'bin', system,
                                           arch, binary_name))
        if os.path.exists(binary_path):
            return binary_path

    # At this point there is not available binary for the working platform
    # Is the binary really in the path?
    if which(binary_name):
        return binary_name

    msg = '{} not found in the path. Please install it to use ngs_crumbs'
    raise MissingBinaryError(msg.format(binary_name)) NUCL = 'nucl' PROT = 'prot' def _text_blasts_in_file(fhand):
def _split_description(string):
def _text_blast_parser(fhand):
    def __init__(self, fhand):
    def __iter__(self):
    def next(self):
def _lines_for_every_tab_blast(fhand, line_format):
def _group_match_parts_by_subject(match_parts):
def _tabular_blast_parser(fhand, line_format):
    def __init__(self, fhand, line_format=DEFAULT_TABBLAST_FORMAT):
    def __iter__(self):
    def next(self):
    def __init__(self, fhand, subj_def_as_accesion=None):
    def __iter__(self):
    def _create_result_structure(self, bio_result):
    def _get_blast_metadata(self):
    def next(self):
    def __init__(self, fhand):
    def __iter__(self):
    def _results_query_from_exonerate(self):
    def _create_structure_result(query_result):
    def next(self):
def _strand_transform(strand):
def _match_num_if_exists_in_struc(subject_name, struct_dict):
def get_alignment_parser(kind):
def get_match_score(match, score_key, query=None, subject=None):
def get_match_scores(match, score_keys, query, subject):
def alignment_results_scores(results, scores, filter_same_query_subject=True):
def build_relations_from_aligment(fhand, query_name, subject_name):
def _get_match_score(match, score_key, query=None, subject=None):
def _score_above_threshold(score, min_score, max_score, log_tolerance,
                           log_best_score):
def _create_scores_mapper_(score_key, score_tolerance=None,
                           max_score=None, min_score=None):
    def map_(alignment):
def _create_best_scores_mapper(score_key, score_tolerance=None,
                              max_score=None, min_score=None):
def _create_scores_mapper(score_key, max_score=None, min_score=None):
def _create_deepcopy_mapper():
    def map_(alignment):
def _create_empty_filter():
    def filter_(alignment):
def _fix_match_scores(match, score_keys):
def _fix_match_start_end(match):
def _fix_matches(alignment, score_keys=None):
def _create_fix_matches_mapper():
def covered_segments_from_match_parts(match_parts, in_query=True,
                                      merge_segments_closer=1):
def elongate_match_part_till_global(match_part, query_length, subject_length,
                                    align_completely):
def elongate_match_parts_till_global(match_parts, query_length,
                                     subject_length, align_completely):
def _match_length(match, length_from_query):
def _match_part_length(match_part, length_in_query):
def _match_long_enough(match_length, total_length, min_num_residues,
                       min_percentage, length_in_query):
def _create_min_length_mapper(length_in_query, min_num_residues=None,
                              min_percentage=None, filter_match_parts=False):
    def map_(alignment):
def filter_alignments(alignments, config):
class TextBlastParser(object):
class TabularBlastParser(object):
class BlastParser(object):
class ExonerateParser(object):
DEFAULT_TABBLAST_FORMAT = ('query', 'subject', 'identity', 'alignment_length',
                           'mismatches', 'gap_open', 'query_start',
                           'query_end', 'subject_start', 'subject_end',
                           'expect', 'score')
MAPPER = 1
FILTER = 2
FILTER_COLLECTION = {'best_scores':
                                 {'funct_factory': _create_best_scores_mapper,
                                  'kind': MAPPER},
                     'score_threshold':
                                 {'funct_factory': _create_scores_mapper,
                                 'kind': MAPPER},
                     'min_length': {'funct_factory': _create_min_length_mapper,
                                    'kind': MAPPER},
                     'deepcopy': {'funct_factory': _create_deepcopy_mapper,
                                  'kind': MAPPER},
                     'fix_matches':
                                 {'funct_factory': _create_fix_matches_mapper,
                                  'kind': MAPPER},
                     'filter_empty':
                                 {'funct_factory': _create_empty_filter,
                                  'kind': FILTER},
                    } 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]
import os.path
import subprocess
import tempfile
from crumbs.utils.optional_modules import NCBIWWW
from crumbs.seq.seqio import seqio, guess_seq_type, write_seqs
from crumbs.utils.bin_utils import (check_process_finishes, popen,
                                    get_binary_path)
from crumbs.utils.tags import NUCL, PROT
from crumbs.seq.alignment_result import (filter_alignments, ELONGATED, QUERY,
                                         covered_segments_from_match_parts,
                                         elongate_match_parts_till_global,
                                         TabularBlastParser, BlastParser)
from crumbs.utils.file_utils import TemporaryDir
from crumbs.settings import get_setting
    'It does a blast'
    if not params:
        params = {}
    evalue, task = _parse_blast_params(params, program)

    if 'outfmt' in params:
        outfmt = params['outfmt']
        del params['outfmt']
    else:
        outfmt = 5

    if program not in ('blastn', 'blastp', 'blastx', 'tblastx', 'tblastn'):
        raise ValueError('The given program is invalid: ' + str(program))
    binary = get_binary_path(program)
    cmd = [binary, '-query', query_fpath, '-db', db_fpath, '-out', out_fpath]
    cmd.extend(['-evalue', str(evalue), '-outfmt', str(outfmt)])
    if task:
        cmd.extend(['-task', task])
    if params:
        for key, value in params.viewitems():
            cmd.extend(('-' + key, str(value)))
    process = popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    check_process_finishes(process, binary=cmd[0])


def _do_blast_2(db_fpath, queries, program, dbtype=None, blast_format=None,
                params=None, remote=False):
    '''It returns an alignment result with the blast.

    It is an alternative interface to the one based on fpaths.
    db_fpath should be a plain sequence file.
    queries should be a SeqRecord list.
    If an alternative blast output format is given it should be tabular, so
    blast_format is a list of fields.
    '''

    query_fhand = write_seqs(queries, file_format='fasta')
    query_fhand.flush()

    if remote:
        blastdb = db_fpath
        fmt = 'XML' if blast_format is None else blast_format.upper()
    else:
        blastdb = get_or_create_blastdb(db_fpath, dbtype=dbtype)
        if blast_format is None:
            blast_format = ['query', 'subject', 'query_length', 'subject_length',
                            'query_start', 'query_end', 'subject_start',
                            'subject_end', 'expect', 'identity', ]
        fmt = generate_tabblast_format(blast_format)

    if params is None:
        params = {}
    params['outfmt'] = fmt

    blast_fhand = tempfile.NamedTemporaryFile(suffix='.blast')
    do_blast(query_fhand.name, blastdb, program, blast_fhand.name, params,
             remote=remote)
    if remote:
        blasts = BlastParser(blast_fhand)
    else:
        blasts = TabularBlastParser(blast_fhand, blast_format)

    return blasts, blast_fhand


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)