monomer libraryとtorsion angle CCP4 monomer libraryにはねじれ角(torsion angle)の理想値・周期・標準偏差が記述されています． 例えばARGの定義は https://github.com/MonomerLibrary/monomers/blob/master/a/ARG.cif にありますが，現時点でのtorsion angle restraintsを眺めてみると，以下のようになっています． ARGの構造と原子名の対応についてはRCSB PDB - ARG Ligand Summary PageでLabelsボタンを押して確認してください． loop_ _chem_comp_tor.comp_id _chem_comp_tor.id _chem_comp_tor.atom_id_1 _chem_comp_tor.atom_id_2 _chem_comp_tor.atom_id_3 _chem_comp_tor.atom_id_4 _chem_comp_tor.value_angle _chem_comp_tor.value_angle_esd _chem_comp_tor.period ARG chi1 N CA CB CG -60.000 10.0 3 ARG chi2 CA CB CG CD 180.000 10.0 3 ARG chi3 CB CG CD NE -60.000 10.0 3 ARG chi4 CG CD NE CZ 180.000 10.0 6 ARG chi5 CD NE CZ NH2 180.000 5.0 2 ARG hh1 NE CZ NH1 HH12 180.000 5.0 2 ARG hh2 NE CZ NH2 HH22 0.000 5.0 2 ARG sp3_sp3_1 C CA N H 180.000 10.0 3 ARG sp2_sp3_1 O C CA N 0.000 10.0 6 右側3つの数字がvalue_angle, value_angle_esd, periodで，理想値・その標準偏差・周期を表しています． 例えばchi1は理想値-60±10°ですが，周期が3なので360/3=120°おきに理想値が存在する，つまり-60, 60, 180°は全て理想値ということになります． ...

昔，ANODEで計算する異常分散差フーリエが，他のプログラム(phenixやRefmac)と比べてもS/Nが良いことが気になり，調査したときのメモを発掘したのでここに清書． ANODEが良いのはその前処理に使っているSHELXCによるoutlier rejectionが大いに効いていることが分かった． SHELXCは以下の2つの条件に基づいてrejectionを行う模様． 条件1: $$ \sigma^2(I^+)>6.25\sigma^2(I^-) {\rm or } \sigma^2(I^-)>6.25\sigma^2(I^+) $$ 条件2: $$ \min(\sqrt{I^+}, \sqrt{I^-}) \le 4\sqrt{ (\sqrt{I^+ +\sigma(I^+)} - \sqrt{I^+})^2 + (\sqrt{I^- +\sigma(I^-)} - \sqrt{I^-})^2 } $$ ただし\(I < 0\)のときは\(I = 0\)とする． いかにもSHELXらしい，empiricalな方法という感じ． このrejection後，SHELXCは\(e^{-2.5s^2/4}\)を掛けて出力している．さらにANODEは\(B=10\)を掛けて計算しているようだ． unmerged intensityを与えた場合の挙動は調査しきれていない．昔は重みなしで平均を計算してしまっていたが，ある時点から重み付けするようにはなったらしい(10年前に聞いた話)． あとは\(\sigma(\langle I \rangle)\)はinternal/external varianceの最大値として計算されているようだ． 上記rejectionを確認するために書いたと思われるコードも発掘したので，動くかも分からないがいちおう載せておく．Python2だ… import iotbx.file_reader from cctbx.array_family import flex import math if __name__ == "__main__": f1 = "xscale.sca" fano = "anode_fa.hkl" i_org = iotbx.file_reader.any_file(f1).file_server.miller_arrays[0] i_org.show_summary() print fa = iotbx.file_reader.any_file(fano+"=amplitudes").file_server.miller_arrays[0] fa.show_summary() fa = fa.customized_copy(crystal_symmetry=i_org.crystal_symmetry()) i_org_fa = i_org.anomalous_differences() matches = i_org_fa.match_indices(fa) flags = flex.bool(i_org.size(), False) flags.set_selected(matches.pairs().column(0), True) ofs = open("sca_common.dat", "w") print >>ofs, "h k l d I sig common" for hkl, r, d, s, f in zip(i_org_fa.indices(), i_org_fa.d_spacings().data(), i_org_fa.data(), i_org_fa.sigmas(), flags): print >>ofs, "% 4d % 4d % 4d" % hkl, print >>ofs, r, d,s,f asu, matches = i_org.match_bijvoet_mates() ip = i_org.select(matches.pairs().column(0)) im = i_org.select(matches.pairs().column(1)) # flag is True when not rejected ofs = open("test_reject.dat", "w") print >>ofs, "h k l Ip Im SA flag" for hkl, Ip, SigIp, Im, SigIm in zip(ip.indices(), ip.data(), ip.sigmas(), im.data(), im.sigmas()): print >>ofs, "% 4d % 4d % 4d" % hkl, print >>ofs, Ip, Im, S, T = SigIp**2, SigIm**2 flag = not (S>6.25*T or T>6.25*S) if flag: Fp, Fm = math.sqrt(max(Ip, 0)), math.sqrt(max(Im, 0)) sa = math.sqrt( (math.sqrt(Fp**2+SigIp)-Fp)**2 + (math.sqrt(Fm**2+SigIm)-Fm)**2 ) flag = min(Fp,Fm) > 4*sa print >>ofs, sa, flag 当時，かなりがっつり反射数が減ることに驚いたような記憶がある． 適切なweightingを考えるよりも経験的な方法でrejectしてしまう方がいいというのはちょっと残念な感じがするが，もっと良い方法にたどり着けたらいいですね． ...

