magattacaのブログ

日付以外誤報

PyMOLをJupyter Notebookと連携させる方法について

前回の記事で、PyMOLのスクリプト機能を利用して3D PSAを求めました。『PyMOL=生物学向けの描画ソフト』と思い込んでいたのですが、まだまだ便利な機能がたくさんあって、化学方面にも応用できそうです。

素敵なPyMOLですが、プログラミングができない私にとってコマンドラインでの操作はとっつきづらいです。コマンドと一緒にコメントや出力結果を残しやすいJupyter notebookで使えた方が後でたどりやすいのになー、と、いうことでJupyterから使う方法を調べました。

PyMOL WikiJupyterという、そのものな項目があったので記載の方法を試していきたいと思います。ほぼそのままです。

検証する操作の確認

Jupyter notebookを試す前に、検証で行う操作を念のためPyMolアプリのGUIで確認します。

  1. ラニンを読み込む(fragmentcmd.fragment)
  2. 描画を画像(PNG)で保存(pngcmd.png)

アミノ酸残基はBuildメニューのResidueから選択して作成することができます。コマンドラインfragment alaとするか、cmd.fragment('ala')としてもOKです。

f:id:magattaca:20210109154316p:plain

画像の保存はFileメニューのExport Image Asあるいは右上のDraw/Rayパネルからできます。コマンドラインで行うには、png testとすると作業ディレクトリの中に「test.png」として保存できます。cmd.png('test')でもOKです。

f:id:magattaca:20210109154344p:plain

詳しい使い方はPyMol Bookが日本語でわかりやすいです。なお、cmdPyMOL Python APIを使うためのモジュールで、PyMol上のコマンドラインを使うのであればわざわざ使わなくても良いと思いますが、以降の内容との関係するので併記しました。

JupyterからPyMOL

では本題、Jupyter notebookからPyMOLを利用する方法です。PyMol Wikiには大きく2つ、直接PyMOL APIを使う方法と、サードパーティーのライブラリを使う方法が書かれています。それぞれ制限があるようです。

概要はこんな感じです。

直接PyMOL APIを使う

同じスレッド

まず、同じスレッドで使う方法です。こちらはGUIの画面をみることはできません(headless)。

操作は簡単。Jupyter notebookを起動し、PyMOLのモジュールを読み込むだけです。

# cmd モジュールのインポート
from pymol import cmd 

# アラニンを構築し、拡大した後に画像を保存
cmd.fragment('ala')
cmd.orient()
cmd.png('./test.png', ray=1)

PyMOL not running, entering library mode (experimental)」とかいう出力がでてきますが問題なく使えます。PyMOL WikiLaunching From a Scriptによると、PyMol 2.1から、cmdを使うと自動的にメインスレッドでバックエンドプロセスをGUI無しで起動するらしいので、その関係だと思います。*1

さて、上記で無事「test.png」が作成されたので中身を確認してみます。

from IPython.display import Image
Image(filename='./test.png')

f:id:magattaca:20210109154444p:plain

うまくいきました!

非同期スレッド

次に、非同期スレッド(asynchronous thread)で使う方法です。こちらはGUIが使えます。

が、残念ながらmacOSではこの操作を行えません(GitHub Issue)。

コマンドだけコピペしておきます。

# PyMOLウィンドウを開く
import sys
import pymol
_stdouterr = sys.stdout, sys.stderr
pymol.finish_launching(['pymol', '-q'])
sys.stdout, sys.stderr = _stdouterr

# PyMOLウィンドウにアラニンを読み込む
from pymol import cmd
cmd.fragment('ala')

Macでこれを実行するとカーネルが死にます。。。

コマンドラインモード

PyMOL Wikiから外れますが、コマンドラインモードであればMacでもJupyterからPyMOLを立ち上げることができます。起動のオプション-cを使うと、GUIなしのコマンドラインのみのモードとなり、バッチプロセスに使うことができます(PyMOLWiki Command Line Options)。

import pymol

# PyMolをコマンドラインモードで立ち上げる
pymol.pymol_argv = ['pymol','-c']
pymol.finish_launching()
# Detected 4 CPU cores.  Enabled multithreaded rendering.

ターミナルにPyMOL起動時のバージョン情報が表示され、起動がうまくいったことがわかります。オプション-qをつければこの表示が出ないようにすることもできます(quietモード)。

注意点として、このPyMOLの立ち上げ方はCPUを100%使うらしいです*2Macが悲鳴をあげるので不要であればcmd.quit()で終了しておいた方が良さそうです。

# コマンドラインモードのPyMOLを終了
cmd.quit()

RPCコミュニケーション

つづいて、Jupyter notebookとは別に立ち上げたPyMol(スタンドアローン)をRPCコミュニケーションで操作する方法です。RPC(Remote procedure call、遠隔手続き呼び出し)は「プログラムから別のアドレス空間にあるサブルーチンや手続きを実行することを可能にする技術」で、「遠隔相互作用の詳細を明示的にコーディングする必要がない」そうです(Wikipedia 遠隔手続き呼び出し)。・・・よくわかりません。

まず、pymol -Rで立ち上げておきます。-RオプションはRPCサーバーを立ち上げるためのオプションです。

pymol -R
# xml-rpc server running on host localhost, port 9123

PyMOL GUIアプリが立ち上がりますが、上記のように表示され、xml-rpc serverlocalhostのポート番号9123で起動していることがわかります。

このRPCサーバーにJupyter notebookからアクセスするにはPythonxmlrpc.clientモジュールを使えば良いそうです。

# Jupyter notebookからPyMOL にRPCコミュニケーション
import xmlrpc.client as xmlrpclib
cmd = xmlrpclib.ServerProxy('http://localhost:9123')
cmd.fragment('ala')

以上を実行すると、PyMOLのGUIウィンドウにアラニンが読み込まれます。

xmlrpc.client.ServerProxy()でRPCコミュニケーションを管理するためのオブジェクトを生成し、これにPyMOL Python APIを扱うためのモジュールcmdと同じ名前をつけておくことで、APIを扱う感覚で外部のPyMolを操作ができているようです。・・たぶん。

以上が、直接PyMol APIを使う方法です。

サードパーティーライブラリ

つづいてサードパーティーのライブラリを使った方法についてです。ライブラリとしてIPyMOLRDKitが言及されていますがPyMOL Wikiには詳細は書かれていません。

ライブラリの説明やGitHubをたよりに試してみましょう。

IPyMOL

まずIPyMOLです。IPythonからPyMolのセッションを制御するためのライブラリで、こちらもRPCコミュニケーションを利用します(MITライセンス)。

「Jupyter Notebookできちんと見せたい時や、PyMOLスクリプトのプロトタイプを簡便に作成したい時に理想的なツール」とのことなので、今回の目的にまさにぴったりなライブラリです。

pipでインストール可能(ipymol - PyPI)ですが、古いバージョンらしく不具合があるらしいです(GitHub issue #21)。GitHubバージョンは修正されているらしいので、直接インストールしてみます。

pip install git+https://github.com/cxhernandez/ipymol.git

IPyMOLでは、さきと同様に、あらかじめpymol -Rで立ち上げておいたPyMOLに接続できますが、IPyMolからstart()PyMOL RPCサーバーを起動することもできます。閉じる時はstop()です。

from ipymol import viewer as pymol
pymol.start()

PyMOLのGUIが立ち上がりました!今回は「xml-rpc server running on host localhost, port 9124」と出たので、ポート番号が先とは異なるようです。

ラニンを読み込むには以下のようにします。

pymol.fragment('ala')

display()を使うことでPyMOL GUIのディスプレイをJupyter Notebook上に描画することができます。

pymol.display()

f:id:magattaca:20210109154735p:plain

IPyMolで利用可能なAPIは以下を打ち込むと確認できます。出力は多いので省略します。

print(pymol._server.system.listMethods())

上記のリストに含まれていないコマンドを使いたい時はpymol.do("実行したいコマンド")とすればOKです。

使い方の実例はこちらのノートブックが参考になります。

終わったらstop()でウィンドウを閉じられます。

# IPyMOLで立ち上げたPyMOLの終了
pymol.start()

RDKit PyMolモジュール

続いてRDKitからPyMolに RPCコミュニケーションする方法です。Chem.PyMolモジュールを利用すればOKで、化学の新しいカタチさんが「RDKitからPyMOLを利用する」で詳しい使い方をご紹介してくださっています。

参考までにこれまでと同じ処理を確認します。

まず、pymol -RXML-RPC サーバーを立ち上げておきます。

RDKitで必要なモジュールをインポートし、MolViewerオブジェクトを作成します。

from rdkit import Chem
from rdkit.Chem import PyMol

v = PyMol.MolViewer()

MolViewerに含まれていない操作を行うには、sever.do("実行したいコマンド")とすれば良いそうです。アラニンフラグメントを読み込み、GetPNG()でJupyter notebook上に描画します。

# アラニンの読み込み
v.server.do("fragment ala")

# ズーム
v.Zoom('ala')

# PyMOL GUIの画面をPNGで取得
v.GetPNG()

f:id:magattaca:20210109154938p:plain

SaveFile()を使うことで画像として保存することもできます。

折角ですので、ついでに前回の記事で取り上げた3D PSAのスクリプトを実行してみましょう!

まずDeleteAll()で一度全て消した後に、LoadFile()でファイルを読み込みます。「conformer-00.xyz」というxzyファイルを「conf1」という名前のオブジェクトとしてロードしています。

# viewewrのオブジェクトを全て消す
v.DeleteAll()

# ファイルの読み込み
v.LoadFile('./conformer-00.xyz', 'conf1')

別途作成したスクリプトpsa3d.py)をつかえるように読み込み(run)、計算を実行します。

# スクリプトの実行
v.server.do('run ./psa3d.py')

# 計算の実行
v.server.do('set solvent_radius=1.4; set dot_density=4; set dot_solvent=1; psa3d')

無事計算が実行され、PyMOL上にSA 3D PSAの値「201.973」が表示されました。

以上、サードパーティライブラリの動作も確認できました。

まとめ

以上、PyMOLをJupyter notebookで利用するための方法の紹介でした。ほとんどPyMOL Wikiそのままですが。。。

PyMOLのスクリプトは色々な方が作成、公開してくださっているようなので、Jupyter notebookとうまく連携させればプログラミング苦手な私でももっと活用できるかもしれません。

Jupyter notebookで操作して出力された結果を、PyMOLから受け取ってJupyter Notebookに表示することもしたかったのですが、やり方がよくわかりませんでした。サブプロセス(?)を使えば良いのかもしれませんが私のキャパシティでは無理でした。。。

ではでは。

*1:どういうことかはさっぱりわかりません。。。

*2:Launching From a Script

PyMolで3D PSAを計算しようとする話

先日の記事で、下記文献を例に3次元構造の指標であるRadius of gyrationについて取り上げました。

論文1) Solution Conformations Shed Light on PROTAC Cell Permeability
Yoseph Atilaw, Vasanthanathan Poongavanam, Caroline Svensson Nilsson, Duy Nguyen, Anja Giese, Daniel Meibom, Mate Erdelyi, and Jan Kihlberg
ACS Medicinal Chemistry Letters Article ASAP DOI: 10.1021/acsmedchemlett.0c00556
https://pubs.acs.org/doi/10.1021/acsmedchemlett.0c00556 (CC-BY 4.0)

今回は同文献で使われていた別の指標SA 3D PSAについて調べました。
合わせて、この記事では以下の内容に触れています。

  • シェルスクリプトによるテキストファイルの処理
  • RDKitに適したxyz形式の変換方法
  • PyMolのコマンドと自作関数の使い方

SA 3D PSAについて

SA 3D PSAの定義や求め方について上記文献には記載がありませんでした。求め方として著者らのグループ(Uppsala大学、 Jan Kihlberg教授)が2020年に発表した下記文献が引用されていましたが、オープンアクセスではなかったため中身を確認できませんでした。

論文2) Solution Conformations Explain the Chameleonic Behaviour of Macrocyclic Drugs
E. Danelius, V. Poongavanam, S. Peintner, L. H. E. Wieske, M. Erdélyi, J. Kihlberg
Chem. Eur. J. 2020, 26, 5231.

