生信图表大全（第五篇）：蛋白结构、变异与表达验证的 40 张图

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Line3DCollection

def parse(path):  # PDB 定宽列解析：只取需要的字段
    atoms = []
    for ln in open(path, encoding="ascii", errors="ignore"):
        if ln[:6] in ("ATOM  ", "HETATM"):
            atoms.append(dict(het=ln[:6] == "HETATM",
                              name=ln[12:16].strip(), resn=ln[17:20].strip(),
                              chain=ln[21], resi=int(ln[22:26]),
                              xyz=(float(ln[30:38]), float(ln[38:46]),
                                   float(ln[46:54])),
                              elem=ln[76:78].strip() or ln[12:16].strip()[0]))
    return atoms

atoms = parse("6a15.pdb")
P = np.array([a["xyz"] for a in atoms
              if a["chain"] == "A" and not a["het"] and a["name"] == "CA"])
fig = plt.figure(figsize=(8.2, 6.6))
ax = fig.add_subplot(projection="3d")
ax.view_init(elev=12, azim=-70)
ax.set_axis_off()
ax.set_box_aspect((1, 1, 1))
pts = np.concatenate([P[:-1, None, :], P[1:, None, :]], axis=1)
cmap = matplotlib.colormaps["viridis"]
ax.add_collection3d(Line3DCollection(
    pts, colors=[cmap(i / (len(P) - 2)) for i in range(len(P) - 1)],
    linewidths=2.2))
ax.text(*P[0], " N", color="0.35", fontsize=9)
ax.text(*P[-1], " C", color="0.35", fontsize=9)
m, M = P.min(0), P.max(0)
ctr, span = (m + M) / 2, (M - m).max() / 2 * 1.05
ax.set_xlim(ctr[0] - span, ctr[0] + span)
ax.set_ylim(ctr[1] - span, ctr[1] + span)
ax.set_zlim(ctr[2] - span, ctr[2] + span)
ax.set_title("C-alpha trace of CYP90B1 (%d residues, PDB 6A15, real data)"
             % len(P))
fig.savefig("134-ca-trace.png", dpi=200, bbox_inches="tight")

ELEMCOL = {"C": "#9aa3ad", "N": "#3369c3", "O": "#dd4b39",
           "S": "#e8b636", "FE": "#e07b39"}

def segs_by_res(atoms, lw=0.6):
    """残基内按键序连线 + 肽键 C->N，返回 Line3DCollection。"""
    from mpl_toolkits.mplot3d.art3d import Line3DCollection
    segs, cols, byres, prev_c = [], [], {}, None
    for a in atoms:
        byres.setdefault((a["chain"], a["resi"]), []).append(a)
    for key in sorted(byres):
        lst = byres[key]
        for i in range(len(lst) - 1):
            segs.append([lst[i]["xyz"], lst[i + 1]["xyz"]])
            cols.append(ELEMCOL.get(lst[i]["elem"], "#9aa3ad"))
        if prev_c is not None and lst[0]["name"] == "N":
            segs.append([prev_c, lst[0]["xyz"]])
            cols.append("#c3c3c3")
        for a in lst:
            if a["name"] == "C":
                prev_c = a["xyz"]
    return Line3DCollection(segs, colors=cols, linewidths=lw)

pa = [a for a in parse("6a15.pdb") if a["chain"] == "A" and not a["het"]]
fig = plt.figure(figsize=(8.2, 6.6))
ax = fig.add_subplot(projection="3d")
ax.view_init(elev=12, azim=-70)
ax.set_axis_off()
ax.set_box_aspect((1, 1, 1))
ax.add_collection3d(segs_by_res(pa, 0.6))
P = np.array([a["xyz"] for a in pa])
m, M = P.min(0), P.max(0)
ctr, span = (m + M) / 2, (M - m).max() / 2 * 1.05
ax.set_xlim(ctr[0] - span, ctr[0] + span)
ax.set_ylim(ctr[1] - span, ctr[1] + span)
ax.set_zlim(ctr[2] - span, ctr[2] + span)
ax.set_title("Wireframe of CYP90B1 (PDB 6A15, real data)")
fig.savefig("135-wireframe.png", dpi=200, bbox_inches="tight")

