機械学習の回帰分析の評価関数

概要

回帰分析の評価関数がいくつかあるので、それらをまとめる。 はそれぞれ、正解値、予測値を表す。今回は、全データが個あるとしている。

R² (決定係数)

機械学習でよく用いられる。次元がなく、比べられる値。

RMSE (Root Mean Squared Error)

直感的にわかりやすいが、他のRMSEと比べられるかどうかはしっかり注意しなければならない。

MAE (Mean Absolute Error)

RMSEと同様に直感的でわかりやすいが、比べられるかどうかは注意しなければならない。

ピアソンの相関係数 (R)

ピアソンの相関係数は次元がなく、比べられる値となっている。

スピアマンの順位相関係数 ()

式の形はピアソンの相関係数と全く同じ形をしている。しかし、ここでの変数は値そのものではなく、順位であることに注意してほしい。

参考文献

精度評価指標と回帰モデルの評価 ピアソンの相関係数、Spearman(スピアマン)の順位相関係数　役に立つ薬の情報　専門薬学

Ornstein-Zernike(OZ)方程式の導出

概要

研究で3D-RISM理論を扱うことになり、3D-RISMはOZ式から導出されるので、OZ式の導出をやってみた。カノニカル分布、グランドカノニカル分布、デルタ関数の概要を理解している前提で進める。

1. 分布関数

粒子の座標を、ポテンシャル関数をと定義すると、N個の粒子が座標となる確率密度は

ある物理量の熱平均を以下のように定義する。

デルタ関数を用いた確率密度

微小体積要素に粒子がある確率はと表される。これを拡張し、粒子まで考慮する。つまり、粒子が、粒子が、、、にある確率を考えると以下のようになる。

これを、熱平均することで確率密度を定義できる。

は熱力学的平均をとっているため、エネルギーを考慮した確率である。上記のは粒子１、粒子２、、、を区別している。これらを区別しないでに粒子が一つ、に粒子が一つ、、、という確率密度を考えると

一様な流体では

粒子を区別しない時の関数による定義、つまり、微小体積要素付近に任意の粒子を少なくとも１つ見出す確率は

これをn個の粒子に拡張すると

2. 近距離にある粒子間の相関

との熱力学的平均を計算してみる。

が十分離れている時、となる。近いときには、二つの位置に相関があるので、その相関をとすると

等方的な流体では、距離を用いて以下のようにかける。

気になったので計算してみた（余談）

物理量を考えると、熱力学平均は以下のように計算できる。

さらに、の熱力学平均を考えてみる。上記と同様にして

3. 密度の揺らぎ

密度の揺らぎについて考える。平均値からのずれをとすると

二次のモーメント、つまり密度揺らぎの相関関数を考えると

ここで、とした。ようやく全相関関数がでてきた。

4. グランドカノニカルの分布関数

とすると、

直接相関関数の定義

式(4.4)の逆関数を使って定義する。

5. 鎖則

上記鎖則(5.1)に式(4.4)、(4.5)をそれぞれを代入するとOZ式が得られる。

、デュラックのデルタ関数の性質より

参考文献

溶媒効果とMO計算

ニューラルネットワークの用語

概要

ニューラルネットワークに関連する用語の意味をまとめる。

勾配降下法

バッチ勾配降下法

全ての訓練データを使用して重みを更新する。

• 時間がかかる
• 局所解に陥るとなかなか脱出できない

確率的勾配降下法（SGD、Stochastic Gradient Descent）

1つの訓練データを使用して重みを更新する。

• 時間がかからない
• 局所解に陥り辛い
• 1つのデータしか使用していないため、勾配の方向が不正確

ミニバッチ勾配降下法

複数の訓練データ（ミニバッチ）を使用して重みを更新する。 バッチ勾配降下法と確率的勾配降下法のいいとこどり。

活性化関数

ステップ関数

単純パーセプトロン

シグモイド関数

tanh

ReLU (Rectified Linear Unit)

Bashスクリプト

概要

Bashスクリプト作成の時にいつも調べてしまうのでまとめておく。

環境

macOS Mojave 10.14.5

変数の扱い

a='aaa'
echo ${a}

計算の仕方

two=$((3-1))
echo $(two)

