MSG = 'A python package to run this executable is required,'
BIO = 'biopython'
BIO_BGZF = 'biopython with Bgzf support'
    NCBIXML = create_fake_class(MSG + BIO)
    NCBIWWW = create_fake_class(MSG + BIO)
def create_fake_class(msg):
        def __init__(self, *args, **kwargs):
def create_fake_funct(msg):
    def FakeRequiredfunct(*args, **kwargs):
    class FakePythonRequiredClass(object): TRIMMING_RECOMMENDATIONS = 'trimming_recommendations' QUALITY = 'quality' OTHER = 'other' VECTOR = 'vector' TRIMMING_KINDS = [VECTOR, QUALITY, OTHER] SEQS_PASSED = 'seqs_passed' ORPHAN_SEQS = 'orphan_seqs' def get_uppercase_segments(string):
    '''It detects the unmasked regions of a sequence

    It returns a list of (start, end) tuples'''
    start = 0
    for is_upper, group in itertools.groupby(string, lambda x: x.isupper()):
        group = list(group)
        end = start + len(group) - 1
        if is_upper:
            yield start, end
        start = end + 1 def copy_seq(seqwrapper, seq=None, name=None):
    seq_class = seqwrapper.kind
    seq_obj = seqwrapper.object
    if seq_class == SEQITEM:
        seq = _copy_seqitem(seqwrapper, seq=seq, name=name)
    elif seq_class == SEQRECORD:
        seq_obj = _copy_seqrecord(seq_obj, seq=seq, name=name, id_=name)
        seq = SeqWrapper(kind=seqwrapper.kind, object=seq_obj,
                         file_format=seqwrapper.file_format)
    return seq def get_str_seq(seq):
    seq_class = seq.kind
    if seq_class == SEQITEM:
        seq = seq.object.lines[1].strip()
    elif seq_class == SEQRECORD:
        seq = str(seq.object.seq)
    return seq.strip() def get_annotations(seq):
    return seq.object.annotations def get_length(seq):
    return len(get_str_seq(seq)) def slice_seq(seq, start=None, stop=None):
    seq_class = seq.kind
    if seq_class == SEQITEM:
        seq_obj = _slice_seqitem(seq, start, stop)
    elif seq_class == SEQRECORD:
        seq_obj = seq.object[start:stop]
    return SeqWrapper(seq.kind, object=seq_obj, file_format=seq.file_format) def get_int_qualities(seq):
    seq_class = seq.kind
    if seq_class == SEQITEM:
        return _get_seqitem_qualities(seq)
    elif seq_class == SEQRECORD:
        try:
            quals = seq.object.letter_annotations['phred_quality']
        except KeyError:
            msg = 'The given SeqRecord has no phred_quality'
            raise AttributeError(msg)
        return quals def get_name(seq):
    if 'SeqRecord' in seq.__class__.__name__:
        seq_class = SEQRECORD
    else:
        seq_class = seq.kind
        seq = seq.object
    if seq_class == SEQITEM:
        name = seq.name
    elif seq_class == SEQRECORD:
        name = seq.id
    return name def get_longest_segment(segments):
    'It returns the longest segment'
    longest = None
    longest_size = None
    for segment in segments:
        size = segment[1] - segment[0]
        if longest is None:
            longest = [segment]
            longest_size = size
        elif size > longest_size:
            longest = [segment]
            longest_size = size
        elif size == longest_size:
            longest.append(segment)

    if longest is None:
        return None
    elif len(longest) == 1:
        return longest[0]
    else:
        return random.choice(longest) def get_all_segments(segments, seq_len):
    '''Given a set of some non overlapping regions it returns all regions.

    input:  --- ----       ----
    output: ---+----+++++++----

    '''
    segments_ = copy.deepcopy(segments)
    if not segments_:
        return [((0, seq_len - 1), False)]

    all_segments = []
    # If the first segment doesn't start at zero we create a new one
    if segments_[0][0] == 0:
        start = segments_[0][1] + 1
        all_segments.append((segments_.pop(0), True))
    else:
        start = 0

    for loc in segments_:
        all_segments.append(((start, loc[0] - 1), False))
        all_segments.append((loc, True))
        start = loc[1] + 1
    else:
        # if the last segment does not ends at the end of the sequence we add
        # an extra one
        end = seq_len - 1
        if start <= end:
            all_segments.append(((start, end), False))
    return all_segments def get_longest_complementary_segment(segments, seq_len):
    'Given a seq and the locations it returns the longest region not covered'

    segments = merge_overlaping_segments(segments)

    # we want the non-matched locations
    segments = get_complementary_segments(segments, seq_len)

    # now we return the longest one
    return get_longest_segment(segments) def merge_overlaping_segments(segments, merge_segments_closer=1):
    '''Given a list of segments it returns the overlapped segments.

       segments 1  -------        ----->    -----------
       segments 2       ------
       returns     -----------    ------    -----------
       merge_segments_closer is an integer. Segments closer than the given
       number of residues will be merged.
    '''

    # we collect all start and ends
    limits = []
    for segment in segments:
        start = segment[0]
        end = segment[1]
        if start > end:  # a reversed item
            start, end = end, start
        limit_1 = (START, start)
        limit_2 = (END, end)
        limits.append(limit_1)
        limits.append(limit_2)

    # sort by secondary key: start before end
    limits.sort(key=itemgetter(0))
    # sort by location (primary key)
    limits.sort(key=itemgetter(1))

    # merge the ends and start that differ in only one base
    filtered_limits = []
    previous_limit = None
    for limit in limits:
        if previous_limit is None:
            previous_limit = limit
            continue
        if (previous_limit[0] == END and limit[0] == START and
            previous_limit[1] >= limit[1] - merge_segments_closer):
            # These limits cancelled each other
            previous_limit = None
            continue
        filtered_limits.append(previous_limit)
        previous_limit = limit
    else:
        filtered_limits.append(limit)
    limits = filtered_limits

    # now we create the merged hsps
    starts = 0
    segments = []
    for limit in limits:
        if limit[0] == START:
            starts += 1
            if starts == 1:
                segment_start = limit[1]
        elif limit[0] == END:
            starts -= 1
            if starts == 0:
                segment = (segment_start, limit[1])
                segments.append(segment)
    return segments SEQRECORD = 'seqrecord' def rolling_window(iterator, window, step=1):
    'It yields lists of items with a window number of elements'
    try:
        length_ = len(iterator)
    except TypeError:
        length_ = None
    if length_ is None:
        return _rolling_window_iter(iterator, window, step)
    else:
        return _rolling_window_serie(iterator, window, length_, step) 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 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 group_pairs_by_name(seqs, all_pairs_same_n_seqs=False):
    paired_seqs = []
    prev_name = None
    n_seqs_per_pair = None
    for seq in iter(seqs):
        try:
            name = _parse_pair_direction_and_name(seq)[0]
        except PairDirectionError:
            name = None
        if name is None or (paired_seqs and name != prev_name):
            if all_pairs_same_n_seqs:
                if n_seqs_per_pair is None:
                    n_seqs_per_pair = len(paired_seqs)
                elif n_seqs_per_pair != len(paired_seqs):
                    msg = 'Pair had different number of reads: '
                    msg += prev_name
                    raise InterleaveError(msg)
            if paired_seqs:
                yield paired_seqs
                paired_seqs = []
        paired_seqs.append(seq)
        prev_name = name
    if paired_seqs:
        if all_pairs_same_n_seqs and n_seqs_per_pair is not None:
            if n_seqs_per_pair != len(paired_seqs):
                msg = 'Pair had different number of reads: '
                msg += prev_name
                raise InterleaveError(msg)
        yield paired_seqs def group_pairs(seqs, n_seqs_in_pair=None, check_all_same_n_seqs=True,
                check_name_matches=True):

    seqs = iter(seqs)
    if n_seqs_in_pair is None:
        first_pair, next_read = _get_first_pair_by_name(seqs)
        if first_pair is None:
            n_seqs_in_pair = None
        else:
            yield first_pair
            n_seqs_in_pair = len(first_pair)
            seqs = chain([next_read], seqs)

    if n_seqs_in_pair == 1:
        # No need to check anything, a pair cannot have less than one read
        # or more than one name
        check_all_same_n_seqs = False
        check_name_matches = False

    if n_seqs_in_pair:
        pairs = group_in_packets_fill_last(seqs, packet_size=n_seqs_in_pair)
        for pair in pairs:
            pair = filter(lambda seq: seq is not None, pair)
            if check_all_same_n_seqs and n_seqs_in_pair != len(pair):
                msg = 'The last pair has fewer reads'
                raise InterleaveError(msg)
            if check_name_matches:
                _check_name_and_direction_match(*pair)
            yield pair def get_setting(key):
    'It returns the value for one setting'
    return _settings[key] def _split_mates(alignments_group):
    # We add a tag to differenciate between mates afterwards. This tag should
    # be probably taken from the file
    forwards = []
    reverses = []
    for alignment_read in alignments_group:
        if alignment_read.is_read1:
            alignment_read.qname += ' 1:N:0:GATCAG'
            forwards.append(alignment_read)
        elif alignment_read.is_read2:
            alignment_read.qname += ' 2:N:0:GATCAG'
            reverses.append(alignment_read)
    return [forwards, reverses] def _get_primary_alignment(alignments_group):
    for alignment in alignments_group:
        if not alignment.is_secondary:
            return alignment def _read_is_totally_mapped(alignments_group, max_clipping):
    for alignment_read in alignments_group:
        if alignment_read.is_unmapped:
            return False
        cigar = alignment_read.cigar
        max_clipping_positions = max_clipping * alignment_read.alen
        if _count_cigar_char(cigar, [4, 5]) <= max_clipping_positions:
            return True
    return False def _get_qstart(aligned_read):
    positions = 0
    for element in aligned_read.cigar:
        if element[0] == 0:
            return positions
        positions += element[1] def _get_qend(aligned_read):
    qstart = _get_qstart(aligned_read)
    end = qstart + aligned_read.alen
    for element in aligned_read.cigar:
        if element[0] == 1:
            end += element[1]
    return end def _5end_mapped(aligned_read, max_clipping):
    if aligned_read.is_unmapped:
        return False
    aligned_read_qend = _get_qend(aligned_read)
    aligned_read_qstart = _get_qstart(aligned_read)
    max_clipping_positions = max_clipping * aligned_read.alen
    if aligned_read.is_reverse:
        return aligned_read_qend > (_get_qlen(aligned_read) -
                                    max_clipping_positions)
    else:
        return aligned_read_qstart < max_clipping_positions def _group_alignments_reads_by_qname(sorted_alignments):
    prev_name = None
    group = None
    for alignment in sorted_alignments:
        name = alignment.qname
        if prev_name is None:
            group = [alignment]
        elif prev_name == name:
            group.append(alignment)
        else:
            yield group
            group = [alignment]
        prev_name = name
    yield group def alignedread_to_seqitem(aligned_read, start_pos=0, end_pos=None):
    if aligned_read is None or aligned_read.seq is None:
        return None
    name = aligned_read.qname
    seq = aligned_read.seq[start_pos: end_pos]
    quals = aligned_read.qual
    if aligned_read.is_reverse:
        seq = _reverse(_complementary(seq))
    if quals is None:
        lines = ['>' + name + '\n', seq + '\n']
        file_format = 'fasta'
    else:
        quals = quals[start_pos: end_pos]
        if aligned_read.is_reverse:
            quals = _reverse(quals)
        lines = ['@' + name + '\n', seq + '\n', '+\n', quals + '\n']
        file_format = 'fastq'
    return SeqWrapper(SEQITEM, SeqItem(name, lines), file_format) def map_with_bwamem(index_fpath, unpaired_fpath=None, paired_fpaths=None,
                    interleave_fpath=None, threads=None, log_fpath=None,
                    extra_params=None, readgroup=None):
    'It maps with bwa mem algorithm'
    interleave = False
    num_called_fpaths = 0
    in_fpaths = []
    if unpaired_fpath is not None:
        num_called_fpaths += 1
        in_fpaths.append(unpaired_fpath)
    if paired_fpaths is not None:
        num_called_fpaths += 1
        in_fpaths.extend(paired_fpaths)
    if interleave_fpath is not None:
        num_called_fpaths += 1
        in_fpaths.append(interleave_fpath)
        interleave = True

    if num_called_fpaths == 0:
        raise RuntimeError('At least one file to map is required')
    if num_called_fpaths > 1:
        msg = 'Bwa can not map unpaired and unpaired reads together'
        raise RuntimeError(msg)

    if extra_params is None:
        extra_params = []

    if '-p' in extra_params:
        extra_params.remove('-p')

    if interleave:
        extra_params.append('-p')

    if readgroup is not None:
        rg_str = r"@RG\tID:{ID}\tSM:{SM}\tPL:{PL}\tLB:{LB}".format(**readgroup)
        extra_params.extend(['-R', rg_str])

    binary = get_binary_path('bwa')
    cmd = [binary, 'mem', '-t', str(get_num_threads(threads)), index_fpath]
    cmd.extend(extra_params)
    cmd.extend(in_fpaths)

    if log_fpath is None:
        stderr = NamedTemporaryFile(suffix='.stderr')
    else:
        stderr = open(log_fpath, 'w')
    # raw_input(' '.join(cmd))
    bwa = popen(cmd, stderr=stderr, stdout=PIPE)
    return bwa def map_process_to_sortedbam(map_process, out_fpath, key='coordinate',
                             stderr_fhand=None, tempdir=None):
    if stderr_fhand is None:
        stderr = NamedTemporaryFile(suffix='.stderr')
    else:
        stderr = stderr_fhand

    if tempdir is None:
        tempdir = tempfile.gettempdir()
    picard_jar = get_setting("PICARD_JAR")
    cmd = ['java', '-jar', picard_jar, 'SortSam', 'I=/dev/stdin',
           'O=' + out_fpath, 'SO=' + key, 'TMP_DIR=' + tempdir,
           'VALIDATION_STRINGENCY=LENIENT']
    sort = popen(cmd, stdin=map_process.stdout, stderr=stderr)
    map_process.stdout.close()
    sort.communicate()

    if map_process.returncode:
        raise RuntimeError('Error in mapping process')

    if sort.returncode:
        raise RuntimeError('Error in Sort process')
