Visualize data on a kegg pathway.


    :param df:
    :param pathway:
    :param a:
    :param b:
    :param ids_from:
    :param cmap:
    :param is_log2:
    :param fillna:
    :param z_score:
    :return:
9426
def kegg_pathway(df, pathway, a, b=None, ids_from="Proteins", cmap=cm.PuOr_r, is_log2=False, fillna=None, z_score=1):
    

    df = df.copy()

    if np.any(df.values < 0) and not is_log2:
        warnings.warn("Input data has negative values. If data is log2 transformed, set is_log2=True.")

    if fillna is not None:
        df = df.fillna(fillna)

    if z_score is None:
        pass
    elif z_score == 0:
        df = (df - df.median(axis=0)) / df.std(axis=0)
    elif z_score == 1:
        df = ((df.T - df.median(axis=1).T) / df.std(axis=1).T).T


    if b is not None:
        # Calculate ratio between two groups
        g1, g2 = df[a].values, df[b].values

        if is_log2:
            dr = np.nanmean(g2, axis=1) - np.nanmean(g1, axis=1)
        else:
            dr = np.log2(np.nanmean(g2, axis=1) / np.nanmean(g1, axis=1))

    else:
        g1 = df[a].values
        dr = np.nanmean(g1, axis=1)


    maxi = np.max(abs(dr))
    norm = mpl.colors.Normalize(vmin=-maxi, vmax=maxi)
    mapper = cm.ScalarMappable(norm=norm, cmap=cm.PuOr_r) # Orange up

    node_colors = {}
    for p, v in zip(df.index.get_level_values(ids_from), dr):
        pid = str(p).split(";")[-1]

        if "_" in pid:
            pid = pid[:pid.index("_")]

        node_colors[pid] = mpl.colors.rgb2hex(mapper.to_rgba(v))


    global uniprot_kegg_cache

    # Only do this once
    upids = list( node_colors.keys() )
    upids = [p for p in upids if p not in uniprot_kegg_cache.keys()]

    if upids:
        new_pairs = get_uniprot_id_mapping_pairs('ACC+ID', 'KEGG_ID', upids)
        uniprot_kegg_cache.update(new_pairs)

        for p in upids:
            if p not in uniprot_kegg_cache:
                uniprot_kegg_cache[p] = None # Not found, don't look again


    with StringIO() as f:
        f.write('#hsa\tData\n')
        for k, c in list(node_colors.items()):
            if k in uniprot_kegg_cache and uniprot_kegg_cache[k] is not None:
                kids = uniprot_kegg_cache[k]
                for kegg_id in kids:
                    f.write('%s\t%s\n' % (kegg_id.split(':')[-1], c ))

        # Reset file
        f.seek(0)

        url = 'https://www.kegg.jp/kegg-bin/mcolor_pathway'
        m = MultipartEncoder(
            fields={
                'map': pathway,
                'mapping_list': ('filename', f),
                'mode': 'color',
                'submit': 'Exec',
                'reference': 'white',
                 }
        )

    r = requests.post(url, data=m, headers={'Content-Type': m.content_type})
    if r.status_code == 200:
        # src="/tmp/mark_pathway154353327948969/hsa04010.1.png"
        ms = re.finditer('src="(/tmp/mark_pathway[^"]*.png)"', r.text)
        m = list(ms).pop()

        # Download image data
        image = Image.open(requests.get('http://www.kegg.jp%s' % m.group(1), stream=True).raw)
        width, height = image.size   # Get dimensions
        image = image.crop((1, 1, width-1, height-1)) # Crop black outline
        print("Scale range: %.2f .. %.2f" % (-maxi, maxi))

        return image
1017856
def load_synapses(path=HOME + "/Downloads/pinky100_final.df",
                  scaling=(1, 1, 1)):
    

    scaling = np.array(list(scaling))

    df = pd.read_csv(path)

    locs = np.array(df[["presyn_x", "centroid_x", "postsyn_x"]])

    mask = ~np.any(np.isnan(locs), axis=1)

    df = df[mask]

    df['pre_pt.position'] = list((np.array(df[['presyn_x', 'presyn_y', 'presyn_z']]) / scaling).astype(np.int))
    df['ctr_pt.position'] = list((np.array(df[['centroid_x', 'centroid_y', 'centroid_z']]) / scaling).astype(np.int))
    df['post_pt.position'] = list((np.array(df[['postsyn_x', 'postsyn_y', 'postsyn_z']]) / scaling).astype(np.int))

    df = df[['pre_pt.position', 'ctr_pt.position', 'post_pt.position', 'size']]

    return df