if

bashの比較演算子は文字列と数値とで大きく異なる

文字列比較

=で比較する

a='aaa'
if [ ${a} = 'aaa' ]; then
    echo 'aaa'
elif [ ${a} = 'bbb' ]; then
    echo 'bbb'
else
    echo 'error'
    exit 0
fi

数値比較

-lt、-gt、-eqで比較する

a=3
if [ ${a} -lt 3 ]; then
    echo '${a} is less than 3'
elif [${a} -eq 3 ]; then
    echo '${a} is 3'
elif [${a} -gt 3 ]; then
    echo '${a} is greater than 3'
else
    echo 'else'
fi

ファイルの存在確認

-eはexistの略？

file='test.txt'
if [ ! -e ${file} ]; then
    echo '${file} not found'
    exit 0
fi

for

bashにリストという概念がない。基本は以下のように指定する。

for a in 1 2 3 4 5; do
    for b in {01..06}; do
        echo ${a} ${b}
    done
done

時間を取得

date %Y/%m/%d/%H:%M:%S

コマンドの出力を変数として取得

スクリプトを書くときはよく使う

time=`date +"%Y/%m/%d/%H:%M:%S"`

関数の定義

引数は${1}、${2}のようにコマンドライン引数と同様に扱う。

echo_dir() {
    dir=${1}
    echo ${dir}
}

echo_dir 'test.txt'

Amberで力場パラメータを作成する

概要

リガンド-タンパク質の複合体の力場パラメータ（拡張子が.prmtop 、.inpcrd）を作成する。

環境

AmberTools19

タンパク質の準備

タンパク質についてはパラメータを作成する前にいくつか注意点がある

• 解像度は十分高いか
• 欠損残基はあるか：キャップを形成するかモデリングする必要を検討
• 欠損原子はあるか：leapで補強してくれるが、補強が原因で原子が重なってしまう可能性がある。
• 金属があるかどうか：新たに力場を準備する必要があるか検討する。
• 変異残基はあるか：leapでエラーが出る。変異残基の力場を準備する必要がある
• ジスルフィド結合はあるか：あればシステインの残基名を変更する必要がある
  • CYS：ジスルフィド結合していないシステイン
  • CYX：ジスルフィド結合しているシステイン

• ヒスチジンのプロトン化：プロトン化によって残基名を変更する必要がある。何も変更をしないと、全てHIEになる。

  

  Amberではtleapを実行する前の準備としてpdb4amberを実行することが推奨されている。pdb4amberを実行すると、ジスルフィド結合を自動的に考慮してくれる。

リガンドの準備

リガンド総電荷を計算する。これをlig_chargeとする。その後、以下を実行する。antechamberでgaff（Generate Amber Force Field）を参考にして各原子の電荷を計算し、parmchk2でgaffで補いきれなかったパラメータを作成する。

#!/bin/bash

lig_mol2=lig.mol2
out_name=lig_out
lig_charge=1
LigResName='Lig'

antechamber -fi mol2 -i ${lig_mol2} -fo mol2 -o ${out_name}.mol2 -at gaff -c bcc -nc ${lig_charge} -rn {LigResName} > antechamber.log

parmchk2 -i -f ${out_name}.mol2 -o ${out_name}.frcmod > parmchk2.log

複合体のパラメータ作成

以下のtest.inを作成する。Amberで用意されている各力場パラメータと、先ほど作成したリガンドのmol2、frcmodファイルを読み込んで複合体のパラメータを作成する。

source leaprc.protein.ff14SB
source leaprc.gaff

loadamberparams ${out_name}.frcmod

prt = loadpdb prt.pdb
lig = loadmol2 lig_out.mol2
sys = combine{prt lig}
charge prt
charge lig

saveamberparam sys sys.prmtop sys.inpcrd
savepdb sys sys.pdb
quit

test.inをtleapコマンドで実行する。

$ tleap -f test.in 

ジスルフィド結合を定義する

タンパク質にジスルフィド結合があれば、定義する必要がある。原子間が結合していることを示す。ここでは二つの方法を紹介する。

