Extracting metrics from unaligned bam format
    Extracting lengths
326810
def process_ubam(bam, **kwargs):
    
    logging.info("Nanoget: Starting to collect statistics from ubam file {}.".format(bam))
    samfile = pysam.AlignmentFile(bam, "rb", check_sq=False)
    if not samfile.has_index():
        pysam.index(bam)
        # Need to reload the samfile after creating index
        samfile = pysam.AlignmentFile(bam, "rb")
        logging.info("Nanoget: No index for bam file could be found, created index.")
    datadf = pd.DataFrame(
        data=[(read.query_name, nanomath.ave_qual(read.query_qualities), read.query_length)
              for read in samfile.fetch(until_eof=True)],
        columns=["readIDs", "quals", "lengths"]) \
        .dropna(axis='columns', how='all') \
        .dropna(axis='index', how='any')
    logging.info("Nanoget: ubam {} contains {} reads.".format(
        bam, datadf["lengths"].size))
    return ut.reduce_memory_usage(datadf)
1367911
def _table_info(self, buf=None):
    
    metrics = [self.count().name('nrows')]
    for col in self.columns:
        metrics.append(self[col].count().name(col))

    metrics = self.aggregate(metrics).execute().loc[0]

    names = ['Column', '------'] + self.columns
    types = ['Type', '----'] + [repr(x) for x in self.schema().types]
    counts = ['Non-null #', '----------'] + [str(x) for x in metrics[1:]]
    col_metrics = util.adjoin(2, names, types, counts)
    result = 'Table rows: {}\n\n{}'.format(metrics[0], col_metrics)

    print(result, file=buf)
299119
def detect_bam_features(bamfile, n_entries=1000):
    '''  '''

    inbam = pysam.Samfile(bamfile)
    inbam = inbam.fetch(until_eof=True)

    tags = ["NH", "X0", "XT"]
    available_tags = {x: 1 for x in tags}

    for n, read in enumerate(inbam):
        if n > n_entries:
            break

        if read.is_unmapped:
            continue

        else:
            for tag in tags:
                if not read.has_tag(tag):
                    available_tags[tag] = 0

    return available_tags
991145
def bamtag(sam):
    ''' 
    '''
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    if is_python3():
        queryalignment = next(track)
    else:
        queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            aln.tags += [('XC', match.group('CB'))]
        if "molecular" in annotations:
            aln.tags += [('RX', match.group('MB'))]
        if "sample" in annotations:
            aln.tags += [('XS', match.group('SB'))]

        out_file.write(aln)

    total_time = time.time() - start_time
    logger.info('BAM tag conversion done - {:.3}s, {:,} alns/min'.format(total_time, int(60. * count / total_time)))
    logger.info("Processed %d alignments." % count)
1368480
def format(file_metrics, build_metrics):
    """"""
    def indent(elem, level=0):
        i = "\n" + level*"  "
        if len(elem):
            if not elem.text or not elem.text.strip():
                elem.text = i + "  "
            if not elem.tail or not elem.tail.strip():
                elem.tail = i
            for elem in elem:
                indent(elem, level+1)
            if not elem.tail or not elem.tail.strip():
                elem.tail = i
        else:
            if level and (not elem.tail or not elem.tail.strip()):
                elem.tail = i

    root = ET.Element('metrics')

    # file_metrics
    files = ET.Element('files')
    root.append(files)

    for key in file_metrics.keys():
        tmp_file = ET.SubElement(files, "file",
                                 {'name': key, 'language': file_metrics[key]['language']})
        for name in file_metrics[key].keys():
            if name == 'language':
                continue
            tmp_metric = ET.SubElement(tmp_file, "metric",
                                       {'name': name, 'value': str(file_metrics[key][name])})

    # build_metrics
    if build_metrics:
        build = ET.Element('build')
        root.append(build)
        # TODO

    indent(root)
    if PY3:
        body = ET.tostring(root, encoding='unicode')
    else:
        body = ET.tostring(root)
    return body
136221
function metrics( clbk ) {
	var m = {},
		mem,
		used,
		total;

	// Get the process id:
	m.pid = process.pid;

	// Get the number of milliseconds the process has been running:
	m.uptime = process.uptime() * 1000; // ms

	// Extract process memory usage:
	mem = process.memoryUsage();
	used = mem.heapUsed;
	total = mem.heapTotal;

	m.mem = {
		'rss': mem.rss / 1000, // KB
		'heapFree': (total-used) / 1000, // KB
		'heapTotal': total / 1000, // KB
		'heapUtilization': used / total
	};

	// Get the average amount of time requests have to wait in Node's event queue before being processed:
	m.lag = toobusy.lag(); // ms

	// Get the CPU and memory consumption:
	pidusage.stat( m.pid, onStats );

	/**
	* FUNCTION: onStats( error, stat )
	*	Callback invoked upon collecting CPU and memory data for the current process.
	*
	* @private
	* @param {Object} error
	* @param {Object} stat - process data
	*/
	function onStats( error, stat ) {
		if ( error ) {
			return clbk( error );
		}
		m.mem.memUsed = stat.memory / 1000; // KB
		m.mem.utilization = stat.memory / os.totalmem();
		m.cpu = {
			'utilization': stat.cpu / 100
		};
		clbk( null, m );
	}
}
1005952
def pretty_memory_info():
    '''
    
    '''
    process = psutil.Process(os.getpid())
    return '{}MB memory usage'.format(int(process.memory_info().rss / 2**20))