import sys
import subprocess
from operator import itemgetter
from tempfile import NamedTemporaryFile
from crumbs.utils.optional_modules import Seq, AlignmentFile
from crumbs.utils.tags import (TRIMMING_RECOMMENDATIONS, QUALITY, OTHER,
                               VECTOR, TRIMMING_KINDS, SEQS_PASSED,
                               ORPHAN_SEQS)
from crumbs.seq.utils.seq_utils import get_uppercase_segments
from crumbs.seq.seq import (copy_seq, get_str_seq, get_annotations, get_length,
                            slice_seq, get_int_qualities, get_name)
from crumbs.utils.segments_utils import (get_longest_segment, get_all_segments,
                                         get_longest_complementary_segment,
                                         merge_overlaping_segments)
from crumbs.utils.tags import SEQRECORD
from crumbs.iterutils import rolling_window
from crumbs.blast import BlasterForFewSubjects
from crumbs.seq.seqio import write_seqs
from crumbs.seq.pairs import group_pairs_by_name, group_pairs
from crumbs.settings import get_setting
from crumbs.seq.mate_chimeras import (_split_mates, _get_primary_alignment,
                                      _read_is_totally_mapped, _get_qstart,
                                      _get_qend, _5end_mapped,
                                      _group_alignments_reads_by_qname)
