网上有许多蛋白质pI预测的工具，本文主要介绍一种提高预测准确性的方法，就以（鸡）溶菌酶为例，演示一下计算过程：
-
ProtParam
在UniPort中定位到氨基酸序列后，右侧下拉菜单中选择“ProtParam”或者“Compute pI/MW”，点击“go”跳转到计算界面，溶菌酶1-18aa为信号肽序列，计算pI时需要去除，预测的结果为9.32。 -
PIPC
从UniPort中获取氨基酸的fasta序列，将序列粘贴到PIPC文本框中（要去除信号肽部分的序列），可以得到计算结果为“Isoelectric point of protein is: 8.34”，需要注意的是PIPC会列出不同计算方法的参考值，一般取“IPC protein”的值（8.34）或者总平均值（8.7），但是溶菌酶分子很小应该去取“IPC peptide”的值（9.2）。
- 溶菌酶的真实等电点
这个很容易搜到，pI接近11具体为10.7，杨黄浩课题组发表的论文（Lin Z, Xia Z, Zheng J, et al. Synthesis of uniformly sized molecularly imprinted polymer-coated silica nanoparticles for selective recognition and enrichment of lysozyme[J]. Journal of Materials Chemistry, 2012, 22(34): 17914-17922.）提到了多种蛋白质的等电点，有兴趣可以多找几个进行试验，本文只拿溶菌酶为例。
| 名称 | 分子量(KD) | pI |
|---|---|---|
| BSA | 66.4 | 4.9 |
| 细胞色素C | 13 | 9.8 |
| 溶菌酶 | 13.4 | 10.7 |
| HRP | 44 | 7.2 |
| 肌红蛋白 | 16.7 | 6.9 |
- 校准预测值
两种预测方法得到的值都明显低于真实值。为何出现这样的偏差呢？之前的文章已经提到过，这些软件都是直接按输入的序列计算的，而在蛋白质三维构象中，并不是所有氨基酸都暴露出来的，但测定pI都是在溶液环境中进行的，对于那些溶剂不可及的氨基酸残基没有统计的必要性，如果全部计算在内就会出现拉底平均值的现象，往往预测出来的碱性蛋白质pI偏低，酸性蛋白质pI偏高。如果知道了酶的三维结构就可以计算出哪些残基是可及的哪些是不可及的，把那些不可及的氨基酸从酶的序列中“剔除”，可能会提高预测的准确性。
在“剔除”之前还应该考虑一个问题：酶的三维结构大多是用X-ray得到的，这与酶在溶液中的构象有多大差异？一定是有差异的，但一般差异不会很大，因为蛋白质的晶体都是在溶液中结晶出来的，在低温条件下，随着饱和度的增加晶核不断增长，这也要求蛋白质必然以离子状态聚集。Schwalbe等（Schwalbe H, Grimshaw S B, Spencer A, et al. A refined solution structure of hen lysozyme determined using residual dipolar coupling data[J]. Protein Science, 2001, 10(4): 677-688.）用NMR测定了溶菌酶的三维结构，通过与X-ray的数据对比，RMSD的值为1.49±0.10Å，两者的偏差并不大。这表明，溶菌酶晶体的三维结构可以用来计算溶剂可及性。但是，当酶与底物或者抑制剂等共结晶时，构象可能会有较大的变化。 - 分析蛋白质分子的溶剂可及性
PISA是EMBL开发的一个专门计算蛋白质晶体表面可及性和界面相互作用的工具，可以指定蛋白质的pdb编号或者自己上传pdb文件，远程服务器完成计算后返回“Interfaces”，“Monomers”，“Assemblies”三个结果，点击“Monomers”和“Assemblies”都可以得到可及性的结果。
结果中用颜色对可及性做了一个划分，对于每一行，全部为深色代表不可及残基，有一列为浅色该残基就可及，黄色为界面残基（单体结果中不会出现）。结果中很多列，我们只关注前四列：序号，多肽链编号和氨基酸名称，氨基酸序号，暴露在蛋白质分子表面的面积。最重要的是第四列，它的值代表该残基侧链暴露出来的面积，当值为0时或者很小时（这里将阈值选择为5）代表该残基被包埋到内部，属于“剔除”对象。溶菌酶可及性分析结果 - 结果处理
由于我的浏览器一直无法打开旁边的XML格式，就直接选中结果复制到文本文档中(各列之间是以制表符分开的)。这里只介绍两种处理方式：
- linux下shell代码
awk '{if($4>5) {print $2}}' lysozyme.txt | cut -d ':' -f 2
输出的结果：
LYS
VAL
PHE
......
LEU
打开：http://www.bio-soft.net/sms/找到“Protein Manipulation”的“三到一”选项将氨基酸三字符转换成单字母，当然你也可以用shell数组完成，这里就不再演示了。转化后序列：
KVFGRCEAAKRHGDNYRGYSNKFENFNQATNRNTDGDYQNRWWNDGRPGSRNLCNIPSALSSD
TANKKSDGNGNAVAWRNRKGTDVQAWIRGCRL
重新提交给PIPC，计算结果为IPC peptide (pI=10.293)，ProtParam计算结果为9.94，都已经比较接近真实值。
- python
# _*_ coding: utf-8 _*_
import re
import os
os.chdir('C:\\Users\\Administrator\\Desktop')
# 构建三到一转换的字典
aaDict = {'ASP':'D', 'GLU':'E', 'CYS':'C', 'TYR':'Y', 'HIS':'H',
'LYS':'K', 'ARG':'R', 'MET':'M', 'PHE':'F', 'LEU':'L',
'VAL':'V', 'ALA':'A', 'GLY':'G', 'GLN':'Q', 'ASN':'N',
'ILE':'I', 'TRP':'W', 'SER':'S', 'THR':'T', 'SEC':'U',
'PRO':'P', 'XAA':'X', 'SEC':'U', 'PYL':'O', 'ASX':'B',
'XLE':'J', }
aa_sequence = ''
with open('lysozyme.txt', 'r') as f:
str = f.readlines()
for line in str:
temp = line.split('\t')[1:3] # 切分列，获得第2-4列
aa = temp[0].strip()[2:5] # 得到氨基酸名称
value = float(temp[1].strip()) # 得到ASA值
if value > 5.0:
aa_sequence += aaDict[aa]
print(aa_sequence)
你也可以通过excel，R或者其他方法来处理，只要把那些ASA≤5 Å2的氨基酸去除就可以了。