residue,dist_to_ligand_A
PHE84,3.85
TYR112,3.60
PRO113,4.54
ILE116,4.37
LEU120,4.57
SER124,4.96
MET125,4.39
LEU126,3.91
MET213,3.93
VAL216,3.55
VAL217,3.74
SER307,3.64
LEU308,3.91
PHE310,3.70
ALA311,4.14
GLU314,4.92
THR315,3.79
VAL381,3.87
PHE383,3.17
LEU384,4.03
HIS385,2.88
PRO497,4.30
PHE498,4.71

import numpy as np

atoms = parse("6a15.pdb")
lig = [a for a in atoms if a["het"] and a["resn"] == "CLR"]
heme = [a for a in atoms if a["het"] and a["resn"] == "HEM"]
L = np.array([a["xyz"] for a in lig])
byres = {}
for a in atoms:
    if a["chain"] == "A" and not a["het"]:
        byres.setdefault(a["resi"], []).append(a)
pocket = []
for ri, lst in byres.items():
    X = np.array([a["xyz"] for a in lst])
    d = np.linalg.norm(X[:, None, :] - L[None, :, :], axis=-1).min()
    if d <= 5.0:
        pocket.append((ri, lst[0]["resn"], d))
# 渲染：pocket 残基棍棒 lw=2.4、CLR 橙 lw=3、HEM 红 lw=3，
# 每个 CA 打一个 "RES123" 标签；视角 azim=-96 elev=16。
# 完整渲染代码见仓库 tools/bio-plots/ch11a_struct.py 的 pocket_stick()。
print("pocket residues:", len(pocket))
print("nearest:", sorted(pocket, key=lambda t: t[2])[:3])

import numpy as np

RAD = {"C": 1.70, "N": 1.55, "O": 1.52, "S": 1.80, "FE": 1.45, "P": 1.80}
atoms = [a for a in parse("6a15.pdb")
         if not a["het"] or a["resn"] in ("HEM", "CLR")]
P = np.array([a["xyz"] for a in atoms])
r = np.array([RAD.get(a["elem"], 1.7) for a in atoms])
c = [ELEMCOL.get(a["elem"], "#9aa3ad") for a in atoms]
fig = plt.figure(figsize=(8.4, 6.8))
ax = fig.add_subplot(projection="3d")
ax.view_init(elev=12, azim=-70)
ax.set_axis_off()
ax.set_box_aspect((1, 1, 1))
ax.scatter(P[:, 0], P[:, 1], P[:, 2], s=(r * 26) ** 1.6 / 6, c=c,
           linewidths=0)
m, M = P.min(0), P.max(0)
ctr, span = (m + M) / 2, (M - m).max() / 2 * 1.08
ax.set_xlim(ctr[0] - span, ctr[0] + span)
ax.set_ylim(ctr[1] - span, ctr[1] + span)
ax.set_zlim(ctr[2] - span, ctr[2] + span)
ax.set_title("Space-filling model of CYP90B1 (PDB 6A15, real data)")
fig.savefig("137-spacefill.png", dpi=200, bbox_inches="tight")

from scipy.interpolate import make_interp_spline
from mpl_toolkits.mplot3d.art3d import Line3DCollection

ca = [a for a in parse("6a15.pdb")
      if a["chain"] == "A" and not a["het"] and a["name"] == "CA"]
P = np.array([a["xyz"] for a in ca])
t = np.zeros(len(P))
t[1:] = np.cumsum(np.linalg.norm(np.diff(P, axis=0), axis=1))
S = make_interp_spline(t, P, k=3)(np.linspace(t[0], t[-1], 800))
fig = plt.figure(figsize=(8.2, 6.6))
ax = fig.add_subplot(projection="3d")
ax.view_init(elev=12, azim=-70)
ax.set_axis_off()
ax.set_box_aspect((1, 1, 1))
pts = np.concatenate([S[:-1, None, :], S[1:, None, :]], axis=1)
cmap = matplotlib.colormaps["cool"]
ax.add_collection3d(Line3DCollection(
    pts, colors=[cmap(i / (len(S) - 2)) for i in range(len(S) - 1)],
    linewidths=6.5))
m, M = P.min(0), P.max(0)
ctr, span = (m + M) / 2, (M - m).max() / 2 * 1.05
ax.set_xlim(ctr[0] - span, ctr[0] + span)
ax.set_ylim(ctr[1] - span, ctr[1] + span)
ax.set_zlim(ctr[2] - span, ctr[2] + span)
ax.set_title("Simplified ribbon of CYP90B1 (PDB 6A15, real data)")
fig.savefig("138-ribbon.png", dpi=200, bbox_inches="tight")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

