Omics Plotting

omics-plotting

Overview

When the user wants a figure, generate it with matplotlib / seaborn, applying the shared style block below. The user can hand-tune colors, fonts, or spines per plot, but unless they ask for something specific, paste the style block and reuse the palette so a whole analysis reads as one figure system at a glance.

This skill is self-contained: everything you need (style, palette, recipes) is in this document.

When to use

• The user asks for a plot / figure / chart / visualization from a results table or an in-memory DataFrame (DEG table, enrichment result, expression matrix, long-form measurements, survival table…).
• You are preparing figures for a report, a paper submission or presentation and want a consistent publication style.

Combining several plots into one multi-panel composite, or assembling user-supplied PNG/PDF panels, is handled by the sibling multipanel skill — use this skill to draw each individual panel.

Do NOT use for

• Interactive dashboards or web charts (this is static matplotlib output).
• 3D molecular structure rendering (that is the structure viewer, not a plot).

Key Concepts

One consistent figure system

The core idea is that every figure from a single analysis should look like it came from the same publication. That is enforced by two shared objects: the PUB_STYLE rcParams block (fonts, spines, DPI, editable vector text) and a fixed PALETTE / directional color set (UP, DOWN, NS). Paste both at the top of every plot script and map the same group or direction to the same color across panels, so a reader can carry meaning from one figure to the next.

Diverging vs sequential colormaps

Color encoding is not free choice. Use the diverging colormap (DIVERGING_CMAP = "RdBu_r", always center=0, vmin=-vmax) for signed quantities where zero is meaningful — z-scores, log2 fold changes, correlations. Use the sequential colormap (SEQUENTIAL_CMAP = "viridis") for unsigned magnitudes — densities, -log10 p, counts. Mixing these (a sequential map on signed data) hides the sign and misleads the reader.

Data shape drives figure type

Each recipe expects a specific table shape: a per-gene DEG table (volcano, MA), a genes × samples matrix (heatmap), a samples × features matrix (PCA/UMAP), or long-form tidy rows (box/violin/bar, ridgeline, Kaplan–Meier). Identifying the shape first — then reading the header to confirm the real column names — is what selects the recipe. The column names in each recipe are defaults to override, not fixed requirements.

Decision Framework

Pick the figure type from what the data represents and what question it answers:

Copy code
What does the table hold?
├─ Per-gene stats (log2FC, padj)
│   ├─ emphasize significance ......... Volcano
│   └─ emphasize expression level ..... MA plot
├─ genes × samples matrix
│   ├─ show patterns/clusters ......... Clustered expression heatmap (z-score)
│   └─ show sample-sample QC .......... Correlation heatmap
├─ Enrichment / gene-set result
│   ├─ signed effect (NES) ............ GSEA bar
│   └─ ratio + size + significance .... GSEA dot plot
├─ Long-form measurements (x, y)
│   ├─ compare distributions .......... Box / Violin
│   ├─ compare means .................. Bar (with error bars)
│   └─ many groups, shape matters ..... Ridgeline
├─ samples × features (high-dim) ...... PCA / UMAP / t-SNE
└─ time-to-event + group ............. Kaplan–Meier

Workflow

1. Identify the data source — a workspace-relative CSV/TSV path or a DataFrame already in memory — and the figure type (pick from the table below). If the required columns are unclear, inspect the table's header first.
2. Write one python script: paste the style block, load the data, draw the plot with the matching recipe, and save to a workspace-relative path under figures/.
3. Report the saved path back to the user (and reference it in any report / deck by that relative path, e.g. ![Volcano](figures/volcano.png)).

Shared style — paste at the top of every plot script

pythonCopy code
import matplotlib.pyplot as plt

# Publication style (colorblind-friendly, editable vector text, no top/right spines)
PUB_STYLE = {
    "figure.dpi": 110, "savefig.dpi": 300, "savefig.bbox": "tight",
    "font.family": "sans-serif",
    "font.sans-serif": ["Arial", "Liberation Sans", "Nimbus Sans", "Helvetica", "DejaVu Sans"],
    "font.size": 11, "axes.titlesize": 13, "axes.titleweight": "bold",
    "figure.titlesize": 13, "figure.titleweight": "bold",
    "axes.labelsize": 12, "axes.linewidth": 1.0,
    "axes.spines.top": False, "axes.spines.right": False,
    "xtick.labelsize": 10, "ytick.labelsize": 10,
    "xtick.direction": "out", "ytick.direction": "out",
    "legend.frameon": False, "legend.fontsize": 9,
    "svg.fonttype": "none", "pdf.fonttype": 42, "ps.fonttype": 42,
}
plt.rcParams.update(PUB_STYLE)   # or: with plt.rc_context(PUB_STYLE): ...

