Rosettaとaccelerated MDシミュレーションで小さなタンパク質のde novo構造予測と構造精緻化を行う方法

![brassica.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/64588c7a-a80b-9e5b-82a1-715c46d3bc8f.png)

## 背景
この記事は https://www.nature.com/articles/s41467-020-18698-w の手法に沿って書かれています。

アブラナ科を含む多くの植物は花のなかに雌ずい（めしべ）と雄ずい（おしべ）を持っており、雄ずいから出てきた花粉が雌ずいの柱頭にくっつくことを**受粉**と呼びます。受粉した花粉からは花粉管が花の内部を通って胚珠に向かって伸長していき、これが到達すると受精が行われて種子が形成され始めます。

ここで多くのアブラナ科植物には**自身の花粉で受粉を行う"自家受粉"が発生しないような仕組み**、いわゆる**自家不和合性**の機構が存在しており、正しく受粉が行われるためには他の個体の出す花粉が必要で、これによって近交弱勢（近親交配によって潜在していた有害な表現形質が現れ、集団中に適応度の低い個体が増えること）を回避するとともに集団の多様性を維持することにつながっていました。この話は[Wikipediaの自家不和合性の記事](https://ja.wikipedia.org/wiki/%E8%87%AA%E5%AE%B6%E4%B8%8D%E5%92%8C%E5%90%88%E6%80%A7_(%E6%A4%8D%E7%89%A9))が詳しいので気になる方は読んでみてください。ちなみにアブラナ科の植物の中でもシロイヌナズナは自家不和合性を失っており、自身の花粉で子孫を残す「自家和合性」の種に変化していることが明らかになっており、そちらの研究は三重大学大学院生物資源学研究科の諏訪部先生たちのグループによってなされているそうです（http://www.mie-u.ac.jp/R-navi/release/cat893/post-30.html ）

この自家不和合性自体は昔からよく知られている機構で、進化生物学者をはじめとして多くの科学的興味を引いてきましたが、これを制御しているタンパク質が柱頭で発現している**SRK**(S receptor kinase)と花粉外被（葯のタペート組織）で発現している**SP11**であることが同定されたのは2000年頃だったそうです。SRK, SP11ともに遺伝子上では*S*遺伝子座(*S*-locus)に含まれており、常にセットで1つのユニットとして子孫に受け継がれていきます。このユニットのことを*S*ハプロタイプと呼び、アブラナ科植物では100種類以上あると推定され、*S*<sub>1</sub>, *S*<sub>2</sub>, ..., *S*<sub>99</sub>のようにナンバリングされています（ナンバリングはやや適当っぽいところがあって、欠番があったり一部重複があったりしています）。**つまり、アブラナ科植物において同じ*S*ハプロタイプ由来のSRKとSP11が柱頭上で出会うと、花粉管の伸長を中止させて受粉をキャンセルさせるような仕組みが存在する**ことがわかりました。

またこのタンパク質の発見によって、進化生物学者たちだけでなく生化学者、構造生物学者たちによるこのタンパク質の結晶構造解析を行おうとする試みが始まりました。[PDB:1UGL](https://www.rcsb.org/structure/1UGL)はハプロタイプ*S*<sub>8</sub>のSP11のNMR構造決定による貴重な構造データです。しかし、SRK側のタンパク質は非常に大きくかつ不安定であるため、なかなか構造決定が進んでいませんでした。そんな中、2016年に[*Brassica rapa*](https://ja.wikipedia.org/wiki/%E3%83%96%E3%83%A9%E3%83%83%E3%82%B7%E3%82%AB%E3%83%BB%E3%83%A9%E3%83%91)のハプロタイプ*S*<sub>9</sub>のSRK/SP11複合体の論文が[Cell Research誌](https://doi.org/10.1038/cr.2016.129)で発表されました（[PDB:5GYY](https://www.rcsb.org/structure/5GYY)）。特にSRK（のエクトドメインと呼ばれる膜外領域）のタンパク質構造が解析されたことは非常に価値がありました。なぜならSRKの配列相同性は100あると言われるハプロタイプ中でも非常に保存度が高いため（9割以上同じでアミノ酸挿入や）、1つ解けるだけで残りは[Modellerを使った比較モデリング法](https://bi.biopapyrus.jp/protein/modeling/modeller/)、すなわちタンパク質の主鎖構造はほぼ同じとしてアミノ酸残基をつけかえるだけのモデリング法、などの手法でとても簡単に他のハプロタイプのSRKの尤もらしいタンパク質構造を推定することができるからです。

しかし一方で、SP11タンパク質の構造予測は単純に比較モデリング法ですべて作れるほど簡単ではありませんでした。

比較モデリング法が使えるための大まかな目安として、モデリングした構造未知のタンパク質のアミノ酸配列とタンパク質構造既知のアミノ酸配列とで配列アライメントをとったとき、だいたい3割くらい相同性があり、かつ大きな挿入や欠損部位が存在しないことが条件となっています。ここでSP11のハプロタイプたちの配列アライメントをまとめてとってみると、以下のようになります（ClustalW 2.1を利用）

![SP11alignment.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/758ccf0f-3d62-f0e9-f6b1-01b3454d0f17.png)

こんな感じでハプロタイプごとにアミノ酸の欠損や挿入の差が大きく、9割以上の位置でアミノ酸が保存されていないように見えます。ただ、一応アミノ酸配列に対する主成分分析（PCA）を使うと、構造既知の*S*<sub>8</sub>には*S*<sub>46</sub>, *S*<sub>47</sub>, *S*<sub>61</sub>が、同様に構造既知の*S*<sub>9</sub>には*S*<sub>32</sub>, *S*<sub>36</sub>が、それぞれ近くにクラスタリングされるため、これらについては比較モデリング法を使える余地がありました。しかし、その他のハプロタイプたちは配列の特徴が*S*<sub>8</sub>と*S*<sub>9</sub>からかけ離れており、特に下4つに示してあるclass 2と呼ばれる*S*<sub>29</sub>, *S*<sub>40</sub>, *S*<sub>44</sub>, *S*<sub>60</sub>の4つのハプロタイプは最も異なるアミノ酸配列を有していることがわかりました。

しかし、モデリングをする上での経験的な手がかりはいくつか残されています。1つは、ほとんどのSP11はハプロタイプのアミノ酸配列には8個のシステイン残基がほぼ同じような位置に現れること（例外は*S*<sub>46</sub>の7個、*S*<sub>32</sub>と*S*<sub>36</sub>は10個）であり、2つめは、2個めのシステイン残基が現れる2つ前の残基にはグリシン残基（G）が現れていることです。また、1番目の配列パターンからSP11の構造ドメインは[defensin](https://pfam.xfam.org/clan/Defensin)様の構造をとると予想されていることがヒントになりそうです。これはつまり、タンパク質の二次構造レベルで言えば、N末端からβストランド1→αヘリックス→βストランド2→βストランド3→C末端となる可能性がとても高いということです。ただ、単純にアミノ酸配列を眺めてもどこが二次構造をとるか、ループ領域になるかは簡単にはわかりません。

そこで、この構造モデリングがやりにくいClass 2ハプロタイプのSP11タンパク質構造を、テンプレート構造なしに計算で作り出せれば、対応するSRK構造と合わせてどのアミノ酸がSRK/SP11複合体形成に効いているかを厳密に議論することができ、将来的には自家不和合性を人為的に誘導もしくは阻害を起こすようなSP11の人工分子設計が可能になり、アブラナ科の野菜の優良品種の開発や生産に繋がることが期待されます。

## 目的
このチュートリアル記事ではアブラナ科植物の自家不和合性を制御するタンパク質SP11（*Brassica oleracea (Wild cabbage)* 由来）のハプロタイプの1つ、**SP11-29**（*S*<sub>29</sub>のSP11）のタンパク質立体構造をモデリングします。アミノ酸配列情報はUniprotデータベースから取得することにします。

## 計算環境
- Linux環境
- [Rosettaをインストールしている](https://qiita.com/Ag_smith/items/0e73d6d10e0b35c3f6d5)
- accelerated MD simulationsが可能な分子動力学ソフトウェアパッケージ([AMBER](https://ambermd.org/)または[NAMD](https://www.ks.uiuc.edu/Research/namd/)など)

## 手順
### 1. アミノ酸配列の取得
SP11-29のアミノ酸配列はUniprotデータベースに登録されていますので、そこから取得することができます。URLは https://www.uniprot.org/uniprot/Q84KT6 です。

ここのページにあるSequenceの見出しに行くとアミノ酸配列が表示されています。このページのFASTAボタンを押すと、アミノ酸配列のfasta形式で表示されます。

![スクリーンショット 2020-09-15 17.27.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/2fc4d056-fcb1-91ef-ccea-684a5b4bf588.png)

```shell
>tr|Q84KT6|Q84KT6_BRAOL S-locus protein 11 (Fragment) OS=Brassica oleracea OX=3712 GN=SP11-29 PE=2 SV=1
VQEVKANLRKRCPEHYSLPGVCGNSGNEECKRRYPYPINKIDLPTTCKCERSKFHKRGLC
KCSRNC
```

ちなみにUniprot内の検索で"brassica rapa sp11"と検索すると、SP11様のアミノ酸配列は86登録されているようです。

https://www.uniprot.org/uniprot/?query=brassica+rapa+SP11&sort=score

![スクリーンショット 2020-09-15 17.28.37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/f2227b0d-3f03-3cbd-61ab-ace5d6ff20b9.png)


### 2. Rosettaを使った*ab initio*モデリング

Rosettaのフラグメントアセンブリ法を使ったタンパク質モデリングを行います。このフラグメントアセンブリ法のコンセプトは、[Protein Data Bank](https://www.rcsb.org/)に登録されているタンパク質の構造を3-mer（アミノ酸の3ペプチド断片）や9-merに分割して"**フラグメント**"としてライブラリ化しておき、構造モデリングするようにユーザーに与えられたクエリのアミノ酸配列に応じて、近いアミノ酸配列を検索し、フラグメントライブラリからフラグメントを呼び出して、タンパク質になるように構築していく（アセンブル）というものです。アミノ酸配列が近いものを検知するアルゴリズムとしては、`psiblast`から得られるposition-specific scoring matrix (PSSM)の類縁配列アライメントから得られる配列プロファイルがもとになっています。

Rosettaプログラムをダウンロード＆インストールしてきた時にはデフォルトで`$ROSETTAHOME/tools/fragment_tools/vall.jul19.2011.gz`にフラグメントライブラリの1つが同梱されており(`$ROSETTAHOME`はRosettaをインストールしたディレクトリ名)、現在でも世界中のRosettaユーザーが発表しているタンパク質設計の多くの論文ではこのライブラリを利用しています。

しかし、このフラグメントライブラリは2011年でライブラリの更新が止まっていることもあり、すなわちフラグメントアセンブリの表現能力も止まっているように思えます。これについては[フラグメントライブラリを自作する](https://qiita.com/Ag_smith/items/58a7e28c75d2a51367f8)ことも可能で、実際に私が作ってみたものもあります。

私が2018年9月26日までのPDBデータで作ったRosettaの自作フラグメントライブラリデータは[こちらからダウンロード](http://vivace.bi.a.u-tokyo.ac.jp:8008/vall.bilab.sep26.2018.gz)できます。（739MBくらい）たぶんこれを使えばマシな構造が出る確率が高くなるんじゃないかな〜と思います。

それではこれを使ってRosettaのab initioモデリングを行ってみます。

#### 2-1. 二次構造の指定
まずは二次構造推定ファイル`SP11-29.ss2`を用意する必要があります。

```shell:SP11-29.ss2
# PSIPRED VFORMAT (PSIPRED V4.0)

   1 A C   0.998  0.001  0.001
   2 R C   0.878  0.079  0.043
   3 N C   0.891  0.061  0.048
   4 C C   0.980  0.009  0.018
   5 P E   0.937  0.071  0.019
   6 E E   0.001  0.001  0.998
   7 G E   0.001  0.001  0.998
   8 I E   0.001  0.001  0.998
   9 A E   0.001  0.001  0.998
  10 K E   0.001  0.001  0.998
  11 S E   0.001  0.001  0.998
  12 N C   0.866  0.113  0.054
  13 D C   0.918  0.041  0.058
  14 V C   0.806  0.037  0.164
  15 I C   0.998  0.001  0.001
  16 G C   0.998  0.001  0.001
  17 T C   0.998  0.001  0.001
  18 C C   0.998  0.001  0.001
  19 L C   0.998  0.001  0.001
  20 N H   0.020  0.599  0.375
  21 T H   0.017  0.694  0.321
  22 K H   0.021  0.920  0.070
  23 S H   0.013  0.969  0.023
  24 R H   0.389  0.625  0.004
  25 D H   0.197  0.778  0.009
  26 C H   0.052  0.949  0.006
  27 Q H   0.036  0.964  0.011
  28 K H   0.070  0.930  0.016
  29 H H   0.131  0.848  0.031
  30 F H   0.267  0.661  0.060
  31 G C   0.857  0.081  0.027
  32 P C   0.930  0.042  0.029
  33 N C   0.944  0.023  0.041
  34 V E   0.001  0.001  0.998
  35 T E   0.001  0.001  0.998
  36 N E   0.001  0.001  0.998
  37 C E   0.001  0.001  0.998
  38 L E   0.001  0.001  0.998
  39 C E   0.001  0.001  0.998
  40 Y E   0.001  0.001  0.998
  41 P E   0.001  0.001  0.998
  42 F E   0.001  0.001  0.998
  43 S E   0.001  0.001  0.998
  44 T C   0.884  0.198  0.025
  45 H C   0.943  0.025  0.033
  46 N C   0.930  0.013  0.060
  47 R C   0.759  0.010  0.233
  48 V E   0.001  0.001  0.998
  49 R E   0.001  0.001  0.998
  50 I E   0.001  0.001  0.998
  51 T E   0.001  0.001  0.998
  52 C E   0.001  0.001  0.998
  53 Y E   0.001  0.001  0.998
  54 C E   0.001  0.001  0.998
  55 C E   0.001  0.001  0.998
  56 K E   0.001  0.001  0.998
  57 V C   0.998  0.001  0.001
  58 K C   0.823  0.012  0.163
  59 S C   0.999  0.001  0.002
```


#### 2-2. Rosettaフラグメントピッキング

次にこの二次構造情報を使って、Rosettaの`fragment_picker`プログラムでSP11-29.faにマッチしそうなフラグメント構造を収集してきます（**フラグメントピッキング**）。このとき、どんな要素を重視してフラグメントを集めてくるかの設定を`simple.wghts`に書いておき、これをもとにして
`simple.wghts`は以下の通り

```shell:simple.wghts
# score name          priority  wght   max_allowed  extras
RamaScore               400     2.0     -       predA
SecondarySimilarity     350     1.5     -       predA
FragmentCrmsd             0     0.0     -
```

priorityとweightの列の値を適宜変更することでフラグメントライブラリからの3-mer, 9-merのフラグメントの選ばれ方が変化します。


`flags`というファイルを用意します。

```shell:flags
# Input databases
-in::file::vall /path/to/fragment_tools/vall.bilab.sep26.2018.gz

# Query-related input files
-in::file::fasta    SP11-29.fa
-in::file::s        SP11-29_ref.pdb
-in::detect_disulf  true
-frags::ss_pred     SP11-29.ss2 predA

# Weights file
-frags::scoring::config         simple.wghts

# we need nine-mers and three-mers
-frags::frag_sizes              9 3

# What should we do?
-frags::bounded_protocol

# Select 200 fragments from 1000 candidates. We need more candidates
#  than fragments for the selector to work properly
-frags::n_candidates            1000
-frags::n_frags                 200

-out::file::frag_prefix         frags1
-frags::describe_fragments      frags1.fsc
```

この時点で必要となるファイルは以下の通りです。

- `SP11-29.fa`: 予測する配列ファイル
- `SP11-29_ref.pdb`: 評価用タンパク質構造ファイル
- `SP11-29.ss2`: 二次構造情報の予測確率ファイル
- `simple.wghts`: フラグメント選択時の構造要素の重み設定ファイル
- `flags`: フラグメントピッキング実行の設定ファイル

これらを揃えて、ターミナルから以下のRosettaのコマンドを実行します

```bash
$ROSETTA3/bin/fragment_picker @flags
```

正常に終了すると、`frags.200.3mers`, `frags.200.9mers`, `frags1.fsc.200.3mers`, `frags1.fsc.200.9mers`という4つのファイルが生成されます。これらのファイルの中にはSP11-29の構造になりそうなタンパク質構造フラグメントの3-merと9-merが200ずつ入っています（ただ、次のフラグメントアセンブリで使うのは`fsc`がついていない方です）。フラグメントアセンブリではこれらをランダムに呼び出しながらタンパク質構造を作り上げていきます。

#### 2-3. Rosettaフラグメントアセンブリ
前項のフラグメントピッキングが終わったら、Ab initio Protein Foldingに移ります。
まず`options`というファイルを作成します。

```
# Make sure all variable names have been replaced with absolute path and that no line begins with a $ or ~s
-in
        -detect_disulf true
        -file
                -native SP11-29_ref.pdb        # native PDB file (optional)
                -fasta SP11-29.fa      # protein sequence in fasta format
                -frag3 frags.200.3mers  # protein 3-residue fragments file
                -frag9 frags.200.9mers  # protein 9-residue fragments file
-run
        -rebuild_disulf true
-constraints
        -cst_weight 1.5
        -cst_file coarse.cst
        -cst_fa_weight 1.5
        -cst_fa_file full.cst
-abinitio
        -detect_disulfide_before_relax true
        -relax
        -increase_cycles 10     # Increase the number of cycles at each stage in AbinitioRelax by this factor
        -rg_reweight 0.5        # Reweight contribution of radius of gyration to total score by this scale factor
        -rsd_wt_helix 0.5       # Reweight env, pair, and cb scores for helix residues by this factor
        -rsd_wt_loop 0.5        # Reweight env, pair, and cb scores for loop residues by this factor
-relax
        -fast   # At the end of the de novo protein_folding, do a relax step of type "FastRelax".  This has been shown to be the best deal for speed and robustness.
-out
        -nstruct 1000  # how many structures do you want to generate?  Usually want to fold at least 1,000.
        -file
                -scorefile score.sc
                -silent silent.out
-overwrite      # overwrite any existing output with the same name you may have generated

-nstruct 1000
```

以下の点はユーザーが特に状況に応じてファイルの場所・指定するべきところです。ファイル名の指定は相対パスでも絶対パスでも構いません。
- `-in:file:native SP11-29_ref.pdb` リファレンスとなる構造ファイル名を指定しておきます。アウトプット時に構造rmsdを計算する比較対象用です。
- `-in:file:fasta SP11-29.fa` インプットとして与えるアミノ酸配列ファイル名です。
- `-in:file:frag3 frags.200.3mers`, `-in:file:frag9 frags.200.9mers` は先程フラグメントピッキングで生成された3-mer, 9-merのライブラリファイル名です。
- `-out:nstruct`, `-nstruct`に指定した数だけモデル構造が生成され続けます（2ヶ所あるんだけど両方指定しないといけないような……）。
- `-run:rebuild_disulf true` **モデル生成の結果、一定距離以内に2つのシステイン残基が近づいた場合、ジスルフィド結合を作らせるようにします**。今回のSP11-29のようにジスルフィド結合がたくさんあるタイプのタンパク質の場合有用です。

また

```
-constraints
        -cst_weight 1.5
        -cst_file coarse.cst
        -cst_fa_weight 1.5
        -cst_fa_file full.cst
```

の部分は今回のSP11-29のモデリングに重要で、ジスルフィド結合を形成するであろうシステイン残基ペアの距離拘束と二次構造の形成に関わっています。`-cst_weight 1.5`, `-cst_fa_weight 1.5`でそれぞれCentroid（いわゆる主鎖原子とCβ原子のみのモデル）での拘束が1.5倍, full-atom model(`fa`の方)も1.5倍にしています（値を強くするとその影響が大きくなります）。
まず距離拘束を規定するファイル`coarse.cst`の中身は以下の通り。

```shell:coarse.cst
AtomPair CB 4 CB 55 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 4.0 1 0.5 #
AtomPair CB 37 CB 54 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 4.0 1 0.5 #
AtomPair CB 18 CB 39 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 4.0 1 0.5 #
AtomPair CB 26 CB 52 SCALARWEIGHTEDFUNC 9.014 BOUNDED 0 4.0 1 0.5 #
AtomPair O  11 N  48 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N  11 O  48 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair O   9 N  50 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N   9 O  50 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair O   7 N  52 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N   7 O  52 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N  24 O  20 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N  25 O  21 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N  26 O  22 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N  27 O  23 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N  28 O  24 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
AtomPair N  29 O  25 SCALARWEIGHTEDFUNC 9.000 BOUNDED 0 2.8 1 0.5 #
```
もう1つの距離拘束ファイル`full.cst`は以下の通り。

```shell:full.cst
AtomPair CB 4 CB 55 SCALARWEIGHTEDFUNC 9.099 SIGMOID 4.0 9.538
AtomPair CB 37 CB 54 SCALARWEIGHTEDFUNC 9.014 SIGMOID 4.0 9.538
AtomPair CB 18 CB 39 SCALARWEIGHTEDFUNC 9.014 SIGMOID 4.0 9.538
AtomPair CB 26 CB 52 SCALARWEIGHTEDFUNC 9.014 SIGMOID 4.0 9.538
AtomPair SG 4 SG 55 SCALARWEIGHTEDFUNC 9.099 SIGMOID 2.2 9.538
AtomPair SG 37 SG 54 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.2 9.538
AtomPair SG 18 SG 39 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.2 9.538
AtomPair SG 26 SG 52 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.2 9.538
AtomPair O  11 N  48 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N  11 O  48 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair O   9 N  50 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N   9 O  50 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair O   7 N  52 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N   7 O  52 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N  24 O  20 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N  25 O  21 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N  26 O  22 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N  27 O  23 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N  28 O  24 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
AtomPair N  29 O  25 SCALARWEIGHTEDFUNC 9.014 SIGMOID 2.8 9.538
```


`coarse.cst`の距離拘束はモデリングの初期段階に、`full.cst`の距離拘束は構造のリファインメント(relax)時に効いてきます（たしか）。Rosettaのフラグメントアセンブリ法は、構造生成時にまず第1段階としてCα、Cβ原子の座標のみで計算して主鎖構造の形成から行い、そこから側鎖全体を"生やして"タンパク質を完成させるという流れになっています。したがって`coarse.cst`の場合は主鎖原子のCα, N, C, OとCβ原子までの距離拘束しか導入することはできませんが、`full.cst`の方ではすべての原子名を使って拘束を入れることができます（よってこちらの方でのみシステインのSG原子間の距離拘束が入っています）。今回は、構造のリファインメントは後述のaccelerated MDでやるので、距離拘束の値は大まかでいいかなと考えていました。

ここまでで必要なファイルは以下の通りです。

- `SP11-29.fa`: 予測する配列ファイル
- `SP11-29_ref.pdb`: 評価用タンパク質構造ファイル
- `coarse.cst`: 距離拘束ファイル（モデリング初期段階時）
- `full.cst`: 距離拘束ファイル（リファインメント時）
- `frags1.200.3mers`: フラグメントファイル（3-mer）
- `frags1.200.9mers`: フラグメントファイル（9-mer）
- `options`: フラグメントアセンブリ実行時の設定ファイル

以上を揃えたら、以下のコマンドでRosettaのフラグメントアセンブリを実行します。

```
mpirun -np ${PBS_NP} $ROSETTA3/bin/AbinitioRelax.mpi.linuxgccrelease @options
```

`mpirun`を利用することでフラグメントアセンブリによるタンパク質モデルを（コア数-1）分だけ倍速に生成できるようになります（16コアマシンなら15倍）。

#### 2-4. モデル構造の取り出し
計算が始まってしばらく経つと、`silent_*.out`という名前のファイルがたくさん生成され始めます。これらのファイルに、Rosettaによって作られた予測タンパク質構造が圧縮されて保存されています。それまでに作り出されたタンパク質構造のエネルギースコアをソートして見て見るために、`sort -k2 score.fsc | less -SN`のコマンドを入れてみましょう。これは、`score.fsc`の2番目のカラム（エネルギースコア）をソートして、さらに`less`コマンドのオプションで右端折返しをしないまま表示してくれます。

```
SCORE:     score     fa_atr     fa_rep     fa_sol .. .. .. description
...
...
SCORE:  -132.108   -293.844     50.634    194.927 .. .. .. S_00000774
SCORE:  -138.378   -294.880     46.939    192.731 .. .. .. S_00000280
```

`score.fsc`の結果から、`S_00000280`が最も低いエネルギースコアを持つことが示されていました。この構造が安定な構造である可能性が高いということです。実際にRosettaだけで構造予測をしたい場合には1万から10万程度のモデルの生成計算が必要になることもあるらしいですが、今回は1000くらいで終えていました。

この`S_00000280`の構造をPDB形式で取り出すために、以下のRosettaのコマンドを入力します。

```bash
extract_pdbs -in:file:silent silent_*.out -tags S_00000280
```

これにより、全`silent_*.out`の中から`S_00000280`の名前を持つ構造をサーチしてPDB形式に変換することができます。

このときの構造をPyMOLで見てみるとこんな感じでした（※ランダムなので、必ずこうなるわけではないです）。

![SP11-29init.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/982823a6-a5af-52fe-aa7a-0959d6adc125.png)


![SP11-29init2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/31358342-ad47-e8b5-e80c-4aee7ad565b1.png)

βシートはいい感じに3本きれいに揃っていますが、なんかαヘリックスがねじれまがっていますね（っていうところが気になるのは熟練者の経験によるところが大きいかもしれませんが）。そこで、続くAccelerated MDシミュレーションによってこのRosettaによって出てきた構造に対して構造精緻化を行い、もっともらしいものになるよう検討していきます。

### Accelerated MDシミュレーションによる精緻化

以前私が書いた[Structure-based drug designとAccelerated MDによるタンパク質構造探索の方法の紹介](https://qiita.com/Ag_smith/items/ccd599eed99af8780c01)という記事に沿って行っていきます。

（現在書きかけです）