1158922
def has_aligned_reads(align_bam, region=None):
    
    import pybedtools
    if region is not None:
        if isinstance(region, six.string_types) and os.path.isfile(region):
            regions = [tuple(r) for r in pybedtools.BedTool(region)]
        else:
            regions = [region]
    with pysam.Samfile(align_bam, "rb") as cur_bam:
        if region is not None:
            for region in regions:
                if isinstance(region, six.string_types):
                    for item in cur_bam.fetch(str(region)):
                        return True
                else:
                    for item in cur_bam.fetch(str(region[0]), int(region[1]), int(region[2])):
                        return True
        else:
            for item in cur_bam:
                if not item.is_unmapped:
                    return True
    return False
790320
def subset_bamfile(sam, barcodes):
    
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)
    barcodes = set(barcode.strip() for barcode in barcodes)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            cb = match.group('CB')
            if cb in barcodes:
                out_file.write(aln)
1146566
def extractTwos(m):
    """"""
    # the problem can be break down to count how many '0's are there in
    # the end of bin(m). This can be done this way: m & a stretch of '1's
    # which can be represent as (2 ** n) - 1.
    assert m >= 0
    i = 0
    while m & (2 ** i) == 0:
        i += 1
    return i, m >> i
684882
def bytes_from_readable_size(C, size, suffix='B'):
        """"""
        s = re.split("^([0-9\.]+)\s*([%s]?)%s?" % (''.join(C.SIZE_UNITS), suffix), size, flags=re.I)
        bytes, unit = round(float(s[1])), s[2].upper()
        while unit in C.SIZE_UNITS and C.SIZE_UNITS.index(unit) > 0:
            bytes *= 1024
            unit = C.SIZE_UNITS[C.SIZE_UNITS.index(unit) - 1]
        return bytes
604298
def _readable_part_size(num_bytes):
    ""
    B = num_bytes
    KB = float(1024)
    MB = float(KB * 1024)
    GB = float(MB * 1024)
    TB = float(GB * 1024)

    if B < KB:
        return '{0} {1}'.format(B, 'bytes' if B != 1 else 'byte')
    elif KB <= B < MB:
        return '{0:.2f} KiB'.format(B/KB)
    elif MB <= B < GB:
        return '{0:.2f} MiB'.format(B/MB)
    elif GB <= B < TB:
        return '{0:.2f} GiB'.format(B/GB)
    elif TB <= B:
        return '{0:.2f} TiB'.format(B/TB)
1184277
def get_seqstarts(bamfile, N):
    
    import pysam
    bamfile = pysam.AlignmentFile(bamfile, "rb")
    seqsize = {}
    for kv in bamfile.header["SQ"]:
        if kv["LN"] < 10 * N:
            continue
        seqsize[kv["SN"]] = kv["LN"] / N + 1

    allseqs = natsorted(seqsize.keys())
    allseqsizes = np.array([seqsize[x] for x in allseqs])
    seqstarts = np.cumsum(allseqsizes)
    seqstarts = np.roll(seqstarts, 1)
    total_bins = seqstarts[0]
    seqstarts[0] = 0
    seqstarts = dict(zip(allseqs, seqstarts))

    return seqstarts, seqsize, total_bins
158073
def process(self, metric):
        
        entry = {'timestamp': metric.timestamp, 'value': metric.value,
                 "tags": {}}
        entry["tags"]["hostname"] = metric.host

        if self.cleanMetrics:
            metric = MetricWrapper(metric, self.log)
            if self.skipAggregates and metric.isAggregate():
                return
            for tagKey in metric.getTags():
                entry["tags"][tagKey] = metric.getTags()[tagKey]

        entry['metric'] = (self.prefix + metric.getCollectorPath() +
                           '.' + metric.getMetricPath())

        for [key, value] in self.tags:
            entry["tags"][key] = value

        self.entrys.append(entry)

        # send data if list is long enough
        if (len(self.entrys) >= self.batch):
            # Compress data
            if self.compression >= 1:
                data = StringIO.StringIO()
                with contextlib.closing(gzip.GzipFile(fileobj=data,
                                        compresslevel=self.compression,
                                        mode="w")) as f:
                    f.write(json.dumps(self.entrys))
                self._send(data.getvalue())
            else:
                # no compression
                data = json.dumps(self.entrys)
                self._send(data)