さらにたどると同グループの以下の文献にあたりました。

論文3) Impact of Dynamically Exposed Polarity on Permeability and Solubility of Chameleonic Drugs Beyond the Rule of 5
Matteo Rossi Sebastiano, Bradley C. Doak, Maria Backlund, Vasanthanathan Poongavanam, Björn Over, Giuseppe Ermondi, Giulia Caron, Pär Matsson, and Jan Kihlberg
Journal of Medicinal Chemistry 2018 61 (9), 4189-4202

こちらも論文の中身自体は確認できませんでしたが、Supporting Informationに3D PSA calcuration methodsという記載がありました。

ざっと以下のような内容です。

  • 3Dで分子の極性原子(polar atom)の占める面積
  • 分子の表面積(molecular surface area)は、プローブの半径をとしたもの
    3D PSA
  • 溶媒接触表面積(solvent accessible surface area)は、プローブの半径を1.4Åとしたもの
    SA 3D PSA

これら2種の3D PSAの計算にはPyMol 1.7.0.1が使われたとのことです。

また、極性原子は

  1. 窒素原子(N)
  2. 酸素原子(O)
  3. それらに結合する水素原子(NH、OH)
  4. 閾値よりも大きな部分電荷(absolute partial chage)を持つ原子

と、なっています。このうち4. 部分電荷PM3 semi-empirical methodで計算した値で、SIから判断する限り閾値0.6のようです。

なお、こちらのJ. Med. Chem.文献(論文3)よりも後に出版された、先のChem. Eur. J.文献(論文2)は、電荷の計算方法が改善しており、その手法をACS. Med. Chem. Lett.文献(論文1)では採用しているようです。

ざっとですが3D PSAの概要がつかめました。

同様な値がオープンソースツールを使って計算できるか?早速やって見ましょう!

入力構造の準備(シェルスクリプトの練習)

先の記事で、論文1 のSupporting Informationからコンフォマーの配座情報を取り出しました。このSIには全部で24個のコンフォマー情報が含まれています。

今回はテキストファイル処理の練習も兼ねて、すべて取り出して見たいと思います。・・・シェルスクリプト苦手。

処理対象のファイルは前回pdfminer.sixでPDFから取り出したテキストのうち、コンフォマー情報を含むテキストファイルです(SI.txt)。

sedコマンドを使えば良いとのことですが、MacsedBSDというもので、LinuxsedGNU)と挙動が異なるようです。

後者の方が使い勝手が良いらしいので、GNUsedをインストールしてみます。*1

Homebrewで簡単に入れられます。

brew install gun-sed   

注意点として、sedの代わりにgsedとコマンドを入力する、とのこと。

ここで、行いたい処理は下図の通り、余分なテキストの削除xyz形式に合わせるためのテキスト追加です。

f:id:magattaca:20210103235736p:plain

ファイルを上書きするオプション-iと、文字列の正規表現を使って、以下のようにしました。*2

(Jupyter notebookで行なっているのでマジックコマンドを含めて!gsedとしています)

# 改ページの文字列(\f)の削除
!gsed -i 's/\f//g' SI.txt

# ページ番号の行を削除。行の先頭が数字となるもの(「^[0-9]」)。
!gsed -i '/^[0-9]/d' SI.txt

# 空行を削除。「^ *$」で改行のみおよび空白のみを指定。
!gsed -i '/^ *$/d' SI.txt

# conformer_xx行の前に、原子数(131)の行を追加。行の先頭がconとなるもの(「^con」)。
!gsed -i '/^con/i 131' SI.txt

上記で得られたファイル(SI.txt)の中身は、24個のコンフォマーの情報が連続する3192行のファイルです。

余談ですが、一番最初の「改ページのメタ文字(\f)」を削除しなければいけないことに気づくまでに2時間くらいかかりました。。。。
何度gsedを実行しても、テキストエディタを変えて確認してもうまくいかない理由がわからず、結局「Pages」で開いて改ページが含まれているに気づきました。*3

皆様もPDFから起こしたテキストにはくれぐれもご注意を!  

さて、気を取り直して!

コンフォマー1つあたり133行なので、stripコマンドを133行ごとに別のファイルに区切ります。

ここでもMacGNU版のstripを使うため、Homebrewでcoreutilsをインストールし、gstripコマンドを使います。

brew install coreutils

分割後のファイルは「conformer-数字.txt」となるようにします。

!gsplit -l 133 -d SI.txt conformer-

「conformer-XX(Xは数字)」というファイルが「XX = 00 ~ 23」の24個できました。

拡張子をxyzに変更するためrenameでファイル名を変更します。これもインストールが必要のようです。

brew install rename

-aオプションで末尾に「.xyz」を追加します。

!rename -a .xyz conformer-*

できました! 拡張子が付与され自動的にAvogadroに紐づけられました。

f:id:magattaca:20210103235849p:plain

XYZ2molでRDKit Molオブジェクトに一発変換

先の記事ではxyzファイルの変換にOpenBabelを使いましたが、出力構造に難がありマニュアルで修正する必要がありました。

今回はxyz2molGitHub)を使ってみたいと思います。
こちらはRSCのワークショップ(open-source-tools-for-chemistry)で紹介されていました。こちらに関連資料があります。

xyz2molの使用には以下が必要とのことなので、インストールします。

  • rdkit (>= 2019.9)
  • networkx
conda install -c anaconda networkx

続いて、githubから作業フォルダにそのままcloneします。

!git clone https://github.com/jensengroup/xyz2mol

READMEには「フォルダを移動してmakeするように」、とかいてありますが、しなくてもimportできました。

import xyz2mol.xyz2mol as xyz2mol

先のワークショップのうち、「Chemistry in the cloud: leveraging Google Colab for quantum chemistry (by Jan Jensen氏)」のGoogle Colab Notebook (URL→ https://bit.ly/37fIYbp ) に利用例があるので、コードをお借りして進めます。

まず、read_xyz_file()でxyzファイルを読み込みます。原子(atom)、電荷(chage)、xyz座標(coordinate)を返すので、それらをxyz2mol()に渡します。するとRDKitのMol objectを含むリストが得られます。

atoms, charge, xyz_coordinates = xyz2mol.read_xyz_file('conformer-00.xyz')
conf_mol_0 = xyz2mol.xyz2mol(atoms, xyz_coordinates, charge=charge)

print(type(conf_mol_0))
print(len(conf_mol_0))
print(type(conf_mol_0[0]))
#  <class 'list'>
#  1
#  <class 'rdkit.Chem.rdchem.Mol'>

無事RDKitのMolオブジェクトが得られました!

少し時間がかかる印象ですが、OpenBabelと異なり一発でxyzファイルから生成できるのは便利ですね。

肝心の構造はどうでしょうか?py3Dmolで描画し、OpenBabelでMolファイル変換した構造と比較してみました

f:id:magattaca:20210104000118p:plain

無事正しい構造が起こせていそうです!!これは便利!

ぶっちゃけxyzファイルを変換しなくても以降の計算は可能なのですが、ファイル変換の方法なんてなんぼあっても良いですからね!

3D PSAの計算

さて、構造の準備と変換方法が分かったので、目的の3D PSA計算に挑戦したいと思います。

冒頭確認した通り、必要な情報は

  1. 各原子のxyz座標
  2. 各原子が極性原子に該当するか(原子タイプと部分電荷)?

でした。1.は手元にあるので、2.を求める方法があるか検討してみました。

部分電荷の計算に失敗

WinMopac

論文3) J. Med. Chem. 2018の記載に従えば、PM3 semi-empirical methodで部分電荷を求める必要があります。

PM3って何でしょう???量子力学さっぱりわかりません。

PC CHEM BASICS.COMさんの記事(MOPACの入手とインストール方法)によれば、とりあえずMOPACがお手軽だそうです。

そこで記載を参考に、Macで実行するため以下をダウンロードしました。

  1. WinMopac 7.21
  2. Wine4 (Mac上でWindowsのソフトを使えるようにする。Wine環境の準備)
  3. Avogadro (GUIでMOPACの入力ファイルを用意できる)

ネット回線が弱いのですごい時間かかりました。。。。

で、いざWinPopacに構造を投げ込んだところ「最大で処理可能な原子数を超えてます」みたいエラーで動きませんでした。

残念。お正月休みを返して欲しい。

xTB

ついで、xTB(extended tight binding)を利用してみることにしました。

xTBはオープンソースで利用可能で、さきにご紹介したRSCのワークショップで使われていたソフトウェアです。

検索した結果、「半経験量子計算法:Tight-Binding(強結合近似)計算」という言葉がでてきました。さっぱりわかりません。。。
とりあえずPM3もsemi-empricalって言ってるしこれを試してみよう!ということでインストールしました。

condaで入ります。

conda install -c conda-forge xtb

xTBのドキュメントによると、single point calculationを行うにはシェルでxtb coord --sccとするだけのようです。

coordは原子座標でxyzファイルを入力として使うこともできます。また、--sccself-consistent charge calculationの略のようです。

以下を実行しました。

!xtb ./conformer-00.xyz --scc

・・・止まりました。

今回も!何の成果も得られませんでした!!!!

あきらめました。

座標のみで計算

諦めが早いので、部分電荷にこだわらず座標原子タイプ(N、O、NH、OH)の情報のみで計算することとしました。

論文3) J. Med. Chem. 2018のSupporting InformationよりPyMolで求められるようなので、こちらの情報に従って試していきましょう。

PyMolで使用するスクリプトの準備

計算用のスクリプトを用意し、PyMolのコマンドラインから使用するだけらしく、スクリプトの情報もありました。

まず、以下を記載したファイルpsa3d.pyを用意します。

from pymol import cmd

def psa3d (arg1=0.6, arg2=-0.6):
    obj_list = cmd.get_names('objects')
    for obj in obj_list:
        cmd.select('noh', '(pc.>'+str(arg1)+'|pc.<'+str(arg2)+'|(e. n OR e. o OR e. h AND(neighbor e. n OR e. o))) AND ' + obj)
        print("3DPSA_noh06 %12s %.3f" % (obj, cmd.get_area('noh')))
    return obj

cmd.extend("psa3d", psa3d)

このスクリプトでは大きく3つのことが行われています。

  1. PyMol Python APIを使う準備(cmdモジュール)
  2. 3D PSAを求めるための関数の定義(psa3d)
  3. 定義した関数をPyMolコマンドライン上で使うための準備(cmd.extend)

1.でimportするcmdモジュールですが、これによりPython Python APIにアクセスすることができ、cmd.<command-name>(argument, ...)という使い方ができます。*4

3.cmd.extend(文字列, 関数)は、PyMol外部の関数をPyMolスクリプトで使うために、新しい名前(文字列)のコマンドに結びつけます。ここでは自作の関数psa3d()の名前をpsa3dと結びつけてるので、PyMolのコマンドライン上でpsa3dとすれば使えるようになります。

それでは関数psa3d (arg1=0.6, arg2=-0.6)の中身をみていきます。

ここで使われているPyMol特有のコマンドは以下の3つです。*5

コマンド 説明
cmd.get_names('objects') PyMol上のオブジェクトの名前を取得(リスト形式)
cmd.select('name', 'selection') 'selection'で選択した原子(群)に、名前('name')をつけて、その名前で扱えるようにする
get_area('selection') 'selection'で選択された原子(群)の表面積(Å2)を計算

これを踏まえると、psa3d()は「PyMolに読み込まれている構造(オブジェクト)について、それぞれ3D PSA算出の対象タイプとなる極性原子を選択し、その原子(群)の表面積を求める」という関数になりそうです。

cmd.select()selectionがややこしいので、もう少しみてみます。

論文3のSIより、現在考えたい極性原子は以下です。(再掲)

  • 窒素原子(N)
  • 酸素原子(O)
  • それらに結合する水素原子(NH、OH)
  • 閾値(今回は0.6)よりも大きな部分電荷(absolute partial chage)を持つ原子

