Return the percent identity of a read.

    based on the NM tag if present,
    if not calculate from MD tag and CIGAR string

    read.query_alignment_length can be zero in the case of ultra long reads aligned with minimap2 -L
326462
def get_pID(read):
    
    try:
        return 100 * (1 - read.get_tag("NM") / read.query_alignment_length)
    except KeyError:
        try:
            return 100 * (1 - (parse_MD(read.get_tag("MD")) + parse_CIGAR(read.cigartuples))
                          / read.query_alignment_length)
        except KeyError:
            return None
    except ZeroDivisionError:
        return None
333306
function getTagLength(tag) {
  var lenTag = 4;
  
  while(lenTag > 1) {
    var tmpTag = tag >>> ((lenTag - 1) * 8);
    
    if ((tmpTag & 0xFF) != 0x00) {
      break;
    } 
    
    lenTag--;
  }
  
  return lenTag;
}
1215464
def actual_query_range(self):
    
    a = self.alignment_ranges
    #return GenomicRange(a[0][1].chr,a[0][1].start,a[-1][1].end,self.get_strand())
    if self.get_strand() == '+':
      return GenomicRange(a[0][1].chr,a[0][1].start,a[-1][1].end,self.get_strand())
    #must be - strand
    return GenomicRange(a[0][1].chr,self.query_sequence_length-a[-1][1].end+1,self.query_sequence_length-a[0][1].start+1,dir=self.strand)
426060
def qn(tag):
    
    prefix, tagroot = tag.split(':')
    uri = nsmap[prefix]
    return '{%s}%s' % (uri, tagroot)
966355
def _compute_signature(parameters, access_key_secret):
    '''
    
    '''

    def percent_encode(line):
        if not isinstance(line, six.string_types):
            return line

        s = line
        if sys.stdin.encoding is None:
            s = line.decode().encode('utf8')
        else:
            s = line.decode(sys.stdin.encoding).encode('utf8')
        res = _quote(s, '')
        res = res.replace('+', '%20')
        res = res.replace('*', '%2A')
        res = res.replace('%7E', '~')
        return res

    sortedParameters = sorted(list(parameters.items()), key=lambda items: items[0])

    canonicalizedQueryString = ''
    for k, v in sortedParameters:
        canonicalizedQueryString += '&' + percent_encode(k) \
            + '=' + percent_encode(v)

    # All aliyun API only support GET method
    stringToSign = 'GET&%2F&' + percent_encode(canonicalizedQueryString[1:])

    h = hmac.new(to_bytes(access_key_secret + "&"), stringToSign, sha1)
    signature = base64.encodestring(h.digest()).strip()
    return signature
55851
def count_uniquelymapping_reads(self, file_name, paired_end):
        
        param = " -c -F256"
        if file_name.endswith("sam"):
            param += " -S"
        return self.samtools_view(file_name, param=param)
251264
def read(self, size=None):
        
        if self._is_mini:
            self.seek(self._position)
        else:
            self.source.seek(self._source_position)
        if not size or size < 0:
            size = self.size - self.tell()

        data = b('')
        while len(data) < size:
            if self.tell() > self.size:
                break
            if self._sector_number == ENDOFCHAIN:
                break

            to_read = size - len(data)
            to_end = self.sector_size - self._position_in_sector
            to_do = min(to_read, to_end)
            data += self.stream.read(to_do)

            self._position += to_do
            self._source_position = self.source.tell()

            if to_read >= to_end:
                self._position_in_sector = 0

                self._sector_number = self.next_sector(self._sector_number)
                position = (self._sector_number + int(not self._is_mini)) \
                    << self.sector_shift
                self.stream.seek(position)
            else:
                self._position_in_sector += to_do

        return data
121593
def estimate_fragment_size(bam_file, nreads=5000):
    
    with open_samfile(bam_file) as bam_handle:
        reads = tz.itertoolz.take(nreads, bam_handle)
        # it would be good to skip spliced paired reads.
        lengths = [x.template_length for x in reads if x.template_length > 0]
    if not lengths:
        return 0
    return int(numpy.median(lengths))
