合成医薬品の変遷をバブルチャートで眺める話
あけましておめでとうございます（？） あれよあれよという間に2020年になりましたね。
オリンピックイヤーな今年はバブリーな一年にしたい！！・・・ということでバブルチャートで遊んでみました。
バブルチャートといえば、やはり書籍「ファクトフルネス」でも有名となったハンス・ロスリング博士らによる「動くバブルチャート」ではないでしょうか？
・・・こういうカッコいい図作りたい !!! と思ったら、好きなデータでグラフをかけるようにしてくださっている記事がありました。*1
そのライブラリ名は・・・Bubbly！！
バブリーさせていただきましょう！（訳：コピペした）
題材として、「医薬品として承認されている合成低分子が時代とともにどのように変化してきたか」眺めてみたいと思います。
ChEMBLからのデータ取得
まずはChEMBLからそれっぽいデータをとってきます。
トップページの「See all visualisations」をクリックすると、 「Drugs by Molecule Type and First Approval」というグラフがありました。 どうやら各医薬品分子が最初に承認された年を横軸にプロットされているようです。こちらのデータで遊んでみます。
画面下の「Browse all Approved Drugs」から遷移した先のページで、「Small molecule」のみにフィルターをかけ、csv形式でデータをダウンロードしました。
Pandasで読み込んで中身を確認します。区切りがセミコロン(;)となっていたのでsepで指定して読み込みました。
import pandas ps pd df = pd.read_csv('CHEMBL25.csv', sep=";") df.shape # (1636, 32)
16363エントリー、32列の情報をもつテーブルとなりました。 欠損値も多くこのままでは使えなさそうなので、必要な列だけを残して絞り込もうと思います。
まず、「Drug Type」という列を確認してみました。
print(df["Drug Type"].unique()) # ['9:Inorganic' '1:Synthetic Small Molecule' '10:Polymer' '7:Natural Product-derived']
無機化合物やポリマーも含まれているようです。合成低分子に絞り込みます。
df_syn = df[df["Drug Type"] == "1:Synthetic Small Molecule"] print(len(df_syn)) # 1149
1149エントリー残りました。構造情報は「Smiles」列にあるようですが、まだNaNとなっている行があるので除きます。
df_syn_smi = df_syn.dropna(subset=['Smiles']) print(len(df_syn_smi)) # 1106
SMILESの情報をもつ合成低分子は1106エントリーでした。今回使用する情報として以下の４列を残しました。
① Parent molecule (化合物のChEMBL ID)
② First Approval （初めて承認された年？）
③ oral （経口投与か否か？）
④ SMILES
各列の情報の内容は私の勝手な推測です。このまま突き進みます・・・すみません。。。
df2 = df_syn_smi[["Parent Molecule", "First Approval","Oral", "Smiles"]] df2.head()
こんな感じになりました。
RDKitで構造の描画と指標の計算
とりあえずDataFrameができたので、RDKitを使ってSMILESをMolオブジェクトに変換し、分子量等の指標を求めます。 RDKitのPandasToolsの使い方は「化学の新しいカタチ」さんの記事を参考にさせていただきました。*2
from rdkit import rdBase, Chem, DataStructs from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools from IPython.display import SVG # SMILESの列を指定して、Molオブジェクトの列を新しく追加 PandasTools.AddMoleculeColumnToFrame(df2, "Smiles") # Molオブジェクトから色々と指標を計算 df2['MW'] = df2.ROMol.map(Descriptors.MolWt) # 分子量 df2['LogP'] = df2.ROMol.map(Descriptors.MolLogP) # LogP df2['NHD'] = df2.ROMol.map(Descriptors.NumHDonors) # 水素結合ドナーの数 df2['NHA'] = df2.ROMol.map(Descriptors.NumHAcceptors) # 水素結合アクセプターの数 df2['NHetero'] = df2.ROMol.map(Descriptors.NumHeteroatoms) #ヘテロ原子の数 df2['NAR'] = df2.ROMol.map(Descriptors.NumAromaticRings) # 芳香環の数
こんな感じです。
注意点ですが、RDKitのPandasToolsは Pandas v0.25.1 以降のバージョンでは構造式の描画がうまくされず文字列になってしまうようです。イシュー (?)に書いてありました。*3
同様の症状がみられた場合は Pandas 0.25.0以前のバージョンにインストールしなおすと良いようです。私もこれで直りました。
最後に、時系列ならべるために「First Approval」でソートしました。
df3 = df2.sort_values('First Approval')
データの準備完了！
Bubbleチャートをプロット
本題のバブルチャートを作成します。
先の記事のパッケージ「bubbly」と「Plotly」のインストールが必要になります。
どちらも「pip install ~~」で大丈夫でした。
まずは必要なものをimport
import plotly from plotly.offline import init_notebook_mode, iplot init_notebook_mode() from bubbly.bubbly import bubbleplot
バブリー.バブリー・・・
以下の図では、各分子を
- X軸：分子量（MW）
- Y軸：LogP
- バブルのサイズ：芳香環の数(NAR)
とした設定でプロットし、継時的な変化を眺めてみます。さあコピペだ！
figure = bubbleplot(dataset=df3, x_column='MW', y_column='LogP', bubble_column='Parent Molecule', time_column='First Approval', size_column='NAR', x_range=[0,1000], y_range=[-5, 10], x_title='Molecular Weight', y_title='LogP', title='Changes of Drug Properties', scale_bubble=1) iplot(figure)
こんなんできました。
1980年代と90年代にポイントの数が一気に増える年があるようです。 また2010年代に入ってからバブルのサイズが大きくなったように見えました。芳香環の数が増えたのでしょうか？
もう少し情報を加えて、「Oralか否か」で色付けしてみたいと思います。 「True/False」のままではうまくいかなかったので、「1/0」に変換してプロットしました。
# 1をかけるとTrue/Falseが1/0に変換できるらしい df3['oral'] = df3['Oral']*1 # color_columnとcolorbar_titleの設定を付け足し figure = bubbleplot(dataset=df3, x_column='MW', y_column='LogP', bubble_column='Parent Molecule', time_column='First Approval', size_column='NAR', color_column='oral', x_range=[0,1000], y_range=[-5, 10], x_title='Molecular Weight', y_title='LogP', title='Changes of Drug Properties', colorbar_title='Oral', scale_bubble=1) iplot(figure)
色付けは基本的に「青色（0; oral = False）」「黄色（1; oral = True）」となっています。
が、どうも不具合があるようで、その年にカテゴリの一方しかない場合は青色・黄色以外の色となっています。
また、最初の再生ではカラーバーもカテゴリーの数に合わせて変化していますが、一度Pauseするとカラーバーがその時点で固定化されてしまいプロットにおける色を正確に反映しなくなってしまいます。
ですので色は「大体そんな感じ」くらいでとらえてください。（私の能力では修正できません。すみません・・・）
なんとなく、2000年代後半から黄色が増えたように思います。
抗体医薬等が承認されるにつれて「低分子は経口投与を重視」というようになったのでしょうか？
雰囲気だけで適当なことを言っています。フェイクニュースです。
化合物空間でプロット
さて、ここまできたらついでに化合物空間における変遷までみてまいましょう。
毎度コピペ記事で申し訳ありませんが@iwatobipen先生のAdvent Calendar記事 からコードを拝借し、フィンガープリントをt-SNEで2次元に落としてプロットします。*4
まずは時系列関係なくエントリー全てをつかって化合物空間を作ります。
import matplotlib.pyplot as plt import seaborn as sns from sklearn.manifold import TSNE # フィンガープリントを格納するための空のリスト fps = [] # DataFrameのMolオブジェクトを取り出して順番にフィンガープリントを計算しリストに追加 mols = df3['ROMol'] for mol in mols: fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2) arr = np.zeros((1,)) DataStructs.ConvertToNumpyArray(fp, arr) fps.append(arr) # numpyのarrayに変換 fps = np.array(fps) # t-SNEで次元圧縮 tsne = TSNE(random_state=123) res = tsne.fit_transform(fps) # プロット plt.clf() plt.figure(figsize=(12, 6)) plt.scatter(res[:,0], res[:,1])
こんな感じで散らばりました。
この空間が時代とともにどのように埋められていったのでしょうか？？
DataFrameにt-SNEでえられた2次元プロットのXとYの値を追加します。
df3['tSNE_X'] = res[:,0] df3['tSNE_Y'] = res[:,1]
X、Y軸はtSNEのプロットそのままに、バブルのサイズを「分子量(MW)」、色は「oral(True/False)」としてバブルチャートを作成してみます。
figure = bubbleplot(dataset=df3, x_column='tSNE_X', y_column='tSNE_Y', bubble_column='Parent Molecule', time_column='First Approval', size_column='MW', color_column='oral', x_range=[-60,60], y_range=[-60, 60], title='Changes in the Chemical Space', colorbar_title='Oral', scale_bubble=1)
あいかわらず色のバグはそのままで申し訳ありませんがこんな感じになりました。
予想以上に満遍なく全体にちらばりました。 時系列に従ってポイントの位置が集まる場所が見えたりしたら面白いかな？と思ったのですが、そんなに綺麗にはいかないですよね・・・
おわり
以上、バブルチャートをつかって合成化合物の変遷を眺めてみました。 これといった中身がなく申し訳ありませんが、なんだかバブリーだからオッケー！！
・・・というのもあれなので真面目な動画のページもご紹介させていただきます。
Norvartis社 Mike Shultz 氏によるウェビナー動画で、Hulinks社のページにありました。
URLはこちら です。StarDropというソフトウェアの使用例の紹介のようです。*5
動画の中で「Drug-likeness」について、過去２０年の承認薬の解析をもとに「Rule of 5」をいかに更新していくか、といった議論がされています。
経口薬の分子量は増大の傾向にあるものの、cLogPと水素結合のドナーの数(HBD)は変化していないそうな・・・
（日本語字幕がなくてよくわからなかったので雰囲気です。すみません。）
・・・なるほど
・・・おしまい！！
色々と間違いがあると思います。ご指摘いただければ幸いです。