以上をふまえてpsa3d(arg1=0.6, arg2=-0.6)の引数の値を代入してselectionを書き直してみます。

(pc.> 0.6 | pc.<-0.6 | (e. n OR e. o OR e. h AND(neighbor e. n OR e. o))) AND オブジェクト)

上記のうち、「縦線(|)」と「OR」は「または」、「AND」は「且つ」と読み替えれば論理式のように読むことができそうです。

以下のSelection Algebraをふまえると、selectionが目的の原子の選択に対応していることがわかると思います。*6

Selection Algebra 説明
pc. partial_charge(部分電荷
pc. <演算子(ex. 大・小・イコール)> <>」
e. elemment(元素記号
e. <元素記号(ex. 窒素 n、酸素 o、水素 h)>」
neighbor 直接結合する近接原子

以上のような条件設定で選んだselectionに対して「noh」という名前を与えて新しいselectionにしているようです。

これでpsa3d.pyの説明は終わりです。・・・思ったより長くなった。

PyMolでスクリプトを利用して計算

用意したスクリプトpsa3d.py)をPyMolで使うには、PyMol起動後、PyMol コマンドラインで以下を打ち込むだけでOKです。

run ~/Desktop/psa3d.py

Pythonファイルのパスは適宜変更してください。今回はわかりやすくするためデスクトップに置いておきました。

次に、計算したい構造ファイルをPyMolに読み込みます。今回は、以前の記事で作成したMolファイルを使いました。

これで準備は完了です!それでは3D PSAを計算しましょう!!

以下のコマンドをコピペすれば2種類の3D PSAをそれぞれ計算することができます。

  1. 分子の表面積(M 3D PSA) 「set dot_density=4; set dot_solvent=0; psa3d
  2. 溶媒接触表面積(SA 3D PSA)「set solvent_radius=1.4; set dot_density=4; set dot_solvent=1; psa3d

コマンドの中身を簡単に説明します。

どちらも最後はpsa3dとなっています。psa3d.pyの最後の行でcmd.extend()を使って、関数psa3d()の名前をpsa3dと結びつけてるので、 PyMolのコマンドライン上でもpsa3dでと打てば関数を使うことができるようになっています。

それよりも前の記述は全て「表面積を求めるコマンド(get_area)」の条件を設定(set)しています。*7

設定条件 説明
dot_solvent 0 :分子の表面積(molecular surface)を計算(デフォルト)
1 : 溶媒接触表面積(SASA)を計算
dot_density 1 ~ 4:サンプリング密度(density)。高密度(数字が大きい)ほど精度は高いが、処理に時間がかかる
solvent_radius 溶媒の半径(デフォルトは1.4

上記とコマンドを見比べると、どちらも最も精度の高いdot_density=4を選択していることがわかります。

また、SA 3D PSAではプローブの半径を1.4Åとしていたので、dot_solvent=1とした上でset solvent_radius=1.4という設定が入っています。

これらの設定が関数psa3dの中のcmd.get_area()に反映されているようですね。

さて、コマンドの中身がわかったので、早速計算を実行してみましょう! コンフォマー1の3D PSAはどのような値を取っているでしょうか?

計算自体は一瞬で終わりました。結果はこんな感じです。

f:id:magattaca:20210104000813p:plain

それぞれ 2132)、2022)となりました。SA 3D PSAの方が小さくなるのは意外でした。

少し見辛いですが、上図の通り、「酸素原子: 赤色のドット」、「窒素原子: 青色のドット」で極性原子の表面積が表されています。

ではこの計算値、論文1)の文献の値とくらべてどうでしょうか? Figure. 5 を引用します。

f:id:magattaca:20210104000842p:plain

Conformer 1のSA 3D PSAは2022)!!同等な値が出ました!今回は部分電荷、あまり関係なかった感じでしょうか?

なお、入力ファイルをmolファイルではなく、座標のみのxyz形式としても同じ結果がでました。

電荷を考えないのであれば3D PSAの計算は原子の座標のみで良さそうです。

芳香属性とか結合の情報はあまり重視しない指標なのでしょうか???

またFig. 5に記載されている残りの2つについても、コンフォマー2のSA 3D PSAは 182.137、コンフォマー3のSA 3D PSAは 194.195と、と同等の結果が得られました。

まとめ

以上、今回は立体構造の指標3D PSAについて調べてみました。論文の本文自体を読めておらず、Supporting Informationの情報のみからの推察ですが、計算結果としてはおおよそ近い値が得られました。

スクリプトが提供されていたお陰ですが、PyMolでさっと計算できてしまうのが個人的にはびっくりでした。描画だけでも盛りだくさんなのに、こんな便利な使い方もできるとは、、、PyMolすごい!

感動したのでPyMolコマンドを調べていたら無駄に長くなってしまいました・・・。

しかし、量子力学がさっぱりわからないです。DFTはおろかHFもよくわかっていません。いったいみなさんどこで勉強されたんでしょうか??

シェルスクリプトの段階で躓いているので、まだまだ道は遠そうです。今回も間違いが多そうなのでご指摘いただければ幸いです。

ではでは!!

「イシューからはじめよ」を頼りに「WHOをゆく-感染症との闘いを超えて」を読んだ

あけましておめでとうございます。本年もよろしくお願いいたします。

さて、今更ですが尾身茂先生のご著書「WHOをゆく」を読みました。

WHOをゆく: 感染症との闘いを超えて

WHOをゆく: 感染症との闘いを超えて

  • 作者:尾身 茂
  • 発売日: 2011/10/21
  • メディア: 単行本

160ページ前後と決して分厚い本ではありませんが、

  • 「どのような経緯を経てWHOで働くに至ったか?」という青春時代を振り返るエッセイ(キャリア論)から、
  • 「WHOで如何に西太平洋地域のポリオ根絶を達成したか?」という戦略・実践のお話、
  • また、効果的なリーダーシップを発揮する秘訣・リーダーの持つべき心構え(リーダーシップ論

と、 多岐にわたる話題が盛りだくさんでとても面白かったです。

なんといっても凄いのはポリオ根絶という疑いようのない業績だと思います。恥ずかしながら私は全く存じ上げておりませんでした。

で、このお話を読んで思ったのは、「あれ?これイシューからはじめよじゃん。」

と、いうわけで「第2章 ポリオ根絶:第2の青春物語」を「イシューからはじめよ」に沿ってご紹介してみたいと思います。

※ 私は医学の専門家でも公衆衛生の専門家でもありません。素人の素朴な感想としてお読みいただけると幸いです。
※ 現在、手元に「イシューからはじめよ」がないので、メモとネットのまとめ記事を参考にかいてます。すみません。 *1

生産性の観点から

さて、安宅和人氏の著書「イシューからはじめよ-知的生産の「シンプルな本質」」は、 元になった氏のブログ記事「圧倒的に生産性の高い人(サイエンティスト)の研究スタイル」 からも分かるように、「生産性の高い仕事を行うにはどうすれば良いか?」が論じられています。

では、まず生産性の観点から「WHOをゆく 第2章」(以降、本書)を読み返してみましょう。

生産性の定義は以下。

f:id:magattaca:20210101231854p:plain
「イシューからはじめよ」における生産性の定義

本書におけるアウトプット(成果)は「西太平洋地域のポリオ根絶」です。・・・もう、これだけで圧倒的。。。

インプットを投下した時間とすると、
* 1990年9月 WHO西太平洋地域事務局(WPRO) 着任とともに着手
* 2000年10月29日 西太平洋地域 ポリオ根絶宣言
ですので、約10年で達成されたことになります。

根絶の確証を得るための期間を3年と設定し、最後の患者発生報告が1997年3月であることを考えると 6年半と捉えることもできるかもしれません。*2

WPRO着任までの尾身先生のご経歴は地域医療や肝炎ウイルスの研究などで、ポリオの専門ではなく、 また、当時オフィスに「地域内のポリオに関する詳細なデータもな」かったそうなので、 ほぼゼロからのスタートだったといっても良いかもしれません。

上の図式に当てはめてみます。

f:id:magattaca:20210101232459p:plain
10年で根絶を達成

参考までにWHO西太平洋地域は地図上のオレンジ色の地域です。*3

面積を見ても人口の多い中国を含むことを考えても、この領域から10年でウイルスを根絶というのは、 まさしく圧倒的な生産性なのではないでしょうか?

イシューの観点から

さて、「イシューからはじめよ」によると「バリューのある仕事」には、 「今本当に答えを出すべき問題 = イシュー」を見極めることが大切となるそうです。

イシュー度が高く、解の質が高い課題こそが解くべき問題、とのこと。

f:id:magattaca:20210101232825p:plain
解くべき問題(あやふやな記憶で再現した図)

ではイシューとは何か?以下の2つを満たすものだそうです。

  1. 2つ以上の集団で決着がついていない問題
  2. 根本に関わる、もしくは白黒がはっきりしていない問題

では、「ポリオの根絶」とはイシューの観点からみてどうだったのでしょうか?

まず、尾身先生のWPRO着任に先立つこと1988年、WHO総会で天然痘に続く次なる感染症根絶の目標として 「2000年までにポリオの根絶」が決議されました。また、日本では1980年以降、野生株によるポリオ患者が確認されていなかったそうです。 *4

つまりポリオは、WHOが根絶すべき、かつ根絶可能とみなしていた感染症で、 当時既に限定的ではあるがポリオの発症を抑えられている地域もあったようです。

ところが、尾身先生が着任1年後(1991年4月)に開催した専門家会議では、呈示した戦略に対して厳しい質問が相次ぎ、 ワクチン購入のための資金提供についても一つも約束が得られませんでした。会議後、資金援助を求めて様々な公的機関を巡り歩いても門前払い同然だったそうです。 

つまり、専門家会議、公的機関ともに「非常に難しい(ムリ)」とみなしていた、ということになると思います。

まとめると、「根絶すべき重要な感染症公衆衛生の根本的な問題ではあるが、WHOとその他の機関では根絶の可否に意見の相違があった」。 イシューの2つの定義を満たしていそうです。

f:id:magattaca:20210101233204p:plain

解の質の観点から

解くべき課題の横軸「イシュー度」について確認しました。 では、縦軸「解の質」について、WHOの目指す解「ポリオの根絶」とは「具体的にはどのような状態をあらわす」のでしょうか?

ここで「根絶」は、単にその地域の「患者の発症数をゼロ」にするだけではなく、 「その病原体(ウイルス)を完全に消滅」させることを意味します。

なぜ発症を抑えるだけでなく、ウイルスをゼロにする必要があるのか?それにはポリオの特徴が関わっています。

ポリオは100~200人が感染し、1人だけ発症する不顕性感染が多い感染症です。 さらに感染後症状がなくても糞便中に数週間にわたりウイルスが排泄され、新たな感染源となります。 *5

つまり「"ゼロ"の証明」にこそ意味がある感染症なのです。

「"ゼロ"の証明」、まるで悪魔の証明のようです。*6

実際、「ポリオウイルスが環境中に存在しない」ことを言う必要条件は、 「良質なサーベイランス」(患者をきちんと発見し、診断できる)の下で、「2ヶ月以上ポリオ患者の報告がない」こと、と述べられてますが、 尾身先生らはさらに「ポリオ"根絶"の証明を確実にする」 ために、3年という非常に長い期間を設定されています。

まとめると、ポリオ根絶プロジェクトの目指す解として、ウイルスの特徴を踏まえた非常に高い目標が設定されていた、ということになりそうです。

f:id:magattaca:20210101233330p:plain

イシュー度、解の質、ともに高く、重要な課題であったことがわかりました。それではこの問題の解決にどのように取り組まれたのでしょうか?
仮説ドリブンアウトプットドリブンメッセージドリブンの3つの観点を頼りにご紹介したいと思います。

1) 仮説ドリブンの観点から

「イシューからはじめよ」では「仮説ドリブン」で課題解決に取り組むことの重要性が言われています。 課題に対してスタンスを取り、 ①イシューをサブイシューに分解したうえで、 ②ストーリーラインを組み立てる、 ことが説かれています。