from crumbs.mapping import (alignedread_to_seqitem, map_with_bwamem,
                            map_process_to_sortedbam)
# This file is part of ngs_crumbs.
# ngs_crumbs is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.

# ngs_crumbs is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR  PURPOSE.  See the
# GNU General Public License for more details.

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


# pylint: disable=R0903


def seq_to_trim_packets(seq_packets, group_paired_reads=False):
    'It yields packets suitable for the filters'

    for packet in seq_packets:
        if group_paired_reads:
            packet = list(group_pairs_by_name(packet))
        else:
            packet = list(group_pairs(packet, n_seqs_in_pair=1))
        yield {SEQS_PASSED: packet, ORPHAN_SEQS: []}


class _BaseTrim(object):
    'Base Trim class'
    def __call__(self, trim_packet):
        'It trims the seqs'
        self._pre_trim(trim_packet)
        trimmed_seqs = []
        for paired_seqs in trim_packet[SEQS_PASSED]:
            trimmed_seqs.append([self._do_trim(s) for s in paired_seqs])
        self._post_trim()
        return {SEQS_PASSED: trimmed_seqs,
                ORPHAN_SEQS: trim_packet[ORPHAN_SEQS]}

    def _do_trim(self, seq):
        raise NotImplementedError()

    def _pre_trim(self, trim_packet):
        pass

    def _post_trim(self):
        pass


