>>> data = get_bed_from_phytozome(["AT5G54690", "AT1G01010"])
    >>> print data.read()  #doctest: +NORMALIZE_WHITESPACE
    Chr1	3631	5899	AT1G01010
    Chr5	22219224	22221840	AT5G54690
    <BLANKLINE>
546249
def get_bed_from_phytozome(genelist):
    
    genelist = ",".join(genelist)
    dataset = get_phytozome_dataset()
    filters = dict(gene_name_filter=genelist)
    attributes = "chr_name1,gene_chrom_start,gene_chrom_end,gene_name1".split(",")

    data = dataset.query(filters=filters, attributes=attributes)
    return data
1019750
def apply(self, read_tuple_name, read_tuple_id=None, synchronize_widths=True):
        
        parts = read_tuple_name.split("__")
        parts[0] = self._fill_right(parts[0], "-", self.prefix_width)
        if read_tuple_id is not None:
            parts[1] = "{:x}".format(read_tuple_id)
        parts[1] = self._fill_left(parts[1], "0", self.read_tuple_id_width)

        if synchronize_widths:
            new_segments = []
            segments = parts[2][1:-1].split("),(")
            for segment in segments:
                values = segment.split(",")
                values[0] = values[0].zfill(self.genome_id_width)
                values[1] = values[1].zfill(self.chr_id_width)
                values[3] = values[3].zfill(self.coor_width)
                values[4] = values[4].zfill(self.coor_width)
                new_segments.append("(" + ",".join(values) + ")")
            parts[2] = ",".join(new_segments)

        return "__".join(parts)
1242986
def _add_genes_to_bed(in_file, gene_file, fai_file, out_file, data, max_distance=10000):
    
    try:
        input_rec = next(iter(pybedtools.BedTool(in_file)))
    except StopIteration:  # empty file
        utils.copy_plus(in_file, out_file)
        return
    # keep everything after standard chrom/start/end, 1-based
    extra_fields = list(range(4, len(input_rec.fields) + 1))
    # keep the new gene annotation
    gene_index = len(input_rec.fields) + 4
    extra_fields.append(gene_index)
    columns = ",".join([str(x) for x in extra_fields])
    max_column = max(extra_fields) + 1
    ops = ",".join(["distinct"] * len(extra_fields))
    # swap over gene name to '.' if beyond maximum distance
    # cut removes the last distance column which can cause issues
    # with bedtools merge: 'ERROR: illegal character '.' found in integer conversion of string'
    distance_filter = (r"""awk -F$'\t' -v OFS='\t' '{if ($NF > %s || $NF < -%s) $%s = "."} {print}'""" %
                       (max_distance, max_distance, gene_index))
    sort_cmd = bedutils.get_sort_cmd(os.path.dirname(out_file))
    cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
    # Ensure gene transcripts match reference genome
    ready_gene_file = os.path.join(os.path.dirname(out_file), "%s-genomeonly.bed" %
                                   (utils.splitext_plus(os.path.basename(gene_file))[0]))
    ready_gene_file = bedutils.subset_to_genome(gene_file, ready_gene_file, data)
    exports = "export TMPDIR=%s && %s" % (os.path.dirname(out_file), utils.local_path_export())
    bcbio_py = sys.executable
    gsort = config_utils.get_program("gsort", data)
    cmd = ("{exports}{cat_cmd} {in_file} | grep -v ^track | grep -v ^browser | grep -v ^# | "
           "{bcbio_py} -c 'from bcbio.variation import bedutils; bedutils.remove_bad()' | "
           "{gsort} - {fai_file} | "
            "bedtools closest -g {fai_file} "
            "-D ref -t first -a - -b <({gsort} {ready_gene_file} {fai_file}) | "
            "{distance_filter} | cut -f 1-{max_column} | "
            "bedtools merge -i - -c {columns} -o {ops} -delim ',' -d -10 > {out_file}")
    do.run(cmd.format(**locals()), "Annotate BED file with gene info")