# Palette — reuse the SAME colors across every panel of an analysis
UP, DOWN, NS = "#d73721", "#204897", "#d9d9d9"   # up / down / not-significant
PALETTE = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4",
           "#008300", "#4a3aa7", "#e34948", "#12a4c0", "#a66a2e"]  # categorical (CVD-safe order)
GROUP_COLORS = {"Group1": "#204897", "Group2": "#e34948", "Group3": "#E7B800"}
DIVERGING_CMAP = "RdBu_r"    # z-score / log2FC heatmaps — set center=0, vmin=-vmax
SEQUENTIAL_CMAP = "viridis"  # magnitude / -log10 p / density

Multi-panel / composite figures

For combining several plots into one multi-panel journal figure (panels A, B, C…), or assembling already-rendered PNG/PDF panels the user supplies, use the sibling multipanel skill — it owns the composition discipline (one subplot_mosaic canvas, per-panel legends, correctly placed panel letters, text-legibility rules, image assembly). Draw each panel with the single-panel recipes below, then compose per that skill. The recipes here each build their own figure, so do not call them directly for a composite — copy the recipe body onto a mosaic axis as multipanel describes.

Plot catalogue

Pick the recipe by figure type. Columns listed are the defaults — override the column-name variables to match the actual table.

Recipes

Each is a full python script body. Adjust column names, thresholds, and the save path. All save under figures/.

Volcano (-log10 p vs log2 fold change):

pythonCopy code
import numpy as np, pandas as pd
df = pd.read_csv("deg_results.csv").dropna(subset=["log2FoldChange", "padj"])
fc, p = df["log2FoldChange"].to_numpy(float), df["padj"].to_numpy(float)
nlp = -np.log10(np.clip(p, 1e-300, None))
fc_t, p_t = 0.58, 0.05
up, down = (fc >= fc_t) & (p < p_t), (fc <= -fc_t) & (p < p_t)
ns = ~(up | down)
fig, ax = plt.subplots(figsize=(7, 6))
ax.scatter(fc[ns], nlp[ns], c=NS, s=12, alpha=0.5, edgecolors="none", rasterized=True, label=f"NS ({ns.sum()})")
ax.scatter(fc[down], nlp[down], c=DOWN, s=18, alpha=0.85, edgecolors="none", label=f"Down ({down.sum()})")
ax.scatter(fc[up], nlp[up], c=UP, s=18, alpha=0.85, edgecolors="none", label=f"Up ({up.sum()})")
for v in (fc_t, -fc_t): ax.axvline(v, ls="--", lw=0.8, color="0.5")
ax.axhline(-np.log10(p_t), ls="--", lw=0.8, color="0.5")
lab = (df["gene"] if "gene" in df else pd.Series(df.index)).astype(str).to_numpy()
sig = np.where(up | down)[0]
top = sig[np.argsort(nlp[sig])[::-1][:10]]      # standalone: top ~10; composite panel: cut to <=5
try:                                            # repel labels so they never overlap
    from adjustText import adjust_text
    texts = [ax.text(fc[i], nlp[i], lab[i], fontsize=7) for i in top]
    adjust_text(texts, ax=ax, expand=(1.3, 1.6),
                arrowprops=dict(arrowstyle="-", color="0.6", lw=0.5))
except ImportError:                             # no adjustText -> label fewer, with an offset
    for i in top[:5]:
        ax.annotate(lab[i], (fc[i], nlp[i]), xytext=(6, 6), textcoords="offset points",
                    fontsize=7, ha="left", va="bottom")
ax.set_xlabel(r"$\log_{2}$ fold change"); ax.set_ylabel(r"$-\log_{10}$ padj")
ax.set_title("Volcano plot"); ax.legend(loc="upper right", markerscale=1.4)
fig.tight_layout(); fig.savefig("figures/volcano.png")

