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
    # 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))
        else:
            packet = list(group_pairs(packet, n_seqs_in_pair=1))
        yield {SEQS_PASSED: packet, SEQS_FILTERED_OUT: []}


def _reverse_bools(bools):
    return (not bool_ for bool_ in bools)


class _BaseFilter(object):
    def __init__(self, failed_drags_pair=True, reverse=False):
        self.reverse = reverse
        self.failed_drags_pair = failed_drags_pair

    def _setup_checks(self, filterpacket):
        pass

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

    def __call__(self, filterpacket):
        self._setup_checks(filterpacket)
        reverse = self.reverse
        failed_drags_pair = self.failed_drags_pair
        seqs_passed = []
        filtered_out = filterpacket[SEQS_FILTERED_OUT][:]
        for paired_seqs in filterpacket[SEQS_PASSED]:
            checks = (self._do_check(seq) for seq in paired_seqs)
            if reverse:
                checks = _reverse_bools(checks)

            if failed_drags_pair:
                pair_passed = True if all(checks) else False
            else:
                pair_passed = True if any(checks) else False

            if pair_passed:
                seqs_passed.append(paired_seqs)
            else:
                filtered_out.append(paired_seqs)

        return {SEQS_PASSED: seqs_passed, SEQS_FILTERED_OUT: filtered_out}


class FilterByFeatureTypes(_BaseFilter):
    'It filters out sequences not annotated with the given feature types'
    def __init__(self, feature_types, reverse=False, failed_drags_pair=True):
        '''The initiator

        feat_types is a list of types of features to use to filter'''

        self._feat_types = feature_types
        super(FilterByFeatureTypes, self).__init__(reverse=reverse,
                                           failed_drags_pair=failed_drags_pair)

    def _do_check(self, seq):
        seq = seq.object
        f_in_seq = [f.type for f in seq.features if f.type in self._feat_types]
        return True if f_in_seq else False


class FilterByRpkm(_BaseFilter):
    def __init__(self, read_counts, min_rpkm, reverse=False,
                 failed_drags_pair=True):
        '''The init

        read_counts its a dict:
            - keys are the sequence names
            - values should be dicts with length, mapped_reads and
            unmmapped_reads, counts
        '''
        self._read_counts = read_counts
        self._min_rpkm = min_rpkm
        total_reads = sum([v['mapped_reads'] + v['unmapped_reads'] for v in read_counts.values()])
        self._million_reads = total_reads / 1e6
        super(FilterByRpkm, self).__init__(reverse=reverse,
                                           failed_drags_pair=failed_drags_pair)

    def _do_check(self, seq):
        count = self._read_counts[get_name(seq)]

        kb_len = count['length'] / 1000
        rpk = count['mapped_reads'] / kb_len  # rpks
        rpkm = rpk / self._million_reads  # rpkms

        return True if rpkm >= self._min_rpkm else False


class FilterByLength(_BaseFilter):
    'It removes the sequences according to their length.'
    def __init__(self, minimum=None, maximum=None, ignore_masked=False,
                 failed_drags_pair=True):
        '''The initiator.

        threshold - minimum length to pass the filter (integer)
        reverse - if True keep the short sequences and discard the long ones
        ignore_masked - If True only uppercase letters will be counted.
        '''
        if min is None and max is None:
            raise ValueError('min or max threshold must be given')
        self.min = minimum
        self.max = maximum
        self.ignore_masked = ignore_masked
        super(FilterByLength, self).__init__(reverse=False,
                                           failed_drags_pair=failed_drags_pair)

    def _do_check(self, seq):
        min_ = self.min
        max_ = self.max

        length = uppercase_length(get_str_seq(seq)) if self.ignore_masked else get_length(seq)