• pdbファイルの最後に結合情報を記す。以下のようにインデックス番号で指定する。pdb4amberを実行するとシステイン残基が2Å以内であればジスルフィド結合と判定され自動的にCONECT情報が追加される。

ATOM    950  CG2 VAL A 124      28.143  10.116  17.555  1.00 22.08           C  
ATOM    951  OXT VAL A 124      27.661  12.956  19.992  1.00 22.31           O  
TER   
CONECT  194  642
CONECT  310  727
END

• tleapを実行時にleapファイル内で指定する。以下のように残基番号で指定する。

bond prt.32.SG prt.35.SG

金属を考慮する

• イオンとして計算する
  amber内で用意されている力場を指定する。この他にもいくつか金属イオンの力場は用意されているので適当なものを選択する。

loadoff atomic_ions.lib
loadamberparams frcmod.ions2341m_126_spce
loadamberparams frcmod.ions1lm_126_spce

• 新しい力場を導入する
  Offファイルとfrcmodファイルの二つを用意し、以下のように読み込む

loadoff test.off
loadamberparams test.frcmod

Ridge回帰の実装

概要

実装はsklearnを利用する。

環境

sklearn 0.21.2
matplotlib 3.1.1

実装

データ

ボストンデータを利用する。

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

data = load_boston()
x_train, x_test, y_train, y_test = train_test_split(data['data'], data['target'], test_size=0.2, random_state=0)
scaler = StandardScaler()
scaler.fit(data['data'])
x_train = scaler.transform(x_train)
x_test = scaler.transform(x_test)

モデル

from sklearn.linear_model import Ridge

model = Ridge()
model

デフォルトのパラメータ

Ridge(alpha=1.0, copy_X=True, fit_intercept=True, max_iter=None,
      normalize=False, random_state=None, solver='auto', tol=0.001)

ハイパーパラメータチューニング

今回のハイパーパラメータはalphaのみ。グリットサーチクロスバリデーションで最適化する。5-fold cross validationをする。

#grid_search cross validation
n_splits = 5
params = {'alpha' : [1.0]}

グリッドサーチクロスバリデーション

グリッドサーチクロスバリデーションを定義する。訓練データでも評価したいのでGridSearchCVを使わずに自分で定義する。

from sklearn.model_selection import train_test_split, KFold, ParameterGrid
from statistics import mean, stdev
import copy

def gscv(datas, grid_params, n_splits=5):
    x_train, y_train = datas
    kf = KFold(n_splits=n_splits, shuffle=False, random_state=0)
    results = {
        'params' : [],
        'validation_score_li' : [],
        'validation_score_mean' : [],
        'validation_score_std' : [],
        'train_score' : []
    }
    for params in ParameterGrid(grid_params):
        results['params'].append(params)
        for k, v in params.items():
            results.setdefault(f'param_{k}', []).append(v)
        score_li = []
        model.set_params(**params)
        for i_train, i_test in kf.split(x_train):
            x_cv_train = x_train[i_train]
            y_cv_train = y_train[i_train]
            x_cv_test = x_train[i_test]
            y_cv_test = y_train[i_test]
            model.fit(x_cv_train, y_cv_train)
            score = model.score(x_cv_test, y_cv_test)
            score_li.append(score)
        results['validation_score_li'].append(score_li)
        results['validation_score_mean'].append(mean(score_li))
        results['validation_score_std'].append(stdev(score_li))
        model.fit(x_train, y_train)
        score = model.score(x_train, y_train)
        results['train_score'].append(score)
    return results

次に、validationスコアの平均と訓練データでのスコアをプロットする関数を定義する。

import matplotlib.pyplot as plt

def ShowLiner(results, x_title):
    x_title = f'param_{x_title}'
    y_title = 'score'
    x_ticks = results[x_title]
    x = range(0, len(x_ticks))
    y_cv = results['validation_score_mean']
    y_train = results['train_score']
    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_subplot(111)
    ax.plot(x, y_cv, label='cv_score')
    ax.plot(x, y_train, label='train_score')
    ax.set_xlabel(x_title)
    ax.set_ylabel(y_title)
    ax.set_xticks(x)
    ax.set_xticklabels(x_ticks)
    plt.grid()
    plt.legend()
    plt.show()