427795
def loss(args):
    
    p = OptionParser(loss.__doc__)
    p.add_option("--bed", default=False, action="store_true",
                 help="Genomic BLAST is in bed format [default: %default]")
    p.add_option("--gdist", default=20, type="int",
                 help="Gene distance [default: %default]")
    p.add_option("--bdist", default=20000, type="int",
                 help="Base pair distance [default: %default]")
    p.set_beds()
    opts, args = p.parse_args(args)

    if len(args) not in (1, 2):
        sys.exit(not p.print_help())

    blocksfile = args[0]
    emptyblast = (len(args) == 1)
    if emptyblast:
        genomicblast = "empty.blast"
        sh("touch {0}".format(genomicblast))
    else:
        genomicblast = args[1]

    gdist, bdist = opts.gdist, opts.bdist
    qbed, sbed, qorder, sorder, is_self = check_beds(blocksfile, p, opts)
    blocks = []
    fp = open(blocksfile)
    genetrack = {}
    proxytrack = {}
    for row in fp:
        a, b = row.split()
        genetrack[a] = b
        blocks.append((a, b))

    data = []
    for key, rows in groupby(blocks, key=lambda x: x[-1]):
        rows = list(rows)
        data.append((key, rows))

    imax = len(data) - 1
    for i, (key, rows) in enumerate(data):
        if i == 0 or i == imax:
            continue
        if key != '.':
            continue

        before, br = data[i - 1]
        after, ar = data[i + 1]
        bi, bx = sorder[before]
        ai, ax = sorder[after]
        dist = abs(bi - ai)
        if bx.seqid != ax.seqid or dist > gdist:
            continue

        start, end = range_minmax(((bx.start, bx.end), (ax.start, ax.end)))
        start, end = max(start - bdist, 1), end + bdist
        proxy = (bx.seqid, start, end)
        for a, b in rows:
            proxytrack[a] = proxy

    tags = {}
    if opts.bed:
        bed = Bed(genomicblast, sorted=False)
        key = lambda x: gene_name(x.accn.rsplit(".", 1)[0])
        for query, bb in groupby(bed, key=key):
            bb = list(bb)
            if query not in proxytrack:
                continue

            proxy = proxytrack[query]
            tag = "NS"
            best_b = bb[0]
            for b in bb:
                hsp = (b.seqid, b.start, b.end)
                if range_overlap(proxy, hsp):
                    tag = "S"
                    best_b = b
                    break

            hsp = (best_b.seqid, best_b.start, best_b.end)
            proxytrack[query] = hsp
            tags[query] = tag

    else:
        blast = Blast(genomicblast)
        for query, bb in blast.iter_hits():
            bb = list(bb)
            query = gene_name(query)
            if query not in proxytrack:
                continue

            proxy = proxytrack[query]
            tag = "NS"
            best_b = bb[0]
            for b in bb:
                hsp = (b.subject, b.sstart, b.sstop)
                if range_overlap(proxy, hsp):
                    tag = "S"
                    best_b = b
                    break

            hsp = (best_b.subject, best_b.sstart, best_b.sstop)
            proxytrack[query] = hsp
            tags[query] = tag

    for b in qbed:
        accn = b.accn
        target_region = genetrack[accn]
        if accn in proxytrack:
            target_region = region_str(proxytrack[accn])
            if accn in tags:
                ptag = "[{0}]".format(tags[accn])
            else:
                ptag = "[NF]"
            target_region = ptag + target_region

        print("\t".join((b.seqid, accn, target_region)))

    if emptyblast:
        sh("rm -f {0}".format(genomicblast))
