parse catalytic RNAs to gff format
194389
def parse_catalytic(insertion, gff):
    
    offset = insertion['offset']
    GeneStrand = insertion['strand']
    if type(insertion['intron']) is not str:
        return gff
    for intron in parse_fasta(insertion['intron'].split('|')):
        ID, annot, strand, pos = intron[0].split('>')[1].split()
        Start, End = [int(i) for i in pos.split('-')]
        if strand != GeneStrand:
            if strand == '+':
                strand = '-'
            else:
                strand = '+'
            Start, End = End - 2, Start - 2
        Start, End = abs(Start + offset) - 1, abs(End + offset) - 1
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append('Rfam')
        gff['feature'].append('Catalytic RNA')
        gff['start'].append(Start)
        gff['end'].append(End)
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('ID=%s; Name=%s' % (ID, annot))
    return gff
823121
def _parse_rnaseq_qualimap_metrics(report_file):
    
    from bs4 import BeautifulSoup
    out = {}
    parsers = ["Reads alignment", "Reads genomic origin", "Transcript coverage profile"]
    with open(report_file) as in_handle:
        root = BeautifulSoup(in_handle.read(), "html.parser")
    for table in root.find_all("div", class_="table-summary"):
        h3 = table.find("h3")
        if h3.text in parsers:
            out.update(_parse_qualimap_rnaseq(table.find("table")))
    return out
750573
def _parse_parallel_sentences(f1, f2):
  """"""
  def _parse_text(path):
    """Returns the sentences from a single text file, which may be gzipped."""
    split_path = path.split(".")

    if split_path[-1] == "gz":
      lang = split_path[-2]
      with tf.io.gfile.GFile(path) as f, gzip.GzipFile(fileobj=f) as g:
        return g.read().split("\n"), lang

    if split_path[-1] == "txt":
      # CWMT
      lang = split_path[-2].split("_")[-1]
      lang = "zh" if lang in ("ch", "cn") else lang
    else:
      lang = split_path[-1]
    with tf.io.gfile.GFile(path) as f:
      return f.read().split("\n"), lang

  def _parse_sgm(path):
    """Returns sentences from a single SGML file."""
    lang = path.split(".")[-2]
    sentences = []
    # Note: We can't use the XML parser since some of the files are badly
    # formatted.
    seg_re = re.compile(r"<seg id=\"\d+\">(.*)</seg>")
    with tf.io.gfile.GFile(path) as f:
      for line in f:
        seg_match = re.match(seg_re, line)
        if seg_match:
          assert len(seg_match.groups()) == 1
          sentences.append(seg_match.groups()[0])
    return sentences, lang

  parse_file = _parse_sgm if f1.endswith(".sgm") else _parse_text

  # Some datasets (e.g., CWMT) contain multiple parallel files specified with
  # a wildcard. We sort both sets to align them and parse them one by one.
  f1_files = tf.io.gfile.glob(f1)
  f2_files = tf.io.gfile.glob(f2)

  assert f1_files and f2_files, "No matching files found: %s, %s." % (f1, f2)
  assert len(f1_files) == len(f2_files), (
      "Number of files do not match: %d vs %d for %s vs %s." % (
          len(f1_files), len(f2_files), f1, f2))

  for f1_i, f2_i in zip(sorted(f1_files), sorted(f2_files)):
    l1_sentences, l1 = parse_file(f1_i)
    l2_sentences, l2 = parse_file(f2_i)

    assert len(l1_sentences) == len(l2_sentences), (
        "Sizes do not match: %d vs %d for %s vs %s." % (
            len(l1_sentences), len(l2_sentences), f1_i, f2_i))

    for s1, s2 in zip(l1_sentences, l2_sentences):
      yield {
          l1: s1,
          l2: s2
      }