325981
def plot(args):
    
    from jcvi.graphics.base import plt, savefig
    from jcvi.graphics.chromosome import ChromosomeMap

    p = OptionParser(plot.__doc__)
    p.add_option("--firstn", type="int", help="Only plot the first N genes")
    p.add_option("--ymax", type="int", help="Y-axis max value")
    p.add_option("--log", action="store_true",
                help="Write plotting data [default: %default]")
    opts, args, iopts = p.set_image_options(args, figsize="6x4")

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

    taggedbed, chr = args
    bed = Bed(taggedbed)
    beds = list(bed.sub_bed(chr))
    old, new = [], []
    i = 0
    for b in beds:
        accn = b.extra[0]
        if "te" in accn:
            continue

        accn, tag = accn.split("|")
        if tag == "OVERLAP":
            continue

        c, r = atg_name(accn)
        if tag == "NEW":
            new.append((i, r))
        else:
            old.append((i, r))
        i += 1

    ngenes = i
    assert ngenes == len(new) + len(old)

    logging.debug("Imported {0} ranks on {1}.".format(ngenes, chr))
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    xstart, xend = .2, .8
    ystart, yend = .2, .8
    pad = .02

    ngenes = opts.firstn or ngenes
    ymax = opts.ymax or 500000

    title = "Assignment of Medtr identifiers"
    if opts.ymax:
        subtitle = "{0}, first {1} genes".format(chr, ngenes)
    else:
        subtitle = "{0}, {1} genes ({2} new)".format(chr, ngenes, len(new))

    chr_map = ChromosomeMap(fig, root, xstart, xend, ystart, yend, pad, 0,
                        ymax, 5, title, subtitle)

    ax = chr_map.axes

    if opts.log:
        from jcvi.utils.table import write_csv
        header = ["x", "y"]
        write_csv(header, new, filename=chr + ".new")
        write_csv(header, old, filename=chr + ".old")

    x, y = zip(*new)
    ax.plot(x, y, "b,")
    x, y = zip(*old)
    ax.plot(x, y, "r,")

    # Legends
    ymid = (ystart + yend) / 2
    y = ymid + pad
    root.plot([.2], [y], "r.", lw=2)
    root.text(.2 + pad, y, "Existing Medtr ids", va="center", size=10)
    y = ymid - pad
    root.plot([.2], [y], "b.", lw=2)
    root.text(.2 + pad, y, "Newly instantiated ids", va="center", size=10)

    ax.set_xlim(0, ngenes)
    ax.set_ylim(0, ymax)
    ax.set_axis_off()

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

    image_name = chr + ".identifiers." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
1305843
def plot_colormap_components(cmap):
    """"""
    from ._helpers import set_ax_labels  # recursive import protection

    plt.figure(figsize=[8, 4])
    gs = grd.GridSpec(3, 1, height_ratios=[1, 10, 1], hspace=0.05)
    # colorbar
    ax = plt.subplot(gs[0])
    gradient = np.linspace(0, 1, 256)
    gradient = np.vstack((gradient, gradient))
    ax.imshow(gradient, aspect="auto", cmap=cmap, vmin=0., vmax=1.)
    ax.set_title(cmap.name, fontsize=20)
    ax.set_axis_off()
    # components
    ax = plt.subplot(gs[1])
    x = np.arange(cmap.N)
    colors = cmap(x)
    r = colors[:, 0]
    g = colors[:, 1]
    b = colors[:, 2]
    RGB_weight = [0.299, 0.587, 0.114]
    k = np.sqrt(np.dot(colors[:, :3] ** 2, RGB_weight))
    r.clip(0, 1, out=r)
    g.clip(0, 1, out=g)
    b.clip(0, 1, out=b)
    xi = np.linspace(0, 1, x.size)
    plt.plot(xi, r, "r", linewidth=5, alpha=0.6)
    plt.plot(xi, g, "g", linewidth=5, alpha=0.6)
    plt.plot(xi, b, "b", linewidth=5, alpha=0.6)
    plt.plot(xi, k, "k", linewidth=5, alpha=0.6)
    ax.set_xlim(0, 1)
    ax.set_ylim(-0.1, 1.1)
    set_ax_labels(ax=ax, xlabel=None, xticks=False, ylabel="intensity")
    # grayified colorbar
    cmap = grayify_cmap(cmap)
    ax = plt.subplot(gs[2])
    gradient = np.linspace(0, 1, 256)
    gradient = np.vstack((gradient, gradient))
    ax.imshow(gradient, aspect="auto", cmap=cmap, vmin=0., vmax=1.)
    ax.set_axis_off()