431769
def parse_tmhmm_long(tmhmm_results):
    
    with open(tmhmm_results) as f:
        lines = f.read().splitlines()

    infodict = defaultdict(dict)

    for l in lines:
        if 'Number of predicted TMHs:' in l:
            gene = l.split(' Number')[0].strip('# ')
            infodict[gene]['num_tm_helices'] = int(l.split(': ')[1])

        if 'WARNING' in l:
            log.warning('{}: no TMHMM predictions'.format(l))
            continue

        # TODO: POSSIBLE N-term signal sequence - parse this

        # Look for the lines without #, these are the TM predicted regions
        if '#' not in l:
            stuff = l.split()
            if stuff[1] == 'TMHMM2.0':

                gene = stuff[0]
                region = stuff[2]
                region_start = stuff[3]
                region_end = stuff[4]

                if 'sequence' in infodict[gene]:
                    tm_seq = infodict[gene]['sequence']
                else:
                    tm_seq = ''

                if region == 'outside':
                    info = 'O'
                elif region == 'inside':
                    info = 'I'
                elif region == 'TMhelix':
                    info = 'T'
                else:
                    log.error('{}: unknown region type'.format(info))
                    info = '-'

                for r in range(int(region_start), int(region_end) + 1):
                    tm_seq += info

                infodict[gene]['sequence'] = tm_seq

    return infodict
556459
def _get_regional_bed_file(data):
    
    variant_regions = bedutils.merge_overlaps(tz.get_in(["config", "algorithm", "variant_regions"], data),
                                              data)
    is_genome = data["config"]["algorithm"].get("coverage_interval", "exome").lower() in ["genome"]
    if variant_regions and utils.file_exists(variant_regions) and not is_genome:
        return variant_regions