1208082
def convertReadAlignment(self, read, readGroupSet, readGroupId):
        
        samFile = self.getFileHandle(self._dataUrl)
        # TODO fill out remaining fields
        # TODO refine in tandem with code in converters module
        ret = protocol.ReadAlignment()
        # ret.fragmentId = 'TODO'
        ret.aligned_quality.extend(read.query_qualities)
        ret.aligned_sequence = read.query_sequence
        if SamFlags.isFlagSet(read.flag, SamFlags.READ_UNMAPPED):
            ret.ClearField("alignment")
        else:
            ret.alignment.CopyFrom(protocol.LinearAlignment())
            ret.alignment.mapping_quality = read.mapping_quality
            ret.alignment.position.CopyFrom(protocol.Position())
            ret.alignment.position.reference_name = samFile.getrname(
                read.reference_id)
            ret.alignment.position.position = read.reference_start
            ret.alignment.position.strand = protocol.POS_STRAND
            if SamFlags.isFlagSet(read.flag, SamFlags.READ_REVERSE_STRAND):
                ret.alignment.position.strand = protocol.NEG_STRAND
            for operation, length in read.cigar:
                gaCigarUnit = ret.alignment.cigar.add()
                gaCigarUnit.operation = SamCigar.int2ga(operation)
                gaCigarUnit.operation_length = length
                gaCigarUnit.reference_sequence = ""  # TODO fix this!
        ret.duplicate_fragment = SamFlags.isFlagSet(
            read.flag, SamFlags.DUPLICATE_READ)
        ret.failed_vendor_quality_checks = SamFlags.isFlagSet(
            read.flag, SamFlags.FAILED_QUALITY_CHECK)
        ret.fragment_length = read.template_length
        ret.fragment_name = read.query_name
        for key, value in read.tags:
            # Useful for inspecting the structure of read tags
            # print("{key} {ktype}: {value}, {vtype}".format(
            #     key=key, ktype=type(key), value=value, vtype=type(value)))
            protocol.setAttribute(ret.attributes.attr[key].values, value)

        if SamFlags.isFlagSet(read.flag, SamFlags.MATE_UNMAPPED):
            ret.next_mate_position.Clear()
        else:
            ret.next_mate_position.Clear()
            if read.next_reference_id != -1:
                ret.next_mate_position.reference_name = samFile.getrname(
                    read.next_reference_id)
            else:
                ret.next_mate_position.reference_name = ""
            ret.next_mate_position.position = read.next_reference_start
            ret.next_mate_position.strand = protocol.POS_STRAND
            if SamFlags.isFlagSet(read.flag, SamFlags.MATE_REVERSE_STRAND):
                ret.next_mate_position.strand = protocol.NEG_STRAND
        if SamFlags.isFlagSet(read.flag, SamFlags.READ_PAIRED):
            ret.number_reads = 2
        else:
            ret.number_reads = 1
        ret.read_number = -1
        if SamFlags.isFlagSet(read.flag, SamFlags.FIRST_IN_PAIR):
            if SamFlags.isFlagSet(read.flag, SamFlags.SECOND_IN_PAIR):
                ret.read_number = 2
            else:
                ret.read_number = 0
        elif SamFlags.isFlagSet(read.flag, SamFlags.SECOND_IN_PAIR):
            ret.read_number = 1
        ret.improper_placement = not SamFlags.isFlagSet(
            read.flag, SamFlags.READ_PROPER_PAIR)
        ret.read_group_id = readGroupId
        ret.secondary_alignment = SamFlags.isFlagSet(
            read.flag, SamFlags.SECONDARY_ALIGNMENT)
        ret.supplementary_alignment = SamFlags.isFlagSet(
            read.flag, SamFlags.SUPPLEMENTARY_ALIGNMENT)
        ret.id = readGroupSet.getReadAlignmentId(ret)
        return ret