N, rng = 520, np.random.default_rng(5143)
i = np.arange(N)[:, None]
j = np.arange(N)[None, :]
sep = np.abs(i - j)
base = 0.8 + 4.0 * (sep / N) ** 1.7
loop = ((i > 300) & (i < 322)) | ((j > 300) & (j < 322))
M = np.clip(base + rng.normal(0, 0.3, (N, N))
            + loop * np.where(sep > 25, 2.6, 0), 0.3, 9.0)
fig, ax = plt.subplots(figsize=(6.6, 5.6))
im = ax.imshow(M, cmap="viridis", origin="lower", vmin=0, vmax=9,
               extent=[1, N, 1, N])
ax.plot([1, N], [1, N], color="white", lw=0.6, alpha=0.7)
ax.set_xlabel("scored residue")
ax.set_ylabel("aligned residue")
fig.colorbar(im, ax=ax, shrink=0.85, label="PAE (A)")
ax.set_title("Predicted aligned error (simulated)")
fig.savefig("143-pae-map.png", dpi=200, bbox_inches="tight")

import numpy as np

ca = [a for a in parse("6a15.pdb")
      if a["chain"] == "A" and not a["het"] and a["name"] == "CA"]
P = np.array([a["xyz"] for a in ca])
D = np.linalg.norm(P[:, None, :] - P[None, :, :], axis=-1)
fig, ax = plt.subplots(figsize=(6.9, 6.2))
im = ax.imshow(D, cmap="viridis_r", vmin=4, vmax=45)
ax.contour(np.arange(len(P)), np.arange(len(P)),
           (D < 8).astype(float), levels=[0.5], colors="black",
           linewidths=0.3)
ax.set_xlabel("residue i")
ax.set_ylabel("residue j")
fig.colorbar(im, ax=ax, shrink=0.85, label="distance (A)")
ax.set_title("C-alpha distance map (PDB 6A15, real data)")
fig.savefig("144-contact-map.png", dpi=200, bbox_inches="tight")

import numpy as np

def kabsch(P, Q):  # 把 P 旋转平移到 Q 上（最小二乘）
    pc, qc = P.mean(0), Q.mean(0)
    U, _, Vt = np.linalg.svd((P - pc).T @ (Q - qc))
    d = np.sign(np.linalg.det(Vt.T @ U.T))
    R = Vt.T @ np.diag([1, 1, d]) @ U.T
    return (P - pc) @ R + qc

def cas(path, chain="A"):
    out = []
    for ln in open(path, encoding="ascii", errors="ignore"):
        if ln[:6] == "ATOM  " and ln[21] == chain and ln[12:16].strip() == "CA":
            out.append((int(ln[22:26]),
                        (float(ln[30:38]), float(ln[38:46]),
                         float(ln[46:54]))))
    return dict(out)

w, m = cas("1bi5.pdb"), cas("1i89.pdb")
common = sorted(set(w) & set(m))
Q = np.array([w[i] for i in common])          # WT
Pt = kabsch(np.array([m[i] for i in common]), Q)
dev = np.linalg.norm(Pt - Q, axis=1)
print("common=%d global RMSD=%.2f G256=%.2f max=%.2f@%d"
      % (len(common), np.sqrt((dev ** 2).mean()),
         dev[common.index(256)], dev.max(), common[dev.argmax()]))

position,score,class
96,0.18,neutral
105,0.35,hypomorphic
124,0.43,hypomorphic
158,0.14,neutral
184,0.81,loss_of_function
198,0.22,neutral
200,0.49,hypomorphic
256,0.93,loss_of_function
258,0.10,neutral
276,0.28,neutral
295,0.27,neutral
300,0.18,neutral
311,0.19,neutral
350,0.35,hypomorphic
385,0.86,loss_of_function
408,0.22,neutral
432,0.57,loss_of_function
448,0.44,hypomorphic
458,0.44,hypomorphic
471,0.52,loss_of_function
474,0.24,neutral
483,0.24,neutral

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

rows = [ln.split(",") for ln in open("146-lollipop.csv").read().splitlines()[1:]]
pos = np.array([int(r[0]) for r in rows])
sc = np.array([float(r[1]) for r in rows])
cols = np.where(sc > 0.5, "#C44E52", np.where(sc > 0.3, "#DD8452", "#9AA0A6"))
fig, ax = plt.subplots(figsize=(9.4, 4.4))
ax.vlines(pos, 0, sc, color=cols, lw=1.6)
ax.scatter(pos, sc, s=44, c=cols, zorder=3)
ax.set_xlabel("residue position")
ax.set_ylabel("deleteriousness score")
ax.set_ylim(0, 1.05)
fig.savefig("146-mutation-lollipop.png", dpi=200, bbox_inches="tight")