class TrimLowercasedLetters(_BaseTrim):
    'It trims the masked segments of the seqrecords.'

    def _do_trim(self, seq):
        str_seq = get_str_seq(seq)
        unmasked_segments = get_uppercase_segments(str_seq)
        segment = get_longest_segment(unmasked_segments)
        if segment is not None:
            segments = []
            if segment[0] != 0:
                segments.append((0, segment[0] - 1))
            len_seq = len(str_seq)
            if segment[1] != len_seq - 1:
                segments.append((segment[1] + 1, len_seq - 1))

            _add_trim_segments(segments, seq, kind=OTHER)

        else:
            segments = [(0, len(seq))]
            _add_trim_segments(segments, seq, kind=OTHER)
        return seq


def _add_trim_segments(segments, sequence, kind):
    'It adds segments to the trimming recommendation in the annotation'
    assert kind in TRIMMING_KINDS
    if not segments:
        return
    annotations = sequence.object.annotations
    if TRIMMING_RECOMMENDATIONS not in annotations:
        annotations[TRIMMING_RECOMMENDATIONS] = {}
        for trim_kind in TRIMMING_KINDS:
            annotations[TRIMMING_RECOMMENDATIONS][trim_kind] = []

    trim_rec = annotations[TRIMMING_RECOMMENDATIONS]
    trim_rec[kind].extend(segments)


class TrimEdges(_BaseTrim):
    'It adds a trimming recommendation a fixed number of bases from the seqs.'
    def __init__(self, left=0, right=0):
        '''The initiator.

        left - number of bases to trim from the left side
        right - number of bases to trim from the right side
        mask - If True the edges will be masked instead of trimmed
        '''
        self.left = left
        self.right = right
        super(TrimEdges, self).__init__()

    def _do_trim(self, seq):
        'It trims the edges of the given seqs.'
        left = self.left
        right = self.right
        segments = [(0, left - 1)] if left else []
        if right:
            seq_len = get_length(seq)
            segments.append((seq_len - right, seq_len - 1))
        _add_trim_segments(segments, seq, kind=OTHER)
        return seq


def _mask_sequence(seq, segments):
    'It masks the given segments of the sequence'

    if not segments:
        return seq
    segments = merge_overlaping_segments(segments)

    segments = get_all_segments(segments, get_length(seq))