ポリオ根絶に向けてWHOでどのような問題の細分化、計画がたてられていたか、詳細を本書の記述だけでは伺い知ることはできません。

ここでは、①1991年の専門家会議で呈示された戦略の内容と、②根絶に至る時系列をたどる、という形をとりたいと思います。

まずは見通しをよくするために②時系列を辿りましょう。本書で書かれている事項をまとめてみました。

f:id:magattaca:20210101233501p:plain
「WHOをゆく」第2章の記述に基づく時系列

上図の通り、いかに資金を確保するか、域内で最も重要な中国でどうやってワクチン投与を成功させるかが鍵となっていたようです。 特に、後者は「中国の成功なしに地域全体の成功はあり得ない」と、言い切られています。

93年に中国で大規模な予防接種を実現していますが、 当時一人っ子政策を取っていた中国で「いかに(政策上、国家が認めていない第2子以降も含めた)すべての子どもを対象とした 予防接種を実現するか?」、は政治的な問題も絡むセンシティブな問題で、 このプロジェクトが公衆衛生に止まらない困難な課題であったことが伺われます。

続いて、①専門家会議で呈示された3つの戦略の内容をみてみましょう。

  1. サーベイランスシステムの確立(途上国含めた地域で感度特異度の向上)
  2. 定期予防接種率の向上(80%)
  3. 特別予防接種デー(ウィーク)

1.は患者をいかに確実に見つけるか?で、これは感染の状況を把握し拡大を防ぐという前向きの重要性に加えて、 "根絶"の確認(ゼロの証明)というゴールにも深く関わっています。

また、2. 3.については本書において「地域全体を「ワクチン漬け」にし、ウイルスの生存の場を奪う」激しい戦略と例えられるほど、 大胆な戦略だったようです。先に見たポリオウイルスの特徴から、なぜこのような戦略が必要だったかがわかりますね。

以上の戦略・時系列から、ポリオ根絶に向けたサブイシューは以下のように分解できるかもしれません。

f:id:magattaca:20210101233844p:plain

2) アウトプットドリブンの観点から

つづいて「イシューから始めよ」では「アウトプットドリブン」ということが言われています。 問題解決では、最終的な結論やストーリの骨格に関わる根源的な課題に重きをおいて検証しアウトプットを得ることが重要です。 そして、得られたアウトプットがロジックに影響するなら、それに合わせて全体を見直し修正していくことが必要になります。

「イシューからはじめよ」では、いきなり問題を解こうとするのではなく、まずは考えることの重要性が繰り返し強調されます。 ですが、ロジックの展開に大きな影響を与える「前提」や「洞察」を考え抜いていても、 実際に取り組み始めると様々なトラブルに出くわします。

よくあるトラブルとして「① 欲しい数字や証明がでない」「② 自分の知識や技では埒があかない」が挙げられています。
これらトラブルの対策として、興味深いことに「足で稼ぐ」や「人に聞きまくる」、ということがあげられています。

私が「イシューからはじめよ」がとても良い本だと思うのは、このように「問題を実際に解決すること」を大前提に設定しており、 目標のない努力を徹底的に批判しつつも、必要な労力は惜しまないスタンスをとっている点です。

少し話が逸れましたが、「WHOをゆく」に話を戻します。

本書で非常に印象的だったのが、尾身先生がためらいなく「足を運び」「人に聞く」方だ、ということです。

WRPOに着任し、ポリオ根絶という壮大な課題を与えられた尾身先生がまず最初にされたこと。。。

  1. WHO天然痘根絶事業で指導的役割を果たした蟻田功先生を、熊本の自宅まで訪問
  2. 当時、ポリオ根絶計画が順調だった汎アメリカ地域事務局の担当者を招き、1週間ぶっ通しで、何をすべきか徹底的に議論

また、先に述べた中国でのワクチン投与では、実現に向けて「直接、中国の保健大臣にお会いし、率直にこの問題を話」しあっています。 さらに<特別予防接種期間>には、自ら中国に赴き、予防接種に参加されています。

で、ここが本当にすごいと思ったのですが、中国での予防接種の際、「関係者には黙って、見学の予定外の接種現場に直接赴き、当の子供やその親に第何子か質問」し、 第2子以降の子供が間違いなく含まれていることを実際に確認されたそうです。

「現場に足を運べ」とはよく言われますが、ここまで実践される方は中々いらっしゃらないのではないでしょうか? ポリオ根絶実現に向けた強い意志を感じるエピソードでした。

f:id:magattaca:20210101234238p:plain

3) メッセージドリブンの観点から

最後に、「イシューからはじめよ」では「いかに人に伝えるか?」という点も重視されています。 どうやって聞き手にメッセージを伝えるか?論拠と構造のしっかりしたストーリーを組み立てることが言われています。 これは交渉力という点にもつながるかもしれません。

この点で本書で非常にとんちの効いた(?)面白いエピソードがありました。

先に、ポリオ根絶プロジェクトにおいて資金獲得に非常に苦労された、ということにふれました。 尾身先生は、「ワクチン購入の金策」のため援助国や公的機関を巡り歩いたことを、 「いわば"営業マン"としての生活」とユーモアを交えて表現されています。

さて、資金を得る上で日本の政府開発援助(ODA)の無償資金協力を得ることが必須となりました。

ですが、「”日本のODAは消耗品には支出しない"」という原則があり、ワクチンは消耗品のためODAの対象にならないと断られてしまいます。

これに対して尾身先生達は次のように答えられたそうです。

「ワクチンは・・・(略)・・・子供に免疫ができることによって、その効果は一生続く”固定資産”になります。 さらにこの戦略が成功すれば、ポリオワクチンそのものが必要になります。橋や建物はいつかは滅びますが、 この財産は永遠に続きます」

こうした議論を重ねることで、最終的にワクチンをODAの対象とすることに理解が得られたそうです。

消耗品のワクチンを、一生の固定資産、永遠の財産としてしまう、この論理とメッセージ。。。私には一生、こんなメッセージをくみたてられるとは思えません。

f:id:magattaca:20210101234606p:plain

まとめ

以上、尾身先生の著作「WHOをゆく」第2章を「イシューからはじめよ」に沿ってご紹介してみました。

できるだけ誤解を生まないように書いたつもりですが、細かい点など端折っているので、ぜひ実際に読んでいただきたいです。

ここでは紹介しきれなかった面白いエピソードがまだまだたくさんありますが、この第2章20ページくらいしかないんですよ!! 本当に圧倒的な内容です。

3000円とちょっとお高いですが、電子書籍もあるそうですので、ステイホームのお正月の読書にぜひぜひ!*7

*1:参考記事1. 『イシューからはじめよ』をまとめてみた
参考記事2. 【図解まとめ】『イシューからはじめよ』を図解で分かりやすく要約

*2:決して根絶確認の検証期間と認定作業を軽視しているわけではありません。

*3:日本WHO協会 HPより WHOはどこにあるの

*4:岡安 裕正 「世界ポリオ根絶イニシアティブの現状と展望

*5:厚生労働省 関西空港検疫所 疾患別解説 ポリオの記載より」

*6:誤用かも

*7:っていうかすべての学校に配りたい。教科書に載せて欲しい。

RDKitで変換可能なペプチド配列と核酸配列について

以前、マクロ分子の表現HELMについて取り上げた記事の中で、ChEMBLに登録されているモノマー表記を調べました*1。その結果、天然のアミノ酸だけでなく様々な修飾アミノ酸が登録されており、略号(ID)で表せるようになっていることがわかりました。

ところでRDKitも同様に略号を使って構造を表すことができます。@iwatobipen先生のこちらの記事に詳しいです。

とても便利な機能ですが、どのような略号が使えるかわからないと、うまく活用できません。そこで今回はRDKitで扱える略号はどのようなものがあるのか調べてみました。

RDKitでの略号からの構造生成例

まずはイメージを掴むために実際に動きを確認してみましょう。

RDKitのドキュメンテーションMolFromHELM()の項目では、「Construct a molecule from a HELM string (currently only supports peptides).」となっています。一方でGitHubのissueで指摘されているように、実際にはペプチド以外にも核酸の配列からMolオブジェクトを作成することができます。

確認してみます。

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole, rdMolDraw2D
from IPython.display import SVG

print(rdBase.rdkitVersion)  
#  2020.03.2

# ペプチドの例
Glutathione = Chem.MolFromHELM("PEPTIDE1{E.C.G}$$$$")

f:id:magattaca:20201231223107p:plain

# RNAの例
RNA_AU = Chem.MolFromHELM("RNA1{P.R(A)P.R(U)P}$$$$")

f:id:magattaca:20201231223151p:plain

# DNAの例
DNA_ACGA = Chem.MolFromHELM("RNA1{[dR](A)P.[dR](C)P.[dR](G)P.[dR](A)}$$$$")

f:id:magattaca:20201231223226p:plain

確かに、アミノ酸の一文字表記に加えてRNA、DNAの核酸の略号にも対応しているようです。

以上はHELMからMolオブジェクトを生成する例です。

ついで、配列(Sequence)からMolオブジェクトを生成する例も確認しておきます。こちらは上でも引用させていただいた@iwatobipen先生の記事で扱われている機能です。

Plan! Do! Check! Act!

# 描画の設定を変更
IPythonConsole.drawOptions.addStereoAnnotation = True
IPythonConsole.drawOptions.annotationFontScale = 1.3

# SequenceからMolオブジェクトを生成
PDCA = Chem.MolFromSequence('PDCA')

f:id:magattaca:20201231223257p:plain

できました!*2

flavor=1という引数を与えるとD-アミノ酸も扱えるようになり、小文字でD体を表せます。*3

# D-アミノ酸の使用例
d_PDCA = Chem.MolFromSequence('pdca', flavor=1)

f:id:magattaca:20201231223348p:plain

くさび表記は見づらいですが、ステレオ表記から立体が反転していることが確認できます。

コードから利用可能な略号を探ろう!

RDKitドキュメンテーションの情報も最新ではない、ということなので、該当と思われる箇所のコードを確認してみましょう。以下がそれっぽい箇所です。*4

どうやらC++名前空間として略号とその表す構造が定義されているようです。

さっぱりわかりません。。。

RDKit UGM_2015 Roger Sayle氏の発表資料*5の記述が、仕組みを理解する上で役に立ちました。

  1. シークエンス情報は、内部でPDBの残基コード(residue number)を使ってエンコードされている
  2. (residue number のような3文字表記は、)1文字表記よりもたくさんの構造を表現できる

これを踏まえてコード(SequenceParsers)を見直すと、確かにPDBで見られる残基、原子の表現を用いながらオブジェクトを作っている感じが伝わってきます。

ここまでHELMとシークエンスの話をごちゃ混ぜに進めてきました。見通しをよくするために、まずはシークエンスからの構造生成に焦点をあてて見ていきたいと思います。

アミノ酸の構築方法

シークエンスを構造に起こすにあたって、1文字表記から各アミノ酸をどうやって構造におこしているのでしょうか?

SequenceParsers 84行目から抜粋します。

具体的な構造構築の情報は下図の真ん中あたりからです。最初のアラニン(ALA)、アルギニン(ARG)の例を見比べていただくと分かりやすいです。

f:id:magattaca:20201231223510p:plain

上図の通り、アミノ酸は以下のように構築しています。

  1. アミノ酸主鎖を構築(CreateAABackbone)
  2. 側鎖を構築(CreateAAAtomで原子を指定し、CreateAABondで結合関係と結合次数を指定)

側鎖の原子の表現にはPDBファイルでよく見かけるγ炭素(CG)、δ炭素(CD)といった表現が使われています。  

注意点としては、ALAの表現から分かるようにβ炭素までCreateAABackboneで構築しています。プログラム実装上の都合と思いますが、通常の主鎖とはズレています。

実際に、グリシン(GLY)は以下のようにCreateAABackboneを使わずに表現されています。

f:id:magattaca:20201231223817p:plain

シークエンスから構築したMolオブジェクトには上で見た原子の情報が含まれています。AtomオブジェクトのGetPDBResidueInfoというメソッドを利用することで取り出せそうです。