MA plot (log2 fold change vs mean expression; columns baseMean, log2FoldChange, optional padj):

pythonCopy code
import numpy as np, pandas as pd
df = pd.read_csv("deg_results.csv").dropna(subset=["baseMean", "log2FoldChange"])
x = np.log10(df["baseMean"].to_numpy(float) + 1)
fc = df["log2FoldChange"].to_numpy(float)
sig = (df["padj"].to_numpy(float) < 0.05) if "padj" in df else np.zeros(len(df), bool)
fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(x[~sig], fc[~sig], c=NS, s=10, alpha=0.5, edgecolors="none", rasterized=True, label="NS")
ax.scatter(x[sig], fc[sig], c=UP, s=14, alpha=0.85, edgecolors="none", label=f"padj<0.05 ({int(sig.sum())})")
ax.axhline(0, color="0.4", lw=0.8)
ax.set_xlabel(r"$\log_{10}$(mean norm. count + 1)"); ax.set_ylabel(r"$\log_{2}$ fold change")
ax.set_title("MA plot"); ax.legend(loc="upper right")
fig.tight_layout(); fig.savefig("figures/ma.png")

Clustered expression heatmap (z-scored, seaborn clustermap):

pythonCopy code
import seaborn as sns, pandas as pd
mat = pd.read_csv("expression.csv", index_col=0).select_dtypes("number")
g = sns.clustermap(mat, z_score=0, cmap=DIVERGING_CMAP, center=0,
                   figsize=(8, 8), xticklabels=True, yticklabels=mat.shape[0] <= 60,
                   cbar_kws={"label": "z-score"})
g.ax_col_dendrogram.set_title("Expression heatmap", pad=12)
g.figure.savefig("figures/heatmap.png")

Correlation heatmap (sample QC):

pythonCopy code
import numpy as np, seaborn as sns, pandas as pd
corr = pd.read_csv("expr.csv", index_col=0).select_dtypes("number").corr(method="pearson")
mask = np.triu(np.ones_like(corr, dtype=bool), k=1)
fig, ax = plt.subplots(figsize=(7, 6))
sns.heatmap(corr, mask=mask, cmap=DIVERGING_CMAP, vmin=-1, vmax=1, center=0,
            annot=True, # if the n>20 samples then set annot=False
            fmt=".2f", annot_kws={"size": 7}, square=True,
            linewidths=0.5, 
            cbar_kws={"label": "Pearson r", "shrink": 0.7}, 
            ax=ax)
ax.set_title("Sample correlation"); fig.tight_layout(); fig.savefig("figures/corr.png")

GSEA bar (top gene sets, colored by direction):

pythonCopy code
import pandas as pd
df = pd.read_csv("gsea.csv").dropna(subset=["Term", "NES"]).copy()
df = df.assign(_a=df["NES"].abs()).nlargest(15, "_a").sort_values("NES")
colors = [UP if v >= 0 else DOWN for v in df["NES"]]
fig, ax = plt.subplots(figsize=(7, 6))
ax.barh(df["Term"].astype(str), df["NES"], color=colors)
ax.axvline(0, color="0.4", lw=0.8); ax.set_xlabel("NES"); ax.set_title("GSEA")
fig.tight_layout(); fig.savefig("figures/gsea_bar.png")

GSEA dot plot (clusterProfiler-style; GeneRatio may be "k/n"):

pythonCopy code
import numpy as np, pandas as pd
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalize
df = pd.read_csv("enrichment.csv").dropna(subset=["Term", "GeneRatio", "Adjusted P-value"]).copy()
df["_ratio"] = df["GeneRatio"].map(lambda v: float(v.split("/")[0]) / float(v.split("/")[1]) if isinstance(v, str) and "/" in v else float(v))
df["_nlp"] = -np.log10(np.clip(df["Adjusted P-value"].astype(float), 1e-300, None))
df = df.nlargest(15, "_nlp").sort_values("_ratio")
cnt = df["Count"].to_numpy(float); sizes = 40 + 220 * (cnt - cnt.min()) / (np.ptp(cnt) + 1e-9)
norm = Normalize(df["_nlp"].min(), df["_nlp"].max())
fig, ax = plt.subplots(figsize=(7, 6))
ax.scatter(df["_ratio"], df["Term"].astype(str), s=sizes, c=df["_nlp"], cmap=SEQUENTIAL_CMAP, norm=norm, edgecolors="0.3", linewidths=0.5, zorder=3)
ax.grid(axis="y", ls=":", color="0.8", zorder=0); ax.set_xlabel("Gene ratio"); ax.set_title("Enrichment")
cb = fig.colorbar(ScalarMappable(norm=norm, cmap=SEQUENTIAL_CMAP), ax=ax, shrink=0.6, pad=0.02)
cb.set_label(r"$-\log_{10}$ adj. $p$"); fig.tight_layout(); fig.savefig("figures/dotplot.png")

