pygamessがTD-DFTとNMR計算をサポートしました

なんとなくこの本の15章をやりたくなったので、朝からゴニョゴニョしてた。色々できたのでバージョン上げておいた。

TD-DFT

データはおなじみのPubchemQCからダウンロードしてきてopenbabelでmolファイルにコンバートしたものを利用。

>>> from pygamess import Gamess
>>> from rdkit import Chem
>>> m = Chem.MolFromMolFile("examples/methyl_yellow.mol", removeHs=False)
>>> g = Gamess()
>>> g.dft_type("b3lyp", tddft=True)
>>> g.basis_sets("6-31G*")
>>> r = g.run(m)
>>> r.uv_spectra # (exitation ev, oscillator strength)

NMR spectra

これは時間がかかった(2020 MBAで30分以上)

>>> from pygamess import Gamess
>>> from rdkit import Chem
>>> m = Chem.MolFromMolFile("examples/C=CCBr.mol", removeHs=False)
>>> g = Gamess(num_cores=1) # PARALLEL EXECUTION IS NOT ENABLED.
>>> g.basis_sets("6-31G*") 
>>> g.run_type("nmr")
>>> r = g.run(m)
>>> r.isotropic_shielding

独習量子化学計算について

わかりやすくて良い本だった。ケミストだったら本書に沿って一通り手を動かしておけば、後々役に立つのではないかと思った。

こっちはいわゆるクックブックみたいな感じ

Free Ligand 1D NMR Conformational Signatures To Enhance Structure Based Drug Design

This was a very interesting paper.

The correlation of kon with free ligand preorganization implies 1D NMR signatures can augment structure−kinetic relationships; when kon is optimized by highly populated bioactive states, potency will be driven by koff.

Storing GAMESS calculation results in RDKit Chem object

Japan is under a state of emergency for the third time, which forced me to do self-quarantine. After all, my five days of GW holidays were spent developing the pygamess library.

This is an example of storing GAMESS calculation results into RDKit's Chem object. I will calculate the two molecules, ethane, and ethanol, to export as an SDF file later.

>>> from pygamess import Gamess
>>> from pygamess.utils import rdkit_optimize
>>> from rdkit import Chem
>>> ethane = rdkit_optimize("CC")
>>> ethanol = rdkit_optimize("CCO")
>>> g = Gamess()
>>> g.run_type("optimize")
>>> ethane_result = g.run(ethane)
>>> optimized_ethane = ethane_result.mol
>>> ethanol_result = g.run(ethanol)
>>> optimized_ethanol = ethanol_result.mol

Orbital energies and dipole moments including HOMO and LUMO energies are stored in the object. Each atom is also assigned a mulliken charge/lowdin charge.

>>> optimized_ethane.GetProp("HOMO")
'-0.46029999999999999'
>>> optimized_ethane.GetProp("LUMO")
'0.65559999999999996'
>>> optimized_ethane.GetProp("dipole_moment")
'1.9999999999999999e-06'
>>> optimized_ethane.GetProp("orbital_energies")
'-11.0355 -11.0352 -0.981 -0.8121 -0.572 -0.572 -0.4717 -0.4603 -0.4603 0.6556 0.6556 0.687 0.737 0.7814 0.7893 0.7893'
>>> for a in optimized_ethane.GetAtoms():
...   print("{}:\t{:.4f}".format(a.GetSymbol(), float(a.GetProp("mulliken_charge"))))
... 
C:  -0.1748
C:  -0.1748
H:  0.0583
H:  0.0583
H:  0.0583
H:  0.0583
H:  0.0583
H:  0.0583

Finally, I'll export the object to SDF.

>>> optimized_mols = [optimized_ethane, optimized_ethanol]
>>> w = Chem.SDWriter('test.sdf')
>>> for m in optimized_mols:
...   w.write(m)
... 
>>> w.close()

Looking at the exported SDF, we can see these results are persisted.

$ cat test.sdf
     RDKit          3D

  8  7  0  0  0  0  0  0  0  0999 V2000
   -0.7687    0.0082   -0.0152 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7687   -0.0082    0.0152 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1742   -0.1016    0.9863 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1352    0.9432   -0.4285 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1499   -0.8048   -0.6260 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1499    0.8048    0.6260 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1742    0.1016   -0.9863 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1352   -0.9432    0.4285 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0
  1  3  1  0
  1  4  1  0
  1  5  1  0
  2  6  1  0
  2  7  1  0
  2  8  1  0
M  END
>  <total_energy>  (1) 
-78.306179626800002

>  <HOMO>  (1) 
-0.46029999999999999

>  <LUMO>  (1) 
0.65559999999999996

>  <nHOMO>  (1) 
-0.46029999999999999

>  <nLUMO>  (1) 
0.65559999999999996

>  <dipole_moment>  (1) 
1.9999999999999999e-06

>  <dx>  (1) 
-0

>  <dy>  (1) 
-1.9999999999999999e-06

>  <dz>  (1) 
0

>  <orbital_energies>  (1) 
-11.0355 -11.0352 -0.981 -0.8121 -0.572 -0.572 -0.4717 -0.4603 -0.4603 0.6556 0.6556 0.687 0.737 0.7814 0.7893 0.7893

>  <atom.dprop.mulliken_charge>  (1) 
-0.174764 -0.174764 0.058255000000000001 0.058255000000000001 0.058254 0.058255000000000001 0.058255000000000001 0.058255000000000001

TODO: Add GAMESS calculation conditions as attributes, like "HF/6-31G**".

getGEOで取得した後、不要なサンプルを除去する