最近たくさん登録作業をしたので，以前学生向けに書いたメモをここに再掲しようと思う． 2021年のPDBj講習会資料も参照． 登録には，how to deposit EM MAPにあるように，PDB or mmCIFファイルのほか，以下の情報が必要です． Primary map along with voxel size and recommended contour level 特にルールは無さそう．モデリングに使った主なマップで良いはず．postprocess_masked.mrcなど Image of the map (500 x 500 pixels in .jpg, .png, etc. format) ウェブでサムネイルとして表示される絵．何でも良い．例 Half maps (as used for FSC calculation; two maps must be uploaded) Mask FSC curves RELIONの場合はPostProcessで出力されるpostprocess_fsc.xmlをアップできる 決めておくべきこと・調べておくべきこと contact author全員のORCiD PDB/EMDBと連絡を取る人+PIのみ．公開されない Entry title & authors Surname, F.M.形式で書く Funding (organisationとnumber) unreleasedに載せるか隠すか 分子ごとのmolecule name 配列（タグの切れ残り等含め，凍結した試料に含まれた配列全部） モデルと配列が一致していることを確認しておく．食い違っていると進めない sourceとexpression host 各マップのrecommended contour level (モデルのvalidationにも使われる) グリッド作製時のpH, 冷媒，defocus範囲 (nm), 電顕と検出器の型，Average electron dose 解析手法, 粒子数, 分解能 よくあること “Still processing. Processing done.“と出てるのにページ遷移しないときがある．リロードしてしまって大丈夫そう 他のエントリから情報をインポートできるのは最初の画面だけ パスワード自動入力が有効だと，Do you want to import information from a previous wwPDB deposition?をNoにしてるのに（ID/PWも空欄にしてるのに），アップロードから先に進むときThe related deposition ID is the same as the current deposition…と出て入力し直しになる現象が起きる．進む前にもう一度確認するようにした方が良い 配列アラインメントがあまり賢くなく，gapがあるときに失敗することが多い (annotatorに連絡してsubmitできるようにしてもらう必要がある) Validationに関して Bond angle outlier bond angleはsigmaが異様に小さい場合があり，そういう場合は無視して良いと思われる． 例えばARGのNE-CZ-NH1角は，PDBの理想値は120.3±0.5となっているが，CCP4 monomer libraryでは120.1±1.5である．このため，例えばモデルの値が117.6だったとき，CCP4的には1.7σの残差だが，PDB的には5.4σも理想値からずれていることになり，怒られる．どちらのsigmaが妥当かは分からないが，0.5°は厳しすぎる． ...