from\to,A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y
A,0.0,1.8,2.2,2.4,0.4,1.1,2.3,1.1,2.5,0.9,0.0,2.4,2.9,2.3,2.8,1.1,1.0,1.0,1.4,1.4
C,1.8,0.0,4.2,4.0,1.3,2.4,3.9,2.1,4.3,1.5,1.3,4.0,4.6,3.7,4.6,2.5,2.7,1.9,3.1,2.9
D,2.2,4.2,0.0,0.0,2.6,1.3,2.5,3.9,2.8,3.0,2.1,0.3,2.3,0.1,2.7,1.4,1.3,3.6,1.0,1.1
E,2.4,4.0,0.0,0.0,3.1,1.3,2.3,3.4,2.0,3.1,2.4,0.0,2.3,0.0,2.8,1.5,1.2,3.5,1.2,1.0
F,0.4,1.3,2.6,3.1,0.0,1.6,2.6,0.9,3.3,0.5,0.5,3.0,3.5,2.8,3.2,1.4,1.5,0.6,1.7,1.4
G,1.1,2.4,1.3,1.3,1.6,0.0,1.1,2.4,1.6,2.0,0.9,1.3,2.1,1.3,2.1,0.4,0.0,1.9,0.4,0.7
H,2.3,3.9,2.5,2.3,2.6,1.1,0.0,3.6,0.1,3.1,2.4,0.0,2.3,0.2,0.8,0.9,1.1,3.6,1.0,0.8
I,1.1,2.1,3.9,3.4,0.9,2.4,3.6,0.0,3.9,0.3,0.9,3.3,4.1,3.7,3.9,2.5,2.2,0.3,2.3,2.9
K,2.5,4.3,2.8,2.0,3.3,1.6,0.1,3.9,0.0,3.2,2.6,0.1,2.6,0.0,0.3,1.4,1.4,3.4,1.3,1.1
L,0.9,1.5,3.0,3.1,0.5,2.0,3.1,0.3,3.2,0.0,1.1,3.4,3.9,3.1,3.8,2.3,2.0,0.4,1.9,2.5
M,0.0,1.3,2.1,2.4,0.5,0.9,2.4,0.9,2.6,1.1,0.0,2.6,3.0,2.2,2.8,1.0,1.5,1.1,1.2,1.7
N,2.4,4.0,0.3,0.0,3.0,1.3,0.0,3.3,0.1,3.4,2.6,0.0,2.3,0.1,0.5,1.4,1.1,3.5,1.3,0.9
P,2.9,4.6,2.3,2.3,3.5,2.1,2.3,4.1,2.6,3.9,3.0,2.3,0.0,2.4,2.7,1.9,2.2,4.0,2.1,1.4
Q,2.3,3.7,0.1,0.0,2.8,1.3,0.2,3.7,0.0,3.1,2.2,0.1,2.4,0.0,0.4,1.1,1.3,3.3,1.5,1.1
R,2.8,4.6,2.7,2.8,3.2,2.1,0.8,3.9,0.3,3.8,2.8,0.5,2.7,0.4,0.0,1.5,1.7,4.1,1.6,1.7
S,1.1,2.5,1.4,1.5,1.4,0.4,0.9,2.5,1.4,2.3,1.0,1.4,1.9,1.1,1.5,0.0,0.3,2.3,0.2,0.1
T,1.0,2.7,1.3,1.2,1.5,0.0,1.1,2.2,1.4,2.0,1.5,1.1,2.2,1.3,1.7,0.3,0.0,2.4,0.0,0.4
V,1.0,1.9,3.6,3.5,0.6,1.9,3.6,0.3,3.4,0.4,1.1,3.5,4.0,3.3,4.1,2.3,2.4,0.0,2.3,2.7
W,1.4,3.1,1.0,1.2,1.7,0.4,1.0,2.3,1.3,1.9,1.2,1.3,2.1,1.5,1.6,0.2,0.0,2.3,0.0,0.0
Y,1.4,2.9,1.1,1.0,1.4,0.7,0.8,2.9,1.1,2.5,1.7,0.9,1.4,1.1,1.7,0.1,0.4,2.7,0.0,0.0

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

