# ---------------------------------------------------------------------------- # Copyright (c) 2013--, scikit-bio development team. # # Distributed under the terms of the Modified BSD License. # # The full license is in the file COPYING.txt, distributed with this software. # ---------------------------------------------------------------------------- import re import warnings import numpy as np from skbio.util import cardinal_to_ordinal _whitespace_regex = re.compile(r'\s') _newline_regex = re.compile(r'\n') def _decode_qual_to_phred(qual_str, variant=None, phred_offset=None): phred_offset, phred_range = _get_phred_offset_and_range( variant, phred_offset, ["Must provide either `variant` or `phred_offset` in order to decode " "quality scores.", "Decoding Solexa quality scores is not currently supported, " "as quality scores are always stored as Phred scores in " "scikit-bio. Please see the following scikit-bio issue to " "track progress on this:\n\t" "https://github.com/biocore/scikit-bio/issues/719"]) qual = np.fromstring(qual_str, dtype=np.uint8) - phred_offset if np.any((qual > phred_range[1]) | (qual < phred_range[0])): raise ValueError("Decoded Phred score is out of range [%d, %d]." % (phred_range[0], phred_range[1])) return qual def _encode_phred_to_qual(phred, variant=None, phred_offset=None): phred_offset, phred_range = _get_phred_offset_and_range( variant, phred_offset, ["Must provide either `variant` or `phred_offset` in order to encode " "Phred scores.", "Encoding Solexa quality scores is not currently supported. " "Please see the following scikit-bio issue to track progress " "on this:\n\t" "https://github.com/biocore/scikit-bio/issues/719"]) qual_chars = [] for score in phred: if score < phred_range[0]: raise ValueError("Phred score %d is out of range [%d, %d]." % (score, phred_range[0], phred_range[1])) if score > phred_range[1]: warnings.warn( "Phred score %d is out of targeted range [%d, %d]. Converting " "to %d." % (score, phred_range[0], phred_range[1], phred_range[1]), UserWarning) score = phred_range[1] qual_chars.append(chr(score + phred_offset)) return ''.join(qual_chars) def _get_phred_offset_and_range(variant, phred_offset, errors): if variant is None and phred_offset is None: raise ValueError(errors[0]) if variant is not None and phred_offset is not None: raise ValueError( "Cannot provide both `variant` and `phred_offset`.") if variant is not None: if variant == 'sanger': phred_offset = 33 phred_range = (0, 93) elif variant == 'illumina1.3': phred_offset = 64 phred_range = (0, 62) elif variant == 'illumina1.8': phred_offset = 33 phred_range = (0, 62) elif variant == 'solexa': phred_offset = 64 phred_range = (-5, 62) raise ValueError(errors[1]) else: raise ValueError("Unrecognized variant %r." % variant) else: if not (33 <= phred_offset <= 126): raise ValueError( "`phred_offset` %d is out of printable ASCII character range." % phred_offset) phred_range = (0, 126 - phred_offset) return phred_offset, phred_range def _get_nth_sequence(generator, seq_num): # i is set to None so that an empty generator will not result in an # undefined variable when compared to seq_num. i = None if seq_num is None or seq_num < 1: raise ValueError('Invalid sequence number (`seq_num`=%s). `seq_num`' ' must be between 1 and the number of sequences in' ' the file.' % str(seq_num)) try: for i, seq in zip(range(1, seq_num + 1), generator): pass finally: generator.close() if i == seq_num: return seq raise ValueError('Reached end of file before finding the %s sequence.' % cardinal_to_ordinal(seq_num)) def _parse_fasta_like_header(line): id_ = '' desc = '' header = line[1:].rstrip() if header: if header[0].isspace(): # no id desc = header.lstrip() else: header_tokens = header.split(None, 1) if len(header_tokens) == 1: # no description id_ = header_tokens[0] else: id_, desc = header_tokens return id_, desc def _format_fasta_like_records(generator, id_whitespace_replacement, description_newline_replacement, require_qual, lowercase=None): if ((id_whitespace_replacement is not None and '\n' in id_whitespace_replacement) or (description_newline_replacement is not None and '\n' in description_newline_replacement)): raise ValueError( "Newline character (\\n) cannot be used to replace whitespace in " "sequence IDs, nor to replace newlines in sequence descriptions.") for idx, seq in enumerate(generator): if len(seq) < 1: raise ValueError( "%s sequence does not contain any characters (i.e., it is an " "empty/blank sequence). Writing empty sequences is not " "supported." % cardinal_to_ordinal(idx + 1)) if 'id' in seq.metadata: id_ = '%s' % seq.metadata['id'] else: id_ = '' if id_whitespace_replacement is not None: id_ = _whitespace_regex.sub(id_whitespace_replacement, id_) if 'description' in seq.metadata: desc = '%s' % seq.metadata['description'] else: desc = '' if description_newline_replacement is not None: desc = _newline_regex.sub(description_newline_replacement, desc) if desc: header = '%s %s' % (id_, desc) else: header = id_ if require_qual and 'quality' not in seq.positional_metadata: raise ValueError( "Cannot write %s sequence because it does not have quality " "scores associated with it." % cardinal_to_ordinal(idx + 1)) qual = None if 'quality' in seq.positional_metadata: qual = seq.positional_metadata['quality'].values if lowercase is not None: seq_str = seq.lowercase(lowercase) else: seq_str = str(seq) yield header, "%s" % seq_str, qual def _line_generator(fh, skip_blanks=False, strip=True): for line in fh: if strip: line = line.strip() skip = False if skip_blanks: skip = line.isspace() or not line if not skip: yield line def _too_many_blanks(fh, max_blanks): count = 0 too_many = False for line in _line_generator(fh, skip_blanks=False): if line: break else: count += 1 if count > max_blanks: too_many = True break fh.seek(0) return too_many