DeepChem源碼分析（一）

本人化學信息學&藥物設計領域小白一枚。最近在學習使用深度學習框架DeepChem開展工作。DeepChem是斯坦福大學Bharath Ramsundar等人所開發的Python工具包，整合了許多深度學習和化學信息學中的常用工具，功能非常強大。但其Tutorial寫得卻稍顯混亂，而且本身DeepChem封裝程度較高，所以對於初學者而言直接看示例難免顯得有些摸不著頭腦。在此斗膽想寫文從源碼入手剖析一下DeepChem的組織結構，同時也是作為讀書筆記進行一下整理。

0 原始數據

數據處理都是從獲取原始數據開始的。藥物設計中的原始數據一般都是以表格形式存儲。本地保存時常見的是csv格式的文件。下圖是一個示例：

DeepChem希望在處理之前原始數據集至少要包含兩個信息，一個是每一個樣本未特徵化前的表示，如上圖中的drug列，二是每個樣本的標籤信息，如上圖中的n1列和n2列。一個樣本可以有多個任務（如一個化合物對應多個蛋白靶標的活性值），因此可以有多個標籤。有幾個任務，標籤列就有幾列。

1 原始數據的特徵化：Featurizer

獲得原始數據後下一步是對原始數據的特徵化。原始數據可以是蛋白序列，基因序列，小分子的字元表示（SMILES)等。特徵化是將原始表示轉換為程序可以識別和處理的特徵表示的過程。從SMILES式到其對應的ECFP指紋就是典型的特徵化過程。如何進行特徵化是一個很大的課題，在此不深入討論。

在DeepChem中，由Featurizer類處理特徵化過程。來看一下Featurizer抽象類的核心部分：

class Featurizer(object):

raise NotImplementedError(Featurizer is not defined.)


可以看出，Featurizer類非常簡單，只有一個可調用方法，即self.featurize. 其輸入是一個列表，其中每個元素是每個樣本的初始表示，如SMILES式。featurize方法對列表中的每個元素循環調用self._featurize處理，並將結果添加到features列表中。如果處理失敗，則添加一個空列表進去。最後將features返回。即使沒有注釋也很容易理解，self._featurize是一個子類需要覆蓋的方法，功能是輸入一個樣本的原始表示，返回其特徵。通過這種方法我們可以定義不同特徵化方法的Featurizer類。下面是一個用分子質量直接作為其特徵的簡單示例（deepchem.weight.basic.MolecularWeight）：

class MolecularWeight(Featurizer):
"""
Molecular weight.
"""
name = [mw, molecular_weight]

def _featurize(self, mol):
"""
Calculate molecular weight.
Parameters
----------
mol : RDKit Mol
Molecule.
"""
from rdkit.Chem import Descriptors
wt = Descriptors.ExactMolWt(mol)
wt = [wt]
return wt


注意倒數第二行。返回的不是wt，而是[wt]，因為wt是一個標量。而對於另一個示例deepchem.weight.basic.RDKitDescriptor:

class RDKitDescriptors(Featurizer):
"""
RDKit descriptors.
See http://rdkit.org/docs/GettingStartedInPython.html
#list-of-available-descriptors.
"""
name = descriptors

def _featurize(self, mol):
"""
Calculate RDKit descriptors.
Parameters
----------
mol : RDKit Mol
Molecule.
"""
rval = []
for desc_name, function in self.descList:
rval.append(function(mol))
return rval


該特徵暴力計算了所有可以用RDKit計算的分子描述符並返回一個由這些分子描述符構成的列表。

2 從序列的特徵化到DataFrame的特徵化

Featurizer可以做到從一串SMILES式構成的列表中提取特徵，但如前所述，我們的原始數據一般是表格，以csv存儲在本地，以pandas的DataFrame格式存儲於內存中，所以我們更關心如何從DataFrame（以下簡稱df）中獲取特徵。deepchem.data.featurize_smiles_df函數提供了這樣的功能：

def featurize_smiles_df(df, featurizer, field):
sample_elems = df[field].tolist()

features = []
#1
from rdkit import Chem
from rdkit.Chem import rdmolfiles
from rdkit.Chem import rdmolops
for ind, elem in enumerate(sample_elems):
mol = Chem.MolFromSmiles(elem)
if mol:
new_order = rdmolfiles.CanonicalRankAtoms(mol)
mol = rdmolops.RenumberAtoms(mol, new_order)
#2
features.append(featurizer.featurize([mol]))
valid_inds = np.array(
[1 if elt.size > 0 else 0 for elt in features], dtype=bool)
features = [elt for (is_valid, elt) in zip(valid_inds, features) if is_valid]
return np.squeeze(np.array(features), axis=1), valid_inds


由於這個函數默認了原始表示是SMILES式，所以會有#1到#2中間的部分。其實我們完全可以對其簡化，使得這個函數更一般化。該函數接收3個輸入，df是原始數據表，featurizer是我們構造好的特徵提取器，field是一個沒見過的參數，它指定了你給的df中哪一列是你所希望特徵化的樣本列，對應的是df中相應列的column名。

這個函數中循環的寫法非常繞。現在假設我們有三個樣本，第一個樣本的特徵是[0,0,1]，第二個樣本的特徵是[1,1,0]，第三個樣本的特徵提取失敗，特徵為空，則我們模擬上述過程應當為:

>>> feature1 = np.array([[0,0,1]]) #featurizer.featurize([mol1])的返回值
>>> feature2 = np.array([[1,1,0]])
>>> feature3 = np.array([np.array([])])
>>> features = 0
>>> features.append(feature1)
>>> features.append(feature2)
>>> features.append(feature3)
>>> features
[array([[0, 0, 1]]), array([[1, 0, 0]]), array([], shape=(1, 0), dtype=float64)]
>>> valid_inds = np.array([1 if elt.size > 0 else 0 for elt in features], dtype=bool)
array([ True, True, False], dtype=bool)


valid_inds保存了哪些樣本計算失敗的信息。

>>> features = [elt for (is_valid, elt) in zip(valid_inds, features) if is_valid]
>>> features
[array([[0, 0, 1]]), array([[1, 1, 0]])]


使用valid_inds去除features中計算失敗的特徵。此時每個樣本特徵的最外層外套就可以脫掉了。每個樣本特徵最外層外套對應np.asarray(features)的axis=1,使用np.squeeze脫掉外套。

>>> np.squeeze(np.array(features), axis=1)
array([[0, 0, 1],
[1, 1, 0]])


最後返回兩個np.ndarray，第一個是特徵向量，第二個是成功計算特徵的樣本標籤。

推薦閱讀：