805578
def from23andme(args):
    
    p = OptionParser(from23andme.__doc__)
    p.set_ref()
    opts, args = p.parse_args(args)

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

    txtfile, seqid = args
    ref_dir = opts.ref
    fastafile = op.join(ref_dir, "hs37d5.fa")
    fasta = Fasta(fastafile)

    pf = txtfile.rsplit(".", 1)[0]
    px = CM[seqid]
    chrvcf = pf + ".{0}.vcf".format(px)
    legend = op.join(ref_dir, "1000GP_Phase3/{0}.rsids".format(px))
    register = read_rsid(seqid, legend)

    fw = open(chrvcf, "w")
    print(get_vcfstanza(fastafile, fasta, txtfile), file=fw)

    fp = open(txtfile)
    seen = set()
    duplicates = skipped = missing = 0
    for row in fp:
        if row[0] == '#':
            continue
        rsid, chr, pos, genotype = row.split()
        if chr != seqid:
            continue
        pos = int(pos)
        if (chr, pos) in seen:
            duplicates += 1
            continue
        seen.add((chr, pos))
        genotype = list(genotype)
        if "-" in genotype:  # missing daa
            missing += 1
            continue

        # Y or MT
        if not register:
            assert len(genotype) == 1
            ref = fasta[chr][pos - 1].seq.upper()
            if "D" in genotype or "I" in genotype:
                skipped += 1
                continue
            genotype = genotype[0]
            code = "0/0" if ref == genotype else "1/1"
            alt = "." if ref == genotype else genotype
            print("\t".join(str(x) for x in \
                    (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)
            continue

        # If rsid is seen in the db, use that
        if rsid in register:
            pos, ref, alt = register[rsid]
        elif pos in register:
            pos, ref, alt = register[pos]
        else:
            skipped += 1  # Not in reference panel
            continue

        assert fasta[chr][pos - 1:pos + len(ref) - 1].seq.upper() == ref
        # Keep it bi-allelic
        not_seen = [x for x in alt if x not in genotype]
        while len(alt) > 1 and not_seen:
            alt.remove(not_seen.pop())
        if len(alt) > 1:
            alt = [alt[0]]
        alleles = [ref] + alt

        if len(genotype) == 1:
            genotype = [genotype[0]] * 2

        alt = ",".join(alt) or "."
        if "D" in genotype or "I" in genotype:
            max_allele = max((len(x), x) for x in alleles)[1]
            alleles = [("I" if x == max_allele else "D") for x in alleles]
            assert "I" in alleles and "D" in alleles
        a, b = genotype
        try:
            ia, ib = alleles.index(a), alleles.index(b)
        except ValueError:  # alleles not seen
            logging.error("{0}: alleles={1}, genotype={2}".\
                            format(rsid, alleles, genotype))
            skipped += 1
            continue
        code = "/".join(str(x) for x in sorted((ia, ib)))

        print("\t".join(str(x) for x in \
                (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)

    logging.debug("duplicates={0} skipped={1} missing={2}".\
                    format(duplicates, skipped, missing))
1026404
def _find_gene_list_from_bed(bed_file, base_file, data):
    
    # Check for a gene list, we can just return that.
    if is_gene_list(bed_file):
        return bed_file
    out_file = "%s-genes.txt" % utils.splitext_plus(base_file)[0]
    if not os.path.exists(out_file):
        genes = set([])
        import pybedtools
        with utils.open_gzipsafe(bed_file) as in_handle:
            for r in pybedtools.BedTool(in_handle):
                if r.name:
                    if not r.name.startswith("{"):
                        genes.add(r.name)
        with file_transaction(data, out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                if len(genes) > 0:
                    out_handle.write("\n".join(sorted(list(genes))) + "\n")
    if utils.file_exists(out_file):
        return out_file
510355
def is_gene_list(bed_file):
    
    with utils.open_gzipsafe(bed_file) as in_handle:
        for line in in_handle:
            if not line.startswith("#"):
                if len(line.split()) == 1:
                    return True
                else:
                    return False
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)
1080368
@Override
  public void readData(BufferedReader in) throws IOException {
    final String LEFT = "left";
    int lineNum = 1;
    // all lines have one rule per line
    boolean doingStop = false;

    for (String line = in.readLine(); line != null && line.length() > 0; line = in.readLine()) {
      try {
        if (line.equals("BEGIN_STOP")) {
          doingStop = true;
          continue;
        }
        String[] fields = StringUtils.splitOnCharWithQuoting(line, ' ', '\"', '\\'); // split on spaces, quote with doublequote, and escape with backslash
        //        System.out.println("fields:\n" + fields[0] + "\n" + fields[1] + "\n" + fields[2] + "\n" + fields[3] + "\n" + fields[4] + "\n" + fields[5]);

        
        short distance = (short)Integer.parseInt(fields[4]);
        IntTaggedWord tempHead = new IntTaggedWord(fields[0], '/', wordIndex, tagIndex);
        IntTaggedWord tempArg = new IntTaggedWord(fields[2], '/', wordIndex, tagIndex);
        IntDependency tempDependency = new IntDependency(tempHead, tempArg, fields[3].equals(LEFT), distance);

        double count = Double.parseDouble(fields[5]);
        if (doingStop) {
          expandStop(tempDependency, distance, count, false);
        } else {
          expandArg(tempDependency, distance, count);
        }
      } catch (Exception e) {
        IOException ioe = new IOException("Error on line " + lineNum + ": " + line);
        ioe.initCause(e);
        throw ioe;
      }
      //      System.out.println("read line " + lineNum + ": " + line);
      lineNum++;
    }
  }
1348160
def wgsim(args):
    
    p = OptionParser(wgsim.__doc__)
    p.add_option("--erate", default=.01, type="float",
                 help="Base error rate of the read [default: %default]")
    p.add_option("--noerrors", default=False, action="store_true",
                 help="Simulate reads with no errors [default: %default]")
    p.add_option("--genomesize", type="int",
                 help="Genome size in Mb [default: estimate from data]")
    add_sim_options(p)
    opts, args = p.parse_args(args)

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

    fastafile, = args
    pf = op.basename(fastafile).split(".")[0]

    genomesize = opts.genomesize
    size = genomesize * 1000000 if genomesize else Fasta(fastafile).totalsize
    depth = opts.depth
    readlen = opts.readlen
    readnum = int(math.ceil(size * depth / (2 * readlen)))

    distance = opts.distance
    stdev = distance / 10

    outpf = opts.outfile or "{0}.{1}bp.{2}x".format(pf, distance, depth)

    logging.debug("Total genome size: {0} bp".format(size))
    logging.debug("Target depth: {0}x".format(depth))
    logging.debug("Number of read pairs (2x{0}): {1}".format(readlen, readnum))

    if opts.noerrors:
        opts.erate = 0

    cmd = "dwgsim -e {0} -E {0}".format(opts.erate)
    if opts.noerrors:
        cmd += " -r 0 -R 0 -X 0 -y 0"

    cmd += " -d {0} -s {1}".format(distance, stdev)
    cmd += " -N {0} -1 {1} -2 {1}".format(readnum, readlen)
    cmd += " {0} {1}".format(fastafile, outpf)
    sh(cmd)
474198
def readstoredmetadata(fname, skiprows=0, linenumber=None, comments='#', 
                       metametadata=None, verbosity=DEFAULT_VERBOSITY):
    
    if is_string_like(fname):
        fh = file(fname, 'rU')
    elif hasattr(fname, 'readline'):
        fh = fname
    else:
        raise ValueError('fname must be a string or file handle')
       
    if not metametadata:
        for _ind in range(skiprows):
            fh.readline()

        phlines = []
        if linenumber is None:
            if comments:
                for line in fh:
                    if not line.startswith(comments):
                        if len(phlines) == 0:
                            phlines = [line]
                        break
                    else:
                        phlines.append(line)
                if len(phlines) == 1:
                    if verbosity >= 10:
                        print '''Looking for metametadata on line 0 
                                 (no comment lines present).'''
                else:
                    if verbosity >= 9:
                        print '''Searching for metametadata lines up to 
                                 and including line %d where comments 
                                 end.''' % (len(phlines) - 1)
            else:
                phlines = [fh.readline()]
                if verbosity >=9:
                    print '''No comments found, looking for metametadata on 
                             line 0.'''
        else:
            for _ind in range(linenumber):
                fh.readline()
            phlines = [fh.readline()]
                               
        metametadata = None
        for (ln, metametaline) in enumerate(phlines):
            s = re.compile(r'metametadata[\s]*=[\s]*{').search(metametaline)
            if s:
                l = s.start()
                if len(uniqify(metametaline[:l])) <= 1:
                    metametaline = metametaline[l:].rstrip()
                    try:
                        X = compiler.parse(metametaline)
                    except SyntaxError:
                        pass
                    else:
                        if IsMetaMetaDict(X):
                            exec metametaline
                            if verbosity > 6:
                                print 'Found valid metametadata at line', ln, 'in file.  Metametadata is:', metametadata
                            break
   
    if metametadata:
        
        metadata = {}
        metadata['metametadata'] = metametadata
        Tval = max([v if isinstance(v,int) else max(v) 
                    for v in metametadata.values()])
        fh = file(fname, 'rU')
        data = [fh.readline() for _ind in range(Tval + 1 + skiprows)][skiprows:]
 
        if (max([v if isinstance(v,int) else max(v) 
                for v in  metametadata.values()]) < len(data)):
            for n in metametadata.keys():
                [s, e] = [metametadata[n], metametadata[n]+1] \
                          if isinstance(metametadata[n],int) \
                          else [metametadata[n][0],metametadata[n][1]]
                metadata[n] = ''.join(data[s:e]).strip('\n')
            processmetadata(metadata, comments=comments, verbosity=verbosity)
        
            return metadata
1369968
def clean_file(in_file, data, prefix="", bedprep_dir=None, simple=None):
    
    # Remove non-ascii characters. Used in coverage analysis, to support JSON code in one column
    #   and be happy with sambamba:
    simple = "iconv -c -f utf-8 -t ascii | sed 's/ //g' |" if simple else ""
    if in_file:
        if not bedprep_dir:
            bedprep_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "bedprep"))
        # Avoid running multiple times with same prefix
        if prefix and os.path.basename(in_file).startswith(prefix):
            return in_file
        out_file = os.path.join(bedprep_dir, "%s%s" % (prefix, os.path.basename(in_file)))
        out_file = out_file.replace(".interval_list", ".bed")
        if out_file.endswith(".gz"):
            out_file = out_file[:-3]
        if not utils.file_uptodate(out_file, in_file):
            check_bed_contigs(in_file, data)
            check_bed_coords(in_file, data)
            with file_transaction(data, out_file) as tx_out_file:
                bcbio_py = sys.executable
                cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
                sort_cmd = get_sort_cmd(os.path.dirname(tx_out_file))
                cmd = ("{cat_cmd} {in_file} | grep -v ^track | grep -v ^browser | grep -v ^@ | "
                       "grep -v ^# | {simple} "
                       "{bcbio_py} -c 'from bcbio.variation import bedutils; bedutils.remove_bad()' | "
                       "{sort_cmd} -k1,1 -k2,2n > {tx_out_file}")
                do.run(cmd.format(**locals()), "Prepare cleaned BED file", data)
        vcfutils.bgzip_and_index(out_file, data.get("config", {}), remove_orig=False)
        return out_file
1278208
def write_last_and_beds(pf, GenePositions, ContigStarts):
    
    qbedfile = pf + "tigs.bed"
    sbedfile = pf + "chr.bed"
    lastfile = "{}tigs.{}chr.last".format(pf, pf)
    qbedfw = open(qbedfile, "w")
    sbedfw = open(sbedfile, "w")
    lastfw = open(lastfile, "w")

    GeneContigs = np.searchsorted(ContigStarts, GenePositions) - 1
    for i, (c, gstart) in enumerate(zip(GeneContigs, GenePositions)):
        gene = "gene{:05d}".format(i)
        tig = "tig{:04d}".format(c)
        start = ContigStarts[c]
        cstart = gstart - start
        print("\t".join(str(x) for x in
                                   (tig, cstart, cstart + 1, gene)), file=qbedfw)
        print("\t".join(str(x) for x in
                                   ("chr1", gstart, gstart + 1, gene)), file=sbedfw)
        lastatoms = [gene, gene, 100] + [0] * 8 + [100]
        print("\t".join(str(x) for x in lastatoms), file=lastfw)

    qbedfw.close()
    sbedfw.close()
    lastfw.close()
460551
def GetPeaksExons(bed,parsedGTF):
    

    bedtool_AB=dfTObedtool(bed)

    exonsGTF=parsedGTF[parsedGTF["feature"]=="exon"]
    exonsGTF.reset_index(inplace=True, drop=True)

    exonsBED=GTFtoBED(exonsGTF, "exon_id")
    exonsBED.columns=['chrom', 'chromStart', 'chromEnd', 'exon_id', 'score', 'strand']
    exonsBEDcols=exonsBED.columns.tolist()

    bedcols=bed.columns.tolist()
    exonsBEDcols_=[]
    for c in exonsBEDcols:
        if c in bedcols:
            exonsBEDcols_.append(c+"_exon")
        else:
            exonsBEDcols_.append(c)

    cols=[bedcols,exonsBEDcols_,["overlap"] ]
    cols=[item for sublist in cols for item in sublist]

    bedtool_exons=dfTObedtool(exonsBED)

    bedtool_target_exons=bedtool_AB.intersect(bedtool_exons, wo=True, s=True)
    dfTargetE=pd.read_table(bedtool_target_exons.fn, names=cols)
    ExonsTransGenes=parsedGTF[["exon_id","transcript_id","gene_id"]].drop_duplicates()
    dfTargets=pd.merge(dfTargetE,ExonsTransGenes,on=["exon_id"],how="left")
    dfTargets["count"]=1

    def getCounts(df,field):
        """
        For each field in a bed narrow peak returns the number or times that field is present,\
        the normalized mean of the '-log10(pValue)' and normalized mean of the signal value.

        :param df: a Pandas dataframe of a bed narrow peak
        :param field: field to analyse, ie. exons or transcripts

        :returns: a Pandas dataframe
        """

        tmp=df[[field,'name',"count"]].drop_duplicates()
        tmp=tmp.drop(["name"],axis=1)
        tmp["count"]=tmp["count"].astype(int)
        tmp.columns=[field,"%s_count" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")

        tmp=df[[field,'name',"-log10(pValue)"]].drop_duplicates()
        tmp=tmp.drop(["name"],axis=1)
        tmp["-log10(pValue)"]=tmp["-log10(pValue)"].astype(float)
        tmp=tmp.groupby(field).apply(lambda l: reduce(lambda x, y: x*y, l["-log10(pValue)"]) )
        tmp=pd.DataFrame(tmp)
        tmp.reset_index(inplace=True,drop=False)
        tmp.columns=[field,"%s norm. mean -log10(pValue)" %str(field)]
        df=pd.merge(df,tmp,on=field,how="left")

        tmp=df[[field,'name',"signalValue"]].drop_duplicates()
        tmp=tmp.drop(["name"],axis=1)
        tmp["signalValue"]=tmp["signalValue"].astype(float)
        tmp=tmp.groupby(field).apply(lambda l: reduce(lambda x, y: x*y, l["signalValue"]) )
        tmp=pd.DataFrame(tmp)
        tmp.reset_index(inplace=True,drop=False)
        tmp.columns=[field,"%s signalValue" %str(field)]
        df=pd.merge(df,tmp,on=field,how="left")

        return df

    for f in ["exon_id","transcript_id"]:
        dfTargets=getCounts(dfTargets,f)

    def getCounts_GeneIDs(df):
        """
        For each gene id in a bed narrow peak returns the number or times that field is present,\
        the normalized mean of the '-log10(pValue)' and normalized mean of the signal value.

        :param df: a Pandas dataframe of a bed narrow peak

        :returns: a Pandas dataframe
        """

        field="gene_id"

        tmp=df[[field,"transcript_id","transcript_id_count"]].drop_duplicates()
        tmp=tmp.drop(["transcript_id"],axis=1)
        tmp["transcript_id_count"]=tmp["transcript_id_count"].astype(int)
        tmp.columns=[field,"%s_count" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")

        tmp=df[[field,'transcript_id',"transcript_id norm. mean -log10(pValue)"]].drop_duplicates()
        tmp=tmp.drop(["transcript_id"],axis=1)
        tmp["transcript_id norm. mean -log10(pValue)"]=tmp["transcript_id norm. mean -log10(pValue)"].astype(float)
        tmp.columns=[field,"%s norm. mean -log10(pValue)" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")



        tmp=df[[field,'transcript_id',"transcript_id signalValue"]].drop_duplicates()
        tmp=tmp.drop(["transcript_id"],axis=1)
        tmp["transcript_id signalValue"]=tmp["transcript_id signalValue"].astype(float)
        tmp.columns=[field,"%s signalValue" %str(field)]
        tmp=tmp.groupby(field, as_index=False).sum()
        df=pd.merge(df,tmp,on=field,how="left")

        return df

    dfTargets=getCounts_GeneIDs(dfTargets)


    dfTargets=dfTargets.drop(["count"],axis=1)
    return dfTargets
959592
def bed12(self, score="0", rgb="."):
        
        if not self.is_gene_pred:
            raise CruzException("can't create bed12 from non genepred feature")
        exons = list(self.exons)
        # go from global start, stop, to relative start, length...
        sizes = ",".join([str(e[1] - e[0]) for e in exons]) + ","
        starts = ",".join([str(e[0] - self.txStart) for e in exons]) + ","
        name = self.name2 + "," + self.name if hasattr(self, "name2") \
                                            else self.name
        return "\t".join(map(str, (
            self.chrom, self.txStart, self.txEnd, name,
            score, self.strand, self.cdsStart, self.cdsEnd, rgb,
            len(exons), sizes, starts)))