1789
def run(self, args, pipeline_command):
        

        # output that must be moved but not renamed
        consistentNaming = ['alignments/normal_dna_fix_pg_sorted.bam',
                           'alignments/normal_dna_fix_pg_sorted.bam.bai',
                           'alignments/rna_genome_sorted.bam',
                           'alignments/rna_genome_sorted.bam.bai',
                           'alignments/rna_transcriptome.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam.bai',
                           'mutations/merged/all_merged.vcf',
                           'rankboost/mhcii_rankboost_concise_results.tsv',
                           'rankboost/mhci_rankboost_concise_results.tsv',
                           ]

        # output that must be renamed as well as moved
        # map of the original name to the final name
        renamingNeeded = {'binding_predictions': 'binding_predictions.tar',
                         'expression': 'expression.tar',
                         'haplotyping': 'haplotyping.tar',
                         'peptides': 'peptides.tar',
                         'rankboost': 'rankboost.tar',
                         'reports': 'reports.tar',
                         'mutations/snpeffed/mutations.vcf': 'all_snpeffed.vcf',
                         'mutations/transgened/mutations.vcf': 'all_transgened.vcf',
                         'mutations/merged': 'merged_perchrom.tar',
                         'mutations/muse': 'muse_perchrom.tar',
                         'mutations/mutect': 'mutect_perchrom.tar',
                         'mutations/radia': 'radia_perchrom.tar',
                         'mutations/somaticsniper': 'somaticsniper_perchrom.tar',
                         'mutations/strelka/snv': 'strelka_snv_perchrom.tar',
                         'mutations/strelka/indel': 'strelka_indel_perchrom.tar'}

        def make_output(output_dir, source_dir):
            """
            :param output_dir: dir to write the output to
            :param source_dir: dir containing the directory structure to be parsed
            :return:
            """
            def make_tar(dir, tar):
                with tarfile.open(tar, "w:gz") as tar:
                    tar.add(dir)
            # the output dir is where the real output directories are written
            protect_outputs = os.listdir(source_dir)
            for protectOut in protect_outputs:
                def getName(fileName):
                    return os.path.join(os.path.join(source_dir, protectOut), fileName)
                # move individual files out
                for fileName in consistentNaming:
                    shutil.copyfile(getName(fileName), os.path.join(output_dir, os.path.basename(fileName)))
                for src, dst in renamingNeeded.iteritems():
                    if dst.endswith('.tar'):
                        make_tar(getName(src), os.path.join(output_dir, dst))
                    else:
                        shutil.copyfile(getName(src), os.path.join(output_dir, dst))
            shutil.rmtree(source_dir)


        # prepare workdir
        mount = self._prepare_mount(args)
        self._workdir = os.path.join(mount, 'Toil-' + self._name)
        # insure the pairs are in the same directory, as protect expects
        # This is made more complicated by the fact CWLTool mounts inputs into random, read-only dirs
        # to get around this we copy all inputs into their own directories that we own
        tumor_dna_dir = os.path.expanduser('~/tumorDNA')
        tumor_rna_dir = os.path.expanduser('~/tumorRNA')
        normal_dna_dir = os.path.expanduser('~/normalDNA')
        os.mkdir(tumor_dna_dir)
        os.mkdir(tumor_rna_dir)
        os.mkdir(normal_dna_dir)
        shutil.copy(args.tumor_dna, tumor_dna_dir)
        shutil.copy(args.tumor_rna, tumor_rna_dir)
        shutil.copy(args.normal_dna, normal_dna_dir)
        shutil.copy(args.tumor_dna2, tumor_dna_dir)
        shutil.copy(args.tumor_rna2, tumor_rna_dir)
        shutil.copy(args.normal_dna2, normal_dna_dir)
        args.tumor_dna = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna))
        args.tumor_dna2 = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna2))
        args.tumor_rna = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna))
        args.tumor_rna2 = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna2))
        args.normal_dna = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna))
        args.normal_dna2 = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna2))

        # prepare config
        args_dict = vars(args)
        args_dict['output_dir'] = mount
        self._config = textwrap.dedent(self._config.format(**args_dict))
        self._sample_name = args_dict["sample_name"]
        config_path = os.path.join(self._workdir, 'config')
        command = self._make_prefix(os.path.join(self._workdir, 'jobStore'),
                                          config_path,
                                          self._workdir) + pipeline_command
        if self._resume and args.resume:
            command.append('--restart')
        self._create_workdir(args)
        with open(config_path, 'w') as f:
            f.write(self._config)

        try:
            subprocess.check_call(command)
        except subprocess.CalledProcessError as e:
            print(e, file=sys.stderr)
        finally:
            log.info('Pipeline terminated, changing ownership of output files from root to user.')
            stat = os.stat(self._mount)
            subprocess.check_call(['chown', '-R', '{}:{}'.format(stat.st_uid, stat.st_gid),
                                   self._mount])
            make_output(self._mount, os.path.join(self._mount, 'output'))
            if self._no_clean and args.no_clean:
                log.info('Flag "--no-clean" was used, therefore %s was not deleted.', self._workdir)
            else:
                log.info('Cleaning up temporary directory: %s', self._workdir)
                shutil.rmtree(self._workdir)