684662
def dlogpdf_df_dtheta(self, f, y, Y_metadata=None):
        
        if self.size > 0:
            if self.not_block_really:
                raise NotImplementedError("Need to make a decorator for this!")
            if isinstance(self.gp_link, link_functions.Identity):
                return self.dlogpdf_dlink_dtheta(f, y, Y_metadata=Y_metadata)
            else:
                inv_link_f = self.gp_link.transf(f)
                dlink_df = self.gp_link.dtransf_df(f)
                dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(inv_link_f, y, Y_metadata=Y_metadata)

                dlogpdf_df_dtheta = np.zeros((self.size, f.shape[0], f.shape[1]))
                #Chain each parameter of hte likelihood seperately
                for p in range(self.size):
                    dlogpdf_df_dtheta[p, :, :] = chain_1(dlogpdf_dlink_dtheta[p,:,:], dlink_df)
                return dlogpdf_df_dtheta
                #return chain_1(dlogpdf_dlink_dtheta, dlink_df)
        else:
            # There are no parameters so return an empty array for derivatives
            return np.zeros((0, f.shape[0], f.shape[1]))
802976
def composite_correlation(df, size=(12, 8)):
    
    fig = plt.figure(1, size)
    ax1 = plt.subplot2grid((2, 2), (0, 0))
    ax2 = plt.subplot2grid((2, 2), (0, 1))
    ax3 = plt.subplot2grid((2, 2), (1, 0))
    ax4 = plt.subplot2grid((2, 2), (1, 1))
    chemistry = ["V1", "V2", "V2.5", float("nan")]
    colors = sns.color_palette("Set2", 8)
    color_map = dict(zip(chemistry, colors))

    age_label = "Chronological age (yr)"
    ax1.scatter(df["hli_calc_age_sample_taken"], df["teloLength"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax1.set_ylim(0, 15)
    ax1.set_ylabel("Telomere length (Kb)")

    ax2.scatter(df["hli_calc_age_sample_taken"], df["ccn.chrX"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax2.set_ylim(1.8, 2.1)
    ax2.set_ylabel("ChrX copy number")

    ax4.scatter(df["hli_calc_age_sample_taken"], df["ccn.chrY"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax4.set_ylim(0.8, 1.1)
    ax4.set_ylabel("ChrY copy number")

    ax3.scatter(df["hli_calc_age_sample_taken"], df["TRA.PPM"],
                s=10, marker='.',
                color=df["Chemistry"].map(color_map))
    ax3.set_ylim(0, 250)
    ax3.set_ylabel("$TCR-\\alpha$ deletions (count per million reads)")

    from matplotlib.lines import Line2D
    legend_elements = [Line2D([0], [0], marker='.', color='w', label=chem,
                          markerfacecolor=color, markersize=16) \
                        for (chem, color) in zip(chemistry, colors)[:3]]
    for ax in (ax1, ax2, ax3, ax4):
        ax.set_xlabel(age_label)
        ax.legend(handles=legend_elements, loc="upper right")

    plt.tight_layout()
    root = fig.add_axes((0, 0, 1, 1))
    labels = ((.02, .98, "A"),
              (.52, .98, "B"),
              (.02, .5, "C"),
              (.52, .5, "D"))
    panel_labels(root, labels)
    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()
913918
def _process_pathway_ko(self, limit):
        
        LOG.info("Processing KEGG pathways to kegg ortholog classes")
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        line_counter = 0

        raw = '/'.join((self.rawdir, self.files['pathway_ko']['file']))
        with open(raw, 'r', encoding="iso-8859-1") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            for row in filereader:
                line_counter += 1
                (ko_id, pathway_id) = row

                if self.test_mode and pathway_id not in self.test_ids['pathway']:
                    continue

                pathway_id = 'KEGG-' + pathway_id
                ko_id = 'KEGG-' + ko_id

                p = Pathway(graph)
                p.addGeneToPathway(ko_id, pathway_id)

                if not self.test_mode and limit is not None and line_counter > limit:
                    break

        return
7395
def create_pep_protein_quant_lookup(fns, pqdb, poolnames, featcolnr, patterns,
                                    storefuns, isobqcolpattern=None,
                                    psmnrpattern=None):
    """"""
    tablefn_map = create_tablefn_map(fns, pqdb, poolnames)
    feat_map = pqdb.get_feature_map()
    for pattern, storefun in zip(patterns, storefuns):
        if pattern is None:
            continue
        colmap = get_colmap(fns, pattern, single_col=True)
        if colmap:
            store_single_col_data(fns, tablefn_map, feat_map,
                                  storefun, featcolnr, colmap)
    if isobqcolpattern is not None:
        isocolmap = get_colmap(fns, isobqcolpattern, antipattern=psmnrpattern)
    else:
        return
    if psmnrpattern is not None:
        psmcolmap = get_colmap(fns, psmnrpattern)
    else:
        psmcolmap = False
    create_isobaric_quant_lookup(fns, tablefn_map,
                                 feat_map, pqdb,
                                 featcolnr,
                                 isocolmap, psmcolmap)
176584
def coda(df, window, level):
    
    # pull out the data
    d = df.values

    # smooth the data and standardize it
    smooth_data = movingaverage(d, df.index, window)[0]
    stand_data = (smooth_data - smooth_data.mean()) / smooth_data.std()

    # scale the data to have unit length
    scale_data = d / np.sqrt(np.sum(d ** 2, axis=0))

    # calculate the "mass chromatographic quality" (MCQ) index
    mcq = np.sum(stand_data * scale_data, axis=0) / np.sqrt(d.shape[0] - 1)

    # filter out ions with an mcq below level
    good_ions = [i for i, q in zip(df.columns, mcq) if q >= level]
    return good_ions
295836
def plot_log_histogram(df, palette, title, histnorm=""):
    
    data = [go.Histogram(x=np.log10(df.loc[df["dataset"] == d, "lengths"]),
                         opacity=0.4,
                         name=d,
                         histnorm=histnorm,
                         marker=dict(color=c))
            for d, c in zip(df["dataset"].unique(), palette)]
    xtickvals = [10**i for i in range(10) if not 10**i > 10 * np.amax(df["lengths"])]
    html = plotly.offline.plot(
        {"data": data,
         "layout": go.Layout(barmode='overlay',
                             title=title,
                             xaxis=dict(tickvals=np.log10(xtickvals),
                                        ticktext=xtickvals))},
        output_type="div",
        show_link=False)
    fig = go.Figure(
        {"data": data,
         "layout": go.Layout(barmode='overlay',
                             title=title,
                             xaxis=dict(tickvals=np.log10(xtickvals),
                                        ticktext=xtickvals))})
    return html, fig
373018
def load_from_args_as_dataframe(args):
    '''
    
    '''
    if not args.variants and not args.single_variant:
        return None

    if args.variant_source_name:
        variant_source_names = util.expand(
            args.variant_source_name,
            'variant_source_name',
            'variant source',
            len(args.variants))
    else:
        variant_source_names = util.drop_prefix(args.variants)

    variant_to_sources = collections.defaultdict(list)

    dfs = []
    for i in range(len(args.variants)):
        name = variant_source_names[i]
        prefix = (
            'metadata:' if len(args.variants) == 1 else "metadata:%s:" % name)
        df = load_as_dataframe(
            args.variants[i],
            name=name,
            genome=args.genome,
            max_variants=args.max_variants_per_source,
            only_passing=not args.include_failing_variants,
            metadata_column_prefix=prefix)

        if df.shape[0] == 0:
            logging.warn("No variants loaded from: %s" % args.variants[i])
        else:
            for variant in df.variant:
                variant_to_sources[variant].append(name)
            dfs.append(df)

    if args.single_variant:
        variants = []
        extra_args = {}
        if args.genome:
            extra_args = {
                'ensembl': varcode.reference.infer_genome(args.genome)
            }
        for (locus_str, ref, alt) in args.single_variant:
            locus = Locus.parse(locus_str)
            variant = varcode.Variant(
                    locus.contig,
                    locus.inclusive_start,
                    ref,
                    alt,
                    **extra_args)
            variants.append(variant)
            variant_to_sources[variant].append("commandline")
        dfs.append(variants_to_dataframe(variants))

    df = dfs.pop(0)
    for other_df in dfs:
        df = pandas.merge(
            df,
            other_df,
            how='outer',
            on=["variant"] + STANDARD_DATAFRAME_COLUMNS)

    genomes = df["genome"].unique()
    if len(genomes) > 1:
        raise ValueError(
                "Mixing references is not supported. "
                "Reference genomes: %s" % (", ".join(genomes)))

    df["sources"] = [" ".join(variant_to_sources[v]) for v in df.variant]

    # Apply filters:
    if args.ref:
        df = df.ix[df.ref.isin(args.ref)]
    if args.alt:
        df = df.ix[df.alt.isin(args.alt)]
    loci = loci_util.load_from_args(
        util.remove_prefix_from_parsed_args(args, "variant"))
    if loci is not None:
        df = df.ix[[
            loci.intersects(pileup_collection.to_locus(v))
            for v in df.variant
        ]]
    return df
342766
def zscore_df(df, axis='row', keep_orig=False):
  '''
  
  '''
  df_z = {}

  for mat_type in df:
    if keep_orig and mat_type == 'mat':
      mat_orig = deepcopy(df[mat_type])

    inst_df = df[mat_type]

    if axis == 'row':
      inst_df = inst_df.transpose()

    df_z[mat_type] = (inst_df - inst_df.mean())/inst_df.std()

    if axis == 'row':
      df_z[mat_type] = df_z[mat_type].transpose()

  if keep_orig:
    df_z['mat_orig'] = mat_orig

  return df_z
717292
def _process_pathway(self, row):
        
        model = Model(self.graph)
        self._check_list_len(row, 4)
        (gene_symbol, gene_id, pathway_name, pathway_id) = row

        if self.test_mode and (int(gene_id) not in self.test_geneids):
            return

        entrez_id = 'NCBIGene:' + gene_id

        pathways_to_scrub = [
            'REACT:REACT_116125',  # disease
            "REACT:REACT_111045",  # developmental biology
            "REACT:REACT_200794",  # Mus musculus biological processes
            "REACT:REACT_13685"]   # neuronal system ?

        if pathway_id in pathways_to_scrub:
            # these are lame "pathways" like generic
            # "disease" and "developmental biology"
            return

        # convert KEGG pathway ids... KEGG:12345 --> KEGG-path:map12345
        if re.match(r'KEGG', pathway_id):
            pathway_id = re.sub(r'KEGG:', 'KEGG-path:map', pathway_id)
        # just in case, add it as a class
        model.addClassToGraph(entrez_id, None)

        self.pathway.addPathway(pathway_id, pathway_name)
        self.pathway.addGeneToPathway(entrez_id, pathway_id)

        return
1073959
def df_pdb_ranking(self):
        """"""
        df = pd.DataFrame()
        for g in self.genes_with_experimental_structures:
            protein_df = g.protein.df_pdb_ranking.copy().reset_index()
            if not protein_df.empty:
                protein_df['gene'] = g.id
                df = df.append(protein_df)
        if df.empty:
            log.warning('Empty dataframe')
            return df
        else:
            return ssbio.utils.clean_df(df.set_index('gene'))
1156483
def _cnvkit_prioritize(sample, genes, allele_file, metrics_file):
    
    mdf = pd.read_table(metrics_file)
    mdf.columns = [x.lower() for x in mdf.columns]
    if len(genes) > 0:
        mdf = mdf[mdf["gene"].str.contains("|".join(genes))]
    mdf = mdf[["chromosome", "start", "end", "gene", "log2", "ci_hi", "ci_lo"]]
    adf = pd.read_table(allele_file)
    if len(genes) > 0:
        adf = adf[adf["gene"].str.contains("|".join(genes))]
    if "cn1" in adf.columns and "cn2" in adf.columns:
        adf = adf[["chromosome", "start", "end", "cn", "cn1", "cn2"]]
    else:
        adf = adf[["chromosome", "start", "end", "cn"]]
    df = pd.merge(mdf, adf, on=["chromosome", "start", "end"])
    df = df[df["cn"] != 2]
    if len(df) > 0:
        def passes(row):
            spread = abs(row["ci_hi"] - row["ci_lo"])
            return spread < 0.25
        df["passes"] = df.apply(passes, axis=1)
    df.insert(0, "sample", [sample] * len(df))
    return df
818397
def write_perseus(f, df):
    

    ### Generate the Perseus like type index

    FIELD_TYPE_MAP = {
        'Amino acid':'C',
        'Charge':'C',
        'Reverse':'C',
        'Potential contaminant':'C',
        'Multiplicity':'C',
        'Localization prob':'N',
        'PEP':'N',
        'Score':'N',
        'Delta score':'N',
        'Score for localization':'N',
        'Mass error [ppm]':'N',
        'Intensity':'N',
        'Position':'N',
        'Proteins':'T',
        'Positions within proteins':'T',
        'Leading proteins':'T',
        'Protein names':'T',
        'Gene names':'T',
        'Sequence window':'T',
        'Unique identifier':'T',
    }

    def map_field_type(n, c):
        try:
            t = FIELD_TYPE_MAP[c]
        except:
            t = "E"

        # In the first element, add type indicator
        if n == 0:
            t = "#!{Type}%s" % t

        return t

    df = df.copy()
    df.columns = pd.MultiIndex.from_tuples([(k, map_field_type(n, k)) for n, k in enumerate(df.columns)], names=["Label","Type"])
    df = df.transpose().reset_index().transpose()
    df.to_csv(f, index=False, header=False)
378299
def score_hist(df, columns=None, groupby=None, threshold=0.7, stacked=True,
               bins=20, percent=True, alpha=0.33, show=True, block=False, save=False):
    
    df = df if columns is None else df[([] if groupby is None else [groupby]) + list(columns)].copy()
    if groupby is not None or threshold is not None:
        df = groups_from_scores(df, groupby=groupby, threshold=threshold)
    percent = 100. if percent else 1.
    if isinstance(df, pd.core.groupby.DataFrameGroupBy):
        df = df_from_groups(df, columns=columns) * percent
    columns = df.columns if columns is None else columns
    if bins is None:
        bins = 20
    if isinstance(bins, int):
        bins = np.linspace(np.min(df.min()), np.max(df.max()), bins)
    log.debug('bins: {}'.format(bins))
    figs = []

    df.plot(kind='hist', alpha=alpha, stacked=stacked, bins=bins)
    # for col in df.columns:
    #     series = df[col] * percent
    #     log.debug('{}'.format(series))
    #     figs.append(plt.hist(series, bins=bins, alpha=alpha,
    #                 weights=percent * np.ones_like(series) / len(series.dropna()),
    #                 label=stringify(col)))
    plt.legend()
    plt.xlabel('Score (%)')
    plt.ylabel('Percent')
    plt.title('{} Scores for {}'.format(np.sum(df.count()), columns))
    plt.draw()
    if save or not show:
        fig = plt.gcf()
        today = datetime.datetime.today()
        fig.savefig(os.path.join(IMAGES_PATH, 'score_hist_{:04d}-{:02d}-{:02d}_{:02d}{:02d}.jpg'.format(*today.timetuple())))
    if show:
        plt.show(block=block)
    return figs