引言：为什么选择Matplotlib进行蛋白质结构可视化
蛋白质结构可视化是生物信息学和结构生物学研究中的重要环节。虽然专业的分子可视化软件如PyMOL、ChimeraX功能强大，但在数据分析流程中，使用Python进行自动化、定制化的可视化同样重要。Matplotlib作为Python生态中最基础且功能最强大的绘图库，结合NumPy等科学计算工具，可以实现从简单的氨基酸序列展示到复杂的三维蛋白质结构渲染。
本文将详细介绍如何使用Matplotlib及其扩展工具来展示蛋白质结构，包括：
基础的氨基酸序列可视化
二维投影的蛋白质骨架绘制
三维结构的展示方法
高级可视化技巧，如表面渲染、相互作用分析等
1. 基础准备：环境配置与数据获取
1.1 安装必要的库
# 核心科学计算库
pip install numpy scipy
# 数据可视化
pip install matplotlib
# 生物信息学工具
pip install biopython
# 3D可视化增强（可选）
pip install mpl_toolkits.mplot3d
# 结构数据处理
pip install mdanalysis # 用于处理轨迹文件
1.2 获取蛋白质结构数据
蛋白质结构数据通常从PDB（Protein Data Bank）获取。我们可以使用Biopython的PDB模块下载和解析PDB文件。
from Bio.PDB import PDBList, PDBParser
import os
# 下载PDB文件
def download_pdb(pdb_id, save_dir=’./pdb_files’):
if not os.path.exists(save_dir):
os.makedirs(save_dir)
pdblist = PDBList()
pdblist.retrieve_pdb_file(pdb_id, pdir=save_dir, file_format=’pdb’)
return os.path.join(save_dir, f”pdb{pdb_id}.ent”)
# 解析PDB文件
def parse_pdb(pdb_file):
parser = PDBParser(QUIET=True)
structure = parser.get_structure(‘protein’, pdb_file)
return structure
# 示例：下载并解析1CRN（胰蛋白酶抑制剂）
pdb_file = download_pdb(‘1CRN’)
structure = parse_pdb(pdb_file)
1.3 提取坐标信息
import numpy as np
def extract_coordinates(structure):
“””从PDB结构中提取所有原子的坐标和信息”””
atoms = []
for model in structure:
for chain in model:
for residue in chain:
for atom in residue:
atoms.append({
‘name’: atom.get_name(),
‘residue’: residue.get_resname(),
‘resnum’: residue.get_id()[1],
‘chain’: chain.get_id(),
‘coord’: atom.get_coord(),
‘element’: atom.element
})
return atoms
atoms = extract_coordinates(structure)
print(f”提取到 {len(atoms)} 个原子”)
print(f”示例原子: {atoms[0]}”)
2. 基础可视化：氨基酸序列展示
2.1 绘制单字母氨基酸序列
import matplotlib.pyplot as plt
from Bio.SeqUtils import seq1
def plot_amino_acid_sequence(structure, chain_id=’A’):
“””绘制氨基酸序列，用颜色区分不同性质的氨基酸”””
# 获取序列
sequence = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0]: # 只取标准氨基酸
sequence.append(seq1(residue.get_resname()))
# 氨基酸性质分类
hydrophobic = [‘A’, ‘V’, ‘L’, ‘I’, ‘M’, ‘F’, ‘W’, ‘Y’]
hydrophilic = [‘R’, ‘K’, ‘D’, ‘E’, ‘N’, ‘Q’, ‘S’, ‘T’, ‘C’, ‘G’, ‘H’, ‘P’]
colors = []
for aa in sequence:
if aa in hydrophobic:
colors.append(‘red’)
elif aa in hydrophilic:
colors.append(‘blue’)
else:
colors.append(‘gray’)
# 创建图形
fig, ax = plt.subplots(figsize=(15, 2))
# 绘制每个氨基酸
for i, (aa, color) in enumerate(zip(sequence, colors)):
ax.text(i, 0, aa, ha=’center’, va=’center’,
fontsize=12, fontweight=’bold’,
bbox=dict(boxstyle=’round’, facecolor=color, alpha=0.7))
ax.set_xlim(-1, len(sequence))
ax.set_ylim(-0.5, 0.5)
ax.set_title(f’Amino Acid Sequence (Chain {chain_id})’, fontsize=14)
ax.axis(‘off’)
# 添加图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor=’red’, label=’Hydrophobic’),
Patch(facecolor=’blue’, label=’Hydrophilic’),
Patch(facecolor=’gray’, label=’Other’)
]
ax.legend(handles=legend_elements, loc=’upper right’, bbox_to_anchor=(1.15, 1))
plt.tight_layout()
plt.show()
# 使用示例
plot_amino_acid_sequence(structure)
2.2 序列特征热图
def plot_sequence_features(structure, chain_id=’A’):
“””绘制序列特征热图，展示不同位置的物理化学性质”””
# 提取序列和特征
residues = []
hydrophobicity = []
b_factors = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0]:
residues.append(residue.get_resname())
# 简化的疏水性值（实际应用中可使用更精确的尺度）
hydrophobicity.append(len(residue.get_resname()) * 0.5)
# B因子（温度因子）
b_factor = np.mean([atom.get_bfactor() for atom in residue])
b_factors.append(b_factor)
# 创建特征矩阵
features = np.array([hydrophobicity, b_factors])
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 4), sharex=True)
# 疏水性曲线
ax1.plot(hydrophobicity, ‘o-‘, color=’green’, linewidth=2, markersize=6)
ax1.set_ylabel(‘Hydrophobicity’)
ax1.grid(True, alpha=0.3)
# B因子热图
im = ax2.imshow([b_factors], aspect=’auto’, cmap=’viridis’)
ax2.set_yticks([0])
ax2.set_yticklabels([‘B-factor’])
ax2.set_xlabel(‘Residue Position’)
ax2.set_xticks(range(0, len(residues), 5))
ax2.set_xticklabels(range(0, len(residues), 5))
# 添加颜色条
cbar = plt.colorbar(im, ax=ax2, orientation=’horizontal’, pad=0.2)
cbar.set_label(‘B-factor (Ų)’)
plt.suptitle(f’Sequence Features Analysis (Chain {chain_id})’, fontsize=14)
plt.tight_layout()
plt.show()
plot_sequence_features(structure)
3. 二维投影：蛋白质骨架绘制
3.1 Cα原子连接图（Ribbon图的简化版）
def plot_protein_backbone_2d(structure, chain_id=’A’, figsize=(10, 8)):
“””绘制蛋白质Cα骨架的二维投影”””
# 提取Cα原子坐标
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0] and ‘CA’ in residue:
ca_coords.append(residue[‘CA’].get_coord())
ca_coords = np.array(ca_coords)
# 创建图形
fig, ax = plt.subplots(figsize=figsize)
# 绘制Cα原子
ax.scatter(ca_coords[:, 0], ca_coords[:, 1],
c=range(len(ca_coords)), cmap=’rainbow’,
s=50, alpha=0.8, edgecolors=’black’, linewidth=0.5)
# 连接Cα原子
ax.plot(ca_coords[:, 0], ca_coords[:, 1],
‘k-‘, linewidth=1, alpha=0.6)
# 添加序列标签（每5个残基）
for i in range(0, len(ca_coords), 5):
ax.annotate(f'{i+1}’,
xy=(ca_coords[i, 0], ca_coords[i, 1]),
xytext=(5, 5), textcoords=’offset points’,
fontsize=8, alpha=0.7)
ax.set_xlabel(‘X (Å)’)
ax.set_ylabel(‘Y (Å)’)
ax.set_title(f’Protein Cα Backbone Projection (Chain {chain_id})’)
ax.grid(True, alpha=0.3)
# 添加颜色条表示残基顺序
sm = plt.cm.ScalarMappable(cmap=’rainbow’,
norm=plt.Normalize(vmin=0, vmax=len(ca_coords)))
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, label=’Residue Number’)
plt.tight_layout()
plt.show()
plot_protein_backbone_2d(structure)
3.2 二级结构标注的2D图
def plot_secondary_structure_2d(structure, chain_id=’A’):
“””绘制带二级结构标注的2D骨架图”””
# 首先需要检测二级结构（这里使用简化方法）
# 实际应用中可使用DSSP等工具
from Bio.PDB import DSSP
# 获取二级结构信息
model = structure[0]
dssp = DSSP(model, pdb_file, acc_method=’SASA’)
ss_elements = []
for key in dssp:
residue = dssp[key]
ss = residue[2] # 二级结构字符
if ss in [‘H’, ‘E’, ‘C’]: # Helix, Sheet, Coil
ss_elements.append((key[1], ss))
# 提取Cα坐标
ca_coords = []
ss_labels = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0] and ‘CA’ in residue:
ca_coords.append(residue[‘CA’].get_coord())
# 查找对应的二级结构
ss = ‘C’
for (res_num, ss_char) in ss_elements:
if res_num == residue.get_id()[1]:
ss = ss_char
break
ss_labels.append(ss)
ca_coords = np.array(ca_coords)
# 绘制
fig, ax = plt.subplots(figsize=(12, 8))
# 颜色映射
ss_colors = {‘H’: ‘red’, ‘E’: ‘blue’, ‘C’: ‘gray’}
colors = [ss_colors[ss] for ss in ss_labels]
# 绘制骨架
for i in range(len(ca_coords)-1):
ax.plot([ca_coords[i, 0], ca_coords[i+1, 0]],
[ca_coords[i, 1], ca_coords[i+1, 1]],
color=colors[i], linewidth=3, alpha=0.7)
# 绘制原子点
ax.scatter(ca_coords[:, 0], ca_coords[:, 1],
c=colors, s=60, alpha=0.9, edgecolors=’black’, linewidth=0.5)
# 添加图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor=’red’, label=’α-Helix’),
Patch(facecolor=’blue’, label=’β-Sheet’),
Patch(facecolor=’gray’, label=’Coil’)
]
ax.legend(handles=legend_elements, loc=’upper right’)
ax.set_xlabel(‘X (Å)’)
ax.set_ylabel(‘Y (Å)’)
ax.set_title(f’2D Backbone with Secondary Structure (Chain {chain_id})’)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plot_protein_backbone_2d(structure)
4. 三维结构可视化：从基础到高级
4.1 使用mpl_toolkits.mplot3d绘制3D Cα骨架
from mpl_toolkits.mplot3d import Axes3D
def plot_protein_3d_basic(structure, chain_id=’A’, figsize=(10, 10)):
“””基础3D蛋白质骨架可视化”””
# 提取Cα原子坐标
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0] and ‘CA’ in residue:
ca_coords.append(residue[‘CA’].get_coord())
ca_coords = np.array(ca_coords)
# 创建3D图形
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111, projection=’3d’)
# 绘制Cα原子
sc = ax.scatter(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
c=range(len(ca_coords)), cmap=’viridis’,
s=80, alpha=0.8, edgecolors=’black’, linewidth=0.5)
# 连接Cα原子
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
‘k-‘, linewidth=1.5, alpha=0.6)
# 添加颜色条
cbar = plt.colorbar(sc, ax=ax, shrink=0.6, pad=0.1)
cbar.set_label(‘Residue Number’)
ax.set_xlabel(‘X (Å)’)
ax.set_ylabel(‘Y (Å)’)
ax.set_zlabel(‘Z (Å)’)
ax.set_title(f’3D Protein Cα Backbone (Chain {chain_id})’)
# 设置视角
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_protein_3d_basic(structure)
4.2 3D球棍模型（Atoms & Bonds）
def plot_3d_ball_and_stick(structure, chain_id=’A’, max_atoms=50):
“””绘制3D球棍模型，展示原子和化学键”””
# 提取原子坐标（限制数量以提高性能）
atoms = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0]:
for atom in residue:
atoms.append({
‘name’: atom.get_name(),
‘coord’: atom.get_coord(),
‘element’: atom.element,
‘residue’: residue.get_resname(),
‘resnum’: residue.get_id()[1]
})
if len(atoms) >= max_atoms:
break
if len(atoms) >= max_atoms:
break
if len(atoms) >= max_atoms:
break
if len(atoms) >= max_atoms:
break
# 提取坐标和元素
coords = np.array([atom[‘coord’] for atom in atoms])
elements = [atom[‘element’] for atom in atoms]
# 元素颜色映射
element_colors = {
‘C’: ‘gray’, ‘N’: ‘blue’, ‘O’: ‘red’, ‘S’: ‘yellow’,
‘H’: ‘white’, ‘P’: ‘orange’, ‘F’: ‘green’, ‘Cl’: ‘green’
}
colors = [element_colors.get(el, ‘purple’) for el in elements]
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection=’3d’)
# 绘制原子
for i, (coord, color, element) in enumerate(zip(coords, colors, elements)):
ax.scatter(coord[0], coord[1], coord[2],
c=color, s=100, alpha=0.9,
edgecolors=’black’, linewidth=0.5,
label=element if i == 0 else “”)
# 简化的化学键检测（基于距离）
# 实际应用中应使用更精确的键长标准
bond_threshold = 1.8 # Å
for i in range(len(coords)):
for j in range(i+1, len(coords)):
distance = np.linalg.norm(coords[i] – coords[j])
if distance < bond_threshold:
ax.plot([coords[i, 0], coords[j, 0]],
[coords[i, 1], coords[j, 1]],
[coords[i, 2], coords[j, 2]],
'k-', linewidth=1, alpha=0.4)
ax.set_xlabel('X (Å)')
Ball and Stick Model (Chain {chain_id}, {len(atoms)} atoms)')
# 创建自定义图例
unique_elements = list(set(elements))
legend_elements = [plt.Line2D([0], [0], marker='o', color='w',
markerfacecolor=element_colors.get(el, 'purple'),
markersize=10, label=el) for el in unique_elements]
ax.legend(handles=legend_elements, loc='upper right')
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_3d_ball_and_stick(structure)
4.3 3D表面渲染（使用Marching Cubes算法）
from scipy.ndimage import gaussian_filter
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure
def generate_surface_grid(coords, grid_resolution=1.0, sigma=2.0):
"""生成用于表面渲染的3D网格"""
# 定义网格边界
min_coords = coords.min(axis=0) - 3.0
max_coords = coords.max(axis=0) + 3.0
# 创建网格
grid_size = int(np.max(max_coords - min_coords) / grid_resolution)
grid = np.zeros((grid_size, grid_size, grid_size))
# 将原子坐标映射到网格
for coord in coords:
x, y, z = coord - min_coords
xi = int(x / grid_resolution)
yi = int(y / grid_resolution)
zi = int(z / grid_resolution)
if 0 <= xi < grid_size and 0 <= yi < grid_size and 0 <= zi < grid_size:
grid[xi, yi, zi] = 1.0
# 应用高斯模糊模拟原子半径
grid = gaussian_filter(grid, sigma=sigma / grid_resolution)
return grid, min_coords, grid_resolution
def plot_3d_surface(structure, chain_id='A', grid_resolution=1.0):
"""绘制蛋白质表面渲染图"""
# 提取所有原子坐标
coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0]:
for atom in residue:
coords.append(atom.get_coord())
coords = np.array(coords)
# 生成表面网格
grid, min_coords, res = generate_surface_grid(coords, grid_resolution)
# 使用Marching Cubes算法提取表面
verts, faces, normals, values = measure.marching_cubes(grid, level=0.3)
# 转换回真实坐标
verts = verts * res + min_coords
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
# 创建表面集合
mesh = Poly3DCollection(verts[faces], alpha=0.7, edgecolor='none')
mesh.set_facecolor('lightblue')
ax.add_collection3d(mesh)
# 添加Cα骨架作为参考
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
ca_coords.append(residue['CA'].get_coord())
ca_coords = np.array(ca_coords)
if len(ca_coords) > 0:
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
‘k-‘, linewidth=2, alpha=0.8, label=’Cα Backbone’)
ax.set_xlabel(‘X (Å)’)
ax.set_ylabel(‘Y (Å)’)
ax.set_zlabel(‘Z (Å)’)
ax.set_title(f’3D Protein Surface Rendering (Chain {chain_id})’)
ax.legend()
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_3d_surface(structure)
5. 高级可视化技巧
5.1 相互作用分析：氢键网络
def detect_hydrogen_bonds(structure, chain_id=’A’, distance_cutoff=3.5):
“””检测氢键（简化版）”””
# 提取供体和受体原子
donors = [] # N, O
acceptors = [] # O, N
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0]:
for atom in residue:
if atom.element in [‘N’, ‘O’]:
if atom.element == ‘N’ and atom.get_name() in [‘N’, ‘ND1’, ‘NE2’]:
donors.append(atom)
elif atom.element == ‘O’:
acceptors.append(atom)
# 检测氢键（距离和角度条件）
hbonds = []
for donor in donors:
for acceptor in acceptors:
if donor.get_parent() == acceptor.get_parent():
continue
distance = np.linalg.norm(donor.get_coord() – acceptor.get_coord())
if distance < distance_cutoff:
# 简化的角度检查（实际需要更复杂）
hbonds.append((donor, acceptor, distance))
return hbonds
def plot_hydrogen_bonds(structure, chain_id='A'):
"""绘制氢键网络"""
hbonds = detect_hydrogen_bonds(structure, chain_id)
# 提取Cα坐标用于参考
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
ca_coords.append(residue['CA'].get_coord())
ca_coords = np.array(ca_coords)
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
# 绘制Cα骨架
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
'k-', linewidth=1, alpha=0.5, label='Cα Backbone')
# 绘制氢键
for donor, acceptor, distance in hbonds:
d_coord = donor.get_coord()
a_coord = acceptor.get_coord()
# 绘制氢键线
ax.plot([d_coord[0], a_coord[0]],
[d_coord[1], a_coord[1]],
[d_coord[2], a_coord[2]],
'r-', linewidth=2, alpha=0.7)
# 绘制原子
ax.scatter(d_coord[0], d_coord[1], d_coord[2],
c='blue', s=50, alpha=0.9)
ax.scatter(a_coord[0], a_coord[1], a_coord[2],
c='red', s=50, alpha=0.9)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'Hydrogen Bond Network (Chain {chain_id}, {len(hbonds)} bonds)')
# 添加图例
from matplotlib.lines import Line2D
legend_elements = [
Line2D([0], [0], color='k', lw=1, label='Cα Backbone'),
Line2D([0], [0], color='r', lw=2, label='Hydrogen Bond'),
Line2D([0], [0], marker='o', color='w', markerfacecolor='blue', markersize=8, label='Donor'),
Line2D([0], [0], marker='o', color='w', markerfacecolor='red', markersize=8, label='Acceptor')
]
ax.legend(handles=legend_elements, loc='upper right')
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_hydrogen_bonds(structure)
5.2 动态展示：轨迹动画
import matplotlib.animation as animation
def create_trajectory_animation(trajectory_file, top_file, chain_id='A', frame_step=10):
"""从分子动力学轨迹创建动画"""
try:
import MDAnalysis as mda
# 加载轨迹
u = mda.Universe(top_file, trajectory_file)
# 选择蛋白质原子
protein = u.select_atoms(f"protein and segid {chain_id}")
# 提取Cα轨迹
ca_atoms = protein.select_atoms("name CA")
n_frames = len(u.trajectory[::frame_step])
# 创建图形
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 初始化
sc = ax.scatter([], [], [], c='blue', s=50, alpha=0.8)
line, = ax.plot([], [], [], 'k-', linewidth=1, alpha=0.5)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title('Protein Trajectory Animation')
# 设置坐标范围
all_coords = []
for ts in u.trajectory[::frame_step]:
coords = ca_atoms.positions
all_coords.append(coords)
all_coords = np.concatenate(all_coords, axis=0)
ax.set_xlim(all_coords[:, 0].min()-5, all_coords[:, 0].max()+5)
ax.set_ylim(all_coords[:, 1].min()-5, all_coords[:, 1].max()+5)
ax.set_zlim(all_coords[:, 2].min()-5, all_coords[:, 2].max()+5)
def update(frame):
ts = u.trajectory[frame * frame_step]
coords = ca_atoms.positions
sc._offsets3d = (coords[:, 0], coords[:, 1], coords[:, 2])
line.set_data(coords[:, 0], coords[:, 1])
line.set_3d_properties(coords[:, 2])
ax.set_title(f'Protein Trajectory - Frame {ts.frame}')
return sc, line
anim = animation.FuncAnimation(fig, update, frames=n_frames,
interval=50, blit=False)
# 保存为GIF（需要imagemagick或pillow）
# anim.save('protein_trajectory.gif', writer='pillow', fps=20)
plt.show()
return anim
except ImportError:
print("MDAnalysis not installed. Skipping trajectory animation.")
return None
# 示例调用（需要轨迹文件）
# create_trajectory_animation('trajectory.xtc', 'topology.pdb')
5.3 交互式3D可视化（使用mplcursors）
def interactive_3d_structure(structure, chain_id='A'):
"""创建交互式3D结构，悬停显示原子信息"""
try:
import mplcursors
# 提取原子信息
atoms = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0]:
for atom in residue:
atoms.append({
'name': atom.get_name(),
'residue': residue.get_resname(),
'resnum': residue.get_id()[1],
'coord': atom.get_coord(),
'element': atom.element
})
coords = np.array([atom['coord'] for atom in atoms])
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
# 绘制原子
sc = ax.scatter(coords[:, 0], coords[:, 1], coords[:, 2],
c='blue', s=50, alpha=0.8)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'Interactive 3D Structure (Chain {chain_id})')
# 添加交互式光标
cursor = mplcursors.cursor(sc, hover=True)
@cursor.connect("add")
def on_add(sel):
idx = sel.index
atom = atoms[idx]
sel.annotation.set_text(
f"{atom['residue']} {atom['resnum']}\n"
f"{atom['name']} ({atom['element']})\n"
f"Coord: {atom['coord']}"
)
sel.annotation.get_bbox_patch().set(fc="yellow", alpha=0.9)
plt.show()
except ImportError:
print("mplcursors not installed. Install with: pip install mplcursors")
interactive_3d_structure(structure)
6. 实用技巧与性能优化
6.1 大分子处理策略
def process_large_protein(structure, chain_id='A', max_residues=100):
"""处理大分子蛋白质，采样显示"""
# 提取Cα坐标，按残基采样
ca_coords = []
residue_indices = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
count = 0
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
if count % max(len(residue) // max_residues, 1) == 0:
ca_coords.append(residue['CA'].get_coord())
residue_indices.append(residue.get_id()[1])
count += 1
ca_coords = np.array(ca_coords)
# 使用简化表示
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 绘制采样点
sc = ax.scatter(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
c=residue_indices, cmap='coolwarm', s=30, alpha=0.8)
# 连接线
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
'k-', linewidth=1, alpha=0.4)
cbar = plt.colorbar(sc, ax=ax, shrink=0.6)
cbar.set_label('Residue Number')
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'Sampled Structure (Chain {chain_id}, {len(ca_coords)} points)')
plt.tight_layout()
plt.show()
# process_large_protein(large_structure)
6.2 导出高质量图像
def export_high_quality_image(fig, filename, dpi=300, format='png'):
"""导出高质量图像"""
# 设置字体和样式
plt.rcParams.update({
'font.size': 10,
'font.family': 'sans-serif',
'axes.labelsize': 12,
'axes.titlesize': 14,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'legend.fontsize': 10,
'figure.titlesize': 16
})
# 保存
fig.savefig(filename, dpi=dpi, format=format,
bbox_inches='tight', facecolor='white', edgecolor='none')
print(f"图像已保存为 {filename} (dpi={dpi})")
# 使用示例
# fig = plt.figure(figsize=(10, 8))
# ... 绘图代码 ...
# export_high_quality_image(fig, 'protein_structure.png', dpi=600)
6.3 批量处理多个结构
def batch_visualize_pdb_files(pdb_ids, output_dir='./output'):
"""批量处理多个PDB文件"""
if not os.path.exists(output_dir):
os.makedirs(output_dir)
for pdb_id in pdb_ids:
try:
print(f"Processing {pdb_id}...")
# 下载和解析
pdb_file = download_pdb(pdb_id)
structure = parse_pdb(pdb_file)
# 创建多种可视化
fig1 = plt.figure(figsize=(10, 8))
ax = fig1.add_subplot(111, projection='3d')
# ... 3D绘图代码 ...
export_high_quality_image(fig1, f"{output_dir}/{pdb_id}_3d.png", dpi=300)
plt.close(fig1)
# 序列图
fig2, ax = plt.subplots(figsize=(15, 2))
# ... 序列绘图代码 ...
export_high_quality_image(fig2, f"{output_dir}/{pdb_id}_sequence.png", dpi=300)
plt.close(fig2)
except Exception as e:
print(f"Error processing {pdb_id}: {e}")
# 批量处理示例
# pdb_ids = ['1CRN', '1TIM', '2HHB']
# batch_visualize_pdb_files(pdb_ids)
7. 总结与最佳实践
7.1 何时使用Matplotlib vs 专业软件
Matplotlib适合：
自动化批量处理
整合到数据分析流程
定制化可视化需求
生成出版级静态图像
与NumPy/SciPy生态集成
专业软件（PyMOL/ChimeraX）适合：
交互式探索
实时渲染
复杂的分子表面和静电势
晶体学分析
7.2 性能优化建议
数据采样
：对于大分子，使用Cα骨架而非全原子
简化表示
：球棍模型只显示关键原子
使用向量图形
：导出PDF/SVG而非位图
并行处理
：批量处理时使用multiprocessing
7.3 扩展工具推荐
NGLview
：Jupyter中的交互式3D可视化
PyMOL
：通过Python API调用
BioPython
：结构分析和操作
MDAnalysis
：轨迹分析
通过本文的指南，您应该能够使用Matplotlib及其扩展工具创建从简单到复杂的蛋白质结构可视化。虽然Matplotlib在3D渲染方面不如专业软件强大，但其灵活性和Python生态集成使其成为生物信息学工作流程中不可或缺的工具。
# 使用Matplotlib绘制蛋白质结构图的实用指南 从基础代码到高级可视化技巧详解 如何用Python展示氨基酸序列与三维结构
引言：为什么选择Matplotlib进行蛋白质结构可视化
蛋白质结构可视化是生物信息学和结构生物学研究中的重要环节。虽然专业的分子可视化软件如PyMOL、ChimeraX功能强大，但在数据分析流程中，使用Python进行自动化、定制化的可视化同样重要。Matplotlib作为Python生态中最基础且功能最强大的绘图库，结合NumPy等科学计算工具，可以实现从简单的氨基酸序列展示到复杂的三维蛋白质结构渲染。
本文将详细介绍如何使用Matplotlib及其扩展工具来展示蛋白质结构，包括：
基础的氨基酸序列可视化
二维投影的蛋白质骨架绘制
三维结构的展示方法
高级可视化技巧，如表面渲染、相互作用分析等
1. 基础准备：环境配置与数据获取
1.1 安装必要的库
# 核心科学计算库
pip install numpy scipy
# 数据可视化
pip install matplotlib
# 生物信息学工具
pip install biopython
# 3D可视化增强（可选）
pip install mpl_toolkits.mplot3d
# 结构数据处理
pip install mdanalysis # 用于处理轨迹文件
1.2 获取蛋白质结构数据
蛋白质结构数据通常从PDB（Protein Data Bank）获取。我们可以使用Biopython的PDB模块下载和解析PDB文件。
from Bio.PDB import PDBList, PDBParser
import os
# 下载PDB文件
def download_pdb(pdb_id, save_dir='./pdb_files'):
if not os.path.exists(save_dir):
os.makedirs(save_dir)
pdblist = PDBList()
pdblist.retrieve_pdb_file(pdb_id, pdir=save_dir, file_format='pdb')
return os.path.join(save_dir, f"pdb{pdb_id}.ent")
# 解析PDB文件
def parse_pdb(pdb_file):
parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', pdb_file)
return structure
# 示例：下载并解析1CRN（胰蛋白酶抑制剂）
pdb_file = download_pdb('1CRN')
structure = parse_pdb(pdb_file)
1.3 提取坐标信息
import numpy as np
def extract_coordinates(structure):
"""从PDB结构中提取所有原子的坐标和信息"""
atoms = []
for model in structure:
for chain in model:
for residue in chain:
for atom in residue:
atoms.append({
'name': atom.get_name(),
'residue': residue.get_resname(),
'resnum': residue.get_id()[1],
'chain': chain.get_id(),
'coord': atom.get_coord(),
'element': atom.element
})
return atoms
atoms = extract_coordinates(structure)
print(f"提取到 {len(atoms)} 个原子")
print(f"示例原子: {atoms[0]}")
2. 基础可视化：氨基酸序列展示
2.1 绘制单字母氨基酸序列
import matplotlib.pyplot as plt
from Bio.SeqUtils import seq1
def plot_amino_acid_sequence(structure, chain_id='A'):
"""绘制氨基酸序列，用颜色区分不同性质的氨基酸"""
# 获取序列
sequence = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0]: # 只取标准氨基酸
sequence.append(seq1(residue.get_resname()))
# 氨基酸性质分类
hydrophobic = ['A', 'V', 'L', 'I', 'M', 'F', 'W', 'Y']
hydrophilic = ['R', 'K', 'D', 'E', 'N', 'Q', 'S', 'T', 'C', 'G', 'H', 'P']
colors = []
for aa in sequence:
if aa in hydrophobic:
colors.append('red')
elif aa in hydrophilic:
colors.append('blue')
else:
colors.append('gray')
# 创建图形
fig, ax = plt.subplots(figsize=(15, 2))
# 绘制每个氨基酸
for i, (aa, color) in enumerate(zip(sequence, colors)):
ax.text(i, 0, aa, ha='center', va='center',
fontsize=12, fontweight='bold',
bbox=dict(boxstyle='round', facecolor=color, alpha=0.7))
ax.set_xlim(-1, len(sequence))
ax.set_ylim(-0.5, 0.5)
ax.set_title(f'Amino Acid Sequence (Chain {chain_id})', fontsize=14)
ax.axis('off')
# 添加图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor='red', label='Hydrophobic'),
Patch(facecolor='blue', label='Hydrophilic'),
Patch(facecolor='gray', label='Other')
]
ax.legend(handles=legend_elements, loc='upper right', bbox_to_anchor=(1.15, 1))
plt.tight_layout()
plt.show()
# 使用示例
plot_amino_acid_sequence(structure)
2.2 序列特征热图
def plot_sequence_features(structure, chain_id='A'):
"""绘制序列特征热图，展示不同位置的物理化学性质"""
# 提取序列和特征
residues = []
hydrophobicity = []
b_factors = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0]:
residues.append(residue.get_resname())
# 简化的疏水性值（实际应用中可使用更精确的尺度）
hydrophobicity.append(len(residue.get_resname()) * 0.5)
# B因子（温度因子）
b_factor = np.mean([atom.get_bfactor() for atom in residue])
b_factors.append(b_factor)
# 创建特征矩阵
features = np.array([hydrophobicity, b_factors])
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 4), sharex=True)
# 疏水性曲线
ax1.plot(hydrophobicity, 'o-', color='green', linewidth=2, markersize=6)
ax1.set_ylabel('Hydrophobicity')
ax1.grid(True, alpha=0.3)
# B因子热图
im = ax2.imshow([b_factors], aspect='auto', cmap='viridis')
ax2.set_yticks([0])
ax2.set_yticklabels(['B-factor'])
ax2.set_xlabel('Residue Position')
ax2.set_xticks(range(0, len(residues), 5))
ax2.set_xticklabels(range(0, len(residues), 5))
# 添加颜色条
cbar = plt.colorbar(im, ax=ax2, orientation='horizontal', pad=0.2)
cbar.set_label('B-factor (Ų)')
plt.suptitle(f'Sequence Features Analysis (Chain {chain_id})', fontsize=14)
plt.tight_layout()
plt.show()
plot_sequence_features(structure)
3. 二维投影：蛋白质骨架绘制
3.1 Cα原子连接图（Ribbon图的简化版）
def plot_protein_backbone_2d(structure, chain_id='A', figsize=(10, 8)):
"""绘制蛋白质Cα骨架的二维投影"""
# 提取Cα原子坐标
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
ca_coords.append(residue['CA'].get_coord())
ca_coords = np.array(ca_coords)
# 创建图形
fig, ax = plt.subplots(figsize=figsize)
# 绘制Cα原子
ax.scatter(ca_coords[:, 0], ca_coords[:, 1],
c=range(len(ca_coords)), cmap='rainbow',
s=50, alpha=0.8, edgecolors='black', linewidth=0.5)
# 连接Cα原子
ax.plot(ca_coords[:, 0], ca_coords[:, 1],
'k-', linewidth=1, alpha=0.6)
# 添加序列标签（每5个残基）
for i in range(0, len(ca_coords), 5):
ax.annotate(f'{i+1}',
xy=(ca_coords[i, 0], ca_coords[i, 1]),
xytext=(5, 5), textcoords='offset points',
fontsize=8, alpha=0.7)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_title(f'Protein Cα Backbone Projection (Chain {chain_id})')
ax.grid(True, alpha=0.3)
# 添加颜色条表示残基顺序
sm = plt.cm.ScalarMappable(cmap='rainbow',
norm=plt.Normalize(vmin=0, vmax=len(ca_coords)))
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, label='Residue Number')
plt.tight_layout()
plt.show()
plot_protein_backbone_2d(structure)
3.2 二级结构标注的2D图
def plot_secondary_structure_2d(structure, chain_id='A'):
"""绘制带二级结构标注的2D骨架图"""
# 首先需要检测二级结构（这里使用简化方法）
# 实际应用中可使用DSSP等工具
from Bio.PDB import DSSP
# 获取二级结构信息
model = structure[0]
dssp = DSSP(model, pdb_file, acc_method='SASA')
ss_elements = []
for key in dssp:
residue = dssp[key]
ss = residue[2] # 二级结构字符
if ss in ['H', 'E', 'C']: # Helix, Sheet, Coil
ss_elements.append((key[1], ss))
# 提取Cα坐标
ca_coords = []
ss_labels = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
ca_coords.append(residue['CA'].get_coord())
# 查找对应的二级结构
ss = 'C'
for (res_num, ss_char) in ss_elements:
if res_num == residue.get_id()[1]:
ss = ss_char
break
ss_labels.append(ss)
ca_coords = np.array(ca_coords)
# 绘制
fig, ax = plt.subplots(figsize=(12, 8))
# 颜色映射
ss_colors = {'H': 'red', 'E': 'blue', 'C': 'gray'}
colors = [ss_colors[ss] for ss in ss_labels]
# 绘制骨架
for i in range(len(ca_coords)-1):
ax.plot([ca_coords[i, 0], ca_coords[i+1, 0]],
[ca_coords[i, 1], ca_coords[i+1, 1]],
color=colors[i], linewidth=3, alpha=0.7)
# 绘制原子点
ax.scatter(ca_coords[:, 0], ca_coords[:, 1],
c=colors, s=60, alpha=0.9, edgecolors='black', linewidth=0.5)
# 添加图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor='red', label='α-Helix'),
Patch(facecolor='blue', label='β-Sheet'),
Patch(facecolor='gray', label='Coil')
]
ax.legend(handles=legend_elements, loc='upper right')
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_title(f'2D Backbone with Secondary Structure (Chain {chain_id})')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
plot_secondary_structure_2d(structure)
4. 三维结构可视化：从基础到高级
4.1 使用mpl_toolkits.mplot3d绘制3D Cα骨架
from mpl_toolkits.mplot3d import Axes3D
def plot_protein_3d_basic(structure, chain_id='A', figsize=(10, 10)):
"""基础3D蛋白质骨架可视化"""
# 提取Cα原子坐标
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
ca_coords.append(residue['CA'].get_coord())
ca_coords = np.array(ca_coords)
# 创建3D图形
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111, projection='3d')
# 绘制Cα原子
sc = ax.scatter(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
c=range(len(ca_coords)), cmap='viridis',
s=80, alpha=0.8, edgecolors='black', linewidth=0.5)
# 连接Cα原子
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
'k-', linewidth=1.5, alpha=0.6)
# 添加颜色条
cbar = plt.colorbar(sc, ax=ax, shrink=0.6, pad=0.1)
cbar.set_label('Residue Number')
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'3D Protein Cα Backbone (Chain {chain_id})')
# 设置视角
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_protein_3d_basic(structure)
4.2 3D球棍模型（Atoms & Bonds）
def plot_3d_ball_and_stick(structure, chain_id='A', max_atoms=50):
"""绘制3D球棍模型，展示原子和化学键"""
# 提取原子坐标（限制数量以提高性能）
atoms = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0]:
for atom in residue:
atoms.append({
'name': atom.get_name(),
'coord': atom.get_coord(),
'element': atom.element,
'residue': residue.get_resname(),
'resnum': residue.get_id()[1]
})
if len(atoms) >= max_atoms:
break
if len(atoms) >= max_atoms:
break
if len(atoms) >= max_atoms:
break
if len(atoms) >= max_atoms:
break
# 提取坐标和元素
coords = np.array([atom[‘coord’] for atom in atoms])
elements = [atom[‘element’] for atom in atoms]
# 元素颜色映射
element_colors = {
‘C’: ‘gray’, ‘N’: ‘blue’, ‘O’: ‘red’, ‘S’: ‘yellow’,
‘H’: ‘white’, ‘P’: ‘orange’, ‘F’: ‘green’, ‘Cl’: ‘green’
}
colors = [element_colors.get(el, ‘purple’) for el in elements]
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection=’3d’)
# 绘制原子
for i, (coord, color, element) in enumerate(zip(coords, colors, elements)):
ax.scatter(coord[0], coord[1], coord[2],
c=color, s=100, alpha=0.9,
edgecolors=’black’, linewidth=0.5,
label=element if i == 0 else “”)
# 简化的化学键检测（基于距离）
# 实际应用中应使用更精确的键长标准
bond_threshold = 1.8 # Å
for i in range(len(coords)):
for j in range(i+1, len(coords)):
distance = np.linalg.norm(coords[i] – coords[j])
if distance < bond_threshold:
ax.plot([coords[i, 0], coords[j, 0]],
[coords[i, 1], coords[j, 1]],
[coords[i, 2], coords[j, 2]],
'k-', linewidth=1, alpha=0.4)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'3D Ball and Stick Model (Chain {chain_id}, {len(atoms)} atoms)')
# 创建自定义图例
unique_elements = list(set(elements))
legend_elements = [plt.Line2D([0], [0], marker='o', color='w',
markerfacecolor=element_colors.get(el, 'purple'),
markersize=10, label=el) for el in unique_elements]
ax.legend(handles=legend_elements, loc='upper right')
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_3d_ball_and_stick(structure)
4.3 3D表面渲染（使用Marching Cubes算法）
from scipy.ndimage import gaussian_filter
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure
def generate_surface_grid(coords, grid_resolution=1.0, sigma=2.0):
"""生成用于表面渲染的3D网格"""
# 定义网格边界
min_coords = coords.min(axis=0) - 3.0
max_coords = coords.max(axis=0) + 3.0
# 创建网格
grid_size = int(np.max(max_coords - min_coords) / grid_resolution)
grid = np.zeros((grid_size, grid_size, grid_size))
# 将原子坐标映射到网格
for coord in coords:
x, y, z = coord - min_coords
xi = int(x / grid_resolution)
yi = int(y / grid_resolution)
zi = int(z / grid_resolution)
if 0 <= xi < grid_size and 0 <= yi < grid_size and 0 <= zi < grid_size:
grid[xi, yi, zi] = 1.0
# 应用高斯模糊模拟原子半径
grid = gaussian_filter(grid, sigma=sigma / grid_resolution)
return grid, min_coords, grid_resolution
def plot_3d_surface(structure, chain_id='A', grid_resolution=1.0):
"""绘制蛋白质表面渲染图"""
# 提取所有原子坐标
coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0]:
for atom in residue:
coords.append(atom.get_coord())
coords = np.array(coords)
# 生成表面网格
grid, min_coords, res = generate_surface_grid(coords, grid_resolution)
# 使用Marching Cubes算法提取表面
verts, faces, normals, values = measure.marching_cubes(grid, level=0.3)
# 转换回真实坐标
verts = verts * res + min_coords
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
# 创建表面集合
mesh = Poly3DCollection(verts[faces], alpha=0.7, edgecolor='none')
mesh.set_facecolor('lightblue')
ax.add_collection3d(mesh)
# 添加Cα骨架作为参考
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
ca_coords.append(residue['CA'].get_coord())
ca_coords = np.array(ca_coords)
if len(ca_coords) > 0:
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
‘k-‘, linewidth=2, alpha=0.8, label=’Cα Backbone’)
ax.set_xlabel(‘X (Å)’)
ax.set_ylabel(‘Y (Å)’)
ax.set_zlabel(‘Z (Å)’)
ax.set_title(f’3D Protein Surface Rendering (Chain {chain_id})’)
ax.legend()
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_3d_surface(structure)
5. 高级可视化技巧
5.1 相互作用分析：氢键网络
def detect_hydrogen_bonds(structure, chain_id=’A’, distance_cutoff=3.5):
“””检测氢键（简化版）”””
# 提取供体和受体原子
donors = [] # N, O
acceptors = [] # O, N
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if ‘H’ in residue.get_id()[0]:
for atom in residue:
if atom.element in [‘N’, ‘O’]:
if atom.element == ‘N’ and atom.get_name() in [‘N’, ‘ND1’, ‘NE2’]:
donors.append(atom)
elif atom.element == ‘O’:
acceptors.append(atom)
# 检测氢键（距离和角度条件）
hbonds = []
for donor in donors:
for acceptor in acceptors:
if donor.get_parent() == acceptor.get_parent():
continue
distance = np.linalg.norm(donor.get_coord() – acceptor.get_coord())
if distance < distance_cutoff:
# 简化的角度检查（实际需要更复杂）
hbonds.append((donor, acceptor, distance))
return hbonds
def plot_hydrogen_bonds(structure, chain_id='A'):
"""绘制氢键网络"""
hbonds = detect_hydrogen_bonds(structure, chain_id)
# 提取Cα坐标用于参考
ca_coords = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
ca_coords.append(residue['CA'].get_coord())
ca_coords = np.array(ca_coords)
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
# 绘制Cα骨架
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
'k-', linewidth=1, alpha=0.5, label='Cα Backbone')
# 绘制氢键
for donor, acceptor, distance in hbonds:
d_coord = donor.get_coord()
a_coord = acceptor.get_coord()
# 绘制氢键线
ax.plot([d_coord[0], a_coord[0]],
[d_coord[1], a_coord[1]],
[d_coord[2], a_coord[2]],
'r-', linewidth=2, alpha=0.7)
# 绘制原子
ax.scatter(d_coord[0], d_coord[1], d_coord[2],
c='blue', s=50, alpha=0.9)
ax.scatter(a_coord[0], a_coord[1], a_coord[2],
c='red', s=50, alpha=0.9)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'Hydrogen Bond Network (Chain {chain_id}, {len(hbonds)} bonds)')
# 添加图例
from matplotlib.lines import Line2D
legend_elements = [
Line2D([0], [0], color='k', lw=1, label='Cα Backbone'),
Line2D([0], [0], color='r', lw=2, label='Hydrogen Bond'),
Line2D([0], [0], marker='o', color='w', markerfacecolor='blue', markersize=8, label='Donor'),
Line2D([0], [0], marker='o', color='w', markerfacecolor='red', markersize=8, label='Acceptor')
]
ax.legend(handles=legend_elements, loc='upper right')
ax.view_init(elev=20, azim=30)
plt.tight_layout()
plt.show()
plot_hydrogen_bonds(structure)
5.2 动态展示：轨迹动画
import matplotlib.animation as animation
def create_trajectory_animation(trajectory_file, top_file, chain_id='A', frame_step=10):
"""从分子动力学轨迹创建动画"""
try:
import MDAnalysis as mda
# 加载轨迹
u = mda.Universe(top_file, trajectory_file)
# 选择蛋白质原子
protein = u.select_atoms(f"protein and segid {chain_id}")
# 提取Cα轨迹
ca_atoms = protein.select_atoms("name CA")
n_frames = len(u.trajectory[::frame_step])
# 创建图形
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 初始化
sc = ax.scatter([], [], [], c='blue', s=50, alpha=0.8)
line, = ax.plot([], [], [], 'k-', linewidth=1, alpha=0.5)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title('Protein Trajectory Animation')
# 设置坐标范围
all_coords = []
for ts in u.trajectory[::frame_step]:
coords = ca_atoms.positions
all_coords.append(coords)
all_coords = np.concatenate(all_coords, axis=0)
ax.set_xlim(all_coords[:, 0].min()-5, all_coords[:, 0].max()+5)
ax.set_ylim(all_coords[:, 1].min()-5, all_coords[:, 1].max()+5)
ax.set_zlim(all_coords[:, 2].min()-5, all_coords[:, 2].max()+5)
def update(frame):
ts = u.trajectory[frame * frame_step]
coords = ca_atoms.positions
sc._offsets3d = (coords[:, 0], coords[:, 1], coords[:, 2])
line.set_data(coords[:, 0], coords[:, 1])
line.set_3d_properties(coords[:, 2])
ax.set_title(f'Protein Trajectory - Frame {ts.frame}')
return sc, line
anim = animation.FuncAnimation(fig, update, frames=n_frames,
interval=50, blit=False)
# 保存为GIF（需要imagemagick或pillow）
# anim.save('protein_trajectory.gif', writer='pillow', fps=20)
plt.show()
return anim
except ImportError:
print("MDAnalysis not installed. Skipping trajectory animation.")
return None
# 示例调用（需要轨迹文件）
# create_trajectory_animation('trajectory.xtc', 'topology.pdb')
5.3 交互式3D可视化（使用mplcursors）
def interactive_3d_structure(structure, chain_id='A'):
"""创建交互式3D结构，悬停显示原子信息"""
try:
import mplcursors
# 提取原子信息
atoms = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
for residue in chain:
if 'H' in residue.get_id()[0]:
for atom in residue:
atoms.append({
'name': atom.get_name(),
'residue': residue.get_resname(),
'resnum': residue.get_id()[1],
'coord': atom.get_coord(),
'element': atom.element
})
coords = np.array([atom['coord'] for atom in atoms])
# 创建3D图形
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
# 绘制原子
sc = ax.scatter(coords[:, 0], coords[:, 1], coords[:, 2],
c='blue', s=50, alpha=0.8)
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'Interactive 3D Structure (Chain {chain_id})')
# 添加交互式光标
cursor = mplcursors.cursor(sc, hover=True)
@cursor.connect("add")
def on_add(sel):
idx = sel.index
atom = atoms[idx]
sel.annotation.set_text(
f"{atom['residue']} {atom['resnum']}\n"
f"{atom['name']} ({atom['element']})\n"
f"Coord: {atom['coord']}"
)
sel.annotation.get_bbox_patch().set(fc="yellow", alpha=0.9)
plt.show()
except ImportError:
print("mplcursors not installed. Install with: pip install mplcursors")
interactive_3d_structure(structure)
6. 实用技巧与性能优化
6.1 大分子处理策略
def process_large_protein(structure, chain_id='A', max_residues=100):
"""处理大分子蛋白质，采样显示"""
# 提取Cα坐标，按残基采样
ca_coords = []
residue_indices = []
for model in structure:
for chain in model:
if chain.get_id() == chain_id:
count = 0
for residue in chain:
if 'H' in residue.get_id()[0] and 'CA' in residue:
if count % max(len(residue) // max_residues, 1) == 0:
ca_coords.append(residue['CA'].get_coord())
residue_indices.append(residue.get_id()[1])
count += 1
ca_coords = np.array(ca_coords)
# 使用简化表示
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 绘制采样点
sc = ax.scatter(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
c=residue_indices, cmap='coolwarm', s=30, alpha=0.8)
# 连接线
ax.plot(ca_coords[:, 0], ca_coords[:, 1], ca_coords[:, 2],
'k-', linewidth=1, alpha=0.4)
cbar = plt.colorbar(sc, ax=ax, shrink=0.6)
cbar.set_label('Residue Number')
ax.set_xlabel('X (Å)')
ax.set_ylabel('Y (Å)')
ax.set_zlabel('Z (Å)')
ax.set_title(f'Sampled Structure (Chain {chain_id}, {len(ca_coords)} points)')
plt.tight_layout()
plt.show()
# process_large_protein(large_structure)
6.2 导出高质量图像
def export_high_quality_image(fig, filename, dpi=300, format='png'):
"""导出高质量图像"""
# 设置字体和样式
plt.rcParams.update({
'font.size': 10,
'font.family': 'sans-serif',
'axes.labelsize': 12,
'axes.titlesize': 14,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'legend.fontsize': 10,
'figure.titlesize': 16
})
# 保存
fig.savefig(filename, dpi=dpi, format=format,
bbox_inches='tight', facecolor='white', edgecolor='none')
print(f"图像已保存为 {filename} (dpi={dpi})")
# 使用示例
# fig = plt.figure(figsize=(10, 8))
# ... 绘图代码 ...
# export_high_quality_image(fig, 'protein_structure.png', dpi=600)
6.3 批量处理多个结构
def batch_visualize_pdb_files(pdb_ids, output_dir='./output'):
"""批量处理多个PDB文件"""
if not os.path.exists(output_dir):
os.makedirs(output_dir)
for pdb_id in pdb_ids:
try:
print(f"Processing {pdb_id}...")
# 下载和解析
pdb_file = download_pdb(pdb_id)
structure = parse_pdb(pdb_file)
# 创建多种可视化
fig1 = plt.figure(figsize=(10, 8))
ax = fig1.add_subplot(111, projection='3d')
# ... 3D绘图代码 ...
export_high_quality_image(fig1, f"{output_dir}/{pdb_id}_3d.png", dpi=300)
plt.close(fig1)
# 序列图
fig2, ax = plt.subplots(figsize=(15, 2))
# ... 序列绘图代码 ...
export_high_quality_image(fig2, f"{output_dir}/{pdb_id}_sequence.png", dpi=300)
plt.close(fig2)
except Exception as e:
print(f"Error processing {pdb_id}: {e}")
# 批量处理示例
# pdb_ids = ['1CRN', '1TIM', '2HHB']
# batch_visualize_pdb_files(pdb_ids)
7. 总结与最佳实践
7.1 何时使用Matplotlib vs 专业软件
Matplotlib适合：
自动化批量处理
整合到数据分析流程
定制化可视化需求
生成出版级静态图像
与NumPy/SciPy生态集成
专业软件（PyMOL/ChimeraX）适合：
交互式探索
实时渲染
复杂的分子表面和静电势
晶体学分析
7.2 性能优化建议
数据采样
：对于大分子，使用Cα骨架而非全原子
简化表示
：球棍模型只显示关键原子
使用向量图形
：导出PDF/SVG而非位图
并行处理
：批量处理时使用multiprocessing
7.3 扩展工具推荐
NGLview
：Jupyter中的交互式3D可视化
PyMOL
：通过Python API调用
BioPython
：结构分析和操作
MDAnalysis
：轨迹分析
通过本文的指南，您应该能够使用Matplotlib及其扩展工具创建从简单到复杂的蛋白质结构可视化。虽然Matplotlib在3D渲染方面不如专业软件强大，但其灵活性和Python生态集成使其成为生物信息学工作流程中不可或缺的工具。