私の興味範囲のみのまとめです． 公式情報: https://www.ccp4.ac.uk/ccp4-8-0-updates/ リリースノートに記載されてないアップデートも実際にはあるのがCCP4 update. 実際にどのファイルが変わったのかは $CCP4/restore/update.log に記録されています． 8.0.019 2024-04-11公開 monomer library 下記の問題のために古いバージョンに戻されてしまった https://github.com/keitaroyam/servalcat/issues/14 8.0.018 2024-04-03公開 新しいServalcatが入るはずでしたが，GEMMIのビルドシステム変更のためにGEMMIのバージョンを上げられず， 最新のGEMMIを必要とするServalcatも次のCCP4メジャーアップデートまで見送りとなりました． このためAcedrgのアップデートも見送りとなりました． 既知のバグが直ってないままになるのは残念だけど仕方ない．早く次のバージョンが出ますように Coot 0.9.8.93 https://www.mail-archive.com/coot@jiscmail.ac.uk/msg05645.html Mutate Residue RangeでDNAがRNAになってしまう問題の修正 TYRのHH原子に関する問題の修正: https://github.com/pemsley/coot/commit/5b52df6509afd137d72c604a6dde1232fe6cff70 monomer library https://github.com/MonomerLibrary/monomers/commits/ccp4-8.0.018 AX/RXを追加．phaserで使われる未知の異常散乱/実原子．なぜか今までのRefmacでは定義無しで動いていたらしい？ #51 CCDにおける原子名の変更などを反映 #44 TYRのhh1 (OHのねじれ角)が間違っていたのを修正 #45 PROのNがsp2になってしまっていた問題の修正．すべてのP-peptideをプロトン化(NH2+)状態に #38 B12の修正 #41 PDBによるペプチド主鎖原子の修正を反映 #40 Robbieによる修正 #35 8.0.017 2024-01-18公開 Refmac 5.8.0425 http://fg.oisin.rc-harwell.ac.uk/scm/loggerhead/refmac/5.8/revision/487 残基数の最大値を50000から100000に symmetry related external bond distanceが正しく機能しないバグを修正 LINKヘッダに書く順序で認識されたりされなかったりするバグが発生しており，直っていない．Refmacatを使えば問題なし． CCP4i2 Refmac5インタフェースでRefmacatが使われるようになり，中性子構造の精密化のインタフェースが追加された． 残念ながらServalcatが古いままで，最近行った様々なバグ修正が反映されてない． Windowsでgemmi 0.6.4が正しく動かないバグが見つかったことが原因． doubleHelix 新プログラム https://doi.org/10.1093/nar/gkad553 8.0.016 2023-10-05公開 Coot 0.9.8.92 https://www.mail-archive.com/coot@jiscmail.ac.uk/msg05575.html 前回のad hocな解決方法ではDNA-SERやdisulfでも問題が起きることが分かったので，ついにちゃんとした修正が入りました． 一部まだ挙動が怪しいようですが． monomer library https://github.com/MonomerLibrary/monomers/commits/ccp4-8.0.016 ホスホジエステル結合(p)にtorsion angle alphaが抜けていたので追加 #33 Acedrg 278で20種のアミノ酸を更新 #32 すでにCCD (PDBの低分子ライブラリ)に存在しないエントリを削除 #30 gemmi 0.6.0以前でエラーになる項目を修正 #29 2.は一部のアミノ酸でねじれ角に不適切な理想値が設定されてしまっていた問題の修正です．特にMETやLEUなど．chi*の名前はREFMACでも使われるので，ちょっと大きな問題でした．同時に，nucleus distanceも最新のテーブルの値になっています． ...