とりあえず、ノーマライズされたデータをサクサク再解釈できればいいので、getGEOから。これでとってきたデータに不適切なサンプルが混じったりすることがあるので、最初にサクッと除きたい

gse <- getGEO("GSE138458")
removes <- c("GSM4109031", "GSM4109192", "GSM4109193", "GSM4109194", "GSM4109223", "GSM4109271")
gse <- gse[,-which (colnames(gse) %in% removes)]

とやれば綺麗なデータセットになり嬉しい。

RでGEOのデータ解析チュートリアルをやってみた

朝からQM関連の論文読みまくってちょっと疲れたので写経してみた。Rは作図が楽だし大抵きれいに描いてくれるから良いですよね。

このチュートリアルは大変わかりやすいですが、dplyrとかlimmaはわかっている前提なので別途他のドキュメントで保管する必要があるかもしれないです。

それから、写経して気づいたんですがvolcano plotのとこ多分こうだと思います。

full_results %>% 
+     mutate(Significant = adj.P.Val < p_cutoff & abs(logFC) > fc_cutoff ) %>% 
+     mutate(Rank = 1:n(), Label = ifelse(Rank < topN, Symbol,"")) %>% 
+     ggplot(aes(x = logFC, y = B, col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")

GSE33126

一通りやったら、YouTubeにもう少し複雑な例を同じ流れでやっている動画があるのでみると更に理解が深まると思います。

このタイミングでまたバイオインフォマティクスにコミットしなきゃならんのかよーとか思っていたけど、実際やってみると色々と新しい発見もあるし生物学はやっぱり面白いですね。

独習 Pythonバイオ情報解析はPythonでシングルセル解析をやりたい人が読めば良いと思う

シングルセル解析をPythonでやりたい人には色々と役に立つことが書いてあると思います。私にとっては7, 10-12章が良かったです。1章から6章はPythonとかPandasの使い方なのでスキップしました。

scanpyの使い方が気になっていたのですが、付録Bに10ページくらいの解説だったけど、もう少しページを割いてもらえると嬉しかったかも。まぁオフィシャルドキュメント読めばいいんだろうけど。

全体的に満足な内容でしたが、

  1. リードのトリミング
  2. ゲノム配列のマッピング
  3. マッピングされたリード数の数え上げ

といった部分は割愛されているのでRNA-Seqデータ解析の内容くらいは理解してある必要があるかもしれません。

独習と銘打ってありますが、Pythonにある程度なれていないとハマるかもしれないなぁと思いました。本書で動作確認したPythonやCondaのバージョンとか各ライブラリのバージョンが明記されていないので、Pyhon3.9.4だとcondaでumap入らなかったりとかconda-forge指定しないといけなかったりとかちょいちょいありました。

それからjupyter notebookを勧めているのにサンプルコードがpyなのは初心者は戸惑うんじゃなかろうかなぁと。ipynb同梱しなかったのはなにか理由があったんですかね?

  • 出力7-24の出力結果間違ってる?
  • 8.1.2 joypy, matplotlib_vennのインストールには -c conda-forgeオプションが必要だった
  • 10.2.3 コンベンショナルな方法で正規化してたが実際はどうするのがいいのか?
  • 11.4 Python3.9.4ではumapをcondaでインストールできなかった。pip混ぜたくないのでスキップした
  • 12.5 SAUCIE後で論文読む

ここで差がつく! 英文ライティングの技術

この本の2章が大変参考になった。

また、この本はまだ半分くらいしか読んでいないが、かなり役に立っている。

シルバー民主主義というよりは民意の高齢化

色々思うところはある。

筆者は現在の状況を政治がマジョリティ(高齢化の世代)によっているだけで主権が高齢者にあるわけではないのでシルバー民主主義という立場をとっていない。というより民意が高齢化しているという表現をしていた。

  • シルバー民主主義とシルバーファースト減少の違いは高齢者に政策決定の主導権があるかどうか
  • シルバー優遇の最たるものが公的年金制度。賦課方式は勤労世代から高齢世代への所得移転
  • 賦課方式が可能なのは経済や人口が右肩上がりの場合

民意の高齢化は当分続くわけだから、ドラスティックに変わらない限り色々難しいことになるんじゃないかなぁ。定年引き上げなんて、生産年齢人口をせこく水増ししてるだけだから根本的な解決にならないだろうし、年金だって払い損だろうしねぇ。

  • 余命比例投票はどうか?
  • シルバー民主主義の下では様々な形で負担の先送りが行われる
  • 「若者も安心できる年金制度を」は根本的に誤っている

AIの話に関しては過剰な期待が含まれているので割り引いて考えないといけないが、まぁ読み物としては面白かった

  • 少子化の要因が6つくらいある
  • 森林税のような少額の課税はこれからも増える
  • 社会保険料は本質的に税金なのに会社員は徴収を実感しづらいから取りやすい
  • AI創薬の話は盛りすぎ(私の感想)
  • 自動車の「所有する価値」は今後どうなるのか?
  • テレワークは少子化対策に有効

期せずしてテレワークが増えたので、出生率増えるんでしょうかね?

ポケモンGO(Lv. 45)

前回から40日くらい。リーダーを50回倒すために下っ端と戦うこと300回。

長かった。

次は卵を孵化させてエクセレントスローを投げまくるだけなので、経験値さえ貯ればレベルアップしそう。

あの子は貴族

読みました。久しぶりに気持ちの良い読後感を得ました。

と同時にBlack mirrorという大変に後味の悪いドラマも見始めてしまった。