AA = list("ACDEFGHIKLMNPQRSTVWY")
kd = dict(zip(AA, [1.8, 2.5, -3.5, -3.5, 2.8, -0.4, -3.2, 4.5, -3.9,
                   3.8, 1.9, -3.5, -1.6, -3.5, -4.5, -0.8, -0.7, 4.2,
                   -0.9, -1.3]))                      # Kyte-Doolittle
rng = np.random.default_rng(5147)
M = np.zeros((20, 20))
for a in range(20):
    for b in range(20):
        if a == b:
            continue
        charge = 2.2 if (("KRH".find(AA[a]) >= 0 and "DE".find(AA[b]) >= 0)
                         or ("DE".find(AA[a]) >= 0 and "KRH".find(AA[b]) >= 0)) else 0
        M[a, b] = 0.45 * abs(kd[AA[a]] - kd[AA[b]]) + charge \
            + 1.5 * ("P" in (AA[a], AA[b])) + 1.2 * ("C" in (AA[a], AA[b])) \
            + rng.normal(0, 0.25)
M = np.clip((M + M.T) / 2, 0, None)
np.fill_diagonal(M, 0.0)
fig, ax = plt.subplots(figsize=(7.0, 6.0))
im = ax.imshow(M, cmap="magma_r")
ax.set_xticks(range(20), AA)
ax.set_yticks(range(20), AA)
fig.colorbar(im, ax=ax, shrink=0.85, label="mean ddG (kcal/mol)")
fig.savefig("147-substitution-matrix.png", dpi=200, bbox_inches="tight")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

rows = [ln.split(",") for ln in open("148-ddg.csv").read().splitlines()[1:]]
mut = [r[0] for r in rows]
v = np.array([float(r[1]) for r in rows])
e = np.array([float(r[2]) for r in rows])
fig, ax = plt.subplots(figsize=(8.6, 4.4))
ax.bar(range(len(mut)), v, yerr=e, capsize=3, width=0.62,
       color=np.where(v > 2.0, "#C44E52", "#4C72B0"))
ax.axhline(2.0, ls="--", color="0.35", lw=1)
ax.set_xticks(range(len(mut)), mut, rotation=45, ha="right", fontsize=8.5)
ax.set_ylabel("ddG (kcal/mol)")
fig.savefig("148-ddg-bars.png", dpi=200, bbox_inches="tight")

start,end,type,label
34,46,helix,A
60,72,helix,B
90,102,helix,B'
110,116,sheet,b1
130,136,helix,C
148,160,helix,D
190,205,helix,E
210,216,sheet,b2
220,236,helix,F
250,262,helix,G
292,330,helix,I
342,352,helix,J
360,368,sheet,b3
380,395,helix,K'
420,432,helix,meander
452,462,helix,heme loop
500,508,sheet,b4

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

N = 520
x = np.arange(1, N + 1)
dis = np.clip(0.12 + 0.5 * np.exp(-((x - 40) / 45.0) ** 2)
              + 0.45 * np.exp(-((x - 311) / 16.0) ** 2)
              + 0.55 * np.exp(-((x - 516) / 26.0) ** 2), 0.01, 0.98)
fig, (ax, ax2) = plt.subplots(2, 1, figsize=(9.6, 4.6), sharex=True)
for s, e, c in [(34, 46, "#4C72B0"), (60, 72, "#4C72B0"),
                (90, 102, "#4C72B0"), (292, 330, "#4C72B0"),
                (250, 262, "#4C72B0"), (380, 395, "#4C72B0"),
                (110, 116, "#DD8452"), (210, 216, "#DD8452")]:
    ax.add_patch(Rectangle((s, 0.32), e - s + 1, 0.36, color=c, lw=0))
for k, (s, e) in enumerate([(96, 106), (208, 218), (248, 262),
                            (316, 326), (388, 398)]):
    ax.add_patch(Rectangle((s, 0.78), e - s + 1, 0.16, color="#55A868", lw=0))
    ax.text((s + e) / 2, 0.70, "SRS%d" % (k + 1), ha="center",
            fontsize=7.5, color="#2d6a4f")
for p, lab in [(256, "G256"), (385, "H385"), (462, "heme Cys")]:
    ax.axvline(p, color="0.3", lw=0.8, ls=":")
