C++/Eigen での RMSD計算

タンパク質の構造変化を定量化する概念として RMSD (Root mean square deviation) というのがある。これは要するに距離の差の二乗平均なので非常に汎用性がありそうなのだが、タンパク質の構造以外の用途に使う人を見たことがない。

Biopython に SVD (特異値分解) を使って RMSD を計算するコードがあったので、これを参考に C++ Eigen における RMSD 計算のコードを書いた。

タンパク質の原子座標には構造変化のほかに並進や回転の自由度もあるので、これをあらかじめ除去し、構造変化による部分だけのRMSDを計算する。手順は以下のとおり。

• 座標を N×3 行列に格納する (N は原子数)
• 両構造の重心を原点に合わせる
• 特異値分解で回転行列を求める
• 回転を実行しRMSDを計算する

という手順。

#include <iostream>
#include <Eigen/SVD>
#include <Eigen/LU>

using Mat = Eigen::Matrix<double, Eigen::Dynamic, 3>;
using SqMat = Eigen::Matrix3d;

int main() {
    Mat ref, tar; // reference, target

    // 座標入力は省略
    
    // 重心を原点へ
    Eigen::Vector3d centroid;
    centroid = ref.colwise().sum() / ref.rows();
    ref.rowwise() -= centroid.transpose();
    centroid = tar.colwise().sum() / tar.rows();
    tar.rowwise() -= centroid.transpose();

    // SVD をとる
    SqMat a, U, V, rot;
    a = tar.transpose() * ref;

    Eigen::JacobiSVD<SqMat> svd(a, Eigen::ComputeFullU | Eigen::ComputeFullV);
    U = svd.matrixU();
    V = svd.matrixV();
    rot = (V*U.transpose()).transpose();

    // 鏡像チェック
    if (rot.determinant() < 0) {
        V.row(2) *= -1;
        rot = (V*U.transpose()).transpose();
    }

    // ref のうえに superimpose する
    Mat tar_on_ref = tar * rot;

    // RMSD を計算
    Mat diff = tar_on_ref - ref;
    double rmsd = std::sqrt((diff.array() * diff.array() / diff.rows()).sum());

    std::cout << "rmsd: " << rmsd << '\n';
}