parses and plots tag length distribution files
344115
def parse_tagLength_dist(self):
        """"""
        # Find and parse homer tag length distribution reports
        for f in self.find_log_files('homer/LengthDistribution', filehandles=True):
            s_name = os.path.basename(f['root'])
            s_name = self.clean_s_name(s_name, f['root'])
            parsed_data = self.parse_length_dist(f)
            if parsed_data is not None:
                if s_name in self.tagdir_data['length']:
                    log.debug("Duplicate Length Distribution sample log found! Overwriting: {}".format(s_name))
                self.add_data_source(f, s_name, section='length')
                self.tagdir_data['length'][s_name] = parsed_data

        self.tagdir_data['length'] = self.ignore_samples(self.tagdir_data['length'])
        if len(self.tagdir_data['length']) > 0:
            self.add_section (
                name = 'Tag Length Distribution',
                anchor = 'homer-tagLength',
                description = 'This plot shows the distribution of tag length.',
                helptext = 'This is a good quality control for tag length inputed into Homer.',
                plot = self.length_dist_chart()
            )
522731
def summary(args):
    
    from jcvi.formats.base import SetFile
    from jcvi.formats.bed import BedSummary
    from jcvi.utils.table import tabulate

    p = OptionParser(summary.__doc__)
    p.add_option("--isoform", default=False, action="store_true",
                 help="Find longest isoform of each id")
    p.add_option("--ids", help="Only include features from certain IDs")
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    gff_file, = args
    ids = opts.ids

    if ids:
        ids = SetFile(ids)
        logging.debug("Total ids loaded: {0}".format(len(ids)))

        if opts.isoform:
            pids = set()
            gff = Gff(gff_file)
            for g in gff:
                if g.type != "mRNA":
                    continue
                if g.parent not in ids:
                    continue
                if "longest" not in g.attributes:
                    pids = set(x + ".1" for x in ids)
                    break
                if g.attributes["longest"][0] == "0":
                    continue
                pids.add(g.id)
            ids = pids
            logging.debug("After checking longest: {0}".format(len(ids)))

        # Collects aliases
        gff = Gff(gff_file)
        for g in gff:
            if g.name in ids:
                ids.add(g.id)
        logging.debug("Total ids including aliases: {0}".format(len(ids)))

    gff = Gff(gff_file)
    beds = defaultdict(list)
    for g in gff:
        if ids and not (g.id in ids or g.name in ids or g.parent in ids):
            continue

        beds[g.type].append(g.bedline)

    table = {}
    for type, bb in sorted(beds.items()):
        bs = BedSummary(bb)
        table[(type, "Features")] = bs.nfeats
        table[(type, "Unique bases")] = bs.unique_bases
        table[(type, "Total bases")] = bs.total_bases

    print(tabulate(table), file=sys.stdout)