# 一文字表記からアルギニンを構築
arg_mol = Chem.MolFromSequence("R")

# 原子の名称を取得しリスト化
atom_names_list = []

for a in arg_mol.GetAtoms():
    a_name = a.GetPDBResidueInfo().GetName()
    atom_names_list.append(a_name)

# 原子名と一緒に描画
view = rdMolDraw2D.MolDraw2DSVG(300,300)
tm = rdMolDraw2D.PrepareMolForDrawing(arg_mol)

# optionで描画にatom labelを追加
option = view.drawOptions()
for i, an in enumerate(atom_names_list):
    option.atomLabels[i] = an

option.padding=0.01
option.legendFontSize=20
    
view.DrawMolecule(tm, legend="arginine with atom labels")
view.FinishDrawing()

svg = view.GetDrawingText()
SVG(svg)

f:id:magattaca:20201231223854p:plain

できました!

C末端のカルボキシル基酸素原子はOXTなんですね。なるほど。。

登録されているアミノ酸

構造構築の裏側がわかってきました。デフォルトではどのようなアミノ酸を含んでいるのでしょうか?

コードにある3文字と1文字略号を表にしてみました。

f:id:magattaca:20201231223932p:plain

約50種のアミノ酸情報がありました。

標準アミノ酸L体、D体の他にもセレノメチオニンやオルニチンといったものもあります。

D体の1文字表記はL体の表記を小文字にしたものとなっています。

コードに不備があるのか、D-セリンD-チロシンは構造生成のための情報が欠けており、利用することができません。
1文字略号syでMolオブジェクトは生成されますが、中身は空です。

d_SER = Chem.MolFromSequence('s', flavor=1)
print(type(d_ser))
#  rdkit.Chem.rdchem.Mol

print([a for a in d_ser.GetAtoms()])
# []

登録済みのヌクレオチド

こちらは種類が少ないので簡単に。

塩基として標準的なACGTUが用意されています。

詳細はコードの中身を見ていただきたいですが、どうやらリン酸-を別に用意し、略号で指定された塩基を結合させてそれぞれのヌクレオチドを構築しているようです。

他の特徴として、末端の5'側3'側、それぞれについてリン酸でキャップした構造(PCap5PCap3)を生成するための情報も組み込まれています。

DNAとRNAのシークエンスの区別は引数flavorで可能です。

シークエンスの種類を区別するflavor

ここまでRDKitで利用可能な略号を確認してきましたが、アミノ酸ヌクレオチドで同じ略号があるため、区別する必要があります。

MolFromSequenceでは、引数flavorの値により、シークエンスがペプチド、DNA、RNAいずれかを区別することができます。

f:id:magattaca:20201231224152p:plain

同じシークエンスを使って違いを確認しましょう。Proteinについては既にみたのでRNA、DNAを試します。

# 5'のみリン酸キャップしたRNA
RNA_test = Chem.MolFromSequence('ACG', flavor=3)

f:id:magattaca:20201231224809p:plain

# 3'のみリン酸キャップしたDNA
DNA_test = Chem.MolFromSequence('ACG', flavor=8)

f:id:magattaca:20201231224827p:plain

図の通り、シークエンスをヌクレオチドとして認識し、それぞれリン酸キャップの違いが反映されたオブジェクトが生成できました。

以上でMolFromSequence()で構築可能な構造の範囲がおおよそわかりました。

先に見た登録済みのアミノ酸の表を再度見直すと、非標準のアミノ酸3文字表記のみで、1文字表記がありません。ですが、MolFromSequnece()は1文字表記の配列が対象なので、これらのアミノ酸を利用できません。

これらのアミノ酸はHELMであれば利用可能です。では続いてHELMからの構造生成を見てみましょう!

HELMで利用可能な略号

HELMではモノマーが複数文字の場合、「角括弧[ ]」で囲むことで一つとして扱えます。したがってシークエンスで扱えなかった3文字表記のアミノ酸も利用可能です。ただし略記(ID)が若干異なるので注意してください。

以下の表にまとめました。L-アミノ酸はシークエンスと同様、一般的な1文字略号でOKなので省略しています。

f:id:magattaca:20201231225024p:plain

末端のキャップを含めて、非標準のアミノ酸を試してみましょう!

peptide_test = Chem.MolFromHELM("PEPTIDE1{[ac].[Orn].[Glp].[seC].[am]}$$$$")

f:id:magattaca:20201231225106p:plain

「N末端: アセチルーオルニチンL-ピログルタミン酸L-セレノメチオニンーNH2:C末端」という構造が生成されました!

HELMは表現の幅を広げやすいのが良いですね!

ヌクレオチドについては冒頭に見た通り、について

  • RNA : リボース「R
  • DNA : デオキシリボース「dR

を使うことに注意すればOKです。

以上で、RDKitに組み込まれている略号は終了です。

メッセージをペプチドにのせて!

さて、RDKitで扱える配列を眺めてきました。短いエントリですがまとめに代えて。。。

2020年は本当に色々なことがありました。来年のことを言うと鬼が笑う、なんて言いますが、最後はこのメッセージを構造式にしてしまいましょう!

#  シークエンスをペプチドに変換
HNY = Chem.MolFromSequence("HAPPYNEWYEAR")

#  描画の設定
view = rdMolDraw2D.MolDraw2DSVG(400,400)
tm = rdMolDraw2D.PrepareMolForDrawing(HNY)

# optionで調整
option = view.drawOptions()
option.legendFontSize=20
    
view.DrawMolecule(tm, legend="Happy New Year")
view.FinishDrawing()

svg = view.GetDrawingText()
SVG(svg)

f:id:magattaca:20201231225218p:plain

良いお年を!

*1:過去記事はこちら

*2:サイクルは回ってませんが

*3:RDKit ドキュメンテーション MolFromSequence

*4:Pistoia Alliance Public Resources HELM Parsers and Writersの記述を参考

*5:p21

コンフォメーションとRadius of gyrationによる評価について

巨大なPROTACがなぜ細胞膜を透過するのか?溶液中におけるコンフォメーションの違いから解析した興味深い文献を読みました。

Solution Conformations Shed Light on PROTAC Cell Permeability