ax.set_yticks([])
ax.set_ylim(-0.12, 1.02)
ax2.plot(x, dis, color="#8172B3", lw=1.1)
ax2.set_ylabel("disorder")
ax2.set_xlabel("residue")
fig.savefig("149-ss-track.png", dpi=200, bbox_inches="tight")

import numpy as np

atoms = parse("6a15.pdb")
lig = np.array([a["xyz"] for a in atoms
                if a["het"] and a["resn"] == "CLR"])
heme = np.array([a["xyz"] for a in atoms
                 if a["het"] and a["resn"] == "HEM"])
pa = [a for a in atoms if a["chain"] == "A" and not a["het"]]
byres = {}
for a in pa:
    byres.setdefault(a["resi"], []).append(a)
pk = []
for ri, lst in byres.items():
    X = np.array([a["xyz"] for a in lst])
    if np.linalg.norm(X[:, None, :] - lig[None, :, :],
                      axis=-1).min() <= 5.0:
        pk += lst
P = np.array([a["xyz"] for a in pa if a["name"] == "CA"])
X = np.array([a["xyz"] for a in pk])
fig = plt.figure(figsize=(8.4, 6.8))
ax = fig.add_subplot(projection="3d")
ax.view_init(elev=16, azim=-96)
ax.set_axis_off()
ax.set_box_aspect((1, 1, 1))
pts = np.concatenate([P[:-1, None, :], P[1:, None, :]], axis=1)
ax.add_collection3d(Line3DCollection(pts, colors="0.78", linewidths=1.2))
ax.scatter(X[:, 0], X[:, 1], X[:, 2], s=26, c="#C44E52", linewidths=0)
ax.scatter(lig[:, 0], lig[:, 1], lig[:, 2], s=64, c="#e08a2e", linewidths=0)
ax.scatter(heme[:, 0], heme[:, 1], heme[:, 2], s=40, c="#8172B3",
           linewidths=0)
fig.savefig("150-pocket-map.png", dpi=200, bbox_inches="tight")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

rows = [ln.split(",") for ln in open("151-pocket-volume.csv").read().splitlines()[1:]]
labs = [r[0] for r in rows]
v = np.array([float(r[1]) for r in rows])
cols = ["#4C72B0" if l == "WT" else ("#C44E52" if x < 700 else "#9AA0A6")
        for l, x in zip(labs, v)]
fig, ax = plt.subplots(figsize=(8.2, 4.4))
ax.bar(range(len(v)), v, color=cols, width=0.6)
ax.axhline(405, ls="--", color="#e08a2e", lw=1.4)
ax.set_xticks(range(len(v)), labs, rotation=30, ha="right")
ax.set_ylabel("pocket volume (A^3)")
ax.set_ylim(0, 1150)
fig.savefig("151-pocket-volume.png", dpi=200, bbox_inches="tight")

pocket_residue,conservation
84,0.88
112,0.94
113,0.93
116,0.92
120,0.90
124,0.92
125,0.89
126,0.95
213,0.90
216,0.89
217,0.93
307,0.93
308,0.92
310,0.92
311,0.89
314,0.88
315,0.91
381,0.90
383,0.91
384,0.91
385,0.90
497,0.92
498,0.89

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

rng = np.random.default_rng(5152)
N = 520
POCKET = [84, 112, 113, 116, 120, 124, 125, 126, 213, 216, 217, 307,
          308, 310, 311, 314, 315, 381, 383, 384, 385, 497, 498]
x = np.arange(1, N + 1)
cons = np.clip(0.58 + rng.normal(0, 0.10, N), 0.25, 0.9)
cons[np.array(POCKET) - 1] = 0.90
cons[(x >= 448) & (x <= 468)] += 0.25
fig, ax = plt.subplots(figsize=(9.6, 4.2))
ax.plot(x, cons, color="#4C72B0", lw=1.0)
ax.scatter(POCKET, cons[np.array(POCKET) - 1], s=30, color="#C44E52",
           zorder=3)
ax.axvspan(448, 468, color="#8172B3", alpha=0.14)
ax.set_xlabel("residue")
ax.set_ylabel("conservation (0-1)")
fig.savefig("152-conservation.png", dpi=200, bbox_inches="tight")

