ゴミのやま

完全自分向け雑記

RDKitで化合物のコンフォマーを作成する

背景

ある分子のコンフォマー分布を調べたくなった。しかしMDをやるほど厳密性が要求されているわけでもないし, まずLAMMPS, AMBER, GROMACSといった分子動力学計算ソフトをわざわざこのためにセットアップする気力がない。

-> RDkitでコンフォマー作成すれば楽なのでは？

結論

素晴らしい先行者様がいらっしゃるのでこちらを完コピして何をしているか学習する。

RDKitを使った低分子のコンホメーション探索 – 香川大学農学部 ケミカルバイオロジー研究室

 

流れ

  1. SMILESなどから構造取得
  2. 水素原子付加
  3. 任意の数のコンフォマー生成 (AllChem.EmbedMultipleConfs)
    (ETKDGでは分子力場による最適化は必須ではないので今回はしない）
  4. AllChem.MMFFGetMoleculePropertiesで水素原子付加語の分子からmolecule propertiesを作成
  5. AllChem.MMFFGetMoleculeForceField()で前述のpropertyを元にコンフォマーごとの力場データを作成（？）。エネルギー計算
  6. AllChem.AlignMolConformers(m)でコンフォマーをアラインメント
  7. エネルギー順にコンフォマーを並び替え
  8. (ここはなくてもよさそう）Butina.ClusterDataでデータをクラスタライズ
    クラスタライズはRMSとTFDfingerprintに基づいて行う。
    GetConformerRMSMatrixは全コンフォマーペア(nC2の組み合わせ)のRMSを含む下三角行列を返す。TorsionFingerprintsは全コンフォマーペアのTFDの情報を含む行列を返す。どちらも直接Butina.ClusterDataに渡すことが可能。
  9. コンフォマーをsdf形式で書き出す。

調べたこと

ETDKGとは

DG法に化学的な常識と二面角データを追加したもの、らしい。 

future-chem.com 

DG法とは

  1. The molecule’s distance bounds matrix is calculated based on the connection table and a set of rules.

  2. The bounds matrix is smoothed using a triangle-bounds smoothing algorithm.

  3. A random distance matrix that satisfies the bounds matrix is generated.

  4. This distance matrix is embedded in 3D dimensions (producing coordinates for each atom).

  5. The resulting coordinates are cleaned up somewhat using a crude force field and the bounds matrix.

https://www.rdkit.org/docs/GettingStartedInPython.html

なんのこっちゃら。分子構造に基づく制限の下でランダムな構造を生成し分子力場でclean upするという理解。

Torsion Fingerprint Deviation (TFD)とは

RMSDのTorsion版のような認識。RMSDと比べて1. 最大値が1で規格化されている, 2. 計算に分子同士の重ね合わせを必要としないなどの利点あり。

NGL viewとは

MDなどのトラジェクトリを表示するpythonプラグイン。

https://ambermd.org/tutorials/analysis/tutorial_notebooks/nglview_movie/index.html

PubChempyとは

PubChemにアクセスするためのPythonプラグイン。PubChemは化合物自身あるいは薬価などの情報を含む。

Modellerでループ構造最適化, 束縛条件追加などを行う

背景

前回の記事でホモロジーモデリングの全体の流れを実施した。この記事ではより詳細なModellerのオプションを調べてみる。

参考：https://salilab.org/modeller/tutorial/advanced.html

オプション

automodelの精度は？

sequence identiy <50%の場合はalignment.check()で残基ごとの妥当性のチェックをしたほうがよいとのこと(DOPEが具体的に何を計算しているかは知らない）。

面白そうな機能

複数のテンプレートを使用

複数のテンプレートを使うことでループ領域の精度（DOPEスコア）が改善されたり、場合によっては改悪されてしまったりするらしい。長いループ構造を再現する際に使ってみようかと思う。

ループ構造最適化

ホモロジーモデリングで得た構造を元にループ構造の最適化を行う。もっとも精度が出る方法は数百個のconformersを作成してクラスタリングし, 代表的なものをピックアップすることらしい。

リガンドを結合させた上でホモロジーモデリング

リガンドを含む構造をホモロジーモデリングで生成できるらしい。DOCKやAuto Dock等のドッキングソフトを使えばいいんじゃないの、と思う気もするがなにか利点があるのだろう。

 

 

HHpredとModellerでホモロジーモデリングをする

前提

  • HHblitsとHHsearchにより類似構造の配列は既に入手済み。ここからスタートする

まずModellerのチュートリアルを実施

1. Searching for structures related to the query

ここは前述のHHsearchによりテンプレート構造が既知なので省略。

2. Selecting a template

ここも省略

3. Aliging the query with the template

ここも省略

4. Model building 
Pirファイル生成

モデル構築のためには.pirファイルをhhsearchの結果から作ってやる必要あり。それにはhhmakemodel.pyを使うらしい。hhmakemodel.pyはhhsuitesの付属品。以下はhhmakemodel.pyの-hオプションのコピペ

Creates a MODELLER alignment (*.pir) from a HHSearch results file (*.hhr).

positional arguments:
FILE results file from HHsearch with hit list and alignment
DIR path to the folder containing cif files
FILE output file (PIR-formatted multiple alignment)
DIR path to the folder where modified cif files should be
written to

optional arguments:
-h, --help show this help message and exit
-v, --verbose verbose mode
-m INT [INT ...] pick hits with specified indices (e.g. -m 2 5)
-e FLOAT maximum E-Value threshold (e.g. -e 0.001)
-r FLOAT residue ratio (filter alignments that have contribute at
least residues according to the specified ratio).
-c convert non-canonical residues (default = True)

 

上記を参考にしたコマンド例： 

hhmakemodel.py input.hhr dir_to_ciffiles output.pir ./

それぞれの意味は前述のヘルプのpositional argumentsを参考にしてください。ちなみにcifはpdbからwgetでとってきました。

Modeller 実行

チュートリアル上のpythonファイルを適宜いじって実行するだけ。当然だがsequentialな処理には専用のスクリプトを書く必要あり。

from modeller import *
from modeller.automodel import *
#from modeller import soap_protein_od

env = environ()
a = automodel(env, alnfile='TvLDH-1bdmA.ali',
knowns='1bdmA', sequence='TvLDH',
assess_methods=(assess.DOPE,
#soap_protein_od.Scorer(),
assess.GA341))
a.starting_model = 1
a.ending_model = 5
a.make()

以下わからなかったので調べたこと

  • Uniclust databaseとは？

    > The Uniclust90, Uniclust50, Uniclust30 databases cluster UniProtKB sequences at the level of 90%, 50% and 30% pairwise sequence identity. The clusterings show a high consistency of functional annotation owing to an optimised clustering pipeline that runs with our MMseqs2 software for fast and sensitive protein sequence searching and clustering.

    UniprotKBにあるタンパク質配列をクラスター化したもの。Uniclust90なら90%の配列一致度を持つ配列同士で, Uniclust30なら30%の配列一致度を持つ配列同士をクラスター化している。（そしておそらくクラスター中心か何かを抽出している）。性質上, クラスターができやすいUniclust30が一番データベースのサイズが小さい。

  • PDB70とは
    >PSI-BLAST alignments produced with sequences of PDB full chain representativies (<70% sequence identity) as queries.
    とのことなので, PDBから（70%以下の配列一致度を持つ配列同士をクラスター化して圧縮した）配列をクエリとしPSI-BLASTをかけた結果と思われる。Uniclustとの違いはこちらはPDBの配列を使っているために, 構造が既知である。これをmodellerの鋳型に使う。

  •  ホモロジーモデリングにおいて前述のデータベースはどのように使われる？
    Uniclustは与えた配列からhhblitsで類似配列群を隠れマルコフモデルを使って探す際に用いられる
    PDB70は前述のhhblitsの結果を元にhhsearchで類似構造を探索する際に使う

  • HHsearchとHHblitsの違いは？
    HHblitsは入力配列に類似する配列を反復検索し, 類似配列群を作成する
    HHsearchはHHblitsの結果を用いて相同な構造を見つけ出す

次にやりたいこと

  • 任意の束縛下におけるホモロジーモデリング
  • Modellerが具体的に何をしているか論文を読む
  • Modellerのオプション検討（デフォルトで十分か？）

NCBIからFASTAを取ってくる

以前UniprotからFastaを取ってくる方法の記事を書いた。UniprotではURLを適宜いじることでFastaを直接ダウンロードできたが, NCBIで同様の方法は使えない（ようだ）。

具体的にはBiopythonのEntrezを使う。詳細は以下参照。

https://www.biostars.org/p/17715/

Uniprotから任意のFASTAをとってくる

やりたいこと

Uniprotから任意のアミノ酸配列をFASTA形式でPythonなどのプログラミング言語を使うことにより自動取得したい。

参考

以下のサイトでまずUniprotとは何かを理解した。アミノ酸配列とその機能情報を掲載しているらしい。

magattaca.hatenablog.com 

やりかた

任意のデータにはアドレスを改変することによりアクセスできる。詳細はUniprotを参照。

Programmatic access - Retrieving entries via queries

 

今回は

  • マウス（分類番号10090）の
  • reviewedの状態にある（Uniprotではcuratorが結果をチェックしたものをreviewedと呼ぶらしい）
  • hsp90ab1という名前の遺伝子に
  • fasta形式で

アクセスしたい。この場合の書式は以下のようになる。

https://www.uniprot.org/uniprot/?query=reviewed:yes+AND+organism:10090+AND+gene:hsp90ab1&format=fasta 

 前述のUniprotでは遺伝子名を指定する際に"genes:"を使えと言っているが, "gene:"を使わないと参照できなかった。

HHsuite覚書

目的

  • 仕事でホモロジーモデリングをする必要が出てきた。それにはhhsuiteなるものが必要らしい。よくわからないので入力／出力および特性を理解する。アルゴの詳細には深入りしない。

  • 同じようなソフトウェアは他にも有るはずなので, それらのソフトウェアとの違いについても後で調べる。

  • とりあえず先に手を動かしてみる 

参考

この記事は開発元のサイトからほぼ引用したものとなる。hhsuiteはミュンヘン大のグループが開発している。余談だが欧州はウェットの実験に関する規制がきついのでバイオインフォが盛んと聞いた（デマかもしれない）。

github.com

 

あとは日本語の文献をとりあえずスキミングして外観を掴んだ。ありがてえ 

qiita.com

要約

導入

インストール

  • Anacondaで専用の環境を作ってインストールした。gitの方針に従うだけなので難しくはない

    github.com

 データベース

  • 色々有るが先例にならってpdb70とuniclustをダウンロードした。ダウンロードは非常に遅い。私はリモートサーバにsshで接続しているため, 自動切断されても良いようにscreenを立ち上げてからlftpでダウンロードした。

  • それぞれのデータベースの特色はわかったら追記する 

速度向上のために

  • HHsuiteは頻繁にデータベースにアクセスするため, ファイルアクセスがボトルネックとなるケースが多い。主に以下の対策推奨されているようだ。私の使っているPCはRAMが128GBあるため, RAMに置く方法を取ろうかと思う。
  1. RAMにデータベースを常に配置しておく。
  2. SSDにデータベースを置く

例：ホモロジーモデリング

詳細は開発元のページを参照。以下には大まかなステップだけ書く

  1. hhblitsかhhsearchでUniclust30からクエリと類似度の高いタンパク質"配列"を探索

  2. 1の結果を元にしてPDB70データベースから類似度の高いタンパク質"構造"を探索

  3. 2の検索結果はE-value（重要度を示すスコアらしい。詳細は調べていない）順にソートされている。一番上の構造をテンプレートとして用いる。

  4. hhmakemodel.pyを使って, 3で見出したテンプレートをMODELLERが読める形式(.pir)にコンバートする

  5. MODELLERでホモロジーモデリングする（これ以降はMODELLERのドキュメントを参照する必要あり）

よくある質問集（興味あるトピックだけ読んだ）

hh-suiteはいつ使ったらいいの？

  • 調べたいタンパク質のorthologが十分研究されておらずBLAST, FASTA, or PSI-BLASTといった従来の手法が上手く機能しない場合, HHsuiteはより遠縁のタンパク質同士の関係性を推測できる

HomologusとAnalogusなProteinの違い

  • Analogus: Convergent evolutionのように一見似た構造を持つが進化的に関係がない場合（例：コウモリと鳥）
  • Homologus: 似た構造を持ちかつ進化的に関係がある場合

Homologusの程度をどうやって判定するのか

probabilityかE-valueを見る：
  • HH-suiteは20%以下のsequence identityでも関係性を見いだせる。逆に言えばこれはhomologusの判定に適切ではない
  • 最も良いのはprobabilityを使うこと（算出方法は知らない）。これが95%を超えるならhomolgyがあるのはほぼ確実といえる。(1)>50%のprobability, または(2)上位3hitが>30%のprobabilityを持つならhomologyを持つ可能性を検討したほうがいい
  • E-valueはqueryに関係ない情報のみから成るデータベースを用いた場合, E-valueが今回の計算結果より良くなる可能性を示す（よくわからない…が低いほうがいいと思われる）。値が1より低い場合は重要な結果と考えられる。probabilityとは異なり, E-valueはsecondary-structureを考慮していない。従って, probabiltyのほうがより正確といえる
hitの妥当性を確認する

hitはありそうな生物種のものか？類似した機能を持つか？などを検討する

secondary-structureは似ているか？
top hitsの関係性を調べる
motifの類似性を調べる
よりカバー範囲の広いデータベースを用いる
パラメータを変えてシミュレーションする
実験検討する（究極）

ドメインが複数あるタンパク質にも使える？

使えるけど偽陽性率が上がる。小さなドメインに分割してサーチするのも手

HHblitsとHHsearchは何が違うのか

Due to its fast prefilter, HHblits runs between 30 and 3000 times faster than HHsearch at the cost of only a few percent lower sensitivity.

とのことで, ケースにもよるが大量のシミュレーションをするときはHHblitsがいいかもしれない。逆に精度が重要な際はHHsearchを使おうと思う。

 

Linux mint 導入時にやったことメモ

Ubuntuくんが激重なのでLinux mint (Cinnamon)に乗り換えたのメモ。完全自分用

とりあえず結論

めっちゃ快適なのでさっさとやればよかった

とりあえずアップデート

インストール後表示されるFirst Stepsに従ってドライバをインストールしたりアップデートしたりする

キーバインド周り

いつものキーバインドじゃないと作業効率が落ちるのでまずここから設定する。

Synapseを入れる

アプリ名で検索できる便利なランチャー。これがないとストレスでしぬ

Xmodmapを入れる

CapsをデリートにしたりVim風のカーソル移動を実現するために。詳細はそのうち書くかも。これによりタイプ速度が2倍になる（当社比）。独自のkeymapを有効化する前にCapsをbackspaceに割当ておく

setxkbmap -option caps:backspace
連打スピードを上げる

キー押しっぱのデフォルトの認識は遅いので。

xset r rate 220 40

ネットワーク・クラウド周り

Dropboxを入れる

いつも使ってるので。

Chromiumを入れる

デフォルトのFirefoxは好かない

Telegramを入れる

仕事の連絡用

Nixnote2を入れる

Evernote用フリークライアント

SSH周りの整備

>>||
apt-get install openssh-server
|

グラフィック周り

Redshiftを有効にする

目がつかれるのでRedshiftをAutostartにする

Plankを入れる

Mac風ランチャー

sshfを入れる

計算機サーバをローカルにマウントしてファイルのやり取りを容易に

エディタ&プログラミング周り

PyCharm & CLion

使い慣れたPythonおよびC用エディタ

Neovimを入れる

GUIのVimエディタがほしかったので

terminatorを入れる

デフォルトのターミナルより高機能

Miniconda&PIPを入れる

Python環境作り

Pymolを入れる

生体分子表示用


こうしてみると入れる物たくさんありますね。。。。その他適宜追加します