Box plot (long-form, jittered points, optional 2-group test):

pythonCopy code
import seaborn as sns, pandas as pd
from scipy import stats
df = pd.read_csv("measurements.csv"); x, y = "group", "value"
fig, ax = plt.subplots(figsize=(6, 5))
sns.boxplot(data=df, x=x, y=y, hue=x, palette=PALETTE[:df[x].nunique()], legend=False, fliersize=0, width=0.6, ax=ax)
sns.stripplot(data=df, x=x, y=y, color="0.25", size=3, alpha=0.6, ax=ax)
lv = list(df[x].dropna().unique())
if len(lv) == 2:  # optional significance star
    a, b = (df.loc[df[x] == l, y].dropna() for l in lv)
    p = stats.mannwhitneyu(a, b, alternative="two-sided").pvalue
    star = "****" if p < 1e-4 else "***" if p < 1e-3 else "**" if p < 1e-2 else "*" if p < 0.05 else "ns"
    ymax = df[y].max(); h = (ymax - df[y].min()) * 0.08
    ax.plot([0, 0, 1, 1], [ymax + h, ymax + 2*h, ymax + 2*h, ymax + h], lw=1.0, c="0.2")
    ax.text(0.5, ymax + 2*h, star, ha="center", va="bottom")
ax.set_title(f"{y} by {x}"); fig.tight_layout(); fig.savefig("figures/box.png")

For violin swap sns.boxplot → sns.violinplot(..., inner="box", cut=0); for bar of the mean use sns.barplot(..., errorbar="se", capsize=0.15).

Grouped embedding scatter (PCA / UMAP / t-SNE / factor scatter): (samples × features + group column):

pythonCopy code
import numpy as np, pandas as pd

def scatter_emb(ax, emb, grp, xi=0, yi=1, ev=None, name="PC", centered=True):
    """emb: (n_samples, n_comp). ev: variance explained % or None. UMAP/t-SNE -> centered=False."""
    emb, grp = np.asarray(emb), pd.Series(grp).astype(str)
    for g in [g for g in GROUP_COLORS if g in set(grp)]:     
        m = (grp == g).to_numpy()
        ax.scatter(emb[m, xi], emb[m, yi], c=GROUP_COLORS[g], s=45, alpha=.9,
                   edgecolors="white", linewidths=.4, label=g)
    if centered:                   # PC/factor are 0-centered and equally spaced
        ax.axhline(0, lw=.5, ls="--", c="0.75"); ax.axvline(0, lw=.5, ls="--", c="0.75")
        ax.set_aspect("equal", adjustable="datalim")
    lab = lambda i: f"{name}{i+1}" + (f" ({ev[i]:.1f}%)" if ev is not None else "")
    ax.set_xlabel(lab(xi)); ax.set_ylabel(lab(yi))
    ax.spines[["top", "right"]].set_visible(False); ax.legend(frameon=False, fontsize=7)

FOR PCA: samples × features + 'group' column

pythonCopy code
from sklearn.decomposition import PCA

df = pd.read_csv("samples_features.csv", index_col=0)    
grp = df.pop("group") if "group" in df else pd.Series("all", index=df.index)
X = df.select_dtypes("number")
p = PCA(n_components=min(5, *X.shape)).fit(X) 
emb, ev = p.transform(X), p.explained_variance_ratio_ * 100
fig, ax = plt.subplots(figsize=(70/25.4, 70/25.4), constrained_layout=True)
scatter_emb(ax, emb, grp, 0, 1, ev, "PC")
fig.savefig("figures/pca.png", dpi=400)