pose,rmsd_A,score_kcal_mol
1,0.83,-9.64
2,0.95,-9.83
3,1.09,-9.55
4,1.19,-9.83
5,1.46,-9.54
6,0.87,-9.33
7,0.91,-9.61
8,0.87,-9.59
9,1.03,-9.87
10,2.82,-6.43
11,7.09,-6.41
12,3.26,-6.12
13,6.77,-6.84
14,4.35,-7.24
15,3.59,-7.12
16,4.59,-7.30
17,5.92,-5.99
18,5.17,-6.39
19,5.84,-6.33
20,6.80,-7.33
21,4.25,-7.85
22,3.71,-7.09
23,5.07,-6.48
24,5.43,-6.31
25,5.96,-6.90
26,3.94,-6.73
27,4.53,-6.97
28,5.83,-7.40
29,3.99,-6.76
30,4.08,-7.53

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

P = np.genfromtxt("153-docking-poses.csv", delimiter=",", names=True)
C = np.genfromtxt("153-docking-compounds.csv", delimiter=",", names=True,
                  dtype=None, encoding="utf8")
fig, (ax, ax2) = plt.subplots(1, 2, figsize=(10.6, 4.4))
cols = np.where(P["rmsd_A"] < 2.5, "#C44E52", "#9AA0A6")
ax.scatter(P["pose"], P["score_kcal_mol"], c=cols, s=46, zorder=3)
ax.axhspan(-11.5, -9.0, color="#C44E52", alpha=0.07)
ax.set_xlabel("pose rank")
ax.set_ylabel("docking score (kcal/mol)")
ax2.barh(range(len(C["compound"])), C["score_kcal_mol"], color="#4C72B0",
         height=0.62)
ax2.set_yticks(range(len(C["compound"])), C["compound"], fontsize=8.5)
ax2.axvline(-7.0, ls="--", color="0.4", lw=1)
ax2.set_xlabel("best score (kcal/mol)")
fig.savefig("153-docking.png", dpi=200, bbox_inches="tight")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

lines = open("154-interface.csv").read().splitlines()
B = lines[0].split(",")[1:]
A = [ln.split(",")[0] for ln in lines[1:]]
M = np.array([[float(v) for v in ln.split(",")[1:]] for ln in lines[1:]])
fig, ax = plt.subplots(figsize=(6.8, 5.8))
im = ax.imshow(M, cmap="YlOrRd", vmin=0, vmax=1)
ax.add_patch(Rectangle((3.5, 2.5), 3, 3, fill=False, ec="#2b5f8a", lw=1.8))
ax.set_xticks(range(len(B)), B)
ax.set_yticks(range(len(A)), A)
fig.colorbar(im, ax=ax, shrink=0.85, label="contact frequency")
fig.savefig("154-interface-contacts.png", dpi=200, bbox_inches="tight")

residue,sasa_A2
PHE84,23.3
TYR112,16.5
PRO113,15.6
ILE116,0.0
LEU120,0.0
SER124,0.0
MET125,0.0
LEU126,1.2
MET213,5.3
VAL216,4.8
VAL217,20.5
SER307,1.2
LEU308,3.9
PHE310,0.0
ALA311,0.0
GLU314,9.3
THR315,0.0
VAL381,6.1
PHE383,0.0
LEU384,0.0
HIS385,0.0
PRO497,6.3
PHE498,15.8
HIS516,107.7
HEM601,3.9
CLR602,10.1

import numpy as np
from scipy.spatial import cKDTree

RAD = {"C": 1.70, "N": 1.55, "O": 1.52, "S": 1.80, "FE": 1.45, "P": 1.80}
keep = [a for a in parse("6a15.pdb")
        if not a["het"] or a["resn"] in ("HEM", "CLR")]
X = np.array([a["xyz"] for a in keep])
r = np.array([RAD.get(a["elem"], 1.7) + 1.4 for a in keep])  # +水探针
n = 92                                    # Shrake-Rupley 采样点数
k = np.arange(n) + 0.5
phi = np.arccos(1 - 2 * k / n)            # 黄金螺旋均匀布点
theta = np.pi * (1 + 5 ** 0.5) * k
SP = np.c_[np.cos(theta) * np.sin(phi), np.sin(theta) * np.sin(phi),
           np.cos(phi)] * r[:, None, None]
pts = (SP + X[:, None, :]).reshape(-1, 3)
tree = cKDTree(X)
cnt = np.array([len(nb) for nb in
                tree.query_ball_point(pts, r.repeat(n), workers=-1)]
               ).reshape(len(keep), n)