947891
def read_alignment(out_sam, loci, seqs, out_file):
    """"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():
            nm = hit[0][0]
            for l in hit:
                if nm == l[0]:
                    print(l[1], file=out_handle)
    return out_file
519593
def _contains(self, hit):
        ''''''
        return (
            hit.qry_name not in self.contigs_to_keep
            and hit.qry_name != hit.ref_name
            and (100 * hit.hit_length_qry / hit.qry_length >= self.min_contig_percent_match)
            and hit.percent_identity >= self.nucmer_min_id
        )
1301697
def read_md5(self, hex=False):
        
        f = self.open('rb')
        try:
            m = hashlib.md5()
            while True:
                d = f.read(8192)
                if not d:
                    break
                m.update(d)
        finally:
            f.close()
        if hex:
            return m.hexdigest()
        else:
            return m.digest()
1030022
def run_tag_from_session_and_metric(session_name, metric_name):
  
  assert isinstance(session_name, six.string_types)
  assert isinstance(metric_name, api_pb2.MetricName)
  # os.path.join() will append a final slash if the group is empty; it seems
  # like multiplexer.Tensors won't recognize paths that end with a '/' so
  # we normalize the result of os.path.join() to remove the final '/' in that
  # case.
  run = os.path.normpath(os.path.join(session_name, metric_name.group))
  tag = metric_name.tag
  return run, tag
710298
def get_spectra_id(self, fn_id, retention_time=None, scan_nr=None):
        """"""
        cursor = self.get_cursor()
        sql = 'SELECT spectra_id FROM mzml WHERE mzmlfile_id=? '
        values = [fn_id]
        if retention_time is not None:
            sql = '{0} AND retention_time=?'.format(sql)
            values.append(retention_time)
        if scan_nr is not None:
            sql = '{0} AND scan_nr=?'.format(sql)
            values.append(scan_nr)
        cursor.execute(sql, tuple(values))
        return cursor.fetchone()[0]
686150
def check_bam(bam, o):
    
    try:
        p = sp.Popen(['samtools', 'view', bam], stdout=sp.PIPE)
        # Count paired alignments
        paired = 0
        read_lengths = defaultdict(int)
        while o > 0:  # Count down number of lines
            line = p.stdout.readline().decode().split("\t")
            flag = int(line[1])
            read_lengths[len(line[9])] += 1
            if 1 & flag:  # check decimal flag contains 1 (paired)
                paired += 1
            o -= 1
        p.kill()
    except OSError:
        reason = "Note (samtools not in path): For NGS inputs, " \
                 "pep needs samtools to auto-populate " \
                 "'read_length' and 'read_type' attributes; " \
                 "these attributes were not populated."
        raise OSError(reason)

    _LOGGER.debug("Read lengths: {}".format(read_lengths))
    _LOGGER.debug("paired: {}".format(paired))
    return read_lengths, paired
637636
def _check_reads_hit(self, alignment_io, min_aligned_fraction):
        ''''''
        to_return = []
        alignment_length = None
        for s in SeqIO.parse(alignment_io, "fasta"):
            if not alignment_length:
                alignment_length = len(s.seq)
                min_length = int(min_aligned_fraction * alignment_length)
                logging.debug("Determined min number of aligned bases to be %s" % min_length)
            elif len(s.seq) != alignment_length:
                raise Exception("Alignment file appears to not be of uniform length")

            num_unaligned = s.seq.count('-')
            num_aligned = alignment_length-num_unaligned
            logging.debug("Sequence %s has %d aligned positions" % (s.name, alignment_length-num_unaligned))
            if num_aligned <= min_length:
                to_return.append(s.name)
        return to_return
115056
def check_tag_data(data):
    ""
    if len(data) < 10:
        raise ValueError("Tag too short")
    if data[0:3] != b"ID3":
        raise ValueError("Missing ID3 identifier")
    if data[3] >= 5 or data[4] != 0:
        raise ValueError("Unknown ID3 version")
    length = stagger.conversion.Syncsafe.decode(data[6:10]) + 10
    if len(data) != length:
        raise ValueError("Tag size mismatch")