FOR UMAP/t-SNE: no ev

pythonCopy code
import umap           # t-SNE: from sklearn.manifold import TSNE

emb = umap.UMAP(n_neighbors=15, min_dist=0.1).fit_transform(X)   # TSNE(perplexity=30).fit_transform(X)
fig, ax = plt.subplots(figsize=(70/25.4, 70/25.4), constrained_layout=True)
scatter_emb(ax, emb, grp, name="UMAP", centered=False)    # 원점·축 스케일에 의미 없음
fig.savefig("figures/umap.png", dpi=400)

FOR factor scatter:

pythonCopy code
Z  = pd.DataFrame(model.get_factors()["group1"], index=meta.index)
r2 = model.get_variance_explained()["r2_per_factor"]["group1"]    # view × factor
ev = r2.sum(axis=0).to_numpy()                           
keep = [i for i, v in enumerate(ev) if v >= 5.0]        
fig, ax = plt.subplots(figsize=(70/25.4, 70/25.4), constrained_layout=True)
scatter_emb(ax, Z, meta["condition"], keep[0], keep[1], ev, "Factor")
fig.savefig("figures/factor_scatter.png", dpi=400)

Kaplan–Meier (survival by group):

pythonCopy code
import pandas as pd
from lifelines import KaplanMeierFitter
df = pd.read_csv("survival.csv"); kmf = KaplanMeierFitter()
fig, ax = plt.subplots(figsize=(6, 5))
for i, g in enumerate(sorted(df["group"].unique())):
    m = df["group"] == g
    kmf.fit(df.loc[m, "time"], df.loc[m, "event"], label=str(g))
    kmf.plot_survival_function(ax=ax, color=PALETTE[i % len(PALETTE)], ci_show=True)
ax.set_xlabel("Time"); ax.set_ylabel("Survival probability"); ax.set_ylim(0, 1.02)
ax.set_title("Kaplan–Meier"); fig.tight_layout(); fig.savefig("figures/km.png")

Manhattan (GWAS summary stats; columns CHR, BP, P):

pythonCopy code
import numpy as np, pandas as pd
df = pd.read_csv("gwas.csv").dropna(subset=["CHR", "BP", "P"]).copy()
df["CHR"] = df["CHR"].astype(str)
chrom_order = sorted(df["CHR"].unique(), key=lambda c: (len(c), c))
df["_nlp"] = -np.log10(np.clip(df["P"].astype(float), 1e-300, None))
offset, xticks, ticklabels, xpos = 0.0, [], [], np.zeros(len(df))  # lay chromosomes end-to-end
for c in chrom_order:
    m = (df["CHR"] == c).to_numpy()
    bp = df.loc[m, "BP"].astype(float)
    xpos[m] = bp + offset
    xticks.append(offset + bp.median()); ticklabels.append(c)
    offset = xpos[m].max() + 1
fig, ax = plt.subplots(figsize=(180/25.4, 70/25.4), constrained_layout=True)
for i, c in enumerate(chrom_order):
    m = (df["CHR"] == c).to_numpy()
    ax.scatter(xpos[m], df["_nlp"].to_numpy()[m], s=6, c=["#2a78d6", "#9fb3c8"][i % 2], edgecolors="none")
ax.axhline(-np.log10(5e-8), color=UP, lw=0.8, ls="--")   # genome-wide significance
ax.set_xticks(xticks); ax.set_xticklabels(ticklabels, fontsize=6)
ax.set_xlabel("Chromosome"); ax.set_ylabel(r"$-\log_{10} p$"); ax.set_title("Manhattan")
fig.savefig("figures/manhattan.png", dpi=400)

QQ plot (p-value calibration; reports genomic inflation λ):

pythonCopy code
import numpy as np, pandas as pd
from scipy import stats
p = pd.read_csv("gwas.csv")["P"].dropna().astype(float).to_numpy()
p = np.clip(np.sort(p), 1e-300, 1.0)
obs = -np.log10(p)
exp = -np.log10((np.arange(1, len(p) + 1) - 0.5) / len(p))
lam = np.median(stats.chi2.isf(p, 1)) / stats.chi2.isf(0.5, 1)   # genomic inflation
fig, ax = plt.subplots(figsize=(70/25.4, 70/25.4), constrained_layout=True)
ax.scatter(exp, obs, s=6, c="#2a78d6", edgecolors="none")
lim = float(max(exp.max(), obs.max()))
ax.plot([0, lim], [0, lim], color="0.4", lw=0.8, ls="--")
ax.set_xlabel(r"Expected $-\log_{10} p$"); ax.set_ylabel(r"Observed $-\log_{10} p$")
ax.set_title(f"QQ (λ = {lam:.3f})"); fig.savefig("figures/qq.png", dpi=400)