frac = (cnt == 0).mean(1)                 # 没被邻居挡住的采样点比例
area = 4 * np.pi * r ** 2 * frac           # 该原子的 SASA
per = {}
for a, s in zip(keep, area):
    per[a["resi"]] = per.get(a["resi"], 0) + s   # 残基求和
print("CLR602 SASA = %.1f A^2" % per[602])

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

rows = [ln.split(",") for ln in open("156-qpcr-std.csv").read().splitlines()[1:]]
x = np.log10(np.array([float(r[0]) for r in rows]))
y = np.array([float(r[1]) for r in rows])
k, b = np.polyfit(x, y, 1)
r2 = 1 - ((y - (k * x + b)) ** 2).sum() / ((y - y.mean()) ** 2).sum()
eff = (10 ** (-1 / k) - 1) * 100
fig, ax = plt.subplots(figsize=(7.2, 4.8))
xx = np.linspace(1.7, 6.3, 50)
ax.plot(xx, k * xx + b, color="#4C72B0", lw=1.4)
ax.scatter(x, y, s=52, color="#C44E52", zorder=3)
ax.text(0.04, 0.10, "slope=%.2f\nintercept=%.1f\nR2=%.4f\nE=%.1f%%"
        % (k, b, r2, eff), transform=ax.transAxes, fontsize=9.5)
ax.set_xlabel("log10(template copies)")
ax.set_ylabel("Cq")
fig.savefig("156-qpcr-stdcurve.png", dpi=200, bbox_inches="tight")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

lanes = ["WT", "mut-1", "mut-2", "comp-1", "comp-2"]
target = np.array([1.00, 0.21, 0.33, 0.88, 0.74])
actin = np.array([0.97, 1.02, 0.95, 0.98, 1.01])
fig = plt.figure(figsize=(8.6, 5.6))
gs = fig.add_gridspec(2, 1, height_ratios=[2.6, 1.0], hspace=0.42)
axg = fig.add_subplot(gs[0])
axg.add_patch(Rectangle((0.25, 0), 5.4, 2.0, fc="0.04", ec="0.3"))
for yy, lab in [(1.86, "500"), (1.50, "250")]:        # marker 道
    axg.add_patch(Rectangle((0.42, yy - 0.045), 0.26, 0.09,
                            fc="0.85", ec="none"))
    axg.text(0.55, yy + 0.07, lab, fontsize=6.5, color="0.8",
             ha="center")
for i, g in enumerate(target):                        # 目标带 + actin 带
    cx = 1.35 + i * 0.85
    axg.add_patch(Rectangle((cx - 0.16, 1.21), 0.32, 0.18,
                            fc=str(min(0.92, 0.28 + 0.62 * g)), ec="none"))
    axg.add_patch(Rectangle((cx - 0.16, 0.89), 0.32, 0.18,
                            fc=str(min(0.92, 0.30 + 0.62 * actin[i])),
                            ec="none"))
axg.set_xticks([0.55] + [1.35 + i * 0.85 for i in range(5)],
               ["M"] + lanes)
axg.set_yticks([])
axb = fig.add_subplot(gs[1])
axb.bar(range(5), target, color=["#4C72B0", "#C44E52", "#C44E52",
                                 "#55A868", "#55A868"], width=0.6)
axb.axhline(1.0, ls="--", color="0.4", lw=0.9)
axb.set_xticks(range(5), lanes)
axb.set_ylabel("target / actin (rel. WT)")
fig.savefig("157-rt-pcr-gel.png", dpi=200, bbox_inches="tight")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

lines = open("158-gene-tissue.csv").read().splitlines()
tis = lines[0].split(",")[1:]
genes = [ln.split(",")[0] for ln in lines[1:]]
M = np.array([[float(v) for v in ln.split(",")[1:]] for ln in lines[1:]])
fig, ax = plt.subplots(figsize=(7.6, 5.2))
im = ax.imshow(M, cmap="RdBu_r", vmin=-2.5, vmax=2.5)
ax.set_xticks(range(len(tis)), tis, rotation=30, ha="right")
ax.set_yticks(range(len(genes)), genes)
for i in range(len(genes)):
    for j in range(len(tis)):
        ax.text(j, i, "%+.1f" % M[i, j], ha="center", va="center",
                fontsize=7.6,
                color="white" if abs(M[i, j]) > 1.4 else "0.15")
fig.colorbar(im, ax=ax, shrink=0.8, label="log2 ratio")
fig.savefig("158-gene-tissue-heatmap.png", dpi=200, bbox_inches="tight")