SEQS_PASSED = 'seqs_passed' SEQS_FILTERED_OUT = 'seqs_filtered_out' SEQITEM = 'seqitem' SEQRECORD = 'seqrecord' def uppercase_length(string):
    'It returns the number of uppercase characters found in the string'
    return len(re.findall("[A-Z]", string)) 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 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_file_format(seq):
    seq_class = seq.kind

    if seq_class == SEQITEM:
        fmt = seq.file_format
    elif seq_class == SEQRECORD:
        fmt = None
    return fmt 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_length(seq):
    return len(get_str_seq(seq)) class WrongFormatError(Exception):
    'Raised when the file have an unexpected format'
    pass 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 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 calculate_dust_score(seq):
    '''It returns the dust score.

    From: "A Fast and Symmetric DUST Implementation to Mask Low-Complexity DNA
    Sequences"
    doi:10.1089/cmb.2006.13.1028

    and re-implemented from PRINSEQ
    '''
    seq = get_str_seq(seq)
    length = len(seq)
    if length == 3:
        return 0
    if length <= 5:
        return None

    windowsize = get_setting('DUST_WINDOWSIZE')
    windowstep = get_setting('DUST_WINDOWSTEP')

    dustscores = []
    if length > windowsize:
        windows = 0
        for seq_in_win in rolling_window(seq, windowsize, windowstep):
            score = _calculate_rawscore(seq_in_win)
            dustscores.append(score / (windowsize - 2))
            windows += 1
        remaining_seq = seq[windows * windowstep:]
    else:
        remaining_seq = seq

    if remaining_seq > 5:
        length = len(remaining_seq)
        score = _calculate_rawscore(remaining_seq)
        dustscore = score / (length - 3) * (windowsize - 2) / (length - 2)
        dustscores.append(dustscore)

    # max score should be 100 not 31
    dustscore = sum(dustscores) / len(dustscores) * 100 / 31
    return dustscore def get_setting(key):
    'It returns the value for one setting'
    return _settings[key] def map_with_bowtie2(index_fpath, paired_fpaths=None,
                     unpaired_fpath=None, readgroup=None, threads=None,
                     log_fpath=None, preset='very-sensitive-local',
                     extra_params=None):
    '''It maps with bowtie2.

    paired_seqs is a list of tuples, in which each tuple are paired seqs
    unpaired_seqs is a list of files
    '''
    if readgroup is None:
        readgroup = {}

    if extra_params is None:
        extra_params = []

    if paired_fpaths is None and unpaired_fpath is None:
        raise RuntimeError('At least one file to map is required')

    binary = get_binary_path('bowtie2')
    cmd = [binary, '-x', index_fpath, '--{0}'.format(preset),
           '-p', str(get_num_threads(threads))]

    cmd.extend(extra_params)
    if unpaired_fpath:
        cmd.extend(['-U', unpaired_fpath])
    if paired_fpaths:
        cmd.extend(['-1', paired_fpaths[0], '-2', paired_fpaths[1]])

    if 'ID' in readgroup.keys():
        for key, value in readgroup.items():
            if key not in ('ID', 'LB', 'SM', 'PL'):
                msg = 'The readgroup header tag is not valid: {}'.format(key)
                raise RuntimeError(msg)
            if key == 'ID':
                cmd.extend(['--rg-id', value])
            else:
                cmd.extend(['--rg', '{0}:{1}'.format(key, value)])

    if log_fpath is None:
        stderr = NamedTemporaryFile(suffix='.stderr')
    else:
        stderr = open(log_fpath, 'w')

    bowtie2 = popen(cmd, stderr=stderr, stdout=PIPE)
    # print bowtie2.stdout.read()
    return bowtie2 def map_process_to_bam(map_process, bam_fpath, log_fpath=None,
                       tempdir=None):
    ''' It receives a mapping process that has a sam file in stdout and
    calling another external process convert the sam file into a bam file.
    Optionally you can fill the readgroup field
    '''
    if log_fpath is None:
        stderr = NamedTemporaryFile(suffix='.stderr')
    else:
        stderr = open(log_fpath, 'w')
    cmd = [get_binary_path('samtools'), 'view', '-h', '-b', '-S', '-',
           '-o', bam_fpath]

    samtools = popen(cmd, stdin=map_process.stdout, stderr=stderr)
    map_process.stdout.close()  # Allow p1 to receive a SIGPIPE if samtools exits.
    samtools.communicate() 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(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 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
from tempfile import NamedTemporaryFile
    from pysam import Samfile
from crumbs.utils.tags import (SEQS_PASSED, SEQS_FILTERED_OUT, SEQITEM,
                               SEQRECORD)
from crumbs.seq.utils.seq_utils import uppercase_length, get_uppercase_segments
from crumbs.seq.seq import get_name, get_file_format, get_str_seq, get_length
from crumbs.exceptions import WrongFormatError
from crumbs.blast import Blaster, BlasterForFewSubjects
from crumbs.statistics import calculate_dust_score
from crumbs.settings import get_setting
from crumbs.mapping import map_with_bowtie2, map_process_to_bam
from crumbs.seq.seqio import write_seqs
from crumbs.seq.pairs import group_pairs, group_pairs_by_name
# Copyright 2012 Jose Blanca, Peio Ziarsolo, COMAV-Univ. Politecnica Valencia
# This file is part of seq_crumbs.
# seq_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.

# seq_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 seq_crumbs. If not, see <http://www.gnu.org/licenses/>.

# pylint: disable=C0111

from __future__ import division

try:
except ImportError:
    # This is an optional requirement
    pass



def seq_to_filterpackets(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))