1047823
def get_linc_rna_genes(
        path_or_buffer,
        remove_duplicates=True,
        **kwargs):
    

    """
    valid_biotypes = set(['lincRNA'])
    
    df = get_genes(path_or_buffer, valid_biotypes,
                   remove_duplicates=remove_duplicates, **kwargs)
    return df
780980
def iTable2GFF(iTable, fa, contig = False):
    
    columns = ['#seqname', 'source', 'feature', 'start', 'end', 'score', 'strand', 'frame', 'attribute']
    gff = {c:[] for c in columns}
    for insertion in iTable.iterrows():
        insertion = insertion[1]
        if insertion['ID'] not in fa:
            continue
        # rRNA strand
        strand = insertion['sequence'].split('strand=', 1)[1].split()[0]
        # set rRNA positions for reporting features on contig or extracted sequence
        if contig is True:
            gene = [int(i) for i in insertion['sequence'].split('pos=', 1)[1].split()[0].split('-')]
            if strand == '-':
                offset = -1 * (gene[1])
            else:
                offset = gene[0]
        else:
            strand = '+'
            gene = [1, int(insertion['sequence'].split('total-len=', 1)[1].split()[0])]
            offset = gene[0]
        insertion['strand'] = strand
        insertion['offset'] = offset
        # source for prediction
        source = insertion['sequence'].split('::model', 1)[0].rsplit(' ', 1)[-1]
        insertion['source'] = source
        # rRNA gene
        geneAnnot = '%s rRNA gene' % (source.split('from', 1)[0])
        geneNum = insertion['sequence'].split('seq=', 1)[1].split()[0]
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append(source)
        gff['feature'].append('Gene')
        gff['start'].append(gene[0])
        gff['end'].append(gene[1])
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('ID=%s; Name=%s' % (geneNum, geneAnnot))
        # rRNA
        gff = parse_rRNA(insertion, fa[insertion['ID']], gff)
        # insertions
        gff = parse_insertion(insertion, gff)
        # orfs
        gff = parse_orf(insertion, gff)
        # catalytic RNAs
        gff = parse_catalytic(insertion, gff)
    return pd.DataFrame(gff)[columns].drop_duplicates()
341093
def from_command_line():
    
    # Parse options
    parser = argparse.ArgumentParser(
        description='Convert Complete Genomics var files to gVCF format.')
    parser.add_argument(
        '-d', '--refseqdir', metavar='REFSEQDIR', required=True,
        dest='refseqdir',
        help='Directory twobit reference genomes files are stored.')
    parser.add_argument(
        '-i', '--input', metavar='INPUTVARFILE',
        dest='cgivarfile',
        help='Path to Complete Genomics var file to convert. If omitted, data '
        ' also be piped in as standard input.')
    parser.add_argument(
        '-o', '--output', metavar='OUTPUTVCFFILE',
        dest='vcfoutfile',
        help='Path to where to save output VCF file.')
    parser.add_argument(
        '-D', '--download', action='store_true', dest='downloadrefseq',
        help='Download the 2bit file from UCSC to REFSEQDIR, if needed.')
    parser.add_argument(
        '-v', '--var-only', action='store_true', dest='varonly',
        help='Only report variant lines (i.e. VCF, but not gVCF)')
    args = parser.parse_args()

    # Get local twobit file from its directory. Download and store if needed.
    twobit_path, twobit_name = get_reference_genome_file(
        args.refseqdir, build='b37')
    # Handle input
    if sys.stdin.isatty():  # false if data is piped in
        var_input = args.cgivarfile
    else:
        var_input = sys.stdin
    # Handle output
    if args.vcfoutfile:
        convert_to_file(var_input,
                        args.vcfoutfile,
                        twobit_path,
                        twobit_name,
                        args.varonly)
    else:
        for line in convert(
                cgi_input=var_input,
                twobit_ref=twobit_path,
                twobit_name=twobit_name,
                var_only=args.varonly):
            print(line)
1064483
function parseCodeChunkToAST(src, options) {
    src = src
    .replace(/@/g, '\uFFDA')
    .replace(/#/g, '\uFFDB')
    ;
    var ast = recast.parse(src);
    return ast;
}
817591
function parse_Ref8U(blob) {
	var rwFirst = blob.read_shift(2);
	var rwLast = blob.read_shift(2);
	var colFirst = blob.read_shift(2);
	var colLast = blob.read_shift(2);
	return {s:{c:colFirst, r:rwFirst}, e:{c:colLast,r:rwLast}};
}
896155
def parse(self, fo):
        
        motifs = []

        pwm = []
        info = {}
        for line in fo.readlines():
            if line.startswith("#"):
                vals =  line.strip()[1:].split(" = ")
                if len(vals) > 1:
                    info[vals[0]] = vals[1]
            elif len(line) > 1:
                pwm.append([float(x) for x in line.strip().split("\t")])
            else:
                motifs.append(Motif())
                motifs[-1].consensus = info["Consensus"]
                motifs[-1].width = info["W"]
                motifs[-1].id = info["ID"]
                motifs[-1].pwm = pwm[:]
                pwm = []
            
        return motifs
985645
def tigrprepare(args):
    
    p = OptionParser(tigrprepare.__doc__)
    opts, args = p.parse_args(args)

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

    fastafile, asmbl_id, db, pasa_db = args
    if asmbl_id == 'all':
        idsfile = fastafile + ".ids"
        if need_update(fastafile, idsfile):
            ids([fastafile, "-o", idsfile])
    else:
        idsfile = asmbl_id

    oneid = open(idsfile).next().strip()

    weightsfile = "weights.txt"
    if need_update(idsfile, weightsfile):
        cmd = "$EVM/TIGR-only/create_sample_weights_file.dbi"
        cmd += " {0} {1} | tee weights.txt".format(db, oneid)
        sh(cmd)

    evs = ["gene_predictions.gff3", "transcript_alignments.gff3",
           "protein_alignments.gff3"]
    if need_update(weightsfile, evs):
        cmd = "$EVM/TIGR-only/write_GFF3_files.dbi"
        cmd += " --db {0} --asmbl_id {1} --weights {2}".\
                format(db, idsfile, weightsfile)
        sh(cmd)

    evs[1] = fix_transcript()

    partition(evs)
    runfile = "run.sh"
    contents = EVMRUN.format(*evs)
    write_file(runfile, contents)
1037633
def parse_orf(insertion, gff):
    
    offset = insertion['offset']
    if type(insertion['orf']) is not str:
        return gff
    for orf in parse_fasta(insertion['orf'].split('|')):
        ID = orf[0].split('>')[1].split()[0]
        Start, End, strand = [int(i) for i in orf[0].split(' # ')[1:4]]
        if strand == 1:
            strand = '+'
        else:
            strand = '-'
        GeneStrand = insertion['strand']
        if strand != GeneStrand:
            if strand == '+':
                strand = '-'
            else:
                strand = '+'
            Start, End = End - 2, Start - 2
        Start, End = abs(Start + offset) - 1, abs(End + offset) - 1
        annot = orf[0].split()[1]
        if annot == 'n/a':
            annot = 'unknown'
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append('Prodigal and Pfam')
        gff['feature'].append('CDS')
        gff['start'].append(Start)
        gff['end'].append(End)
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('ID=%s; Name=%s' % (ID, annot))
    return gff
676387
def g_to_n(self, g_interval):
        """"""

        grs, gre = g_interval.start.base - 1 - self.gc_offset, g_interval.end.base - 1 - self.gc_offset
        # frs, fre = (f)orward (r)na (s)tart & (e)nd; forward w.r.t. genome
        frs, frs_offset, frs_cigar = self._map(
            from_pos=self.ref_pos, to_pos=self.tgt_pos, pos=grs, base="start")
        fre, fre_offset, fre_cigar = self._map(
            from_pos=self.ref_pos, to_pos=self.tgt_pos, pos=gre, base="end")

        if self.strand == -1:
            frs, fre = self.tgt_len - fre - 1, self.tgt_len - frs - 1
            frs_offset, fre_offset = -fre_offset, -frs_offset

        # The returned interval would be uncertain when locating at alignment gaps
        return hgvs.location.BaseOffsetInterval(
            start=hgvs.location.BaseOffsetPosition(
                base=frs + 1, offset=frs_offset, datum=Datum.SEQ_START),
            end=hgvs.location.BaseOffsetPosition(
                base=fre + 1, offset=fre_offset, datum=Datum.SEQ_START),
            uncertain=frs_cigar in 'DI' or fre_cigar in 'DI')
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)
589841
def _parse_epo(cls, fname):
        """"""

        data = {}
        with open(fname) as fd:
            for el in (cls._strfactory(_) for _ in fd):
                if el:
                    data.setdefault(el.gabid, []).append( el )
        log.info("parsed %d elements from %s" % (len(data), fname))
        return data
303752
def fromgtf(args):
    
    p = OptionParser(fromgtf.__doc__)
    p.add_option("--transcript_id", default="transcript_id",
                 help="Field name for transcript [default: %default]")
    p.add_option("--gene_id", default="gene_id",
                 help="Field name for gene [default: %default]")
    p.add_option("--augustus", default=False, action="store_true",
                 help="Input is AUGUSTUS gtf [default: %default]")
    p.set_home("augustus")
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    gtffile, = args
    outfile = opts.outfile
    if opts.augustus:
        ahome = opts.augustus_home
        s = op.join(ahome, "scripts/gtf2gff.pl")
        cmd = "{0} --gff3 < {1} --out={2}".format(s, gtffile, outfile)
        sh(cmd)
        return

    gff = Gff(gtffile)
    fw = must_open(outfile, "w")
    transcript_id = opts.transcript_id
    gene_id = opts.gene_id
    nfeats = 0
    for g in gff:
        if g.type in ("transcript", "mRNA"):
            g.type = "mRNA"
            g.update_tag(transcript_id, "ID")
            g.update_tag("mRNA", "ID")
            g.update_tag(gene_id, "Parent")
            g.update_tag("Gene", "Parent")
        elif g.type in ("exon", "CDS") or "UTR" in g.type:
            g.update_tag("transcript_id", "Parent")
            g.update_tag(g.type, "Parent")
        elif g.type == "gene":
            g.update_tag(gene_id, "ID")
            g.update_tag("Gene", "ID")
        else:
            assert 0, "Don't know how to deal with {0}".format(g.type)

        g.update_attributes()
        print(g, file=fw)
        nfeats += 1

    logging.debug("A total of {0} features written.".format(nfeats))
3718
def gff(args):
    
    align_choices = ("left", "center", "right")
    p = OptionParser(gff.__doc__)
    p.add_option("--align", default="left", choices=align_choices,
                 help="Horizontal alignment [default: %default]")
    p.add_option("--noUTR", default=False, action="store_true",
                 help="Do not plot UTRs [default: %default]")
    opts, args = p.parse_args(args)

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

    fig = plt.figure(1, (8, 5))
    root = fig.add_axes([0, 0, 1, 1])

    gffiles = args
    ngenes = len(gffiles)

    canvas = .6
    setups, ratio = get_setups(gffiles, canvas=canvas, noUTR=opts.noUTR)
    align = opts.align
    xs = .2 if align == "left" else .8
    yinterval = canvas / ngenes
    ys = .8
    tip = .01
    for genename, mrnabed, cdsbeds in setups:
        ExonGlyph(root, xs, ys, mrnabed, cdsbeds, ratio=ratio, align=align)
        if align == "left":
            root.text(xs - tip, ys, genename, ha="right", va="center")
        elif align == "right":
            root.text(xs + tip, ys, genename, ha="left", va="center")
        ys -= yinterval

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    figname = "exons.pdf"
    savefig(figname, dpi=300)