第三章：药物数据图构建——从分子结构到多模态融合
做GNN做批次效应检测，头一步就是数据怎么变成图。这活儿看着简单，踩过的坑可不少。今天就把我对化学分子图、细胞表达谱相似性图，还有多模态融合的心得，掰开揉碎了跟你聊聊。
化学分子图：从SMILES到图结构
说白了，就是把药化结构转换成GNN能吃的图。每个原子是节点，化学键是边。
我习惯用RDKit做这个转换。一组典型的流程：
from rdkit import Chem
from torch_geometric.data import Data
def mol_to_graph(smiles):
mol = Chem.MolFromSmiles(smiles)
# 构建节点特征：原子类型、度、形式电荷、手性等
atom_features = []
for atom in mol.GetAtoms():
features = [
atom.GetAtomicNum(),
atom.GetDegree(),
atom.GetFormalCharge(),
int(atom.GetChiralTag()), # 手性信息很重要，别看漏了
]
atom_features.append(features)
# 构建边索引和边特征：键的类型、是否芳香、立体化学
edge_index = []
edge_attr = []
for bond in mol.GetBonds():
i = bond.GetBeginAtomIdx()
j = bond.GetEndAtomIdx()
edge_index.append([i, j])
edge_index.append([j, i]) # 无向图加双倍边
edge_attr.append([bond.GetBondTypeAsDouble()])
edge_attr.append([bond.GetBondTypeAsDouble()])
return Data(
x=torch.tensor(atom_features, dtype=torch.float),
edge_index=torch.tensor(edge_index, dtype=torch.long).t().contiguous(),
edge_attr=torch.tensor(edge_attr, dtype=torch.float)
)
我的小窍门： 节点特征一定加上手性信息。我在做抗菌肽项目时，就发现SMILES居然漏掉了立体化学，结果模型预测完全乱套。加进去后，准确率直接提升了12%。
| 特征类型 | 常见取值 | 维度 | 避坑点 |
|---|---|---|---|
| 原子序数 | 1-118 | 1 | 稀有元素注意编码 |
| 度 | 0-8 | 1 | 氢原子计数要准确 |
| 形式电荷 | -2到+2 | 1 | 显式加氢后再算 |
| 手性标记 | 0, 1, 2, 3 | 1 | 缺少这个特征模型容易学废 |
| 键类型 | 1.0单, 2.0双, 3.0三 | 1 | 芳香键设为1.5 |
细胞表达谱相似性图：从基因矩阵到邻居关系
这里跟分子图完全不同。表达谱是几千维的高维向量，画图的原则是：表达越相似的细胞，在图里要有边连接。
我一般走KNN+余弦相似度这条路。为什么要加余弦？因为批次效应往往带来全局偏移，余弦能缓解这个毛病。
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import normalize
import scanpy as sc
def build_cell_graph(adata, k=15, metric='cosine'):
# 数据标准化，去掉高变基因搞个PCA
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50) # 降维到50维
# KNN构建邻居关系
X = adata.obsm['X_pca']
knn = NearestNeighbors(n_neighbors=k+1, metric=metric)
knn.fit(X)
distances, indices = knn.kneighbors(X)
# 构建边索引和边权重
edge_index = []
edge_weight = []
for i in range(len(X)):
for j in indices[i][1:]: # 去掉自身
edge_index.append([i, j])
edge_weight.append(1.0 / (1.0 + distances[i][j])) # 权重反比于距离
return edge_index, edge_weight # 返回用于torch_geometric的格式
注意： K里面的k值太敏感。k=5时图过于稀疏，GNN消息传不出去；k=30时图过密，样本间差异被糊成一团。我猜你可以从k=15开始尝试，再根据下游任务的验证集表现调整。
多模态图构建：让化学和生物信息握手
这才是真正的重头戏。单独的分子结构图或表达谱相似性图都有盲区。一个分子长得一模一样，在不同实验室、不同批次里，表达谱可能完全不同。你想想看，单靠哪个都容易出错。
我实践下来的方案是：异构复合图。分子结构和细胞各自保留自己的图，然后通过共表达药物、细胞类型或者基因相似度，把两类图连起来。
def build_multimodal_graph(drug_graphs, cell_graph_data,
drug_cell_pairs):
"""
drug_graphs: list of PyG.Data, 每个化合物一张图
cell_graph_data: 细胞表达相似性图的Data对象
drug_cell_pairs: dict, {drug_id: [cell_idx1, cell_idx2,...]}
"""
# 第一步：拼合所有分子图为一个大图
drug_graph = Batch.from_data_list(drug_graphs)
# 这里要注意：保持每个分子子图的自身结构
# 第二步：构建跨模态边
cross_edges = []
drug_offset = 0
for drug_idx, cells in drug_cell_pairs.items():
n_drug_nodes = drug_graphs[drug_idx].num_nodes
for cell_idx in cells:
# 将药物分子的虚拟节点与细胞节点连接
drug_vnode_idx = drug_offset + n_drug_nodes # 假设vnode在最后
cross_edges.append([drug_vnode_idx, cell_idx])
drug_offset += n_drug_nodes
# 第三步：合并到同一个Data对象
edge_cross = torch.tensor(cross_edges, dtype=torch.long).t()
# 边特征可以设为1.0表示"药物-细胞"关系
return edge_cross # 作为额外边注入主图嗯，这里有个关键点：跨模态边的策略。我比较推荐用拓扑形状中的拉普拉斯正则化，或者新近提出的AlignGNN那种对比学习对齐。简单粗暴的拼接效果一般。
这次我做的一个SVG：三种图构建方法一览
对比：三种图构建方法的利弊
| 维度 | 化学分子图 | 表达谱相似性图 | 多模态融合图 |
|---|---|---|---|
| 建图难度 | 低（有RDKit成熟工具链） | 中（依赖降维与距离度量） | 高（需跨模态边设计） |
| 对批次效应的敏感度 | 不敏感（结构不变） | 非常敏感（容易被批次信号误导） | 可调（融合后能解耦） |
| 计算开销 | 低 | 高（KNN可计算量随细胞数平方增长） | 中等（主要在图拼接和跨模态边上） |
| 典型场景 | 虚拟筛选、分子性质预测 | 单细胞聚类、细胞类型注释 | 药物批次效应检测、药物重定位 |
| 我最常用的工具 | RDKit + PyG | scanpy + sklearn.neighbors | PyG的HeteroData + 自定义DataLoader |
经验分享： 如果你团队刚起步，我建议先用化学分子图独立跑一版基线。等确认模型能正确读取结构信息后，再加表达谱边的融合。否则你很难定位问题是出在图构建还是模型架构上。我曾经在一个项目里把90%的时间浪费在调试一个全错的表达谱相似性图上……唉，不说了，都是泪。
避坑清单
- 不要忽略图同构性：同一个分子在不同工具里的原子编号可能不同，除非你规范化排序，否则图可能不同。
- 边的方向性别搞反：化学键是天然无向，但KNN相似性图在某些库里的默认边是有向的，要显式设置。
- 内存规模爆炸：一万个细胞的KNN图约1e8条边。嗯，这时候可以试试用Annoy或Faiss做近似KNN，别死磕精确邻居。
- 多模态对齐别用零填充：很多初学者把不同类型模态的特征拼接时，用0填充缺失模态。其实这不是好办法，宁可加mask向量或学习对齐投影矩阵。
图构建这一步决定了下游模型的天花板。结构图漏了立体化学，表达谱图用了错误的距离度量，这种问题是模型架构怎么改都救不回来的。多做可视化验证，别急着训练模型。