【1.1.2】计算残基暴露程度（freesasa）
FreeSASA是用于计算蛋白质，RNA和DNA结构的SASA的高效库。 与其他常用工具相比，它的主要优点是它是开源的，易于配置，并且可以同时用作命令行工具，C库和Python模块。 上面的测试表明，在给定的分辨率下，它的运行速度与某些流行的工具一样或更快，并且可以通过并行化计算来进一步提高。
- github软件下载：http://freesasa.github.io/
- 需要结构
计算溶剂可及表面积（SASA， solvent accessible surface areas）是结构生物学中的常规计算。 尽管有许多程序可用于此计算，但没有为简化工具链集成而设计的独立式开源工具。 FreeSASA是一个用于SASA计算的开源C库，除了其C API外，还提供了命令行和Python接口。 该库同时实现了Lee和Richards以及Shrake和Rupley的近似值，并且高度可配置以允许用户控制分子参数，准确性和输出粒度。 它仅取决于标准C库，因此应易于编译和安装在任何平台上。 该库文件齐全，稳定高效。 命令行界面可以轻松地替换封闭源的旧版程序，具有可比或更高的准确性和速度，并具有某些附加功能。
一、下载和安装
wget -c http://freesasa.github.io/freesasa-2.0.3.tar.gz
tar -xzf freesasa-2.0.3.tar.gz
cd freesasa-2.0.3
./configure --disable-json --prefix=/data/software/freesasa
make
make install
二、使用说明
2.1 输入参数
例子：
[sam@g03 protein]$ freesasa input/6F1H.pdb
## FreeSASA 2.0.3 ##
PARAMETERS
algorithm : Lee & Richards
probe-radius : 1.400
threads : 2
slices : 20
INPUT
source : input/6F1H.pdb
chains : ACDB
model : 1
atoms : 8988
RESULTS (A^2)
Total : 57622.06
Apolar : 31206.08
Polar : 26415.97
CHAIN A : 14632.95
CHAIN C : 14434.11
CHAIN D : 14612.67
CHAIN B : 13942.34
如果需要更高的精度（use 100 slices per atom instead of the default 20）：
freesasa input/6F1H.pdb -n 100
计算SASA用Shrake & Rupley算法， 200 test points, a probe radius of 1.2 Å, 4个线程
freesasa --shrake-rupley -n 200 --probe-radius 1.2 --n-threads 4 3wbm.pdb
使用自定义的配置文件
freesasa --config-file <file> 3wbm.pdb
To use the atomic radii from NACCESS call
freesasa --radii=naccess 3wbm.pdb
Another way to specify a custom set of atomic radii is to store them as occupancies in the input PDB file
freesasa --radius-from-occupancy 3wbm.pdb
2.2 输出例子
- 除了上次的格式，FreeSASA可以通过 –format 输出每个残基的JSON, XML, PDB, RSA, SASA格式，也可以输出每个残基的SASA。
- JSON和XML输出结果的内容，可以通过 –output-depth= 选择在atom, residue, chain and structure输出。如果为atom, SASA values are shown for all levels of the structure, including individual atoms. With chain, only structure and chain SASA values are printed (this is the default).
- The output can include relative SASA values for each residues. To calculate these a reference SASA value is needed, calculated using the same atomic radii. At the moment such values are only available for the ProtOr and NACCESS radii (selected using the option –radii), if other radii are used relative SASA will be excluded (in RSA output all REL columns will have the value ‘N/A’).
- The reference SASA values for residue X are calculated from Ala-X-Ala peptides in a stretched out configuration. The reference configurations are supplied for reference in the directory rsa. Since these are not always the most exposed possible configuration, and because bond lengths and bond angles vary, the relative SASA values will sometimes be larger than 100 %. At the moment there is no interface to supply user-defined reference values.
（原理暂时不理解，后续再来回顾）
2.2.1 json
freesasa --format=xml --output-depth=residue 3wbm.pdb
输出：
{
"source":"FreeSASA 2.0",
"length-unit":"Ångström",
"results":[
{
"input":"3wbm.pdb",
"classifier":"ProtOr",
"parameters":{
"algorithm":"Lee & Richards",
"probe-radius":1.3999999999999999,
"resolution":20
},
"structures":[
{
"chain-labels":"ABCDXY",
"area":{
"total":25190.768387067546,
"polar":13638.391677017404,
"apolar":11552.376710050148,
"main-chain":3337.1622502425053,
"side-chain":21853.606136825045
},
"chains":[
{
"label":"A",
"n-residues":86,
"area":{
"total":3785.4864049452635,
"polar":1733.8560208488598,
"apolar":2051.6303840964056,
"main-chain":723.34358684348558,
"side-chain":3062.1428181017791
}
"residues":[
{
"name":"THR",
"number":"5",
"area":{
"total":138.48216994006549,
"polar":56.887951514571867,
"apolar":81.594218425493622,
"main-chain":38.898190013033592,
"side-chain":99.583979927031905
},
"relative-area":{
"total":104.05152148175331,
"polar":113.98106895325961,
"apolar":98.093554250413092,
"main-chain":96.330336832673567,
"side-chain":107.414496739329
},
"n-atoms":7
},
...
},
...
]
}
]
}
]
2.2.2 XML
示例：
freesasa --format=xml 3wbm.pdb
结果：
<?xml version="1.0" encoding="UTF-8"?>
<results xmlns="http://freesasa.github.io/" source="FreeSASA 2.0" lengthUnit="Ångström">
<result classifier="ProtOr" input="3wbm.pdb">
<parameters algorithm="Lee & Richards" probeRadius="1.400000" resolution="20"/>
<structure chains="ABCDXY">
<area total="25190.768" polar="13638.392" apolar="11552.377" mainChain="3337.162" sideChain="21853.606"/>
<chain label="A" nResidues="86">
<area total="3785.486" polar="1733.856" apolar="2051.630" mainChain="723.344" sideChain="3062.143"/>
</chain>
<chain label="B" nResidues="84">
<area total="4342.334" polar="1957.114" apolar="2385.220" mainChain="853.707" sideChain="3488.627"/>
</chain>
<chain label="C" nResidues="86">
<area total="3961.119" polar="1838.724" apolar="2122.395" mainChain="782.652" sideChain="3178.468"/>
</chain>
<chain label="D" nResidues="89">
<area total="4904.298" polar="2332.306" apolar="2571.991" mainChain="977.459" sideChain="3926.838"/>
</chain>
<chain label="X" nResidues="25">
<area total="4156.455" polar="2919.576" apolar="1236.879" mainChain="0.000" sideChain="4156.455"/>
</chain>
<chain label="Y" nResidues="25">
<area total="4041.076" polar="2856.815" apolar="1184.261" mainChain="0.000" sideChain="4041.076"/>
</chain>
</structure>
</result>
</results>
2.2.3 PDB
…
后续在补充吧。。
三、算法原理介绍
3.1 前言
将非极性分子暴露于水是非常不利的，并且使疏水性自由能最小化是大分子折叠的重要驱动力（Finkelstein＆Ptitsyn，2002）。分子的溶剂可及表面积（SASA）可以衡量分子与溶剂之间的接触面积。尽管表面积和自由能之间的确切定量关系是难以捉摸的，但SASA可用于比较不同分子或同一分子的不同构象，或例如测量由于低聚而被掩埋的表面。
为了定义SASA，让球形探针代表溶剂分子。将探针滚动到较大分子的表面上。探针中心所追踪的表面积是较大分子的SASA（Lee＆Richards，1971）。通常使用两种经典近似方法来计算SASA：
- 一种是Lee和Richards（L＆R）的方法，其中表面通过一组切片的轮廓近似（1971），
- 另一种是Shrake＆Rupley（1973）的方法（S＆R），其中每个球体的表面由一组测试点近似。
通过代表两者的分辨率，可以以任意精度计算SASA。表面积也可以通过分析来计算（Fraczkiewicz＆Braun，1998），当需要梯度时，或通过针对不同目的量身定制的各种其他近似值，这很有用（Cavallo等，2003； Drechsel等，2014； Nordon等人，2014）。 Sanner等，1996； Weiser等，1999； Xu＆Zhang，2009）。
有许多工具可用于计算SASA:
- 最流行的命令行程序可能是NACCESS（Hubbard＆Thornton，1993）（可免费用于学术用途），这是L＆R近似的有效Fortran实现。
- 另一个著名的命令行工具是DSSP，它主要计算蛋白质结构的二级结构和氢键，但也提供SASA（Touw等人，2015）（使用S＆R，开源）。
- 也有一些网络服务可用，例如Getarea（可解析地计算表面）（Fraczkiewicz＆Braun，1998），Triforce使用半分析镶嵌方法（Drechsel等，2014）（也可用于命令行） ）。
- 此外，大多数分子动力学模拟程序包都包括从轨迹分析SASA的工具。
FreeSASA的目的是与NACCESS以及其他许多类似程序保持相同的地位：一个简单，快速的命令行工具，“做一件事情，做得很好”，并且可以轻松地集成到工具链中。 FreeSASA的优点是它是开源的（GNU通用公共许可证3），并且除了命令行界面外还提供C和Python API。它具有合理的默认参数，对临时用户没有强制性的配置（唯一需要的输入是结构），但还可以完全控制所有计算参数。依赖性已降至最低：编译仅需要标准C和GNU库。该库是线程安全的，并且已经付出了一些努力来优雅地处理各种错误。该代码附带了详尽的文档，也可以从 http://freesasa.github.io/doxygen/ 在线获得。尽管功能和可用性是编写该库的主要动机，但性能测试表明，在单个CPU内核上运行时，FreeSASA的速度与传统程序一样快。另外，很大一部分计算已经并行化，在多核处理器上运行时，速度明显提高。
3.2 方法
计算
S＆R和L＆R都非常易于实现：
- 都需要首先确定哪些原子处于接触状态
- 然后计算每个原子与其相邻原子之间的重叠度。
查找contacts是使用单元格列表完成的，这意味着contacts计算是O（N）操作。然后，两种算法都分别处理每个原子，这也是计算O（N）的第二部分。另外，这第二部分是可并行化的。
对于L＆R，不是一次性将整个蛋白质切成（slicing）薄片，而是将每个原子分别切成薄片。因此，L＆R计算由每个原子的切片数来参数化，即小原子的切片比大原子薄。
斐波那契螺旋线（Fibonacci spiral）很好地逼近了球上点的均匀分布（Swinbank＆Purser，2006年），从而可以高效生成任意数量的S＆R测试点。单元格列表为该算法提供了双立方晶格优化中两个晶格的第一个（Eisenhaber等，1995），第二个晶格（用于测试点）目前尚未在FreeSASA中实现。
通过目视检查表面来测试实现的正确性。在两个原子的情况下，将根据分析计算验证结果。另一个验证来自比较使用两种独立算法的高精度SASA计算结果。此外，当使用相同的分辨率和原子半径时，使用L＆R算法可得到与NACCESS相同的结果。
半径分配
计算的重要步骤是为每个原子分配一个半径。 FreeSASA中的默认设置是使用Tsai等人的ProtOr半径 （1999）。 该文库可识别20种标准氨基酸（加上Sec和Pyl）和标准核苷酸（加上一些非标准氨基酸）。 蔡不要提及磷和硒（phosphorus and selenium），这些原子的半径分别为1.8和1.9Å。
默认情况下，蛋白质数据库（PDB）文件中将忽略氢原子和HETATM记录。 如果包含该库，它将识别出三个常见的HETATM条目：乙酰基和NH2封端基团以及水，并为其分配ProtOr半径。 否则，使用元素的范德华半径，取自Mantina等人的论文 （2009）。 对于Mantina等人处理过的44个主要族元素之外的元素，或者如果需要完全不同的半径，用户可以提供自己的配置。
用户可以通过API或通过提供配置文件来指定自己的原子半径.该库附带了一些示例配置文件，包括一个提供NACCESS参数化子集的文件，以及一个具有默认ProtOr参数的文件。 此外，还提供了脚本以根据PDB CONECT条目自动生成ProtOr配置，例如化学成分词典中的那些（Westbrook等，2015）。 然后可以将它们附加到默认配置。
四、我的案例
freesasa input/12E8.pdb -n 100 --depth=residue --format=seq --n-threads=4
结果
SEQ P 206 VAL : 41.73
SEQ P 207 ASP : 97.35
SEQ P 208 LYS : 42.78
SEQ P 209 LYS : 115.82
SEQ P 210 ILE : 2.54
SEQ P 211 VAL : 81.71
SEQ P 212 PRO : 68.03
SEQ P 213 ARG : 157.27
SEQ P 214 ASP : 22.14
第一列代表用的序列模式，第二列是链，第三列是UID，第四列代表残基，最后一列是对应的ASA值
案例2 ：
[sam@g03 protein]$ freesasa input/12E8.pdb --format=rsa --n-threads=4 --radii=naccess|head -n 20
REM FreeSASA 2.0.3
REM Absolute and relative SASAs for input/12E8.pdb
REM Atomic radii and reference values for relative SASA: NACCESS
REM Chains: LHMP
REM Algorithm: Lee & Richards
REM Probe-radius: 1.40
REM Slices: 20
REM RES _ NUM All-atoms Total-Side Main-Chain Non-polar All polar
REM ABS REL ABS REL ABS REL ABS REL ABS REL
RES ASP L 1 131.82 93.8 74.56 75.5 57.26 137.1 27.11 55.0 104.70 114.8
RES ILE L 2 12.14 6.9 9.85 7.4 2.29 5.6 9.85 7.1 2.29 6.3
RES VAL L 3 88.71 58.6 84.43 76.6 4.27 10.4 88.29 76.6 0.42 1.2
RES MET L 4 8.99 4.6 0.00 0.0 8.99 21.5 0.00 0.0 8.99 24.7
RES THR L 5 64.64 46.4 58.92 60.4 5.72 13.7 46.85 61.8 17.79 28.0
RES GLN L 6 7.63 4.3 0.08 0.1 7.56 18.1 0.14 0.3 7.49 5.9
RES SER L 7 82.71 71.0 53.96 73.7 28.75 66.3 47.02 96.9 35.69 52.5
RES GLN L 8 93.00 52.2 91.12 66.7 1.88 4.5 27.93 53.8 65.07 51.5
RES LYS L 9 113.28 56.6 105.29 66.4 7.99 19.1 93.78 80.9 19.50 23.1
RES PHE L 10 86.51 43.4 75.89 47.1 10.62 27.9 76.71 46.5 9.80 28.6
RES MET L 11 47.61 24.6 40.99 27.0 6.62 15.9 47.61 30.3 0.00 0.0
案例3：代码里面采用
def get_residue_exposure(input_pdb,result_dir):
result_fp = '%s/freesasa_out_tmp.tsv' % result_dir
result_fp_2 = '%s/freesasa_out.xlsx' % result_dir
cmd = """ freesasa %s --format=rsa --n-threads=5 --radii=naccess >%s """ % (input_pdb,result_fp)
os.system(cmd)
cmd_2 = """ sed -i 's/\s\+/\t/g' %s """ % result_fp
os.system(cmd_2)
df_1 = pd.read_csv(result_fp,sep='\t',header=None,skiprows=[0,1,2,3,4,5,6,7,8,],skipfooter=5,engine='python').drop(columns=[0,1]).fillna(0)
df_1.columns = ['chain','uid','All_atoms_abs','All_atoms_rel', 'Total_side_abs', 'Total_side_rel','Main_chain_abs','Main_chain_rel', 'Non_polar_abs', 'Non_polar_rel','All_polar_abs','All_polar_rel']
df_1.to_excel(result_fp_2,index=None)
df_2 = df_1[['chain','uid','All_atoms_rel']].rename(columns= {'All_atoms_rel':'Exposure_pct'})
参考资料
