%prog gff *.gff

    Draw exons for genes based on gff files. Each gff file should contain only
    one gene, and only the "mRNA" and "CDS" feature will be drawn on the canvas.
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)
1291224
def fixpartials(args):
    
    p = OptionParser(fixpartials.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    gffile, gfasta, partials, = args

    gff = make_index(gffile)
    genome = Fasta(gfasta, index=True)
    partials = LineFile(partials, load=True).lines

    # all_transcripts = [f.id for f in gff.features_of_type("mRNA", \
    #     order_by=("seqid", "start"))]
    seen = set()
    fw = must_open(opts.outfile, "w")
    for gene in gff.features_of_type('gene',  order_by=('seqid', 'start')):
        children = AutoVivification()
        cflag = False
        transcripts = list(gff.children(gene, level=1, order_by=('start')))
        for transcript in transcripts:
            trid, seqid, strand = transcript.id, transcript.seqid, transcript.strand

            for child in gff.children(transcript, order_by=('start')):
                ftype = child.featuretype
                if ftype not in children[trid]: children[trid][ftype] = []
                children[trid][ftype].append(child)

            five_prime, three_prime = True, True
            nstart, nstop = (None, None), (None, None)
            cds_span = [children[trid]['CDS'][0].start, \
                children[trid]['CDS'][-1].stop]
            new_cds_span = [x for x in cds_span]

            start_codon = (cds_span[0], cds_span[0] + 2)
            stop_codon = (cds_span[1] - 2, cds_span[1])
            if strand == '-': start_codon, stop_codon = stop_codon, start_codon

            if trid in partials:
                seen.add(trid)
                start_codon_fasta = _fasta_slice(genome, seqid, \
                    start_codon[0], start_codon[1], strand)
                stop_codon_fasta = _fasta_slice(genome, seqid, \
                    stop_codon[0], stop_codon[1], strand)

                if not is_valid_codon(start_codon_fasta, type='start'):
                    five_prime = False
                    nstart = scan_for_valid_codon(start_codon, strand, \
                        seqid, genome, type='start')

                if not is_valid_codon(stop_codon_fasta, type='stop'):
                    three_prime = False
                    nstop = scan_for_valid_codon(stop_codon, strand, \
                        seqid, genome, type='stop')

                logging.debug("feature={0} ({1})".format(trid, strand) + \
                ", 5'={0}, 3'={1}".format(five_prime, three_prime) + \
                ", {0} <== {1} ==> {2}".format(nstart if strand == '+' else nstop, \
                    cds_span, nstop if strand == '+' else nstart))

            if not five_prime or not three_prime:
                if nstart != (None, None) and (start_codon != nstart):
                    i = 0 if strand == '+' else 1
                    new_cds_span[i] = nstart[i]
                if nstop != (None, None) and (stop_codon != nstop):
                    i = 1 if strand == '+' else 0
                    new_cds_span[i] = nstop[i]
                new_cds_span.sort()

            if set(cds_span) != set(new_cds_span):
                cflag = True
                # if CDS has been extended, appropriately adjust all relevent
                # child feature (CDS, exon, UTR) coordinates
                for ftype in children[trid]:
                    for idx in range(len(children[trid][ftype])):
                        child_span = (children[trid][ftype][idx].start, \
                            children[trid][ftype][idx].stop)
                        if ftype in ('exon', 'CDS'):
                            # if exons/CDSs, adjust start and stop according to
                            # new CDS start and stop, respectively
                            if child_span[0] == cds_span[0]:
                                children[trid][ftype][idx].start = new_cds_span[0]
                            if child_span[1] == cds_span[1]:
                                children[trid][ftype][idx].stop = new_cds_span[1]
                        elif ftype.endswith('UTR'):
                            # if *_prime_UTR, adjust stop according to new CDS start and
                            #                 adjust start according to new CDS stop
                            if child_span[1] == cds_span[0]:
                                children[trid][ftype][idx].stop = new_cds_span[0]
                            if child_span[0] == cds_span[1]:
                                children[trid][ftype][idx].start = new_cds_span[1]

                transcript.start, transcript.stop = \
                    children[trid]['exon'][0].start, children[trid]['exon'][-1].stop

        if cflag:
            _gene_span = range_minmax([(tr.start, tr.stop) for tr in transcripts])
            gene.start, gene.stop = _gene_span[0], _gene_span[1]

        # print gff file
        print(gene, file=fw)
        for transcript in transcripts:
            trid = transcript.id
            print(transcript, file=fw)
            for cftype in children[trid]:
                for child in children[trid][cftype]:
                    print(child, file=fw)

    fw.close()
753112
def main(argv=None):
    """"""
    if argv is None:
        argv = sys.argv

    log.setLevel(logging.INFO)

    import optparse  # TODO: use argparse

    parser = optparse.OptionParser(
        usage='usage: %prog [options] path',
        description='Display image data in TIFF files.',
        version='%%prog %s' % __version__, prog='tifffile')
    opt = parser.add_option
    opt('-p', '--page', dest='page', type='int', default=-1,
        help='display single page')
    opt('-s', '--series', dest='series', type='int', default=-1,
        help='display series of pages of same shape')
    opt('--nomultifile', dest='nomultifile', action='store_true',
        default=False, help='do not read OME series from multiple files')
    opt('--noplots', dest='noplots', type='int', default=10,
        help='maximum number of plots')
    opt('--interpol', dest='interpol', metavar='INTERPOL', default=None,
        help='image interpolation method')
    opt('--dpi', dest='dpi', type='int', default=96,
        help='plot resolution')
    opt('--vmin', dest='vmin', type='int', default=None,
        help='minimum value for colormapping')
    opt('--vmax', dest='vmax', type='int', default=None,
        help='maximum value for colormapping')
    opt('--debug', dest='debug', action='store_true', default=False,
        help='raise exception on failures')
    opt('--doctest', dest='doctest', action='store_true', default=False,
        help='runs the docstring examples')
    opt('-v', '--detail', dest='detail', type='int', default=2)
    opt('-q', '--quiet', dest='quiet', action='store_true')

    settings, path = parser.parse_args()
    path = ' '.join(path)

    if settings.doctest:
        import doctest
        if sys.version_info < (3, 6):
            print('Doctests work with Python >=3.6 only')
            return 0
        doctest.testmod(optionflags=doctest.ELLIPSIS)
        return 0
    if not path:
        path = askopenfilename(title='Select a TIFF file',
                               filetypes=TIFF.FILEOPEN_FILTER)
        if not path:
            parser.error('No file specified')

    if any(i in path for i in '?*'):
        path = glob.glob(path)
        if not path:
            print('No files match the pattern')
            return 0
        # TODO: handle image sequences
        path = path[0]

    if not settings.quiet:
        print_('\nReading TIFF header:', end=' ', flush=True)
    start = time.time()
    try:
        tif = TiffFile(path, multifile=not settings.nomultifile)
    except Exception as exc:
        if settings.debug:
            raise
        print('\n\n%s: %s' % (exc.__class__.__name__, exc))
        sys.exit(0)

    if not settings.quiet:
        print('%.3f ms' % ((time.time()-start) * 1e3))

    if tif.is_ome:
        settings.norgb = True

    images = []
    if settings.noplots > 0:
        if not settings.quiet:
            print_('Reading image data: ', end=' ', flush=True)

        def notnone(x):
            return next(i for i in x if i is not None)

        start = time.time()
        try:
            if settings.page >= 0:
                images = [(tif.asarray(key=settings.page),
                           tif[settings.page], None)]
            elif settings.series >= 0:
                images = [(tif.asarray(series=settings.series),
                           notnone(tif.series[settings.series]._pages),
                           tif.series[settings.series])]
            else:
                for i, s in enumerate(tif.series[:settings.noplots]):
                    try:
                        images.append((tif.asarray(series=i),
                                       notnone(s._pages),
                                       tif.series[i]))
                    except Exception as exc:
                        images.append((None, notnone(s.pages), None))
                        if settings.debug:
                            raise
                        print('\nSeries %i failed with %s: %s... ' %
                              (i, exc.__class__.__name__, exc), end='')
        except Exception as exc:
            if settings.debug:
                raise
            print('%s: %s' % (exc.__class__.__name__, exc))

        if not settings.quiet:
            print('%.3f ms' % ((time.time()-start) * 1e3))

    if not settings.quiet:
        print_('Generating printout:', end=' ', flush=True)
        start = time.time()
        info = TiffFile.__str__(tif, detail=int(settings.detail))
        print('%.3f ms' % ((time.time()-start) * 1e3))
        print()
        print(info)
        print()
    tif.close()

    if images and settings.noplots > 0:
        try:
            import matplotlib
            matplotlib.use('TkAgg')
            from matplotlib import pyplot
        except ImportError as exc:
            log.warning('tifffile.main: %s: %s', exc.__class__.__name__, exc)
        else:
            for img, page, series in images:
                if img is None:
                    continue
                vmin, vmax = settings.vmin, settings.vmax
                if 'GDAL_NODATA' in page.tags:
                    try:
                        vmin = numpy.min(
                            img[img > float(page.tags['GDAL_NODATA'].value)])
                    except ValueError:
                        pass
                if tif.is_stk:
                    try:
                        vmin = tif.stk_metadata['MinScale']
                        vmax = tif.stk_metadata['MaxScale']
                    except KeyError:
                        pass
                    else:
                        if vmax <= vmin:
                            vmin, vmax = settings.vmin, settings.vmax
                if series:
                    title = '%s\n%s\n%s' % (str(tif), str(page), str(series))
                else:
                    title = '%s\n %s' % (str(tif), str(page))
                photometric = 'MINISBLACK'
                if page.photometric not in (3,):
                    photometric = TIFF.PHOTOMETRIC(page.photometric).name
                imshow(img, title=title, vmin=vmin, vmax=vmax,
                       bitspersample=page.bitspersample,
                       photometric=photometric,
                       interpolation=settings.interpol,
                       dpi=settings.dpi)
            pyplot.show()
    return 0
347685
def canonical_transcripts(gtf, out_file):
    
    if file_exists(out_file):
        return out_file
    db = get_gtf_db(gtf)
    with file_transaction(out_file) as tx_out_file:
        with open(tx_out_file, "w") as out_handle:
            for gene in db.features_of_type('gene'):
                exon_list = []
                for ti, transcript in enumerate(db.children(gene, level=1)):
                    cds_len = 0
                    total_len = 0
                    exons = list(db.children(transcript, level=1))
                    for exon in exons:
                        exon_length = len(exon)
                        if exon.featuretype == 'CDS':
                            cds_len += exon_length
                        total_len += exon_length

                    exon_list.append((cds_len, total_len, transcript, exons))

                # If we have CDS, then use the longest coding transcript
                if max(i[0] for i in exon_list) > 0:
                    best = sorted(exon_list)[0]
                # Otherwise, just choose the longest
                else:
                    best = sorted(exon_list, key=lambda x: x[1])[0]
                for exon in db.children(best[2], level=1):
                    out_handle.write(str(exon) + "\n")
    return out_file
163942
def parents(args):
    
    p = OptionParser(parents.__doc__)

    opts, args = p.parse_args(args)

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

    gff_file, idsfile = args
    g = make_index(gff_file)
    fp = open(idsfile)
    for row in fp:
        cid = row.strip()
        b = next(g.parents(cid, 1))
        print("\t".join((cid, b.id)))
1099835
def run_gvcfgenotyper(data, orig_region, vrn_files, out_file):
    
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            regions = _find_gvcf_blocks(vrn_files[0], bamprep.region_to_gatk(orig_region),
                                        os.path.dirname(tx_out_file))
            if len(regions) == 1:
                _run_gvcfgenotyper(data, regions[0], vrn_files, tx_out_file)
            else:
                split_outs = [_run_gvcfgenotyper(data, r, vrn_files,
                                                 "%s-%s.vcf.gz" % (utils.splitext_plus(out_file)[0],
                                                                   r.replace(":", "_").replace("-", "_")))
                              for r in regions]
                vcfutils.concat_variant_files(split_outs, tx_out_file, regions,
                                              dd.get_ref_file(data), data["config"])
    return vcfutils.bgzip_and_index(out_file, data["config"])
298556
def CDSReader( fh, format='gff' ):
    

    known_formats = ( 'gff', 'gtf', 'bed')
    if format not in known_formats: 
        print('%s format not in %s' % (format, ",".join( known_formats )), file=sys.stderr)
        raise Exception('?')

    if format == 'bed':
        for line in fh:    
            f = line.strip().split()
            chrom = f[0]
            chrom_start = int(f[1])
            name = f[4]
            strand = f[5]
            cdsStart = int(f[6])
            cdsEnd = int(f[7])
            blockCount = int(f[9])
            blockSizes = [ int(i) for i in f[10].strip(',').split(',') ]
            blockStarts = [ chrom_start + int(i) for i in f[11].strip(',').split(',') ]

            # grab cdsStart - cdsEnd
            cds_exons = []
            cds_seq = ''
            genome_seq_index = []
            for base,offset in zip( blockStarts, blockSizes ):
                if (base + offset) < cdsStart: continue
                if base > cdsEnd: continue
                exon_start = max( base, cdsStart )
                exon_end = min( base+offset, cdsEnd ) 
                cds_exons.append( (exon_start, exon_end) )
            

    genelist = {}
    grouplist = []
    if format == 'gff' or format == 'gtf':
        for line in fh:
            if line.startswith('#'): continue
            fields = line.strip().split('\t')
            if len( fields ) < 9: continue
            if fields[2] not in ('CDS', 'stop_codon', 'start_codon'): continue

            # fields

            chrom = fields[0]
            ex_st = int( fields[3] ) - 1 # make zero-centered
            ex_end = int( fields[4] ) #+ 1 # make exclusive
            strand = fields[6]

            if format == 'gtf':
                group = fields[8].split(';')[0]
            else:
                group = fields[8]

            if group not in grouplist: grouplist.append( group )
            if group not in genelist:
                genelist[group] = (chrom, strand, [])
            
            genelist[group][2].append( ( ex_st, ex_end ) )

        sp = lambda a,b: cmp( a[0], b[0] )

        #for gene in genelist.values():
        for gene in grouplist:
            chrom, strand, cds_exons = genelist[ gene ]
            seqlen = sum([ a[1]-a[0] for a in cds_exons ])
            overhang = seqlen % 3
            if overhang > 0:
                #print >>sys.stderr, "adjusting ", gene  
                if strand == '+': 
                    cds_exons[-1] = ( cds_exons[-1][0], cds_exons[-1][1] - overhang )
                else:
                    cds_exons[0] = ( cds_exons[0][0] + overhang, cds_exons[0][1] )
            cds_exons = bitset_union( cds_exons )
            yield chrom, strand, cds_exons, gene
966991
def main(args):

  
  sys.stderr.write("reading reference genome\n")
  ref = FASTAData(open(args.genome).read())
  sys.stderr.write("read in "+str(len(ref.keys()))+" chromosomes\n")

  """ Next make the transcriptome """
  txome = {}
  sys.stderr.write("write the transcriptome\n")
  inf = None
  if args.gpd[-3:] == '.gz':
    inf = gzip.open(args.gpd)
  else:
    inf = open(args.gpd)
  stream = GPDStream(inf)
  tof = open(args.tempdir+'/transcriptome.fa','w')
  z = 0
  for gpd in stream:
    z += 1
    if gpd.transcript_name in txome:
      sys.stderr.write("WARNING already have a transcript "+gpd.transcript_name+" ignoring line "+str(z)+" of the gpd\n")
      continue
    txome[gpd.transcript_name] = gpd.gene_name
    tof.write('>'+gpd.transcript_name+"\n"+str(gpd.get_sequence(ref))+"\n")
  tof.close()
  inf.close()
  sys.stderr.write("wrote "+str(len(txome.keys()))+" transcripts\n")

  """Build the salmon index"""
  sys.stderr.write("building a salmon index\n")
  cmd = 'salmon index -p '+str(args.numThreads)+' -t '+args.tempdir+'/transcriptome.fa -i '+args.tempdir+'/salmon_index'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished building the index\n")

  """Use the index to quanitfy"""
  sys.stderr.write("quanitfy reads\n")
  reads = ''
  if args.rU:
    reads = '-r '+args.rU
  else:
    reads = '-1 '+args.r1+' -2 '+args.r2
  cmd = 'salmon quant -p '+str(args.numThreads)+' -i '+args.tempdir+'/salmon_index -l A '+reads+' -o '+args.tempdir+'/output_quant'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished quanitfying\n")

  """Now parse the salmon output to add gene name"""
  salmon = {}
  with open(args.tempdir+'/output_quant/quant.sf') as inf:
    header = inf.readline()
    for line in inf:
      f = line.rstrip().split("\t")
      # by each transcript name hold a data strcture of the information
      salmon[f[0]] = {'name':f[0],'length':int(f[1]),'EffectiveLength':float(f[2]),'TPM':float(f[3]),'NumReads':float(f[4])}
  genes = {}
  for name in salmon:
    gene = txome[name]
    if gene not in genes: genes[gene] = []
    genes[gene].append(salmon[name])
  genetot = {}
  for gene in genes:
    tot = sum([x['TPM'] for x in genes[gene]])
    genetot[gene] = tot
  ordered_gene_names = sorted(genetot.keys(), key=lambda x: genetot[x],reverse=True)
      
  """Collected enough information to make output"""
  sys.stderr.write("generating output\n")
  of = sys.stdout
  if args.output != '-':
    of = open(args.output,'w')
  of.write("geneName\ttranscriptName\tlength\tEffectiveLength\ttxTPM\tNumReads\tgeneTPM\n")
  for gene in ordered_gene_names:
    txs = sorted(genes[gene],key=lambda x: x['TPM'],reverse=True)
    for tx in txs:
      of.write(gene+"\t"+tx['name']+"\t"+str(tx['length'])+"\t"+str(tx['EffectiveLength'])+"\t"+str(tx['TPM'])+"\t"+str(tx['NumReads'])+"\t"+str(genetot[gene])+"\n")
  of.close()

  # Temporary working directory step 3 of 3 - Cleanup
  if not args.specific_tempdir:
    rmtree(args.tempdir)
1004054
def combine_files(samples):
    
    data = samples[0][0]
    # prefer the supplied transcriptome gtf file
    gtf_file = dd.get_transcriptome_gtf(data, None)
    if not gtf_file:
        gtf_file = dd.get_gtf_file(data, None)
    dexseq_gff = dd.get_dexseq_gff(data)

    # combine featureCount files
    count_files = filter_missing([dd.get_count_file(x[0]) for x in samples])
    combined = count.combine_count_files(count_files, ext=".counts")
    annotated = count.annotate_combined_count_file(combined, gtf_file)

    # add tx2gene file
    tx2gene_file = os.path.join(dd.get_work_dir(data), "annotation", "tx2gene.csv")
    if gtf_file:
        tx2gene_file = sailfish.create_combined_tx2gene(data)

    # combine eXpress files
    express_counts_combined = combine_express(samples, combined)

    # combine Cufflinks files
    fpkm_files = filter_missing([dd.get_fpkm(x[0]) for x in samples])
    if fpkm_files and combined:
        fpkm_combined_file = os.path.splitext(combined)[0] + ".fpkm"
        fpkm_combined = count.combine_count_files(fpkm_files, fpkm_combined_file)
    else:
        fpkm_combined = None
    isoform_files = filter_missing([dd.get_fpkm_isoform(x[0]) for x in samples])
    if isoform_files and combined:
        fpkm_isoform_combined_file = os.path.splitext(combined)[0] + ".isoform.fpkm"
        fpkm_isoform_combined = count.combine_count_files(isoform_files,
                                                          fpkm_isoform_combined_file,
                                                          ".isoform.fpkm")
    else:
        fpkm_isoform_combined = None
    # combine DEXseq files
    to_combine_dexseq = filter_missing([dd.get_dexseq_counts(data[0]) for data
                                        in samples])
    if to_combine_dexseq and combined:
        dexseq_combined_file = os.path.splitext(combined)[0] + ".dexseq"
        dexseq_combined = count.combine_count_files(to_combine_dexseq,
                                                    dexseq_combined_file, ".dexseq")
        if dexseq_combined:
            dexseq.create_dexseq_annotation(dexseq_gff, dexseq_combined)
    else:
        dexseq_combined = None
    samples = spikein.combine_spikein(samples)
    updated_samples = []
    for data in dd.sample_data_iterator(samples):
        if combined:
            data = dd.set_combined_counts(data, combined)
        if annotated:
            data = dd.set_annotated_combined_counts(data, annotated)
        if fpkm_combined:
            data = dd.set_combined_fpkm(data, fpkm_combined)
        if fpkm_isoform_combined:
            data = dd.set_combined_fpkm_isoform(data, fpkm_isoform_combined)
        if express_counts_combined:
            data = dd.set_express_counts(data, express_counts_combined['counts'])
            data = dd.set_express_tpm(data, express_counts_combined['tpm'])
            data = dd.set_express_fpkm(data, express_counts_combined['fpkm'])
            data = dd.set_isoform_to_gene(data, express_counts_combined['isoform_to_gene'])
        if dexseq_combined:
            data = dd.set_dexseq_counts(data, dexseq_combined_file)
        if gtf_file:
            data = dd.set_tx2gene(data, tx2gene_file)
        updated_samples.append([data])
    return updated_samples
1108010
def vcf_to_df_worker(arg):
    
    canvasvcf, exonbed, i = arg
    logging.debug("Working on job {}: {}".format(i, canvasvcf))
    samplekey = op.basename(canvasvcf).split(".")[0].rsplit('_', 1)[0]
    d = {'SampleKey': samplekey}

    exons = BedTool(exonbed)
    cn = parse_segments(canvasvcf)
    overlaps = exons.intersect(cn, wao=True)
    gcn_store = {}
    for ov in overlaps:
        # Example of ov.fields:
        # [u'chr1', u'11868', u'12227', u'ENSG00000223972.5',
        # u'ENST00000456328.2', u'transcribed_unprocessed_pseudogene',
        # u'DDX11L1', u'.', u'-1', u'-1', u'.', u'0']
        gene_name = "|".join((ov.fields[6], ov.fields[3], ov.fields[5]))
        if gene_name not in gcn_store:
            gcn_store[gene_name] = defaultdict(int)

        cn = ov.fields[-2]
        if cn == ".":
            continue
        cn = int(cn)
        if cn > 10:
            cn = 10
        amt = int(ov.fields[-1])
        gcn_store[gene_name][cn] += amt

    for k, v in sorted(gcn_store.items()):
        v_mean, v_median = counter_mean_and_median(v)
        d[k + ".avgcn"] = v_mean
        d[k + ".medcn"] = v_median
    cleanup()
    return d
90783
def guess_infer_extent(gtf_file):
    
    _, ext = os.path.splitext(gtf_file)
    tmp_out = tempfile.NamedTemporaryFile(suffix=".gtf", delete=False).name
    with open(tmp_out, "w") as out_handle:
        count = 0
        in_handle = utils.open_gzipsafe(gtf_file)
        for line in in_handle:
            if count > 1000:
                break
            out_handle.write(line)
            count += 1
        in_handle.close()
    db = gffutils.create_db(tmp_out, dbfn=":memory:", infer_gene_extent=False)
    os.remove(tmp_out)
    features = [x for x in db.featuretypes()]
    if "gene" in features and "transcript" in features:
        return False
    else:
        return True
282948
def run_vep(in_file, data):
    
    if not vcfutils.vcf_has_variants(in_file):
        return None
    out_file = utils.append_stem(in_file, "-vepeffects")
    assert in_file.endswith(".gz") and out_file.endswith(".gz")
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            vep_dir, ensembl_name = prep_vep_cache(data["genome_build"],
                                                   tz.get_in(["reference", "fasta", "base"], data))
            if vep_dir:
                cores = tz.get_in(("config", "algorithm", "num_cores"), data, 1)
                fork_args = ["--fork", str(cores)] if cores > 1 else []
                vep = config_utils.get_program("vep", data["config"])
                # HGVS requires a bgzip compressed, faidx indexed input file or is unusable slow
                if dd.get_ref_file_compressed(data):
                    hgvs_compatible = True
                    config_args = ["--fasta", dd.get_ref_file_compressed(data)]
                else:
                    hgvs_compatible = False
                    config_args = ["--fasta", dd.get_ref_file(data)]
                if vcfanno.is_human(data):
                    plugin_fns = {"loftee": _get_loftee, "maxentscan": _get_maxentscan,"genesplicer": _get_genesplicer,
                                  "spliceregion": _get_spliceregion, "G2P": _get_G2P}
                    plugins = ["loftee", "G2P"]
                    if "vep_splicesite_annotations" in dd.get_tools_on(data):
                        # "genesplicer" too unstable so currently removed
                        plugins += ["maxentscan", "spliceregion"]
                    for plugin in plugins:
                        plugin_args = plugin_fns[plugin](data)
                        config_args += plugin_args
                    config_args += ["--sift", "b", "--polyphen", "b"]
                    if hgvs_compatible:
                        config_args += ["--hgvsg","--hgvs", "--shift_hgvs", "1"]
                if (dd.get_effects_transcripts(data).startswith("canonical")
                      or tz.get_in(("config", "algorithm", "clinical_reporting"), data)):
                    config_args += ["--most_severe"]
                else:
                    config_args += ["--flag_pick_allele_gene"]
                if ensembl_name.endswith("_merged"):
                    config_args += ["--merged"]
                    ensembl_name = ensembl_name.replace("_merged", "")
                resources = config_utils.get_resources("vep", data["config"])
                extra_args = [str(x) for x in resources.get("options", [])]
                cmd = [vep, "--vcf", "-o", "stdout", "-i", in_file] + fork_args + extra_args + \
                      ["--species", ensembl_name,
                       "--no_stats", "--cache",
                        "--offline", "--dir", vep_dir,
                       "--symbol", "--numbers", "--biotype", "--total_length", "--canonical",
                       "--gene_phenotype", "--ccds", "--uniprot", "--domains", "--regulatory",
                       "--protein", "--tsl", "--appris", "--af", "--max_af", "--af_1kg", "--af_esp", "--af_gnomad",
                       "--pubmed", "--variant_class", "--allele_number"] + config_args
                perl_exports = utils.get_perl_exports()
                # Remove empty fields (';;') which can cause parsing errors downstream
                cmd = "%s && %s | sed '/^#/! s/;;/;/g' | bgzip -c > %s" % (perl_exports, " ".join(cmd), tx_out_file)
                do.run(cmd, "Ensembl variant effect predictor", data)
    if utils.file_exists(out_file):
        return vcfutils.bgzip_and_index(out_file, data["config"])
83532
def getquals(args):
    
    p = OptionParser(getquals.__doc__)
    p.add_option("--types", default="gene,mRNA,CDS",
                type="str", dest="quals_ftypes",
                help="Feature types from which to extract qualifiers")
    p.add_option("--ignore", default="locus_tag,product,codon_start,translation",
                type="str", dest="quals_ignore",
                help="Qualifiers to exclude from parsing")
    opts, args = p.parse_args(args)

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

    gbkfile, = args
    quals_ftypes = opts.quals_ftypes.split(",")
    quals_ignore = opts.quals_ignore.split(",")

    locus = dict()
    locus_tag = None
    for rec in SeqIO.parse(gbkfile, "gb"):
        for f in rec.features:
            if f.type in quals_ftypes:
                locus_tag = f.qualifiers[LT][0]
                if locus_tag not in locus:
                    locus[locus_tag] = dict()
                    for ftype in quals_ftypes:
                        if ftype not in locus[locus_tag]:
                            locus[locus_tag][ftype] = []
                        if ftype == "CDS":  # store the CDS protein_id
                            locus[locus_tag]["protein_id"] = []

                quals = []
                for qual in f.qualifiers:
                    if qual in quals_ignore:
                        continue
                    for qval in f.qualifiers[qual]:
                        quals.append((locus_tag, qual, qval))
                        if qual == "protein_id":
                            locus[locus_tag]["protein_id"].append(qval)
                if len(quals) > 0:
                    locus[locus_tag][f.type].append(quals)

    for locus_tag in locus:
        print_locus_quals(locus_tag, locus, quals_ftypes)
1089778
def filter(args):
    
    p = OptionParser(filter.__doc__)
    p.add_option("--type", default="mRNA",
             help="The feature to scan for the attributes [default: %default]")
    g1 = OptionGroup(p, "Filter by identity/coverage attribute values")
    g1.add_option("--id", default=95, type="float",
                 help="Minimum identity [default: %default]")
    g1.add_option("--coverage", default=90, type="float",
                 help="Minimum coverage [default: %default]")
    g1.add_option("--nocase", default=False, action="store_true",
                 help="Case insensitive lookup of attribute names [default: %default]")
    p.add_option_group(g1)
    g2 = OptionGroup(p, "Filter by child feature bp length")
    g2.add_option("--child_ftype", default=None, type="str",
                 help="Child featuretype to consider")
    g2.add_option("--child_bp", default=None, type="int",
                 help="Filter by total bp of children of chosen ftype")
    p.add_option_group(g2)
    p.set_outfile()

    opts, args = p.parse_args(args)
    otype, oid, ocov = opts.type, opts.id, opts.coverage
    cftype, clenbp = opts.child_ftype, opts.child_bp

    id_attr, cov_attr = "Identity", "Coverage"
    if opts.nocase:
        id_attr, cov_attr = id_attr.lower(), cov_attr.lower()

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

    gffile, = args

    gffdb = make_index(gffile)
    bad = set()
    ptype = None
    for g in gffdb.features_of_type(otype, order_by=('seqid', 'start')):
        if not ptype:
            parent = list(gffdb.parents(g))
            ptype = parent[0].featuretype \
                if len(parent) > 0 else otype
        if cftype and clenbp:
            if gffdb.children_bp(g, child_featuretype=cftype) < clenbp:
                bad.add(g.id)
        elif oid and ocov:
            identity = float(g.attributes[id_attr][0])
            coverage = float(g.attributes[cov_attr][0])
            if identity < oid or coverage < ocov:
                bad.add(g.id)

    logging.debug("{0} bad accns marked.".format(len(bad)))

    fw = must_open(opts.outfile, "w")
    for g in gffdb.features_of_type(ptype, order_by=('seqid', 'start')):
        if ptype != otype:
            feats = list(gffdb.children(g, featuretype=otype, order_by=('start')))
            ok_feats = [f for f in feats if f.id not in bad]
            if len(ok_feats) > 0:
                print(g, file=fw)
                for feat in ok_feats:
                    print(feat, file=fw)
                    for child in gffdb.children(feat, order_by=('start')):
                        print(child, file=fw)
        else:
            if g.id not in bad:
                print(g, file=fw)
                for child in gffdb.children(g, order_by=('start')):
                    print(child, file=fw)
    fw.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
885499
def demo(args):
    
    p = OptionParser(demo.__doc__)
    opts, args = p.parse_args(args)

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

    panel_space = .23
    dup_space = .025
    # Draw a gene and two regulatory elements at these arbitrary locations
    locs = [(.5, .9), # ancestral gene
            (.5, .9 - panel_space + dup_space), # identical copies
            (.5, .9 - panel_space - dup_space),
            (.5, .9 - 2 * panel_space + dup_space), # degenerate copies
            (.5, .9 - 2 * panel_space - dup_space),
            (.2, .9 - 3 * panel_space + dup_space), # sub-functionalization
            (.2, .9 - 3 * panel_space - dup_space),
            (.5, .9 - 3 * panel_space + dup_space), # neo-functionalization
            (.5, .9 - 3 * panel_space - dup_space),
            (.8, .9 - 3 * panel_space + dup_space), # non-functionalization
            (.8, .9 - 3 * panel_space - dup_space),
            ]

    default_regulator = "gm"
    regulators = [default_regulator,
            default_regulator, default_regulator,
            "wm", default_regulator,
            "wm", "gw",
            "wb", default_regulator,
            "ww", default_regulator,
            ]

    width = .24
    for i, (xx, yy) in enumerate(locs):
        regulator = regulators[i]
        x1, x2 = xx - .5 * width, xx + .5 * width
        Glyph(root, x1, x2, yy)
        if i == 9:  # upper copy for non-functionalization
            continue

        # coding region
        x1, x2 = xx - .16 * width, xx + .45 * width
        Glyph(root, x1, x2, yy, fc="k")

        # two regulatory elements
        x1, x2 = xx - .4 * width, xx - .28 * width
        for xx, fc in zip((x1, x2), regulator):
            if fc == 'w':
                continue

            DoubleCircle(root, xx, yy, fc=fc)

        rotation = 30
        tip = .02
        if i == 0:
            ya = yy + tip
            root.text(x1, ya, "Flower", rotation=rotation, va="bottom")
            root.text(x2, ya, "Root", rotation=rotation, va="bottom")
        elif i == 7:
            ya = yy + tip
            root.text(x2, ya, "Leaf", rotation=rotation, va="bottom")

    # Draw arrows between panels (center)
    arrow_dist = .08
    ar_xpos = .5
    for ar_ypos in (.3, .53, .76):
        root.annotate(" ", (ar_xpos, ar_ypos),
                (ar_xpos, ar_ypos + arrow_dist),
                arrowprops=arrowprops)

    ar_ypos = .3
    for ar_xpos in (.2, .8):
        root.annotate(" ", (ar_xpos, ar_ypos),
                (.5, ar_ypos + arrow_dist),
                arrowprops=arrowprops)

    # Duplication, Degeneration
    xx = .6
    ys = (.76, .53)
    processes = ("Duplication", "Degeneration")
    for yy, process in zip(ys, processes):
        root.text(xx, yy + .02, process, fontweight="bold")

    # Label of fates
    xs = (.2, .5, .8)
    fates = ("Subfunctionalization", "Neofunctionalization",
            "Nonfunctionalization")
    yy = .05
    for xx, fate in zip(xs, fates):
        RoundLabel(root, xx, yy, fate)

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

    figname = "demo.pdf"
    savefig(figname, dpi=300)
103326
def gc(args):
    
    p = OptionParser(gc.__doc__)
    p.add_option("--binsize", default=500, type="int",
                 help="Bin size to use")
    opts, args = p.parse_args(args)

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

    fastafile, = args
    binsize = opts.binsize
    allbins = []
    for name, seq in parse_fasta(fastafile):
        for i in range(len(seq) / binsize):
            atcnt = gccnt = 0
            for c in seq[i * binsize: (i + 1) * binsize].upper():
                if c in "AT":
                    atcnt += 1
                elif c in "GC":
                    gccnt += 1
            totalcnt = atcnt + gccnt
            if totalcnt == 0:
                continue
            gcpct = gccnt * 100 / totalcnt
            allbins.append(gcpct)

    from jcvi.graphics.base import asciiplot
    from collections import Counter

    title = "Total number of bins={}".format(len(allbins))
    c = Counter(allbins)
    x, y = zip(*sorted(c.items()))
    asciiplot(x, y, title=title)