计算化学公社
标题: 使用VMD结合freeSASA计算分子的溶剂可及表面区域(SASA) [打印本页]
作者Author: ene 时间: 2018-7-17 09:25
标题: 使用VMD结合freeSASA计算分子的溶剂可及表面区域(SASA)
本帖最后由 ene 于 2020-10-30 22:00 编辑
前几天在群里看到有群友问怎么表征蛋白质二聚体的接触面积，正好赶上最近工作中也需要计算分子的SASA（溶剂可及表面区域），因此就写了两个VMD脚本，用来分析动力学轨迹中SASA的变化情况。
SASA的计算不算困难，Gromacs程序本身就支持通过轨迹直接计算SASA。而VMD也内置了直接计算SASA的命令，使用一些免费程序，例如freeSASA同样可以进行计算。由于个人没使用过gmx，因此今天只讨论如何使用VMD和 freeSASA计算溶剂可及表面区域。在VMD中，使用measure sasa 命令可以直接对指定原子的SASA进行计算，具体的命令从VMD手册中可以获得：
(, 下载次数 Times of downloads: 186)
简单对轨迹中的frames进行循环，就可以把每一帧指定区域的SASA值输出到文件中。脚本如下：
(, 下载次数 Times of downloads: 232)
开源软件freeSASA也可以用来进行溶剂可及表面区域计算，可以在http://freesasa.github.io/免费下载安装。编译安装过程很简单，一般的报错根据屏幕输出的提示，禁用json和xml之后一般都能成功编译。和VMD以及gmx不同，freeSASA支持两种SASA计算方法：Shrake & Rupley算法和Lee & Richards算法，后者为默认算法。具体的原理细节请参考官网说明文档。同时，freeSASA的原子选择表达式应用的是pymol规则，具体可以看pymolwiki: https://pymolwiki.org/index.php/Selection_Algebra上面的说明。pymol的规则虽然不算非常热门，但和VMD的Tcl语言也有一定相似的地方。可偏偏freeSASA使用pymol的选择语句也没有模仿到家，最小的可选择结构只能是残基，也就是没法指定原子来进行计算（也很有可能是我没有研究透彻，还希望各位老师不吝赐教）——这就直接导致freeSASA不能够用于亲水区域和憎水区域面积的计算，因为没法通过电荷指定原子（至于残基平均电荷能不能代替原子电荷来进行亲水/憎水区域的判断，这个还有待讨论）。抛除这一点，freeSASA在单纯的计算溶剂可及表面区域的问题上表现不错（记得好像有一篇Nature就使用了gmx+freeSASA表征靶点口袋的开闭状态），结合VMD的Tcl语言同样可以用于轨迹中各帧的SASA计算，脚本如下：
(, 下载次数 Times of downloads: 91)
简单说一下脚本如何使用，首先freeSASA的脚本只能在linux 下使用，而且在使用前必须要安装freeSASA。准备就绪后，首先在VMD中载入轨迹和结构文件并且将脚本放在当前目录下，随后在VMD的Tcl命令行中输入source freeSASA.tcl。首先程序会提示你输入需要进行计算的SASA的区域所处的结构，至于为什么会有这个要求后面会提到。下一步输入你需要计算的结构。这两步输入均使用标准的Tcl正则表达式，且不能包含”atomselect top”关键字。之后程序会要求你输入计算SASA的频率，对于非常长时间的模拟来说，往往没有必要对每一帧都计算SASA。接下来是算法和相关参数的选择。如果没有特殊计算需要，一路回车下来也可以。最后一个输入选项是计算所使用的CPU核心数，freeSASA在默认算法下不能超过16个核心。再一次确认后，即可开始计算。
(, 下载次数 Times of downloads: 184)
(, 下载次数 Times of downloads: 161)
计算时屏幕会不断滚动如下信息：
(, 下载次数 Times of downloads: 142)
这是正常现象，不用担心。当计算完成之后，程序会输出 AllDone！以提示。之后就会在当前目录下找到一个以SASA开头的dat文件，其中就是你选定的每一帧的SASA值。
下面是真正有意思的地方，如果对比freeSASA与VMD的SASA插件计算同一结构所得到的数值，会发现他们之间可能具有相当大的差距。例如，计算G蛋白(PDBID: 5vai)的A链时，freeSASA所给出的结果大概是11534，而VMD所给出的结果则是15607，两者相差几乎达到36%。这显然是错误的结果，即使是算法不同也不可能造成如此大的误差。
(, 下载次数 Times of downloads: 143)
(, 下载次数 Times of downloads: 147)
最后，楼主发现这是两种程序所采用的计算“环境”不同所导致的。freeSASA在进行计算时会默认将所选择结构周边的结构一同考虑进去进行SASA的计算，而VMD会将所选择的结构从整体结构中“抽出“，放在单纯的溶剂环境下进行计算。就好比是G蛋白的A链，freeSASA计算时认为A链与其余B、G和N三条链处于结合的状态，而VMD则认为A链裸露在溶剂环境下。为了证实这一点，可以用VMD从整个蛋白结构中提取出单独的A链，使用freeSASA进行计算。（图中红色为A链）
(, 下载次数 Times of downloads: 151)
(, 下载次数 Times of downloads: 155)
(, 下载次数 Times of downloads: 167)
可以看到，结果在15322左右，与VMD的15607之间误差不到2%，基本可以认为是允许范围内由算法不同引起的误差了。因此，这就是为什么在freeSASA的脚本中必须定义“环境”这个概念的原因了。楼主编写的VMD脚本首先要将每一帧的结构写成临时PDB文件，再通过freeSASA进行运算。如果将每一帧结构（包括水分子，脂膜等等）都写入临时PDB，只会徒增计算量。因此计必须首先明确哪些结构会对最终的结果造成影响。例如计算水中的蛋白质复合物中某条链的SASA，那么第一条输入的就应该是简单的”protein”，若要考察膜蛋白整体的SASA，那么第一条输入写为“protein or lipid”，然后第二条输入写为“protein”似乎更加合理。从某种角度上来说，VMD计算的是SASA可能的“极大值“，而freeSASA计算的是"真实情况”下的SASA。但VMD能够根据电荷计算不同原子的SASA，而且速度高很多。两者各有利弊，实际研究中选用时，还要结合自身研究特点，谨慎选用。以上说法多有不严谨之处，还请各位老师指正
作者Author: Chenglong_li 时间: 2018-7-17 09:55
多谢
作者Author: sobereva 时间: 2018-7-17 13:36
值得一提的是，利用-restrict，也可以计算出整个蛋白中某个链裸露在溶剂区的面积，例
(, 下载次数 Times of downloads: 170)
作者Author: ene 时间: 2018-7-28 13:54
原来如此，是我粗心大意了。谢谢老师
作者Author: 1130240115 时间: 2018-7-29 09:05
谢谢分享
作者Author: lijiayisjtu 时间: 2018-8-1 19:26
您好，我把freeSASA.tcl 脚本和 prmtop文件，以及轨迹文件放在一起，在VMD的tcl 输入却找不到，这该如何解决呢，谢谢
作者Author: ene 时间: 2018-8-1 21:12
正常来说，不可能存在这种情况。不行的话，可以试试在脚本前面加上绝对路径
作者Author: zdwssg123 时间: 2018-8-23 19:10
sob老师，这个vmd除了蛋白质的疏水和亲水的sasa，能计算其他的分子类型吗，为什么我用这种方法，算出来的疏水面积是0呢？还请sob老师指教一下，感激不尽!
作者Author: sobereva 时间: 2018-8-24 17:10
只对蛋白质能用。VMD对标准氨基酸划分成了疏水和非疏水，所以才能直接那么写。其它的分子，你得自行划定哪些原子/残基是疏水哪些是亲水，然后在-restrict后面用相应的选择语句选择
作者Author: zlyuan 时间: 2019-8-19 16:44
sob老师，我想利用-restrict，计算半胱氨酸残基中巯基—SH基团的溶剂可及性，怎么实现？
作者Author: sobereva 时间: 2019-8-20 22:21
-restrict后面加上基团的选择语句就完了
参考
VMD里原子选择语句的语法和例子
http://sobereva.com/504（http://bbs.keinsci.com/thread-14267-1-1.html）
作者Author: 臧臧 时间: 2022-6-6 16:29
请问老师，我计算按照vmd的tcl文件计算得到的结果均为0，这该怎么解决呢？
作者Author: Eridium 时间: 2022-10-16 22:35
请问怎么用VMD计算蛋白质侧链中半径R以内（比如5埃）每一个原子的sasa呢？
| 欢迎光临 计算化学公社 (http://ccc.keinsci.com/) |
Powered by Discuz! X3.3 |