最後に、ハイパーパラメータをチューニングするパイプラインを関数で定義する。ハイパーパラメータをチューニングし、グラフを表示し、もっともよかったパラメータを自動的に更新する。

def search_param1(param_name, param_values, params):
    new_params = params.copy()
    new_params[param_name] = param_values
    print(f'set new_params: {new_params}')
    results = gscv([x_train, y_train], new_params)
    cv_score_li = results['validation_score_mean']
    best_index = cv_score_li.index(max(cv_score_li))
    print(f"cv score: {results['validation_score_mean'][best_index]:.2f}+-{results['validation_score_std'][best_index]:.3f}")
    print(f"best params: {results['params'][best_index]}")
    ShowLiner(results, param_name)
    params[param_name] = [results['params'][best_index][param_name]]

実行する

search_param1('alpha', [0, 1, 10, 100, 1000], params)

set new_params: {'alpha': [0, 1, 10, 100, 1000]}
cv score: 0.50+-0.211
best params: {'alpha': 10}

青が検証データ、黄色が訓練データである。alphaが10のときに一番スコアがよかった。

訓練と予測

#train
model.set_params(**params)
model.fit(x_train, y_train)

#predict
y_pred = model.predict(x_test)

結果と分析

決定係数は組み込み関数で計算できる。

r2 = model.score(x_test, y_test)
print(r2)

0.5796111714164923

予測した値をプロットしてみる。

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.scatter(y_test, y_pred)
ax.set_xlabel('y_test')
ax.set_ylabel('y_pred')
plt.show()

切片と係数を取得し、係数の大きさを確認してみる。

w_vec = model.coef_
b = model.intercept_

import matplotlib.pyplot as plt
import numpy as np

w_vec = model.coef_
names = data['feature_names']
mask = np.where(w_vec > 0, True, False)

w_vec = np.abs(w_vec)
order = np.argsort(w_vec)
mask = mask[order]

fig = plt.figure(figsize=(6, 10))
ax = fig.add_subplot(111)
line_n = ax.barh(names[order], w_vec[order], color='orange')
line_p = ax.barh(names[order][mask], w_vec[order][mask], color='blue')
line_p.set_label('Positive number')
line_n.set_label('Negative number')
ax.legend()
plt.show()

参考文献

https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html

VMDで水素結合を解析する (GUI, CUI)

概要

vmdで水素結合の解析をする。GUIとCUIでの方法を紹介する。

環境

macOS Mojave 10.14.5
vmd 1.9.3

GUIでの方法

まずは構造ファイルとトラジェクトリ ファイルを読み込む。読み込み方はこちら → VMDの使い方の基本 (Mac)

水素結合の数をプロット

1. VMD Main ==> Extensions ==> Analysis ==> Hydrogen Bonds

   

2. 以下を参考に詳細を設定する。

   

   • Selection 1 : protein
   • Selection 2 : resname G39
   • Update selections every frame? : チェックを外すと計算時間が短くなる。
   • Only polar atoms (N, O, S, F)? : チェックする
   • Calculate detailed info for : Residue pairsをチェックする
   • Output directory : お好み
   • Frame/bond data? : お好み
   • Detailed hbond data? : お好み
   • Find hydrogen bonds! : 設定が終わればここをクリックし、計算を始める

3. プロットが現れる。FileからExport to PostScriptを選択すると画像を保存できる。また、Export to ASCII matrix...を選択すると計算結果データを保存できる。

   

アウトプットファイルの説明

アウトプットファイルは2つある。Frame/bond data?で指定された結果ファイル(test.dat)とDetailed hbond data?で指定された詳細ファイル(test-detailed.dat)が出力される。

• test.dat : 計算結果が書かれている。各フレームの水素結合の本数が記述されている。

#frame number
0 2
1 0
2 0

• test-detailed.dat : 計算結果の詳細が書かれている。見つかった水素結合の組み合わせと、出現頻度が書かれている。Calculate detailed info forオプションの設定によって詳細ファイルの内容が変わる。以下、出現頻度について補足を加える。

Residue pairs open 下の例では、ARG290-SideとG39389-Sideでの水素結合出現頻度が0.30%となっている。今回は全部で1000フレームあったので1000*0.30/100=3フレームで水素結合があったことになる。