389100
def tagcount(sam, out, genemap, output_evidence_table, positional, minevidence,
             cb_histogram, cb_cutoff, no_scale_evidence, subsample, sparse,
             parse_tags, gene_tags):
    ''' 
    '''
    from pysam import AlignmentFile

    from io import StringIO
    import pandas as pd

    from utils import weigh_evidence

    logger.info('Reading optional files')

    gene_map = None
    if genemap:
        with open(genemap) as fh:
            try:
                gene_map = dict(p.strip().split() for p in fh)
            except ValueError:
                logger.error('Incorrectly formatted gene_map, need to be tsv.')
                sys.exit()

    if positional:
        tuple_template = '{0},{1},{2},{3}'
    else:
        tuple_template = '{0},{1},{3}'

    if not cb_cutoff:
        cb_cutoff = 0

    if cb_histogram and cb_cutoff == "auto":
        cb_cutoff = guess_depth_cutoff(cb_histogram)

    cb_cutoff = int(cb_cutoff)

    cb_hist = None
    filter_cb = False
    if cb_histogram:
        cb_hist = pd.read_csv(cb_histogram, index_col=0, header=-1, squeeze=True, sep="\t")
        total_num_cbs = cb_hist.shape[0]
        cb_hist = cb_hist[cb_hist > cb_cutoff]
        logger.info('Keeping {} out of {} cellular barcodes.'.format(cb_hist.shape[0], total_num_cbs))
        filter_cb = True

    parser_re = re.compile('.*:CELL_(?P<CB>.*):UMI_(?P<MB>.*)')

    if subsample:
        logger.info('Creating reservoir of subsampled reads ({} per cell)'.format(subsample))
        start_sampling  = time.time()

        reservoir = collections.defaultdict(list)
        cb_hist_sampled = 0 * cb_hist
        cb_obs = 0 * cb_hist

        track = stream_bamfile(sam)
        current_read = 'none_observed_yet'
        for i, aln in enumerate(track):
            if aln.qname == current_read:
                continue

            current_read = aln.qname

            if parse_tags:
                CB = aln.get_tag('CR')
            else:
                match = parser_re.match(aln.qname)
                CB = match.group('CB')

            if CB not in cb_hist.index:
                continue

            cb_obs[CB] += 1
            if len(reservoir[CB]) < subsample:
                reservoir[CB].append(i)
                cb_hist_sampled[CB] += 1
            else:
                s = pd.np.random.randint(0, cb_obs[CB])
                if s < subsample:
                    reservoir[CB][s] = i

        index_filter = set(itertools.chain.from_iterable(reservoir.values()))
        sam_file.close()
        sampling_time = time.time() - start_sampling
        logger.info('Sampling done - {:.3}s'.format(sampling_time))

    evidence = collections.defaultdict(int)

    logger.info('Tallying evidence')
    start_tally = time.time()

    sam_mode = 'r' if sam.endswith(".sam") else 'rb'
    sam_file = AlignmentFile(sam, mode=sam_mode)
    targets = [x["SN"] for x in sam_file.header["SQ"]]
    track = sam_file.fetch(until_eof=True)
    count = 0
    unmapped = 0
    kept = 0
    nomatchcb = 0
    current_read = 'none_observed_yet'
    count_this_read = True
    missing_transcripts = set()
    for i, aln in enumerate(track):
        if count and not count % 1000000:
            logger.info("Processed %d alignments, kept %d." % (count, kept))
            logger.info("%d were filtered for being unmapped." % unmapped)
            if filter_cb:
                logger.info("%d were filtered for not matching known barcodes."
                            % nomatchcb)
        count += 1

        if aln.is_unmapped:
            unmapped += 1
            continue

        if gene_tags and not aln.has_tag('GX'):
            unmapped += 1
            continue

        if aln.qname != current_read:
            current_read = aln.qname
            if subsample and i not in index_filter:
                count_this_read = False
                continue
            else:
                count_this_read = True
        else:
            if not count_this_read:
                continue

        if parse_tags:
            CB = aln.get_tag('CR')
        else:
            match = parser_re.match(aln.qname)
            CB = match.group('CB')

        if filter_cb:
            if CB not in cb_hist.index:
                nomatchcb += 1
                continue

        if parse_tags:
            MB = aln.get_tag('UM')
        else:
            MB = match.group('MB')

        if gene_tags:
            target_name = aln.get_tag('GX').split(',')[0]
        else:
            txid = sam_file.getrname(aln.reference_id)
            if gene_map:
                if txid in gene_map:
                    target_name = gene_map[txid]
                else:
                    missing_transcripts.add(txid)
                    target_name = txid
            else:
                target_name = txid

        e_tuple = tuple_template.format(CB, target_name, aln.pos, MB)

        # Scale evidence by number of hits
        if no_scale_evidence:
            evidence[e_tuple] += 1.0
        else:
            evidence[e_tuple] += weigh_evidence(aln.tags)

        kept += 1

    tally_time = time.time() - start_tally
    if missing_transcripts:
        logger.warn('The following transcripts were missing gene_ids, so we added them as the transcript ids: %s' % str(missing_transcripts))
    logger.info('Tally done - {:.3}s, {:,} alns/min'.format(tally_time, int(60. * count / tally_time)))
    logger.info('Collapsing evidence')

    logger.info('Writing evidence')
    with tempfile.NamedTemporaryFile('w+t') as out_handle:
        for key in evidence:
            line = '{},{}\n'.format(key, evidence[key])
            out_handle.write(line)

        out_handle.flush()
        out_handle.seek(0)

        evidence_table = pd.read_csv(out_handle, header=None)

    del evidence

    evidence_query = 'evidence >= %f' % minevidence
    if positional:
        evidence_table.columns=['cell', 'gene', 'umi', 'pos', 'evidence']
        collapsed = evidence_table.query(evidence_query).groupby(['cell', 'gene'])['umi', 'pos'].size()

    else:
        evidence_table.columns=['cell', 'gene', 'umi', 'evidence']
        collapsed = evidence_table.query(evidence_query).groupby(['cell', 'gene'])['umi'].size()

    expanded = collapsed.unstack().T

    if gene_map:
        # This Series is just for sorting the index
        genes = pd.Series(index=set(gene_map.values()))
        genes = genes.sort_index()
        # Now genes is assigned to a DataFrame
        genes = expanded.ix[genes.index]

    elif gene_tags:
        expanded.sort_index()
        genes = expanded

    else:
        # make data frame have a complete accounting of transcripts
        targets = pd.Series(index=set(targets))
        targets = targets.sort_index()
        expanded = expanded.reindex(targets.index.values, fill_value=0)
        genes = expanded

    genes.fillna(0, inplace=True)
    genes = genes.astype(int)
    genes.index.name = "gene"

    logger.info('Output results')

    if subsample:
        cb_hist_sampled.to_csv('ss_{}_'.format(subsample) + os.path.basename(cb_histogram), sep='\t')

    if output_evidence_table:
        import shutil
        buf.seek(0)
        with open(output_evidence_table, 'w') as etab_fh:
            shutil.copyfileobj(buf, etab_fh)

    if sparse:
        pd.Series(genes.index).to_csv(out + ".rownames", index=False, header=False)
        pd.Series(genes.columns.values).to_csv(out + ".colnames", index=False, header=False)
        with open(out, "w+b") as out_handle:
            scipy.io.mmwrite(out_handle, scipy.sparse.csr_matrix(genes))
    else:
        genes.to_csv(out)