とりさんから一年前に約束した当該機能まだですかと聞かれ，全く身に覚えが無かったものの，実装したらDVD買ってあげるということで書きました． Coot 0.9.8.7で動作確認． def change_elem_ext(): def change_elem(atom_spec, elem): imol, chain, resi, ins, name, alt = atom_spec[1:] set_atom_string_attribute(imol, chain, resi, ins, name, alt, "element", elem) def f(elem): add_status_bar_text("Click an atom") user_defined_click(1, lambda x: change_elem(x, elem)) generic_single_entry("Element?", "", "OK", lambda text: f(text)) coot_toolbar_button("Change elem", change_elem_ext, None) 上記pythonコードを適当な名前(.py)で保存して Calculate - Run Script から実行するか，~/.coot-preferences 以下に置いておくと起動時に読まれます． 成功すればtoolbarに"Change elem"のボタンが追加されます． ちなみにこういったスクリプトを書くときにAPIを調べる方法ですが，体系的には把握してないので毎回ソースコードから探してます． Python側で定義されてる場合は.pyのファイルにあるのでそこから探し，今回は見つからなかったので src/coot_wrap_python.cc の方を見てatomでgrepしたらset_atom_attributeを見つけました． これは実数値を変更する関数だったのですが，下の方にset_atom_string_attributeを発見し，attribute_nameとしてelementを受け付けてくれることが分かりました． https://github.com/pemsley/coot/blob/abe4c5d9e372b91ba43d4bd9a7e924042c7f00ff/src/molecule-class-info-other.cc#L550 文字列をユーザに入力させるとか，クリックした原子の情報を取得するとかいった部分は，類似の他の機能のコードを探して真似します． 特に拡張的な機能はたいていpythonで書かれています．

順次追記予定 注: Coot 0.9.8.92以上 (CCP4 8.0.016以上)を使ってください．0.9.8.6は核酸単体，0.9.8.7-8は核酸-タンパク複合体の精密化に深刻なバグがあります．0.9.8.91でもDNA-SERやDNA-TYRに由来する問題あり． base pair & stacking restraints LIBGを使う CCP4に含まれているプログラムLIBGで作製したrestraintをCootに読ませて使うことができる． LIBGの文献は以下を参照 Brown et al. “Tools for macromolecular model building and refinement into electron cryo-microscopy reconstructions” Acta Cryst. (2015). D71, 136-153 Kovalevskiy et al. “Overview of refinement procedures within REFMAC5: utilizing data from different sources” Acta Cryst. (2018). D74, 215-227 注: mmCIFの読み込みバグがCCP4 8.0.011で修正されているので，当バージョンの使用を推奨． 使い方は libg と打てば表示される．PDBファイルを使う場合は libg -p input.pdb -o libg.txt とするとRefmacのexternal restraintsが記述されたlibg.txtが作られる． 以下は2zm5を使った例 (抜粋): exte dist first chain C resi 23 ins . atom N6 second chain C resi 12 ins . atom O4 value 2.94 sigma 0.15 type 1 exte dist first chain C resi 23 ins . atom N1 second chain C resi 12 ins . atom N3 value 2.84 sigma 0.1 type 1 exte torsion first chain C resi 23 ins . atom C2 second chain C resi 23 ins . atom N1 third chain C resi 12 ins . atom N3 fourth chain C resi 12 ins . atom C4 value 180 sigma 15 type 1 ... exte stac plan 1 firs resi 10 ins . chai C atoms { C1' N9 C8 N7 C5 C6 O6 N1 C2 N2 N3 C4 } plan 2 firs resi 11 ins . chai C atoms { C1' N1 C2 O2 N3 C4 N4 C5 C6 } dist 3.4 sddi 0.2 sdan 6.0 type 1 ... 塩基対の水素結合に対して結合距離制約(exte dist)が，塩基のスタッキングに対して平面間角度と距離の制約(exte stac)が作られる． 記法についてはrefmac keywords (version 5.5.0026 and later)を参照． ...

