19. Putting it together — from raw counts to one finished figure · Dr. Seonggu Sohn's Data Visualization, First Steps

The last exercise. Starting from the raw counts the paper published, you build the equivalent of one figure from beginning to end.

1️⃣

Read

counts, 24,419 × 7

2️⃣

Normalise

CPM → log2. The log axis of lesson 5

3️⃣

Test

WT 4 versus KD 3, a t-test

4️⃣

Draw

volcano + heatmap + dots. Lessons 9, 10, 14

Library sizes differ between samples. Skip this and your conclusion becomes 'genes are more expressed in deeply sequenced samples'.

Python Runs in your browser

# Step 1 - read the data and look at its size
import viz

c = viz.counts()
print("raw counts :", c.shape[0], "genes x", c.shape[1], "samples")
print("samples    :", ", ".join(c.columns))
print()
print("Each sample was sequenced to a different depth - which is why we normalise")
for s in c.columns:
    print("   %-12s %12s" % (s, format(int(c[s].sum()), ",")))

Most values sit near zero. Expression is a quantity that moves by multiplication, so taking a log is the standard move.

Python Runs in your browser

# Step 2 - equalise the totals with CPM, then take log2
import viz, numpy as np

lg = viz.logcpm()
print("genes left after dropping the barely detected :", len(lg))
print()
print("distribution of log2(CPM+1)")
v = lg.values.ravel()
for q in [0, 25, 50, 75, 100]:
    print("   %3d pct : %6.2f" % (q, np.percentile(v, q)))

These numbers are the raw material for the figure. Look at the numbers before you draw — that is the order.

Python Runs in your browser

# Step 3 - differential expression, WT versus KD
import viz

d = viz.de()
up = ((d.log2FC >= 1) & (d.p < .05)).sum()
dn = ((d.log2FC <= -1) & (d.p < .05)).sum()
print("genes tested : %d" % len(d))
print("up in KD   (>= 2x, p<0.05) : %d" % up)
print("down in KD                 : %d" % dn)
print()
print("the ten that moved most")
top = d.reindex((d.log2FC.abs() * d.mlog10p).sort_values(ascending=False).index).head(10)
print(top[["log2FC", "p"]].round(3).to_string())

Step through the six stages. Watch how the chart reads differently as each element is added.

Loading…

This is the skeleton of a paper figure — show the whole (A), show the shortlist (B), show the few that are the point (C). Most figures are built from those three.

Python Runs in your browser

# Step 4 - one figure with three panels
import matplotlib.pyplot as plt
import viz, numpy as np

d  = viz.de()
lg = viz.logcpm()
wt = [c for c in lg.columns if c.startswith("WT")]
kd = [c for c in lg.columns if c.startswith("KD")]

fig = plt.figure(figsize=(12, 4.2))
gs  = fig.add_gridspec(1, 3, width_ratios=[1.25, 1, .85], wspace=.42)

# (A) volcano
axA = fig.add_subplot(gs[0])
up = (d.log2FC >= 1) & (d.p < .05)
dn = (d.log2FC <= -1) & (d.p < .05)
axA.scatter(d.log2FC[~(up|dn)], d.mlog10p[~(up|dn)], s=3, color=viz.GREY, alpha=.35)
axA.scatter(d.log2FC[up], d.mlog10p[up], s=7, color=viz.PAL[1])
axA.scatter(d.log2FC[dn], d.mlog10p[dn], s=7, color=viz.PAL[0])
for g in ["Ddx54", "Dicer1", "Cd47", "Myc", "Cd38"]:
    if g in d.index:
        r = d.loc[g]
        axA.annotate(g, (r.log2FC, r.mlog10p), fontsize=7.5,
                     xytext=(4, 2), textcoords="offset points")
axA.axhline(-np.log10(.05), ls="--", lw=.9, color="#aaa")
for x in (-1, 1):
    axA.axvline(x, ls="--", lw=.9, color="#aaa")
axA.set_xlabel("log2 FC (KD - WT)", fontsize=8.5)
axA.set_ylabel("-log10(p)", fontsize=8.5)
axA.set_title("A   differential expression", fontsize=9.5, loc="left")
viz.bare(axA)

# (B) heatmap
axB = fig.add_subplot(gs[1])
sel = (d.log2FC.abs() * d.mlog10p).sort_values(ascending=False).index[:20]
Z = viz.zscore(lg.loc[sel].values)
im = axB.imshow(Z, cmap="RdBu_r", vmin=-1.8, vmax=1.8, aspect="auto")
axB.set_yticks(range(len(sel))); axB.set_yticklabels(sel, fontsize=6)
axB.set_xticks(range(lg.shape[1]))
axB.set_xticklabels([s.replace("_DDX54", "") for s in lg.columns], rotation=90, fontsize=6)
axB.set_title("B   top 20, row z-score", fontsize=9.5, loc="left")
fig.colorbar(im, ax=axB, shrink=.7, label="z")

# (C) dot panel
axC = fig.add_subplot(gs[2])
show = [g for g in ["Dicer1", "Drosha", "Myc", "Cd47", "Cd38"] if g in lg.index]
for i, g in enumerate(show):
    for k, cols in enumerate([wt, kd]):
        v = lg.loc[g, cols].values
        axC.plot([i + (k-.5)*.26]*len(v), v, "o", ms=4, color=viz.PAL[k], alpha=.85)
        axC.plot([i + (k-.5)*.26 - .08, i + (k-.5)*.26 + .08], [v.mean()]*2,
                 lw=2, color=viz.PAL[k])
axC.set_xticks(range(len(show))); axC.set_xticklabels(show, rotation=45, fontsize=7.5)
axC.set_ylabel("log2(CPM+1)", fontsize=8.5)
axC.set_title("C   WT (blue) vs KD (orange)", fontsize=9.5, loc="left")
viz.bare(axC)

🤔 Check yourself

What is the A (volcano) → B (heatmap) → C (dot panel) arrangement doing?

🤔 Check yourself

Several panels of one figure describe the same conditions. What must be consistent?