328044
def distributions_impl(self, tag, run):
    """"""
    (histograms, mime_type) = self._histograms_plugin.histograms_impl(
        tag, run, downsample_to=self.SAMPLE_SIZE)
    return ([self._compress(histogram) for histogram in histograms],
            mime_type)
1139479
function fparse(fpath) {
  if (fs.existsSync(!fpath)) {
    console.error(config_file + ': No such file.');
    return [false, false];
  }

  var fd        = fs.openSync(fpath, 'r');
  var read      = '';
  var pos       = 0;
  var start     = 0;
  var Ns        = {};
  var currentNs = null;
  var currentN  = null;
  var remnant   = '';
  var result    = {};
  var buffsize  = 65535;
  var summary   = null;
  var length    = 0;
  var emptyline = 0;

  function setAtEnd(_summary, _length, _emptyline) {
    // result, and fpath is in outer scope.
    _summary.length = _length - _emptyline;
    result[_summary.id] = new FASTA(_summary, fpath);
  }


  do {
    read = fs.readSync(fd, buffsize, pos);

    var lines = (remnant + read[0]).split('\n');
    remnant = lines.pop();

    lines.forEach(function(line, i) {

      if (line.match('N')) {
        if (!currentN) {
          currentN = {
            line: i,
            start: start + length,
            end : start + length+ line.length
          };
        }
        else {
          currentN.line++;
          currentN.end += line.length;
        }
      }
      else if (currentN) {
        currentNs.push(currentN);
        currentN = null;
      }

      if (line == '') {
        emptyline++;
      }
      else if (line.charAt(0) == '>') {
        // register a previous summary
        if (summary) {
          setAtEnd(summary, length, emptyline);
        }
        start += length;
        emptyline = 0;
        length    = 0;
        // make a new summary
        var id_desc = line.slice(1).split(' ');
        var id = id_desc.shift();
        var desc = (id_desc.length) ? id_desc.join(' ') : null;
        summary = {id: id, desc: desc, start: start, linelen: 0};
        // make a new Ns
        Ns[summary.id] = [];
        currentNs = Ns[summary.id];
      }
      else {
        if (!summary) {
          console.error(fpath +' does not seem to be FASTA format.');
          process.exit();
        }
        if (!summary.linelen) {
          summary.linelen = line.length;
        }
      }
      length += line.length + 1;
    });
    pos = pos + read[1];
  }
  while (read[1] > 0);

  // end
  length += remnant.replace('\n', '').length;
  setAtEnd(summary, length, emptyline);

  fs.closeSync(fd);
  return [result, Ns];
}
446134
def add_readlen_dist_plot(self):
        
        pconfig = {
            'id': 'skewer_read_length_histogram',
            'title': 'Skewer: Read Length Distribution after trimming',
            'xDecimals': False,
            'ylab': '% of Reads',
            'xlab': 'Read Length',
            'xmin': 0,
            'ymin': 0,
            'tt_label': '<b>{point.x}</b>: {point.y:.1f}%',
        }
        self.add_section(
            plot = linegraph.plot(self.skewer_readlen_dist, pconfig)
        )