https://doi.org/10.1107/S2059798323002413 Servalcat論文第二弾出ました．去年のCCP4 study weekendのproceedingsの一部です．同じvirtual issueにCCP4の新しい論文も載ります． REFMAC5は20年以上前に書かれたプログラムで，特にその中のMAKECIFと呼ばれる機能がモデル(pdb/mmcif)を読んでrestraintを作るんですが，歴史的経緯からロジックが必要以上に複雑になっており，メンテナンスがしづらくなっていました． そこでいま開発が活発に進んでいるGEMMIのライブラリを利用して置き換えようというのが今回の内容です． これによって，例えば今までREFMAC5が扱えなかったmicroheterogeneityなどが扱えるようになりました． また，REFMAC5は最初にmonomer library (今や35298個ものmonomerが収録されています)の中身を全部読んでから処理する仕組みだったりして精密化がスタートするまで時間がかかっていたんですが，この部分が一瞬で終わるようになりました．精密化計算の方がずっと時間はかかるので，全体の中ではわずかな改善ですが． 詳細は論文にありますが，この機能はServalcatに実装された"refmacat"コマンドを通じて使えるようになっています． 先日出たCCP4 8.0.011にこっそりとServalcat 0.4.0が入ったので，CCP4をアップデートすれば使えるようになります． 今まで refmac5 xyzin ... のようにして使っていたのを， refmacat xyzin ...のようにコマンド名だけ置き換えればそのまま動くというわけです． モデルファイルの読み込みはGEMMIライブラリを通じて行われるので，GEMMIがサポートしている形式なら何でも読めることになります． Gzipのまま読めるのは地味に嬉しいポイントです． 出力モデルもGEMMIを使って再出力されるので，hybrid-36によりPDB形式でもほとんどの巨大分子が扱えます． 結局各種プログラムのmmCIFへの対応状況が悪く，hybrid-36 PDBの方が勝手が良いというのが現状かと思われます． 未対応の問題として，4文字以上の残基名はREFMAC5ではまだ使えません． そう遠くない未来に3文字IDが枯渇するので，各種プログラムは対応を迫られています． それまでにREFMAC5が直るか，あるいはServalcatがREFMAC5をすっかり置き換えるようになっているかもしれません(？) またこうなるとPDB形式は本当に使えなくなるので，ついにmmCIFへの移行も必要です． とはいえ手元で何らかの3文字IDで作り直してしまえば動かせるので，結局みんな何とかしてPDB形式を使い続けるかもしれませんが．

Servalcat 0.4.6からC++コードを含むようになり，少々配布方法が複雑になってしまった． C++側からもPython側からもGEMMIが必要なため，(Pythonライブラリとしての) GEMMIと同じコンパイラでビルドしなければGEMMIオブジェクトの互換性が失われる． PyPI上で配布している(pipでインストールされる)バイナリは本家GEMMIと同一方法でビルドしているので，パッケージ版は問題なく動くはず． ただし特殊なプラットフォームでは自力でのビルドが必要になる．またPython 3.6以前を使う場合も自力でのビルドが必要である． パッケージ版 上述のように， python -mpip install -U servalcat で入るはず．ただしここでビルドが始まった場合(= 環境に合う配布バイナリが無かった)は，一見うまくいってもGEMMIを配布バイナリからインストールしてる場合は動かない． 下記の「自分でビルドする場合」を参照． Github最新版 Github最新版 (パッケージ版に未収録)をインストールする場合，Github Actionsで自動ビルドされたバイナリを使うことができる． Githubにログインした状態でないとダウンロードできないので注意． https://github.com/keitaroyam/servalcat/actions/workflows/ci.yml のworkflow runsのうち，最新のrunをクリック Artifactsから"wheels2"をダウンロード zipを解凍し，以下のコマンドを実行． wheels2はzipの中身が展開された場所(ディレクトリ)の名前 python -mpip install -U --find-links=wheels2/ servalcat 自分でビルドする場合 まずパッケージ版のGEMMIを自分でビルドする． Github最新版のGEMMIは動かない可能性があるので注意(Servalcatはパッケージ版の最新のGEMMIで動くように作っている)． どのようにしても良いが，pipに --no-binary オプションがあるのでそれが使えるかも（未検証） python -mpip install --force-reinstall --no-binary gemmi gemmi それからGithub最新版のServalcatをソースからビルド python -mpip install -U git+https://github.com/keitaroyam/servalcat.git これで動くはず(未検証)

とりさんのブログを参考に，Hugoを使ったblogを始めてみることにした． Twilogの停止により，twitterに有用なメモを残していくことができなくなったため． テーマにはPaperModを選択．カスタマイズメモ Search page Add menu to site Math typesetting Use Lastmod with PaperMod Shortcodes (tweetやyoutubeの埋め込み方) 数式のテスト: $$ F(s) = \int \rho(x) e^{2\pi i x^T s} dv $$