Found 16 hbonds.
donor        acceptor    occupancy
ARG290-Side      G39389-Side     0.30%
ARG212-Side      G39389-Side     0.10%

close

Unique hbondopen 下の例では、残基名は同じ（ARG290-SideとG39389）だが、NH1、NH2、O1B、O1Aという原子名の違いで出現頻度が計算されている。

ARG290-Side-NH1      G39389-Side-O1B     78.22%
ARG290-Side-NH2      G39389-Side-O1A     76.22%
ARG290-Side-NH1      G39389-Side-O1A     0.10%

close

All hbonds open 下の例では、ARG290-SideとG39389-Side間での水素結合出現頻度が154.55%となっており、100%を超えている。この理由は上記Unique hbondの例からわかる通り、78.22+76.22+0.10=154.54%となる。よってAll hbondsでは出現頻度は原子間で計算しているが、組み合わせが残基間で表示されていることがわかる。

ARG37-Side   G39389-Side     0.10%
ARG290-Side      G39389-Side     154.55%

close

Hydrogen Bondsメニューオプション

• Input options
  • Molecule：解析したい分子の指定
  • Selection 1：解析したい原子の指定
  • Selection 2：解析したい原子の指定
  • Frame：解析したいフレームの指定
  • Update selections every frame?：フレームごとにSelection 1,2で指定した原子の選択を毎回更新するかどうかの指定。原子を時間依存しないように指定（ex 原子名など）していればチェックを外すことで計算時間を節約できる。原子を時間依存するように指定（ex 座標での指定、またはprotein within 5 of (resname G39)のように、時間経過によって選択される原子が変化する）していれば毎回更新しなければならない。
  • Only polar atoms：解析対象原子をN, O, S, Fに限定する
  • Selection 1 is the：Selection 1で選んだ原子をドナー、アクセプター、またはその両方として解析するかどうか。Selection 1でドナーとすると、Selection 2で選んだ原子はアクセプターとなる。
  • Donor-Acceptor distance：水素結合を定義するときの距離
  • Angle cutoff：水素結合を定義するときの角度
  • Calculate detailed info for：解析結果の詳細情報を定義する。ここが一番理解しにくい。具体例はtest-detailed.dat。
    • None：何もしない
    • All hbonds：原子間で水素結合があると１とカウントする。しかし、集計自体は残基間で行われる。よって、出現頻度が100%を超えるものもある。
    • Residue pairs：残基間で水素結合があると１とカウントする。残基間で水素結合が2本以上あったとしても、１とカウントする。
    • Unique hbond：原子間で水素結合があると１とカウントする。
• Output options
  • Output directory：出力ファイルのディレクトリ
  • Log file：ログファイル名
  • Frame/bond data：出力ファイル名
  • Detailed hbond data：解析結果の詳細ファイル名

CUIでの方法

vmdのスクリプト作成が初めてという方はこちら→VMDのスクリプトの書き方。

measure hbonds

水素結合情報はmeasure hbondsで取得できる。このコマンドで取得できるのは水素結合を形成しているドナー、アクセプター、水素のインデックスを要素とするリスト。リストの長さは３つ（ドナー、アクセプター、水素）なので、ここから水素結合の本数、角度、距離など様々な情報を取得できる。ここでは、例としてGUIで出力した結果ファイル（test.dat）と同様の出力をする。

mol load pdb test.pdb dcd test.dcd

set out_file [open test.dat w]
set num_frames [molinfo top get numframes]

for {set i 1} {$i < $num_frames} {incr i} {
    set donor [atomselect top "protein" frame $i]
    set acceptor [atomselect top "resname G39" frame $i]
    set hbond_li [measure hbonds 3 20 $donor $acceptor]
    set num_hbonds [llength [lindex $hbond_li 0]]
    puts $out_file "$i\t$num_hbonds"
}
exit

水素結合の本数はmeasure hbondsで取得したリストのいづれかの長さと等しくなる。結合の距離を計算したいときはmeasure bond {3 5}のようにすればインデックスが３と５の原子間距離を計算できる。詳しくは参考文献を参照。

参照文献

VMD Hbonds Plugin
VMD User's Guide: measure