835954
function parseFile(fs, file, type) {
    var filepath = file.getPath();
    var parser = file.getParser();

    if (!parser) {
        return Promise.reject(
            error.FileNotParsableError({
                filename: filepath
            })
        );
    }

    return fs.readAsString(filepath)
    .then(function(content) {
        if (type === 'readme') {
            return parser.parseReadme(content);
        } else if (type === 'glossary') {
            return parser.parseGlossary(content);
        } else if (type === 'summary') {
            return parser.parseSummary(content);
        } else if (type === 'langs') {
            return parser.parseLanguages(content);
        } else {
            throw new Error('Parsing invalid type "' + type + '"');
        }
    })
    .then(function(result) {
        return [
            file,
            result
        ];
    });
}
646291
def mirtrace_length_plot(self):
        

        data = dict()
        for s_name in self.length_data:
            try:
                data[s_name] = {int(d): int(self.length_data[s_name][d]) for d in self.length_data[s_name]}
            except KeyError:
                pass
        if len(data) == 0:
            log.debug('No valid data for read length distribution')
            return None

        config = {
            'id': 'mirtrace_length_plot',
            'title': 'miRTrace: Read Length Distribution',
            'ylab': 'Read Count',
            'xlab': 'Read Lenth (bp)',
            'ymin': 0,
            'xmin': 0,
            'xDecimals': False,
            'tt_label': '<b>Read Length (bp) {point.x}</b>: {point.y} Read Count',
            'xPlotBands': [
                {'from': 40, 'to': 50, 'color': '#ffebd1'},
                {'from': 26, 'to': 40, 'color': '#e2f5ff'},
                {'from': 18, 'to': 26, 'color': '#e5fce0'},
                {'from': 0, 'to': 18, 'color': '#ffffe2'},
            ]
        }

        return linegraph.plot(data, config)