Forest (effect size ± 95% CI per study/variant; columns label, estimate, ci_low, ci_high):

pythonCopy code
import pandas as pd
df = pd.read_csv("effects.csv").iloc[::-1].reset_index(drop=True)  # first row at top
y = list(range(len(df)))
fig, ax = plt.subplots(figsize=(90/25.4, max(60, 14 * len(df)) / 25.4), constrained_layout=True)
ax.errorbar(df["estimate"], y,
            xerr=[df["estimate"] - df["ci_low"], df["ci_high"] - df["estimate"]],
            fmt="s", color="#204897", ecolor="0.4", capsize=2, markersize=4, lw=0.8)
ax.axvline(1.0, color=UP, lw=0.8, ls="--")   # null: OR/HR=1 (use 0.0 for log-scale beta)
ax.set_yticks(y); ax.set_yticklabels(df["label"])
ax.set_xlabel("Effect size (95% CI)"); ax.set_title("Forest")
fig.savefig("figures/forest.png", dpi=400)

Best Practices

• Only plot data that exists. Never invent columns, groups, or values; if a needed column is missing, inspect the header and ask or adapt — do not fabricate.
• Workspace-relative paths only. Save under figures/ (create it if needed); never write to absolute paths like /tmp or /home/....
• Reuse the palette across a figure set so related panels share colors for the same group / direction. Use the diverging colormap (centered at 0) for z-scores and log2FC; the sequential colormap for magnitudes and -log10 p. However, if the user explicitly requests a different color, use it.
• Label axes and give a real title. Include units, group n, and thresholds where relevant (e.g. volcano cutoff lines).
• For vector / print output, also save a .pdf (fig.savefig("figures/x.pdf")); text stays editable because pdf.fonttype=42 / svg.fonttype="none".

Common Pitfalls

• Wrong colormap family for the data. A sequential map on signed values (log2FC, z-score) hides the sign. How to avoid: use the diverging colormap with center=0, vmin=-vmax for signed data, and reserve the sequential colormap for unsigned magnitudes only.
• Columns don't match the recipe (KeyError). The recipe's default column names rarely match the real table. How to avoid: read the header first and override the *_col / x / y variables to the real names; never invent a missing column — pick a figure the data supports instead.
• Overcrowded / overlapping labels. Annotating every volcano point, or cramming long pathway names on an x-axis, produces unreadable overlap. How to avoid: label only the top ~10 by |log2FC| × -log10 p (≤5 inside a composite panel) and repel with adjustText; move long category names to a horizontal y-axis, or rotate 90° and wrap to ≤26 chars; skip labels rather than dumping colliding text.
• Dot / bubble markers all one size or too small. A size encoding that maps to invisible dots carries no information. How to avoid: floor the size mapping (s in ~[25, 220]) so the smallest stays visible, and confirm the size column actually varies.
• Legend or clustermap title collides. A legend sits on top of the data, or a clustermap title is cropped. How to avoid: move the legend to a free corner or outside the axes (or shrink its font); put a clustermap title on g.ax_col_dendrogram, not a figure suptitle.
• Empty / all-NS plot. Nothing passes the threshold because the columns or scale are wrong. How to avoid: confirm the p-value and fold-change columns are numeric and that thresholds match the data scale (adjusted vs raw p).

For multi-panel composition issues (panel letters misplaced, legends floating in margins, panels colliding), see the multipanel skill.

Further Reading

• Matplotlib documentation — https://matplotlib.org/stable/
• Seaborn documentation — https://seaborn.pydata.org/
• Wong, "Points of view: Color blindness", Nature Methods (2011) — https://www.nature.com/articles/nmeth.1618
• ColorBrewer diverging/sequential palettes — https://colorbrewer2.org/