先日の記事で、下記文献を例に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に適したｘｙz形式の変換方法
PyMolのコマンドと自作関数の使い方

SA 3D PSAについて
入力構造の準備（シェルスクリプトの練習）
XYZ2molでRDKit Molオブジェクトに一発変換
3D PSAの計算
- 部分電荷の計算に失敗
  - WinMopac
  - xTB
- 座標のみで計算
  - PyMolで使用するスクリプトの準備
  - 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)は、プローブの半径を0Åとしたもの
→ 3D PSA
溶媒接触表面積(solvent accessible surface area)は、プローブの半径を1.4Åとしたもの
→ SA 3D PSA

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

また、極性原子は

窒素原子(N)
酸素原子(O)
それらに結合する水素原子(NH、OH)
閾値よりも大きな部分電荷(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コマンドを使えば良いとのことですが、MacのsedはBSD版というもので、Linuxのsed（GNU版）と挙動が異なるようです。

後者の方が使い勝手が良いらしいので、GNU版sedをインストールしてみます。*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)」を削除しなければいけないことに気づくまでに２時間くらいかかりました。。。。
何度gsedを実行しても、テキストエディタを変えて確認してもうまくいかない理由がわからず、結局「Pages」で開いて改ページが含まれているに気づきました。*3

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

さて、気を取り直して！

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

ここでもMacでGNU版の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を使いましたが、出力構造に難がありマニュアルで修正する必要がありました。

今回はxyz2mol（GitHub）を使ってみたいと思います。
こちらは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計算に挑戦したいと思います。

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

各原子のxyz座標
各原子が極性原子に該当するか（原子タイプと部分電荷）？

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

部分電荷の計算に失敗

WinMopac

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

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

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

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

WinMopac 7.21
Wine4 (Mac上でWindowsのソフトを使えるようにする。Wine環境の準備)
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ファイルを入力として使うこともできます。また、--sccはself-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つのことが行われています。

PyMol Python APIを使う準備(cmdモジュール)
3D PSAを求めるための関数の定義(psa3d)
定義した関数を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'で選択された原子(群)の表面積（Å²）を計算

これを踏まえると、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ファイルを使いました。

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

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

分子の表面積（M 3D PSA）「set dot_density=4; set dot_solvent=0; psa3d」
溶媒接触表面積(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）