Yoseph Atilaw, Vasanthanathan Poongavanam, Caroline Svensson Nilsson, Duy Nguyen, Anja Giese, Daniel Meibom, Mate Erdelyi, and Jan Kihlberg ACS Medicinal Chemistry Letters Article ASAP DOI: 10.1021/acsmedchemlett.0c00556 (https://pubs.acs.org/doi/10.1021/acsmedchemlett.0c00556)

CC-BY 4.0の元で公開されています。

この論文中で出てきた Radius of gyration (Rgyr) というものがよくわからなかったので調べてみました。

合わせて、この記事では以下の内容を扱っています。

  • PDFファイルからのテキスト情報抽出
  • 学情報のファイル形式xyz形式とその変換

論文の概要

まずは論文の中身についてざっと紹介します。

PROTACは概して分子量が大きく、LipinskiのRule of 5からは大きく外れた構造をとっています。このような巨大な分子が細胞膜を透過する、経口投与で薬効を示す、というのは非常に驚きです。

では何故膜透過が可能なのでしょうか? 

その一つの解として、こちらの論文では、溶液中でのコンフォメーションの違いが寄与しているのではないか?と示唆しています。つまり細胞膜中のような疎水性環境下では、丸まったコンフォメーションをとっており、実質的な体積(大きさ)は小さくなっているのではないか?、ということのようです。

解析されている構造例は以下。グラフィカルアブストラクトがわかりやすいですね。

f:id:magattaca:20201230182328p:plain

以下の解析手法でコンフォメーションを求めているようです。

  • コンフォメーション解析の元となる実験データはNMR
  • NMRデータ解析手法はNAMFIS(NMR analysis of molecular flexibility in solutionアルゴリズム
    • NAMFISでは様々なコンフォメーションの時間平均であるNMRデータを、各コンフォーメションの占有率に分解して解析可能 *1
  • NAMFISに必要なコンフォメーションのアンサンブルはMCMM (Monte Carlo molecular mechanics) conformational searchで生成*2

得られたコンフォメーションを評価する指標として、以下を用いています。

  • SA 3D PSA (solvent accessible 3D polar surface area)
  • Rgyr (radius of gyration)

SA 3D PSAの具体的な計算方法は理解できていませんが、化学の新しいカタチさんの記事「トポロジカル極性表面積は計算コストの低い推定PSA」が参考になると思います。

PSAはオクタノール/水分配比率(logP)などと同様に、分子の極性・脂溶性を示す記述子になります。PSAと膜透過性に関する実験データとの間にはよい相関があることが判明しています。例えばあまりに大きなPSA(140 Å2以上)を有する化合物は膜透過性が低いことが経験的に知られています。(化学の新しいカタチ 上記記事より引用)

得られたコンフォメーション解析結果から、PROTAC 1は疎水性環境中では丸まった構造をとっていることが示唆されました。

このコンフォメーションを安定化させる構造的特徴として以下が指摘されています。

  • 分子内水素結合(IMHB、intramolecular hydrogen bond)
  • 芳香族へテロ環(チアゾール、ピリミジン)間のπ-π interaction
  • ヒドロキシ基とピリミジン間の非古典的水素結合(nonclasical hydrogen bond

IMHBやnonclasical hydrogen bondは、アミドNHやヒドロキシ基を溶媒接触面からマスクすることでもSA 3D PSAを低下させ、膜透過の向上に関与しているのではないか?ということが指摘されています。

この解析結果が非常に興味深いと思ったのは、一見すると2つの大事なユニットをつなぐだけに見えるflexibleなPEGリンカーが、分子内水素結合という形で3次元構造の安定化に関わり、膜透過の向上という物性変化に積極的に寄与している可能性が示唆されている点です。

回転可能な結合数(Nubmer of rotatablebond、NRotB)は極性表面積(PSA)と同様に、一般に膜透過と逆相関することが指摘されている指標のように思います。NRotBを増加させる自由なPEGリンカーが、実は膜透過に大きな役割を果たしているかもしれない、カメレオン的特性のキモ、というのはとても面白い示唆だと思います。

Radius of gyration (Rgyr)

で本題、Radius of gyrationとは何なのか?

Wikipediaによると断面回転半径と訳されるらしく、断面の性質を表すパラメータの一つだそうです。

f:id:magattaca:20201230183314p:plain

何のこっちゃ????

英語版のWikipediaの方がもう少しわかりやすいです(Radius of gyration (wikipedia))。

曰く、「慣性モーメント(moment of inertia)を持つ点までの半径方向距離(radial distance)」です、と。

数学的には「重心、あるいは与えられた軸からの物体の各パーツまでの距離の二乗平均平方根」となるそうです。

・・・なるほどわからん。ポリマー科学におけるIUPACの定義も記載されているので、ポリマー界では一般的に使われる指標なのでしょうか???

Wikipediaの表式は(axis)についてのものっぽいです。重心(center of mass)についての式はこちらのページにありました。

f:id:magattaca:20201230183352p:plain

Rgyrが小さいほど原子が重心に近く、大きいほど遠くなるそうです。つまりRgyrが小さい方が3次元的に小さい(丸まっている?)ということのようです。

とりあえず座標と質量があれば計算できるようです。

ほな、やってみましょうか!!

PDFから座標を抜き出す

先のACS Med. Chem. Lett.論文ではコンフォマーの座標がSupporting Informationで提供されています。・・・PDFで。。。

とりあえずpdfminer.sixを使えばテキストを抽出することができるそうです*3 。pdfminer.sixのドキュメントはこちらで、GitHubはこちらです。

pip install pdfminer.six

一緒にダウンロードされるpdf2txt.pyを使えば簡単らしいです。

私はプログラミングが苦手なので、以降Jupyter notebookで実行しています。

まずPython3でpythonファイルを実行するため、エクスクラメーションマーク(!)をつけて!python3としています。pdf2txt.pyはAnacondaのbinの中にインストールされているようなので相対Pathを使ってpythonファイルを指定しました。

入出力は以下の通りです。

  • 入力ファイル : ACSからダウンロードしたSupporting Information (「ml0c00556_si_001.pdf」)
  • 出力ファイル(-o):「SI.txt

引数--page-numbersを使うと、抜き出したいページを指定できます(ページ番号をスペースを挟んで並べる)が、今回は目的のページ数が多く面倒なのですべて処理します。

!python3 ../../.pyenv/versions/anaconda3-5.3.1/bin/pdf2txt.py -o SI.txt  ml0c00556_si_001.pdf

無事「SI.txt」ファイルができました。該当箇所はこんな感じ。

f:id:magattaca:20201230183724p:plain

ただし、 ページ番号も同時に抽出されるので、座標の合間に改行やベージ番号が挟まっているためそのままでは使えません。 とりあえずコンフォマー1つ分座標があれば良いので手作業で取り除きました。

このままでもRgyrを計算することはできそうですが、折角ですので化学構造の取り扱いに適したファイル形式への変換を試みたいと思います。

座標のXYZ形式とPyMolでの確認

今回の座標の記載方法はXYZ形式に相当すると思います。化学の新しいカタチさんの記事「XYZ形式とZ-マトリックスは分子の立体構造を表す入力フォーマット」に詳しいのでご参照下さい。

上記記事によると、

原子数
コメント
元素記号 z座標 y座標 z座標
元素記号 z座標 y座標 z座標

という記載になっていれば良いようです。

論文の記載より今回の化合物1の組成式はC50H58N9O10F3S となっていますので原子の数は全部で131となります。 この情報を加えて「conformer_1.xyz」というファイルを作成しました。

f:id:magattaca:20201230183857p:plain

xyzファイルが正しくできているか? 確認のためPyMolで描画してみます。

論文本文中のFig. 3とおおよそ同じ向きになるようにして並べてみました。

f:id:magattaca:20201230183934p:plain

どうやらうまく情報を抽出できていそうです。  

xyz形式は座標の情報しか入っていないようなので、結合次数の情報が無く、すべて単結合で表されているようです。・・・なるほど。

Open Babelによるファイル形式の変換と修正

欲しい情報を得られていそうなことがわかりました。ところで私はRDKitが好きなのですが、RDKitではxyzファイルを直接読み込むことができません。

そこで、Open Babelを用いてmolファイル形式に変換してみようと思います。ここではOpen Babelを扱えるPythonライブラリpybelを利用します。Pybelのドキュメントはこちらです。使い方は化学の新しいカタチさんの記事 「Open Babelを使って化学情報のファイル形式を変換」がとても参考になります。

扱えるファイル形式とその拡張子はドキュメントのこちら(Supported File Formats and Options)で確認できます。Jupyter notebook上でprint(pybel.informats)あるいはprint(pybel.outformats)として確認することもできます(入力 or 出力)。

早速、xyzファイルを読み込んでmolファイルに変換してみましょう!

readfile()でファイルを読み込めます。ただしこちらはジェネレーターオブジェクトです。今回は分子一つだけなので、__next__()で中身にアクセスして取り出しています。複数分子の時はforループで取り出せば良いらしいです。

import pybel

single_obmol = pybel.readfile("xyz", "conformer_1.xyz").__next__()

出力はwrite()を使います。今回は1分子なのでそのままでOKです。

single_obmol.write('mol', './conformer_1.mol')

無事、conformer_1.molというファイルができました。

確認のためPyMolで描画してみましょう。先のxyzファイルの描画と並べてみます。

f:id:magattaca:20201230184427p:plain

Open Babelで変換したMolファイルは単結合だけでなく2重結合も見えている一方で、ベンゼン環の結合が切れてしまっています。

どうしてこうなった??

おかしくなっている原子とその周りにAtom Identifiers: indexでラベルをつけてみました。*4。合わせてMolファイルの該当箇所を確認してみます。

f:id:magattaca:20201230184508p:plain

図のように、バグのある原子(Atom index: 47)は、molファイルのアトムリスト元素記号が「*」となっていました。また、結合リストでは47を中心とした結合関係が書かれていませんでした。結合次数が原子価に足りないAtom index 46、48、54はラジカルとして扱われてしまっています。

以上を踏まえて次のように修正してみました(conformer_1mod.mol)。

f:id:magattaca:20201230184552p:plain

修正したmolファイルを再度PyMolで読み込み、修正前と比較してみました。
下図の通り、構造のバグが修正できているようです!!

f:id:magattaca:20201230184619p:plain

修正したmolファイルがRDKitでも読み込めるか試してみましょう。

Molオブジェクトに変換したのち、py3DMolで3次元描画を行います。使い方は化学の新しいカタチさんの記事「py3Dmolを使って化学構造をJupyter上で美しく表示する」に詳しいです。

from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole
import py3Dmol

conf1 = Chem.MolFromMolFile('./conformer_1mod.mol', removeHs=False) # 水素原子も必要

IPythonConsole.drawMol3D(conf1)

f:id:magattaca:20201230184708p:plain

修正したmolファイルは無事RDKitで読み込むことができ、正しい構造を構築できていそうです。

ダメ押しで組成式も確認してみます。

from rdkit.Chem.rdMolDescriptors import CalcMolFormula
formula = CalcMolFormula(conf1)
print(formula)
#    C50H58F3N9O10S

Supporting Informationに記載のデータと同じ組成式が得られました。良かった良かった!!

なお、RDKitでxyzファイルを読み込めるようにするためのプログラムもあるようです。

RDKitによるRadius of gyrationの計算

かなり遠回りしてしまいましたが、ここまでで論文からコンフォマーの座標の抽出と一般的な化学構造ファイルへの変換が終わりました。

本題、Radius of gyrationの計算を行ってみましょう!

PyMolWikiのRadius of gyrationを参考に、RDKitのMolオブジェクトから直接求められる形に修正を試みたいと思います。

まずは重心の座標(r_CM)を求めます。以下の手順でいいはず。。。

  1. Molオブジェクトから各原子の座標の取り出す
  2. 各原子の座標と原子量を掛け算して総和をとる
  3. 分子量で割り算

座標の取得は化学の新しいカタチさんの記事「RDKitによるコンフォマーの生成」を参考にさせていただきました。

こんな感じの関数を作ってみました。

def calc_r_CM(mol):
    from rdkit import Chem
    from rdkit.Chem import Descriptors
    import numpy as np
    
    # 各原子の原子量と座標を準備
    mass_list = []
    coordinate_list = []
    
    mol_conf = mol.GetConformer(0)
    for i, (x, y, z) in enumerate(mol_conf.GetPositions()):
        # 原子量
        atom = mol.GetAtomWithIdx(i)
        mass_list.append(atom.GetMass())
        # 座標
        coordinate_list.append((x,y,z))
    
    # (x,y,z)それぞれに原子量x座標を計算し総和をとる
    r_CM_fraction = np.zeros(3)
    
    for m, (x, y, z) in zip(mass_list, coordinate_list):
        r_CM_fraction[0] += m*x
        r_CM_fraction[1] += m*y
        r_CM_fraction[2] += m*z
    
    # 分子量
    MW = Descriptors.MolWt(mol)
    
    # 重心の座標
    r_CM = r_CM_fraction / MW
    
    return r_CM

コンフォマー1のMolオブジェクトを投げ込んで重心の座標を求めてみます。

print(calc_r_CM(conf1))
# [-0.07220759  0.07596426 -0.02222654]

求まりました。ほぼ原点ですね。

続いて目的のRadius of gyration(RG)を求めます。先の関数と重複が多くなってしまいましたがとりあえずこんな感じ?

def calc_RG(mol):
    from rdkit import Chem
    from rdkit.Chem import Descriptors
    import numpy as np
    
    # 各原子の原子量と座標を準備
    mass_list = []
    coordinate_list = []
    
    mol_conf = mol.GetConformer(0)
    for i, (x, y, z) in enumerate(mol_conf.GetPositions()):
        # 原子量
        atom = mol.GetAtomWithIdx(i)
        mass_list.append(atom.GetMass())
        # 座標
        coordinate_list.append((x,y,z))
    
    # 重心の座標
    r_CM = calc_r_CM(mol)
    
    # 分母は分子量
    MW = Descriptors.MolWt(mol)
    
    # 分子は(原子量 x (原子座標 - 重心の座標)^2 )の総和
    r2_RG_fraction = 0
    
    for m, (x, y, z) in zip(mass_list, coordinate_list):
        r2_RG_fraction += m * ((x - r_CM[0])**2 + (y - r_CM[1])**2 + (z - r_CM[2])**2)
    
    # 回転半径の2乗
    r2_RG = r2_RG_fraction / MW
    
    # 回転半径
    r_RG = np.sqrt(r2_RG)
    
    return r_RG

作成した関数でRadius of gyrationを求めてみます。

print(calc_RG(conf1))
#  5.481507768388569

5.5Åと求まりました。さて正しい計算結果となっているでしょうか??

論文のFig.5 を引用します。

f:id:magattaca:20201230184903p:plain

Conformation 1のRgyrは 5.5Åとなっています。どうやら合っていそうです。良かった!!

PyMolでradius of gyrationを描画

正しそうな結果は出ましたが、実際この数値の意味するものがうまくつかめません。

せっかく座標があるのでPyMolでコンフォマーと一緒に描画してみましょう!

Psudoatomを使えばいけそうです。

先に求めた重心の座標([-0.072, 0.076, -0.022] )とRgyr(5.48)を使って、PyMolのコマンドラインに以下を打ち込みました。

pseudoatom Rgyr, vdw=5.48, pos=[-0.072, 0.076, -0.022]

このままでは座標の点だけなのでShowas meshとしてみました。

f:id:magattaca:20201230192236p:plain

おおっ!!! 

コンフォマーと重なるように球体が描かれました!!!!

確かにこれなら3次元的な構造を表す指標と言われても納得です!!!すげー!!!

まとめ

以上、今回は立体的な構造の指標、Radius of gyrationについて調べてみました。知らない指標だったのでPyMolWikiにも載っていてびっくりしました。

生物界隈でもよく使われる指標なのでしょうか??詳しい方教えていただけると幸いです。・・・あとNAFMISとSA 3D PSAについても教えてください。。。

また、実際に求める過程でxyz形式とその変換も試しました。予想外にOpen Babelで行った変換結果にはバグがありました。 ファイル形式の変換によく使われるソフトと聞いていたので、あれっ!!!ちゃうやん!!!って感じでした。

xyz形式は座標と原子の種類の情報しかないため、情報量が少ないことがバグが生じた原因なのかもしれません。

情報なんてなんぼ合ってもいいですからね!

いずれにせよファイル形式の変換は注意しないといけないということがわかって良かったです。

今回もいろいろと間違いが多そうなのでご指摘いただけると幸いです。

ではでは!

*1:どういうことかさっぱりわかりません。 ref. J. Am. Chem. Soc. 1995(117)1027

*2: Supportgint Informationより、ソフトウェアはシュレーディンガー社のMacromodelのようです

*3:【Python】pdfファイルからテキストを超簡単に抽出する方法

*4:PyMolのAtomラベルのサイズと色はPyMol Wiki:Labelを参考に変えています

「20歳の自分に受けさせたい文章講義」を読んだ

私には2つ悩みがあります。

 

1つ目は「脳内の一人会話が止まらない」こと。お仕事中でも車の運転中でも、場所も時間もわきまえず、くだらないおしゃべりが突然始まり止まらないので危険です。

 

2つ目は読書を趣味と言いつつ「へー、ほー、オモシロ・・・」と、小学生並みの感想ですぐに内容を忘れてしまうことです。

 

で、そんな私が「文章術」に関する本を読んでみました。

 

 

 

どんな本?

こちらはフリーランスライターの古賀史健氏*1が、かつての自分、「20歳の自分」に宛てて書かれた「文書の書き方、組み立て方」に関する本です。

 

「嫌われる勇気」でも有名な著者で、amazonのレビューも評価が高いので読まれた方もたくさんいらっしゃるのでは、と思います。

 

講義形式で全4講となっています。 

f:id:magattaca:20201229110003p:plain

同書「はじめに」をもとに作成

 

各講義のテーマと問いを見ればわかるように、文体や構成といった書く側の技術だけではなく、読者・編集という受け取り手との関わりについても扱われています。

 

これは同書の出発点が学校の国語教育への疑問に端を発していることと関係しています。

 

  •  「"名作の品評会"」ばかりで「どうして誰も"文章の授業"をしてくれないのか?」
  • 「作文がつまらない」のは添削する「先生の目」を基準に「いいこと」を書こうとしてしまうからではないか?

 

書き方の技術を教わらない子ども達の「感じたまま」「思ったまま」の作文は、「先生にほめられる」ためのものとなってしまい、これでは文章の指導ではなく「形を変えた"生活指導"」だ、と著者は指摘しています。

 

では、書く技術はなぜ必要なのか?それは「”考える技術”を身につける」ことだからだ、と著者は述べます。書くアウトプットによって、ものの見方、考え方が変わる。世界の見方が変われば、困難を解決する糸口がみえてくるのではないか?、以上が著者の主張です。

 

お前の悩み文章術関係ないよね?

で、いったいこれが冒頭の悩みとどう関係するのよ??、というお話です。でもその前にもう少し中身について。。。

文章とは何か?

ガイダンス その気持ちを「翻訳」しよう」で、著者は文章について以下のような魅力的な定義をします。

 

  • 頭のなかの「ぐるぐる」を、伝わる言葉に"翻訳"したものが文章なのである

 

書くのではなく”翻訳"する、そう発想することで文章が書けるようになる、というのです。

どういうことでしょうか?

われわれが文書を書く上でぶつかる問題を以下の2点に整理し、それぞれこう説明します。

 

f:id:magattaca:20201229120507p:plain

同書 「ガイダンス」をもとに作成

頭の中を整理してから書こうとするのではなく、「ぐるぐる」をうまく翻訳するのが文章を書くということ。うまく文章にできないのは翻訳の技術が足りていないだけ。まずは翻訳の意識を身につけよう!というのが著者の主張です。

 

ではなぜ翻訳なのか?

 

それは翻訳が、「受け止める相手」を前提とし、「伝えること」を目的とした技術だからです。「ぐるぐる」は言葉にできない。ならば翻訳して伝えるしかない。そう捉え直すのが文章を書く出発点となります。

 

書くことは考えること

ところで、先に著者は「"書く技術"を身に付けることは、そのまま"考える技術"を身に付けることにつながる」と述べていると書きました。この点について「「あー、面白かった」しか言えない人」というタイトルの節中で、読書感想文を例に解説されています。

 

本を読んだ後、読書感想文を書くためには、内容を「自分の頭で整理・再構築し、アウトプット」する必要があります。あいまいとした感想に「言葉を与えなければ」なりません。

 

この面倒な作業を行うことで、読書の理解度は全く異なるものになる。例え内容を理解できていなかったとしても、「自分の言葉に"翻訳"する」ことで自分なりの解を見つけられるようになる。「考えてから書きなさい」ではなく「考えるために書きなさい」。

 

以上が著者の主張です。

 

私のお悩み解決

以上を読んで、浅はかな私はこう考えました。

 

この脳内の独り言を止めるためには、とりあえず文章に書き出せば良いのではないか?、と。ぐるぐる思考を翻訳して吐き出してしまえば、まともな考えとして追い出すことができるのではないか?、と。

 

 私の脳内おしゃべりは、最近読んだ本についてのことが多いです。「なんか面白かったけど、周りに同じ本読んでるおしゃべり相手がいない。もう独り言でいいや。」と、独り言言ってニヤニヤしてます。ヤバいですね。いよいよマスクが手放せません。

 

と、いうわけで、小学生の感想文を抜け出し、脳内おしゃべりを止めるためにも、ちょっとずつ色々な本の感想文をブログに書いていこうと思います。アウトプットすれば理解も深まるそうですし、、、

 

・・・オシマイ

 

・・・では、あまりに内容がなさすぎる!!!ので、もう少し中身について、少し視点を変えて紹介して見たいと思います。

 

リズム・構成・読者・編集

 

文章を書くとはどういうことか?、なんのために書く技術を身につけ文章を書くのか?、という大きな問いについては「ガイダンス」で明らかとなりました。

 

では書くための翻訳の技術はとは具体的にどのようなものか?本書では残りの4つの講義を使って解説されています。ざっと以下の通りです。

f:id:magattaca:20201229134605p:plain

同書 各講のまとめより抜粋して作成

 

以上は、ほぼ本書の目次を抜粋しただけですので、詳細については実際に読んでいただきたいと思います。文章の技術に関してなのに、なぜ映画的なカメラの目線や編集の目線が必要なのか?、といった興味深い内容がたくさんあります。

 

補足として、論理と接続詞については、哲学者 野矢茂樹氏の一連の著作が参考になると思います。著作のエピソードの中で、「学生の時によくできる友人からノートを借りたところ、その友人は板書には書かれていなかった接続詞を自分で書き加えており、そのことでとてもわかりやすい内容となっていた」と書かれていたのがとても印象に残っています。*2

 

とりあえずこの辺が面白かった気がします。。。

 

新版 論理トレーニング (哲学教科書シリーズ)

新版 論理トレーニング (哲学教科書シリーズ)

  • 作者:野矢 茂樹
  • 発売日: 2006/11/01
  • メディア: 単行本
 

 

 

書く技術で書籍の構成を捉え直す

 

文章を書く技術について、本書では4つの視点(リズム、構成、読者、編集)に焦点を当てて解説されていることを見ました。とても勉強になる視点の整理です。この視点、文章単独だけではなく、書籍の構成全体を理解するうえでも役立つのではないでしょうか?

 

早速、この視点を応用して本書で実践されている工夫を捉え直してみたいと思います。

 

手がかりとして、本書が星海社新書に所収されているというから出発します。星海社新書は「次世代による次世代のための武器としての教養」をキャッチコピーに掲げるレーベルです。私は故 瀧本哲史氏による一連の書籍のイメージが強いのですが、恐らく本書と読者層は被っているのではないかと思います。

 

 

武器としての決断思考 (星海社 e-SHINSHO)

武器としての決断思考 (星海社 e-SHINSHO)

 

 

 

武器としての交渉思考 (星海社 e-SHINSHO)

武器としての交渉思考 (星海社 e-SHINSHO)

 

 

 

星海社新書は、次世代武器、といった言葉から予想されるように、若い世代を動かそうとする"熱い"文章が多い印象です。また大きな特徴として、各章の最後に章の内容のまとめがあり、読解を助けるようになっています。さらに、まとめは本文と異なる字体が使われており、本書では板書風の字体、瀧本氏の書籍ではフォントサイズを大きくし、あえて1ページに収まらないようにする、といった他の新書レーベルではみられないアプローチをとっています。

 

そんな"熱い"星海社新書において、本書では読者に届けるためにどのような工夫・仕掛けがなされているでしょうか?

 

先の文章を書く4つの技術(リズム、構成、読者、編集)を使って、以下のように捉え返してみました。

 

f:id:magattaca:20201229154333p:plain

書く技術から書籍の構成をとらえ直した

 

実際のところはわかりませんが、”予備校講師的”な書き方をしているのは、想定の読者層として大学受験を終えたエリート層を設定しているのかもしれませんね。

 

・・・あと何故か著者の一人称がぼく。

 

おふざけっぽくなってしまいましたが、「書く技術の4つの視点(リズム、構成、読者、編集)」は、より大きな書籍の構成を理解するという点においても、見通しよく整理する手段を与えてくれそうです。

 

まとめ

以上、「20歳の自分に受けさせたい文章講義」の読書感想文でした。

 

記事の導入と、後半の話題の一貫性の無さに呆れますね!ひどい文章ですが、吐き出してるとちょっとだけ脳内おしゃべりが収まったので効果があったと思いたい。

 

はじめに」 で著者が書かれているように、私も文章を書く技術についてほとんど指導をうけたことが無いように思います。「20歳」をとうに過ぎてしまいましたが、非常に勉強になることが多くオススメの書籍です。電子書籍にもなっているようなので年末年始のステイホームにいかがでしょうか?

 

あと、一応この記事では著者の主張を踏まえて以下を実践することを目指していました・・・

 

1. 書くことには書くこと以上の効果がある

    → とりとめない(ぐるぐる)思考、脳内おしゃべりを止めるのにも役立つか?

 2. 文書の構成は図解・可視化して「眼」で考えよ

 → 講義らしくスライド形式の図にまとめ、「→」を利用する

3. 仮説を提示し、一緒に”検証する”

 → 書く技術の射程を書籍の構成の理解に広げる

4. 読者の椅子に座る

 → 脳内おしゃべりがとまらなくて、一つの話題だけで文章を書けない人

  ・・・そんな人、私以外にいるのでしょうか???

 

ではでは!*3

*1:こが ふみたけ さんです。コガシ ケンさんだと思ってました。

*2:記憶違いでしたらすみません

*3:色々と怒られそうな内容ですがアフィリエイトじゃ無いから許してください。

ChEMBLのHELMモノマーライブラリーを解析した話 〜XMLをDataFrameに変換して部分構造検索〜

マクロ分子の表現HELMについて調べています。HELMは階層的な構造をとっており、モノマー (ex. アミノ酸)を組み合わせてポリマー (ex.ペプチド)を表現します。

HELMの特徴はその表現の拡張性にあり、モノマーライブラリーにオリジナルなモノマーを追加することで非天然な構造を表現することもできます。

一方で、HELM表現は略記(ID)を用いるので、モノマーライブラリーが共有されていない場合、同一のIDで異なるモノマーを指定する危険性があるため、モノマーライブラリを把握しておくことが重要となりそうです。

そこで今回はHELMも扱っているおなじみのデータベース、ChEMBLでどのようなモノマー情報が格納されているか、探ってみたいと思います。

具体的にはXMLファイルで提供されているモノマーライブラリーを読み込んでPandasDataFrameに変換することを目指します。

つまりXML良くわからない私の単なるお勉強メモ!!

ChEMBL HELM モノマーライブラリーの取得

ChEMBL のHELM モノマーライブラリーをはHELM表記のある化合物のCompound Report Cardから取得することができます。

例えば、GnRH(性腺刺激ホルモン放出ホルモン)アンタゴニストであるアバレリクス(Abarelix)の場合、Compound Report CardにおけるHELM表記は以下のようになっています。

f:id:magattaca:20201125000223p:plain

図の真ん中あたりの「here」 というところをクリックするとxmlファイルのダウンロードが始まります。

2020/11/23現在「chembl_27_monomer_library.xml」というファイルが提供されています。

XMLの中身

XMLについてよくわかっていないので、順番にファイルの中身を見ていきたいと思います。

話がそれますがこちらの記事(Macs in ChemistryDetermining the Amino Acids in a collection of peptides」)でChEMBLモノマーライブラリーXMLの解析事例が紹介されているので、正しい情報が欲しい方はご参照ください。

閑話休題

とりあえずテキストエディタで開くと中身はこんな感じの構造になっています。

f:id:magattaca:20201125000250p:plain

HELMの記載方法同様、モノマーライブラリも階層的な構造となっているようです。

  1. <PolymerList>
  2. <Polymer>
  3. <monomer>

という順にタグが入れ子構造(ネスト)になっており、<monomer> ~ </monomer>ごとに各モノマーの情報が入っているようです。

monomerに記載されている情報は、Pistoia Alliance HPのHELMの解説におけるMonomerと見比べてみると分かりやすいですね。

f:id:magattaca:20201125002237p:plain

上図のテーブルに該当するタグがXMLに含まれているようです。

それでは具体的にどのようなモノマーが含まれているのか、その総数や各モノマーの特徴解析を進めていきたいと思います。

・・・・が、早速冒頭2行から意味がわかりませんでした。

XML宣言と名前空間

ファイルの冒頭2行、以下のように記載されています。

<?xml version="1.0" encoding="UTF-8"?>
<MonomerDB xmlns="lmr" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

1行目はXML宣言というもので、使用するXMLのバージョンや文字コードについて記述してあるそうです。*1

この場合バージョン1.0文字コードUTF-8となっているようです。

2行目、タグの中にある「xmlns="lmr"」という箇所は、XML名前空間を表しているようです。*2

XMLマスターというこちらのページの説明が分かりやすかったです。

XML名前空間とは、XML文書内にある「定義内容は異なるが同じ名前の要素名または属性名」を区別し、名前の衝突を回避するための仕組みです。(XMLマスターポイントレッスン 〜ベーシック編〜 第11回 XML名前空間 より引用)

ここでは接頭辞(prefix)を使わずに名前空間宣言が行われているため、デフォルトの名前空間に相当する記載となっています。

従って、MonomerDBという要素に対して、lmrをデフォルトの名前空間識別子(URI: Uniform Resource Indentifier)として宣言していることになり、 基本的には全てのタグ{lmr}タグという形に展開されることになります。

一方で、2行目後半「xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"」 という箇所はXMLスキーマ定義に関わっているようです。xsiという接頭辞にXML Schema instanceをURIとして宣言しています。*3

どういうこっちゃ?

ElementTreeでXMLを解析

「入門 Python3」によると、XMLの解析にはとりあえずElementTreeを使えば良いとのことなので、試してみます。

ライブラリをインポートした後、まずはparseでファイルを読み込みます。その後getroot()でルートとなるノードを取得して処理を行えば良いそうです。*4

import xml.etree.ElementTree as ET

tree = ET.parse('./chembl_27_monomer_library.xml')
root = tree.getroot()

タグはtag、属性はattribに入っている、とのことなので確認してみます。

print(root.tag)
 #  {lmr}MonomerDB

print(root.attrib)
#  {}

先に見た通り、デフォルトの名前空間lmrとなっているため、タグのURIとして{lmr}が付与されています。

また、MonomerDBタグは名前空間の宣言以外の属性がないのでattribは空となっているようです。

子ノードを取得してみます。*5イテレート可能、とのことなのでforループでタグを表示させてみます。

for child in root:
    print('タグ: ', child.tag, '属性: ', child.attrib)
#  タグ:  {lmr}PolymerList 属性:  {}
#  タグ:  {lmr}AttachmentList 属性:  {}

PolymerListの他にAttachmentListというタグがあるようです。

PolmerLystのもう一つ下の階層にアクセスします。インデックスを利用してみます。

polymerlist = root[0]
print(len(polymerlist))
for child in polymerlist:
    print(child.tag, child.attrib)
#   1
#  {lmr}Polymer {'polymerType': 'PEPTIDE'}

PolymerTypeはPEPTIDEのみでした。HELMには他に RNASACCHEMといったPolymerTypeがありますが、ChEMBLではPEPTIDEのみしか含まれていないようです。

モノマーライブラリーがペプチドのみ、ということがわかったので含まれるモノマーは全てアミノ酸(誘導体)となっていると思われます。

いくつのモノマーが含まれているでしょうか?名前空間URIに気をつけて{lmr}Monomerの全ての要素をfindall()で取得します。

polymer = root[0][0]
monomer_elements = polymer.findall('{lmr}Monomer')
print(len(monomer_elements))
# 2851

全部で2851個のモノマーが登録されているようです。PEPTIDEのモノマーってそんなに種類あるんですか???

上記はindexで階層を辿っていますがXPath機能を使うこともできるそうです。.//を使うことで木全体から要素を選択することができます。*6

monomer_elements_2 = root.findall(".//{lmr}Monomer")
print(monomer_elements == monomer_elements_2)
#  True

先に取り出したものと一致することが確認できました。

モノマー情報の取得

ElementとしてMonomerを取り出すことができました。その中身を表示して見ます。

タグで囲まれた内容は.textで確認することができます。

mon = monomer_elements[0]
monomer_tag_list = []

for m in mon:
    monomer_tag_list.append(m.tag)
    print(m.tag, m.text)

#  {lmr}MonomerID X2484
#  {lmr}MonomerSmiles CC(=O)SCCCC[C@H](N[*])C([*])=O |$;;;;;;;;;;_R1;;_R2;$|
#  {lmr}MonomerMolFile H4sIAAAAAAAAAI2TTU7DMBCF9z7FSHTLaMb/swbUFQW1iAtUCPX+F2BsN7ZTUYgVRe/Fn1/G48QAnM6Xj6/vy5mYLQcmS/bZGGAP7ACoXlzvY4gIfFoiMsUxOpVNSWiKkMqzJ2gRNxxh+oMjtJFiUxzITdzbzFmMtEoZ6jRzDmXTewPae3krLmLYlJcwlUb9mydInbuTp5tdxFgTc+uCDZSnNYc5OyOFhcte/MQdHwbHhG5LrYN75NturfIYvZPryQWR9PsZct0XN6PKtX024/uMqtCNqtiNqtSNqtyNLCZWw7TMqGKeMO4VqCqf/DD+Ws4rwHH/3mplqEQBwZaZl4P+Ljsd5gd2o4VdTQMAAA==
#  {lmr}MonomerType Backbone
#  {lmr}PolymerType PEPTIDE
#  {lmr}NaturalAnalog X
#  {lmr}MonomerName X2484
#  {lmr}Attachments 

モノマーを表す特徴として8つの情報が含まれています。先に図で示したXMLファイルの中身と同様の内容がtextで取り出されているようです。

最後、Attachmentは更に深い階層に情報が書かれています。一つ一つ階層を辿るのが面倒なので手っ取り早くiter()での処理を試してみます。

Elementオブジェクトはiter()で現在の要素を根とする木のイテレータが作成できるそうです。*7

Att = mon.find('{lmr}Attachments')

for Att_el in Att.iter():
    print(Att_el.tag, Att_el.text)

#  {lmr}Attachments 
                        
#  {lmr}Attachment 
                            
#  {lmr}AttachmentID R1-H
#  {lmr}AttachmentLabel R1
#  {lmr}CapGroupName H
#  {lmr}CapGroupSmiles [*][H] |$_R1;$|
#  {lmr}Attachment 
                            
#  {lmr}AttachmentID R2-OH
#  {lmr}AttachmentLabel R2
#  {lmr}CapGroupName OH
#  {lmr}CapGroupSmiles O[*] |$;_R2$|

2箇所のAttachmentについてそれぞれの下の階層の情報まで取得することができました。

IDLabelCapGroupの情報が含まれているようです。ここでもSMILESはCXSMILESのようですね。

DataFrameへの変換

Monomer情報の取得ができるようになったので、より扱いやすくするためモノマーライブラリーをDataFrameに変換してみたいと思います。

具体的にはRDKitのPandasTools*8を利用して、各行に各Monomerの情報Molオブジェクトを含むDataFrameを作成します。作成したDataFrameは部分構造検索などが利用できるので、簡易的なデータベースとして使えるはず!

まず、8つの特徴のうち、MonomerMolFileとAttachmentsを除いた以下の6つを対象として、DataFrameのカラムとします。

  1. MonomerID
  2. MonomerSmiles
  3. MonomerType
  4. PolymerType
  5. NaturalAnalog
  6. MonomerName
monomer_tag_list.remove('{lmr}MonomerMolFile')
monomer_tag_list.remove('{lmr}Attachments')

# 見やすさのため{lmr}を取り除いて列名とする
column_names = [s.replace('{lmr}', '') for s in monomer_tag_list]

必要なライブラリを読み込みます。

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, rdDepictor, PandasTools
from rdkit.Chem.Draw import IPythonConsole
import pandas as pd

# DataFrameの準備
df = pd.DataFrame(columns = column_names )

モノマーごとに取得した情報をリストに格納後、pd.Seriesに変換して、DataFrameに追加していきます。

for m in monomer_elements:
    data = []
    
    for tag in monomer_tag_list:
        e = m.find(tag)
        data.append(e.text)
    
    s = pd.Series(data, index=column_names)
    df = df.append(s, ignore_index=True)

print(df.shape)
#  (2851, 6)

2851個のモノマーを含むDataFrameが作成できました!

中身を確認してみます。

df.head()

f:id:magattaca:20201125001505p:plain

きちんと情報が入っていそうです。

RDKitのPandasToolsを使って各行のMonomerSMILESからROMolオブジェクトを作成します。

PandasTools.AddMoleculeColumnToFrame(df, 'MonomerSmiles')
df.head()

f:id:magattaca:20201125000925p:plain

構造式データを含むDataFrameが無事作成できました!

DataFrame内の構造式検索

ではPandasToolsで作成したDataFrameをが構造式データベースとして利用可能か? 部分構造検索等、試していきましょう!

例としてチロシンの誘導体にどのような構造が含まれているか検索してみます。

まずは簡単にNaturalAnalogとしてYが設定されているものを抜き出してみます。

Y_derivatives = df.query('NaturalAnalog == "Y"')

print(Y_derivatives.shape)
# (176, 7)

176個のチロシン関連構造が含まれているようです。L-チロシンそのものの情報を確認してみます。

Y_derivatives.query('MonomerName == "Tyrosine"')

f:id:magattaca:20201125001811p:plain

L-Tyrosineがどのように格納されているかがわかったので、この構造情報をつかって部分構造検索を行ってみます。

元のDataFrameから「>= (ge比較演算子)」を使って検索します。

# L-tyrosineのROMolオブジェクトをDataFrameから取得
L_ty = df.iat[183, 6]

# ge比較演算子で部分構造マッチ
L_ty_sub = df[df.ROMol >= L_ty]

print(L_ty_sub.shape)
# (77, 7)

L_ty_sub.head()

f:id:magattaca:20201125001844p:plain

L-tyrosineを部分構造として含むものは77個見つかりました。先にNaturalAnalog == Yとしたものは176個でしたので、チロシンのアナログのなかには元の部分構造自体にも変化が加わっているものがあるようです。

では逆にL-tyrosineと部分構造が一致したもののNaturalAnalogは全てYとなっているのでしょうか?

Analog_vc = L_ty_sub['NaturalAnalog'].value_counts()
Analog_vc

# Y    72
# X     5
# Name: NaturalAnalog, dtype: int64

NaturalAnalogXとなっているものも5つあるようです。

L_ty_sub_X = L_ty_sub.query('NaturalAnalog == "X"')

L_ty_sub_X.head()

f:id:magattaca:20201125003153p:plain

NaturalAnalogXとなっていたものは、L-tyrosineと一致する構造をもちつつ、ビシクロ環構造となっているもののようです。

こんな構造もモノマーに入ってるんですね。予想以上に非天然です。。。

試しに一番下のモノマー「MonomerID X316」を含む構造にどのようなものがあるかChEMBLで検索してみました。

f:id:magattaca:20201125003334p:plain

左から2番目、化合物「CHEMBL375661」の情報元の文献は以下でした。

Spatial Conformation and Topography of the Tyrosine Aromatic Ring in Substrate Recognition by Protein Tyrosine Kinases
J. Med. Chem. 2006, 49, 6, 1916–1924

この部分構造が含まれているのも納得なタイトルですね!

まとめ

以上、今回はChEMBLデータベースが提供しているHELMのモノマーライブラリーの中身を解析してみました。

アミノ酸の誘導体だけで3000個近くものモノマーがあるというのはびっくりです。

また、XMLの解析の練習をかねてDataFrameへの変換を行ってみました。RDKitのPandasToolsを使うことで部分構造検索もできるので、うまく使うことができれば簡単なデータベース感覚で利用できそうです。比較演算子が使えることがここまで便利とは思いませんでした。もうちょっと良い使い方があるか考えてみよう。

今回も色々と間違いが多そうです。ご指摘いただければ幸いです。

*1:Let'sプログラミング 「XML宣言を記述する

*2:ELementTree XML API名前空間のあるXMLの解析

*3:都元ダイスケ IT-PRESS 「XMLのスキーマの定義の仕方(2)

*4:参考 LIFE WITH PYTHONライブラリ:ElementTree

*5:ElementTree XML API

*6:ELementTree XML APIXPath サポート

*7:ELementTree XML APIElement オブジェクト

*8:化学の新しいカタチ 「RDKitのPandasToolsでデータ分析を加速する