使用VMD结合freeSASA计算分子的溶剂可及表面区域(SASA)

ene

前几天在群里看到有群友问怎么表征蛋白质二聚体的接触面积，正好赶上最近工作中也需要计算分子的SASA（溶剂可及表面区域），因此就写了两个VMD脚本，用来分析动力学轨迹中SASA的变化情况。

SASA的计算不算困难，Gromacs程序本身就支持通过轨迹直接计算SASA。而VMD也内置了直接计算SASA的命令，使用一些免费程序，例如freeSASA同样可以进行计算。由于个人没使用过gmx，因此今天只讨论如何使用VMD和 freeSASA计算溶剂可及表面区域。在VMD中，使用measure sasa 命令可以直接对指定原子的SASA进行计算，具体的命令从VMD手册中可以获得：



简单对轨迹中的frames进行循环，就可以把每一帧指定区域的SASA值输出到文件中。脚本如下：
       vmd_SASA.tcl (1 KB, 下载次数 Times of downloads: 208)

2018-7-17 09:10 上传 Uploaded

点击下载
Click to download

      开源软件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计算，脚本如下：
       freeSASA.tcl (3.47 KB, 下载次数 Times of downloads: 79)

2018-7-17 09:11 上传 Uploaded

点击下载
Click to download

  简单说一下脚本如何使用，首先freeSASA的脚本只能在linux 下使用，而且在使用前必须要安装freeSASA。准备就绪后，首先在VMD中载入轨迹和结构文件并且将脚本放在当前目录下，随后在VMD的Tcl命令行中输入source freeSASA.tcl。首先程序会提示你输入需要进行计算的SASA的区域所处的结构，至于为什么会有这个要求后面会提到。下一步输入你需要计算的结构。这两步输入均使用标准的Tcl正则表达式，且不能包含”atomselect top”关键字。之后程序会要求你输入计算SASA的频率，对于非常长时间的模拟来说，往往没有必要对每一帧都计算SASA。接下来是算法和相关参数的选择。如果没有特殊计算需要，一路回车下来也可以。最后一个输入选项是计算所使用的CPU核心数，freeSASA在默认算法下不能超过16个核心。再一次确认后，即可开始计算。

计算时屏幕会不断滚动如下信息：

这是正常现象，不用担心。当计算完成之后，程序会输出 AllDone！以提示。之后就会在当前目录下找到一个以SASA开头的dat文件，其中就是你选定的每一帧的SASA值。

下面是真正有意思的地方，如果对比freeSASA与VMD的SASA插件计算同一结构所得到的数值，会发现他们之间可能具有相当大的差距。例如，计算G蛋白(PDBID: 5vai)的A链时，freeSASA所给出的结果大概是11534，而VMD所给出的结果则是15607，两者相差几乎达到36%。这显然是错误的结果，即使是算法不同也不可能造成如此大的误差。

      最后，楼主发现这是两种程序所采用的计算“环境”不同所导致的。freeSASA在进行计算时会默认将所选择结构周边的结构一同考虑进去进行SASA的计算，而VMD会将所选择的结构从整体结构中“抽出“，放在单纯的溶剂环境下进行计算。就好比是G蛋白的A链，freeSASA计算时认为A链与其余B、G和N三条链处于结合的状态，而VMD则认为A链裸露在溶剂环境下。为了证实这一点，可以用VMD从整个蛋白结构中提取出单独的A链，使用freeSASA进行计算。（图中红色为A链）

   

可以看到，结果在15322左右，与VMD的15607之间误差不到2%，基本可以认为是允许范围内由算法不同引起的误差了。因此，这就是为什么在freeSASA的脚本中必须定义“环境”这个概念的原因了。楼主编写的VMD脚本首先要将每一帧的结构写成临时PDB文件，再通过freeSASA进行运算。如果将每一帧结构（包括水分子，脂膜等等）都写入临时PDB，只会徒增计算量。因此计必须首先明确哪些结构会对最终的结果造成影响。例如计算水中的蛋白质复合物中某条链的SASA，那么第一条输入的就应该是简单的”protein”，若要考察膜蛋白整体的SASA，那么第一条输入写为“protein or lipid”，然后第二条输入写为“protein”似乎更加合理。从某种角度上来说，VMD计算的是SASA可能的“极大值“，而freeSASA计算的是"真实情况”下的SASA。但VMD能够根据电荷计算不同原子的SASA，而且速度高很多。两者各有利弊，实际研究中选用时，还要结合自身研究特点，谨慎选用。以上说法多有不严谨之处，还请各位老师指正

Chenglong_li

多谢

sobereva

值得一提的是，利用-restrict，也可以计算出整个蛋白中某个链裸露在溶剂区的面积，例

ene

sobereva 发表于 2018-7-17 13:36
值得一提的是，利用-restrict，也可以计算出整个蛋白中某个链裸露在溶剂区的面积，例

原来如此，是我粗心大意了。谢谢老师

1130240115

谢谢分享

lijiayisjtu

ene 发表于 2018-7-28 13:54
原来如此，是我粗心大意了。谢谢老师

您好，我把freeSASA.tcl 脚本和 prmtop文件，以及轨迹文件放在一起，在VMD的tcl 输入却找不到，这该如何解决呢，谢谢

ene

lijiayisjtu 发表于 2018-8-1 19:26
您好，我把freeSASA.tcl 脚本和 prmtop文件，以及轨迹文件放在一起，在VMD的tcl 输入却找不到，这该如何 ...

正常来说，不可能存在这种情况。不行的话，可以试试在脚本前面加上绝对路径

zdwssg123

sobereva 发表于 2018-7-17 13:36
值得一提的是，利用-restrict，也可以计算出整个蛋白中某个链裸露在溶剂区的面积，例

sob老师，这个vmd除了蛋白质的疏水和亲水的sasa，能计算其他的分子类型吗，为什么我用这种方法，算出来的疏水面积是0呢？还请sob老师指教一下，感激不尽!

sobereva

zdwssg123 发表于 2018-8-23 19:10
sob老师，这个vmd除了蛋白质的疏水和亲水的sasa，能计算其他的分子类型吗，为什么我用这种方法，算出来的 ...

只对蛋白质能用。VMD对标准氨基酸划分成了疏水和非疏水，所以才能直接那么写。其它的分子，你得自行划定哪些原子/残基是疏水哪些是亲水，然后在-restrict后面用相应的选择语句选择

zlyuan

sobereva 发表于 2018-7-17 13:36
值得一提的是，利用-restrict，也可以计算出整个蛋白中某个链裸露在溶剂区的面积，例

sob老师，我想利用-restrict，计算半胱氨酸残基中巯基—SH基团的溶剂可及性，怎么实现？

sobereva

zlyuan 发表于 2019-8-19 16:44
sob老师，我想利用-restrict，计算半胱氨酸残基中巯基—SH基团的溶剂可及性，怎么实现？

-restrict后面加上基团的选择语句就完了
参考
VMD里原子选择语句的语法和例子
http://sobereva.com/504（http://bbs.keinsci.com/thread-14267-1-1.html）

臧臧

请问老师，我计算按照vmd的tcl文件计算得到的结果均为0，这该怎么解决呢？

Eridium

请问怎么用VMD计算蛋白质侧链中半径R以内（比如5埃）每一个原子的sasa呢？