874207
def _extract_metric_from_map(self, _m, messages, obj_map, obj_help):
        
        _obj = metrics_pb2.MetricFamily()
        _obj.name = _m
        _obj.type = self.METRIC_TYPES.index(obj_map[_m])
        if _m in obj_help:
            _obj.help = obj_help[_m]
        # trick for histograms
        _newlbl = _m
        if obj_map[_m] == 'histogram':
            _newlbl = '{}_bucket'.format(_m)
        # Loop through the array of metrics ({labels, value}) built earlier
        for _metric in messages[_newlbl]:
            # in the case of quantiles and buckets, they need to be grouped by labels
            if obj_map[_m] in ['summary', 'histogram'] and len(_obj.metric) > 0:
                _label_exists = False
                _metric_minus = {k: v for k, v in list(iteritems(_metric['labels'])) if k not in ['quantile', 'le']}
                _metric_idx = 0
                for mls in _obj.metric:
                    _tmp_lbl = {idx.name: idx.value for idx in mls.label}
                    if _metric_minus == _tmp_lbl:
                        _label_exists = True
                        break
                    _metric_idx = _metric_idx + 1
                if _label_exists:
                    _g = _obj.metric[_metric_idx]
                else:
                    _g = _obj.metric.add()
            else:
                _g = _obj.metric.add()
            if obj_map[_m] == 'counter':
                _g.counter.value = float(_metric['value'])
            elif obj_map[_m] == 'gauge':
                _g.gauge.value = float(_metric['value'])
            elif obj_map[_m] == 'summary':
                if '{}_count'.format(_m) in messages:
                    _g.summary.sample_count = long(self.get_metric_value_by_labels(messages, _metric, _m, 'count'))
                if '{}_sum'.format(_m) in messages:
                    _g.summary.sample_sum = self.get_metric_value_by_labels(messages, _metric, _m, 'sum')
            # TODO: see what can be done with the untyped metrics
            elif obj_map[_m] == 'histogram':
                if '{}_count'.format(_m) in messages:
                    _g.histogram.sample_count = long(self.get_metric_value_by_labels(messages, _metric, _m, 'count'))
                if '{}_sum'.format(_m) in messages:
                    _g.histogram.sample_sum = self.get_metric_value_by_labels(messages, _metric, _m, 'sum')
            # last_metric = len(_obj.metric) - 1
            # if last_metric >= 0:
            for lbl in _metric['labels']:
                # In the string format, the quantiles are in the labels
                if lbl == 'quantile':
                    # _q = _obj.metric[last_metric].summary.quantile.add()
                    _q = _g.summary.quantile.add()
                    _q.quantile = float(_metric['labels'][lbl])
                    _q.value = float(_metric['value'])
                # The upper_bounds are stored as "le" labels on string format
                elif obj_map[_m] == 'histogram' and lbl == 'le':
                    # _q = _obj.metric[last_metric].histogram.bucket.add()
                    _q = _g.histogram.bucket.add()
                    _q.upper_bound = float(_metric['labels'][lbl])
                    _q.cumulative_count = long(float(_metric['value']))
                else:
                    # labels deduplication
                    is_in_labels = False
                    for _existing_lbl in _g.label:
                        if lbl == _existing_lbl.name:
                            is_in_labels = True
                    if not is_in_labels:
                        _l = _g.label.add()
                        _l.name = lbl
                        _l.value = _metric['labels'][lbl]
        return _obj
462205
def _unpack_from(buf, format_s, count=None, little_endian=True):
    """"""

    if count is not None:
        assert count > 0
        format_s = '%i%s' % (count, format_s)

    if little_endian is True:
        format_s = '<' + format_s
    else:
        format_s = '>' + format_s

    size = struct.calcsize(format_s)
    res = struct.unpack_from(format_s, buf.read(size))

    if count is not None:
        return res
    else:
        return res[0]
693360
byte[] pickleMetrics(List<MetricTuple> metrics) throws IOException {
        // Extremely rough estimate of 75 bytes per message
        ByteArrayOutputStream out = new ByteArrayOutputStream(metrics.size() * 75);
        Writer pickled = new OutputStreamWriter(out, charset);

        pickled.append(MARK);
        pickled.append(LIST);

        for (MetricTuple tuple : metrics) {
            // start the outer tuple
            pickled.append(MARK);

            // the metric name is a string.
            pickled.append(STRING);
            // the single quotes are to match python's repr("abcd")
            pickled.append(QUOTE);
            pickled.append(tuple.name);
            pickled.append(QUOTE);
            pickled.append(LF);

            // start the inner tuple
            pickled.append(MARK);

            // timestamp is a long
            pickled.append(LONG);
            pickled.append(Long.toString(tuple.timestamp));
            // the trailing L is to match python's repr(long(1234))
            pickled.append(LONG);
            pickled.append(LF);

            // and the value is a string.
            pickled.append(STRING);
            pickled.append(QUOTE);
            pickled.append(tuple.value);
            pickled.append(QUOTE);
            pickled.append(LF);

            pickled.append(TUPLE); // inner close
            pickled.append(TUPLE); // outer close

            pickled.append(APPEND);
        }

        // every pickle ends with STOP
        pickled.append(STOP);

        pickled.flush();

        return out.toByteArray();
    }