1364850
private static Object[] extractExtendedFeatures(Configuration configuration, int length)
    {
        Object[] featureMap = new Object[length];

        State state = configuration.state;
        Sentence sentence = configuration.sentence;

        int b0Position = 0;
        int b1Position = 0;
        int b2Position = 0;
        int s0Position = 0;

        long svr = 0; // stack right valency
        long svl = 0; // stack left valency
        long bvl = 0; // buffer left valency

        long b0w = 0;
        long b0p = 0;

        long b1w = 0;
        long b1p = 0;

        long b2w = 0;
        long b2p = 0;

        long s0w = 0;
        long s0p = 0;
        long s0l = 0;

        long bl0p = 0;
        long bl0w = 0;
        long bl0l = 0;

        long bl1w = 0;
        long bl1p = 0;
        long bl1l = 0;

        long sr0p = 0;
        long sr0w = 0;
        long sr0l = 0;

        long sh0w = 0;
        long sh0p = 0;
        long sh0l = 0;

        long sl0p = 0;
        long sl0w = 0;
        long sl0l = 0;

        long sr1w = 0;
        long sr1p = 0;
        long sr1l = 0;

        long sh1w = 0;
        long sh1p = 0;

        long sl1w = 0;
        long sl1p = 0;
        long sl1l = 0;

        long sdl = 0;
        long sdr = 0;
        long bdl = 0;

        int[] words = sentence.getWords();
        int[] tags = sentence.getTags();

        if (0 < state.bufferSize())
        {
            b0Position = state.bufferHead();
            b0w = b0Position == 0 ? 0 : words[b0Position - 1];
            b0w += 2;
            b0p = b0Position == 0 ? 0 : tags[b0Position - 1];
            b0p += 2;
            bvl = state.leftValency(b0Position);

            int leftMost = state.leftMostModifier(state.getBufferItem(0));
            if (leftMost >= 0)
            {
                bl0p = leftMost == 0 ? 0 : tags[leftMost - 1];
                bl0p += 2;
                bl0w = leftMost == 0 ? 0 : words[leftMost - 1];
                bl0w += 2;
                bl0l = state.getDependent(leftMost);
                bl0l += 2;

                int l2 = state.leftMostModifier(leftMost);
                if (l2 >= 0)
                {
                    bl1w = l2 == 0 ? 0 : words[l2 - 1];
                    bl1w += 2;
                    bl1p = l2 == 0 ? 0 : tags[l2 - 1];
                    bl1p += 2;
                    bl1l = state.getDependent(l2);
                    bl1l += 2;
                }
            }

            if (1 < state.bufferSize())
            {
                b1Position = state.getBufferItem(1);
                b1w = b1Position == 0 ? 0 : words[b1Position - 1];
                b1w += 2;
                b1p = b1Position == 0 ? 0 : tags[b1Position - 1];
                b1p += 2;

                if (2 < state.bufferSize())
                {
                    b2Position = state.getBufferItem(2);

                    b2w = b2Position == 0 ? 0 : words[b2Position - 1];
                    b2w += 2;
                    b2p = b2Position == 0 ? 0 : tags[b2Position - 1];
                    b2p += 2;
                }
            }
        }

        if (0 < state.stackSize())
        {
            s0Position = state.stackTop();
            s0w = s0Position == 0 ? 0 : words[s0Position - 1];
            s0w += 2;
            s0p = s0Position == 0 ? 0 : tags[s0Position - 1];
            s0p += 2;
            s0l = state.getDependent(s0Position);
            s0l += 2;

            svl = state.leftValency(s0Position);
            svr = state.rightValency(s0Position);

            int leftMost = state.leftMostModifier(s0Position);
            if (leftMost >= 0)
            {
                sl0p = leftMost == 0 ? 0 : tags[leftMost - 1];
                sl0p += 2;
                sl0w = leftMost == 0 ? 0 : words[leftMost - 1];
                sl0w += 2;
                sl0l = state.getDependent(leftMost);
                sl0l += 2;
            }

            int rightMost = state.rightMostModifier(s0Position);
            if (rightMost >= 0)
            {
                sr0p = rightMost == 0 ? 0 : tags[rightMost - 1];
                sr0p += 2;
                sr0w = rightMost == 0 ? 0 : words[rightMost - 1];
                sr0w += 2;
                sr0l = state.getDependent(rightMost);
                sr0l += 2;
            }

            int headIndex = state.getHead(s0Position);
            if (headIndex >= 0)
            {
                sh0w = headIndex == 0 ? 0 : words[headIndex - 1];
                sh0w += 2;
                sh0p = headIndex == 0 ? 0 : tags[headIndex - 1];
                sh0p += 2;
                sh0l = state.getDependent(headIndex);
                sh0l += 2;
            }

            if (leftMost >= 0)
            {
                int l2 = state.leftMostModifier(leftMost);
                if (l2 >= 0)
                {
                    sl1w = l2 == 0 ? 0 : words[l2 - 1];
                    sl1w += 2;
                    sl1p = l2 == 0 ? 0 : tags[l2 - 1];
                    sl1p += 2;
                    sl1l = state.getDependent(l2);
                    sl1l += 2;
                }
            }
            if (headIndex >= 0)
            {
                if (state.hasHead(headIndex))
                {
                    int h2 = state.getHead(headIndex);
                    sh1w = h2 == 0 ? 0 : words[h2 - 1];
                    sh1w += 2;
                    sh1p = h2 == 0 ? 0 : tags[h2 - 1];
                    sh1p += 2;
                }
            }
            if (rightMost >= 0)
            {
                int r2 = state.rightMostModifier(rightMost);
                if (r2 >= 0)
                {
                    sr1w = r2 == 0 ? 0 : words[r2 - 1];
                    sr1w += 2;
                    sr1p = r2 == 0 ? 0 : tags[r2 - 1];
                    sr1p += 2;
                    sr1l = state.getDependent(r2);
                    sr1l += 2;
                }
            }
        }
        int index = 0;

        long b0wp = b0p;
        b0wp |= (b0w << 8);
        long b1wp = b1p;
        b1wp |= (b1w << 8);
        long s0wp = s0p;
        s0wp |= (s0w << 8);
        long b2wp = b2p;
        b2wp |= (b2w << 8);

        /**
         * From single words
         */
        if (s0w != 1)
        {
            featureMap[index++] = s0wp;
            featureMap[index++] = s0w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p;

        if (b0w != 1)
        {
            featureMap[index++] = b0wp;
            featureMap[index++] = b0w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = b0p;

        if (b1w != 1)
        {
            featureMap[index++] = b1wp;
            featureMap[index++] = b1w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = b1p;

        if (b2w != 1)
        {
            featureMap[index++] = b2wp;
            featureMap[index++] = b2w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = b2p;

        /**
         * from word pairs
         */
        if (s0w != 1 && b0w != 1)
        {
            featureMap[index++] = (s0wp << 28) | b0wp;
            featureMap[index++] = (s0wp << 20) | b0w;
            featureMap[index++] = (s0w << 28) | b0wp;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
            featureMap[index++] = null;
        }

        if (s0w != 1)
        {
            featureMap[index++] = (s0wp << 8) | b0p;
        }
        else
        {
            featureMap[index++] = null;
        }

        if (b0w != 1)
        {
            featureMap[index++] = (s0p << 28) | b0wp;
        }
        else
        {
            featureMap[index++] = null;
        }

        if (s0w != 1 && b0w != 1)
        {
            featureMap[index++] = (s0w << 20) | b0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (s0p << 8) | b0p;
        featureMap[index++] = (b0p << 8) | b1p;

        /**
         * from three words
         */
        featureMap[index++] = (b0p << 16) | (b1p << 8) | b2p;
        featureMap[index++] = (s0p << 16) | (b0p << 8) | b1p;
        featureMap[index++] = (sh0p << 16) | (s0p << 8) | b0p;
        featureMap[index++] = (s0p << 16) | (sl0p << 8) | b0p;
        featureMap[index++] = (s0p << 16) | (sr0p << 8) | b0p;
        featureMap[index++] = (s0p << 16) | (b0p << 8) | bl0p;

        /**
         * distance
         */
        long distance = 0;
        if (s0Position > 0 && b0Position > 0)
            distance = Math.abs(b0Position - s0Position);
        if (s0w != 1)
        {
            featureMap[index++] = s0w | (distance << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (distance << 8);
        if (b0w != 1)
        {
            featureMap[index++] = b0w | (distance << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = b0p | (distance << 8);
        if (s0w != 1 && b0w != 1)
        {
            featureMap[index++] = s0w | (b0w << 20) | (distance << 40);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (b0p << 8) | (distance << 28);

        /**
         * Valency information
         */
        if (s0w != 1)
        {
            featureMap[index++] = s0w | (svr << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (svr << 8);
        if (s0w != 1)
        {
            featureMap[index++] = s0w | (svl << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (svl << 8);
        if (b0w != 1)
        {
            featureMap[index++] = b0w | (bvl << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = b0p | (bvl << 8);

        /**
         * Unigrams
         */
        if (sh0w != 1)
        {
            featureMap[index++] = sh0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sh0p;
        featureMap[index++] = s0l;
        if (sl0w != 1)
        {
            featureMap[index++] = sl0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sl0p;
        featureMap[index++] = sl0l;
        if (sr0w != 1)
        {
            featureMap[index++] = sr0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sr0p;
        featureMap[index++] = sr0l;
        if (bl0w != 1)
        {
            featureMap[index++] = bl0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = bl0p;
        featureMap[index++] = bl0l;

        /**
         * From third order features
         */
        if (sh1w != 1)
        {
            featureMap[index++] = sh1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sh1p;
        featureMap[index++] = sh0l;
        if (sl1w != 1)
        {
            featureMap[index++] = sl1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sl1p;
        featureMap[index++] = sl1l;
        if (sr1w != 1)
        {
            featureMap[index++] = sr1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sr1p;
        featureMap[index++] = sr1l;
        if (bl1w != 1)
        {
            featureMap[index++] = bl1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = bl1p;
        featureMap[index++] = bl1l;
        featureMap[index++] = s0p | (sl0p << 8) | (sl1p << 16);
        featureMap[index++] = s0p | (sr0p << 8) | (sr1p << 16);
        featureMap[index++] = s0p | (sh0p << 8) | (sh1p << 16);
        featureMap[index++] = b0p | (bl0p << 8) | (bl1p << 16);

        /**
         * label set
         */
        if (s0Position >= 0)
        {
            sdl = state.leftDependentLabels(s0Position);
            sdr = state.rightDependentLabels(s0Position);
        }

        if (b0Position >= 0)
        {
            bdl = state.leftDependentLabels(b0Position);
        }

        if (s0w != 1)
        {
            featureMap[index++] = (s0w + "|" + sdr);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (s0p + "|" + sdr);
        if (s0w != 1)
        {
            featureMap[index++] = s0w + "|" + sdl;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (s0p + "|" + sdl);
        if (b0w != 1)
        {
            featureMap[index++] = (b0w + "|" + bdl);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (b0p + "|" + bdl);
        return featureMap;
    }
528416
def parse_args():
    ''''''
    parser = argparse.ArgumentParser(
        description='Morphology fit distribution extractor',
        epilog='Note: Outputs json of the optimal distribution \
                and corresponding parameters.')

    parser.add_argument('datapaths',
                        nargs='+',
                        help='Morphology data directory paths')

    parser.add_argument('--mtype', choices=MTYPE_GETTERS.keys(),
                        help='Get mtype from filename or parent directory')

    return parser.parse_args()
1366224
def multihistogram(args):
    
    p = OptionParser(multihistogram.__doc__)
    p.add_option("--kmin", default=15, type="int",
            help="Minimum K-mer size, inclusive")
    p.add_option("--kmax", default=30, type="int",
            help="Maximum K-mer size, inclusive")
    p.add_option("--vmin", default=2, type="int",
            help="Minimum value, inclusive")
    p.add_option("--vmax", default=100, type="int",
            help="Maximum value, inclusive")
    opts, args, iopts = p.set_image_options(args, figsize="10x5", dpi=300)

    if len(args) < 1:
        sys.exit(not p.print_help())

    histfiles = args[:-1]
    species = args[-1]
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])
    A = fig.add_axes([.08, .12, .38, .76])
    B = fig.add_axes([.58, .12, .38, .76])

    lines = []
    legends = []
    genomesizes = []
    for histfile in histfiles:
        ks = KmerSpectrum(histfile)
        x, y = ks.get_xy(opts.vmin, opts.vmax)
        K = get_number(op.basename(histfile).split(".")[0].split("-")[-1])
        if not opts.kmin <= K <= opts.kmax:
            continue

        line, = A.plot(x, y, '-', lw=1)
        lines.append(line)
        legends.append("K = {0}".format(K))
        ks.analyze(K=K)
        genomesizes.append((K, ks.genomesize / 1e6))

    leg = A.legend(lines, legends, shadow=True, fancybox=True)
    leg.get_frame().set_alpha(.5)

    title = "{0} genome K-mer histogram".format(species)
    A.set_title(markup(title))
    xlabel, ylabel = "Coverage (X)", "Counts"
    A.set_xlabel(xlabel)
    A.set_ylabel(ylabel)
    set_human_axis(A)

    title = "{0} genome size estimate".format(species)
    B.set_title(markup(title))
    x, y = zip(*genomesizes)
    B.plot(x, y, "ko", mfc='w')
    t = np.linspace(opts.kmin - .5, opts.kmax + .5, 100)
    p = np.poly1d(np.polyfit(x, y, 2))
    B.plot(t, p(t), "r:")

    xlabel, ylabel = "K-mer size", "Estimated genome size (Mb)"
    B.set_xlabel(xlabel)
    B.set_ylabel(ylabel)
    set_ticklabels_helvetica(B)

    labels = ((.04, .96, 'A'), (.54, .96, 'B'))
    panel_labels(root, labels)

    normalize_axes(root)
    imagename = species + ".multiK.pdf"
    savefig(imagename, dpi=iopts.dpi, iopts=iopts)
956898
function parse_PtgName(blob, length, opts) {
	var type = (blob.read_shift(1) >>> 5) & 0x03;
	var w = (!opts || (opts.biff >= 8)) ? 4 : 2;
	var nameindex = blob.read_shift(w);
	switch(opts.biff) {
		case 2: blob.l += 5; break;
		case 3: case 4: blob.l += 8; break;
		case 5: blob.l += 12; break;
	}
	return [type, 0, nameindex];
}
904854
def parse_file(self, cu, analysis):
        """"""
        if hasattr(analysis, 'parser'):
            filename = cu.file_locator.relative_filename(cu.filename)
            source_lines = analysis.parser.lines
            with cu.source_file() as source_file:
                source = source_file.read()

            try:
                if sys.version_info < (3, 0):
                    encoding = source_encoding(source)
                    if encoding != 'utf-8':
                        source = source.decode(encoding).encode('utf-8')
            except UnicodeDecodeError:
                log.warning(
                    'Source file %s can not be properly decoded, skipping. '
                    'Please check if encoding declaration is ok',
                    os.path.basename(cu.filename))
                return
        else:
            if hasattr(cu, 'relative_filename'):
                filename = cu.relative_filename()
            else:
                filename = analysis.coverage.file_locator.relative_filename(
                    cu.filename)

            token_lines = analysis.file_reporter.source_token_lines()
            source_lines = list(enumerate(token_lines))
            source = analysis.file_reporter.source()

        coverage_lines = [self.get_hits(i, analysis)
                          for i in range(1, len(source_lines) + 1)]

        # ensure results are properly merged between platforms
        posix_filename = filename.replace(os.path.sep, '/')

        results = {
            'name': posix_filename,
            'source': source,
            'coverage': coverage_lines,
        }

        branches = self.get_arcs(analysis)
        if branches:
            results['branches'] = branches

        self.source_files.append(results)
172114
def plots_from_files(imspaths, figsize=(10,5), rows=1, titles=None, maintitle=None):
    
    f = plt.figure(figsize=figsize)
    if maintitle is not None: plt.suptitle(maintitle, fontsize=16)
    for i in range(len(imspaths)):
        sp = f.add_subplot(rows, ceildiv(len(imspaths), rows), i+1)
        sp.axis('Off')
        if titles is not None: sp.set_title(titles[i], fontsize=16)
        img = plt.imread(imspaths[i])
        plt.imshow(img)
25971
def gather_parsing_stats(self):
        """"""
        if self.verbose:
            start_time = time.clock()
            try:
                yield
            finally:
                elapsed_time = time.clock() - start_time
                printerr("Time while parsing:", elapsed_time, "seconds")
                if packrat_cache:
                    hits, misses = ParserElement.packrat_cache_stats
                    printerr("Packrat parsing stats:", hits, "hits;", misses, "misses")
        else:
            yield
352460
def _parse_tags(tag_file):
    

    tag_name = None
    tag_value = None

    # Line folding is handled by yielding values only after we encounter
    # the start of a new tag, or if we pass the EOF.
    for num, line in enumerate(tag_file):
        # If byte-order mark ignore it for now.
        if num == 0:
            if line.startswith(BOM):
                line = line.lstrip(BOM)
        # Skip over any empty or blank lines.
        if len(line) == 0 or line.isspace():
            continue
        elif line[0].isspace() and tag_value is not None:  # folded line
            tag_value += line
        else:
            # Starting a new tag; yield the last one.
            if tag_name:
                yield (tag_name, tag_value.strip())

            if ':' not in line:
                raise BagValidationError("invalid line '%s' in %s" % (line.strip(),
                                                                      os.path.basename(tag_file.name)))

            parts = line.strip().split(':', 1)
            tag_name = parts[0].strip()
            tag_value = parts[1]

    # Passed the EOF.  All done after this.
    if tag_name:
        yield (tag_name, tag_value.strip())
953076
def stats(args):
    
    from jcvi.utils.table import tabulate

    p = OptionParser(stats.__doc__)
    p.add_option("--warn", default=False, action="store_true",
                 help="Warnings on small component spans [default: %default]")
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(p.print_help())

    agpfile, = args

    agp = AGP(agpfile)
    gap_lengths = []
    component_lengths = []
    for a in agp:
        span = a.object_span
        if a.is_gap:
            label = a.gap_type
            gap_lengths.append((span, label))
        else:
            label = "{0}:{1}-{2}".format(a.component_id, a.component_beg, \
                   a.component_end)
            component_lengths.append((span, label))
            if opts.warn and span < 50:
                logging.error("component span too small ({0}):\n{1}".\
                    format(span, a))

    table = dict()
    for label, lengths in zip(("Gaps", "Components"),
            (gap_lengths, component_lengths)):

        if not lengths:
            table[(label, "Min")] = table[(label, "Max")] \
                                  = table[(label, "Sum")] = "n.a."
            continue

        table[(label, "Min")] = "{0} ({1})".format(*min(lengths))
        table[(label, "Max")] = "{0} ({1})".format(*max(lengths))
        table[(label, "Sum")] = sum(x[0] for x in lengths)

    print(tabulate(table), file=sys.stderr)