|
基于UNIPROT ID对蛋白的结构和序列进行差异分析
### 基于UNIPROT ID对蛋白的结构和序列进行差异分析
#### 摘要
随着晶体结构解析技术的发展,同一个UniProt的蛋白越来越多的结构被解析出来。
本文将教你如何巧妙适当的这些出这些结构的差异。
#### 案例 human BACE1 蛋白
在PDB 数据库中很容易就找到对应的UniProt 编号
Find proteins for P56817 (Homo sapiens)
截止2018-10-27号,一共有373个PDB 编号
http://www.rcsb.org/pdb/results/results.do?tabtoshow=Current&qrid=D8C70A10
#### 进入UniProt 数据库进行查看
Go to UniProtKB: P56817
在Sequence 中可以看到该蛋白一共有6种亚型(isoforms)的序列,
```
Sequence statusi: Complete.
Sequence processingi: The displayed sequence is further processed into a mature form.
This entry describes 6 isoformsi produced by alternative splicing.
This entry has 6 described isoforms and 3 potential isoforms that are computationally mapped
```
下载这6种亚型,并进行alignment (网站上自带align 按钮);
https://www.uniprot.org/align/A20181027F725F458AC8690F874DD868E4ED79B88013073E.aln
最长的序列是第一个序列,一共有501个氨基酸残基,其中差异主要发生在1-215个前半部分,后半部分这6个亚型序列都一样。
从进化树来看 1234是一类 56是一类;
相同的位置有317个*;相似的位置有4个点.:号。
#### 对373个蛋白质结构的序列进行分类
首先找一个模板文件,要求是序列比较长,质量比较好,没有插入序列等PyMOL 显示错乱情况,
1FKN PyMOL 显示错乱,因为里面有插入残基编码。
http://www.rcsb.org/pdb/protein/P56817
选择第一个代码4joo 作为模板分析,发现其和第一个序列亚型是一致的。
那么我们现在需要矫正4joo的氨基酸编号。
其GSF编号是从45开始,对应第一个亚型的58号GSF-des453.
取45GSF-DES440区域基于alginment 调整编号。
45+13
1fkn 需要删除或者调整一些杂乱的残基,否则PyMOL显示会错乱,这里我直接删除
1fkn从EMVDL开始编号是1到GYN编号是385 对应第一个亚型的62-446
+61
PyMOL 命令
alter (all),resi=str(int(resi)+61);
sort;
必须使用sort命令更新才会产生结果。
1m4h 从EMVDL 到GYN
1tqf 序列上存在2个突变K136A,E138A,结构上还存在一段缺失372DVATS376.
372-376缺失。
从EMVDL 到GYN
1W51 序列上存在2个突变R56K R57K
而结构是从60FVEMVDL开始,缺失219GFPLNQSEVL228序列;
alter (all),resi=str(int(resi)+61);
1XN2 序列上和第一个序列亚型是一致的,
alter (all),resi=str(int(resi)+61);
从EMVDL 开始的
1XN3 同 1XN2
1XS7序列上和第一个序列亚型是一致的,
alter (all),resi=str(int(resi)+61);
alter (all),resi=str(int(resi)+61);
从EMVDL 开始的
1YM2 序列上和第一个序列亚型是一致的,结构上缺少一段219GFPLNQSEVLAS230 的loop
alter (all),resi=str(int(resi)+61)
1YM4 序列上和第一个序列亚型是一致的,结构上缺少一段219GFPLNQSEVL228 的loop
alter (all),resi=str(int(resi)+61)
2b8L 序列上有2个突变 K136A E138A 同1tqf
结构上还存在一段缺失,219GFPLNQSEVLAS230
2B8V.pdb 同2K8L 有2个突变 K136A E138A
结构上还存在一段缺失,220FPLNQSEVLASVG232,372DVATS376
2F3E 缺少一段loop 219GFPLNQSEVLA229
alter (all),resi=str(int(resi)+61)
2F3F 序列一致
### 操作总结
上述的操作总结起来,无非就是删除插入残基,并将序号加上61
### 批量处理
1. 删除插入残基, 根据插入残基代码删除插入残基,文本处理 ,HHC 给我的PDB文件中已经处理好了。
2. 删除HTEATOM,python 文本处理即可
3. 保存fasta,pymol 批量处理,save 1.fasta,obj
4. 序列比对,schrodinger 中Multiple Sequence Viewer
导入需要比对的Sequece fasta文件,建议一次比对不超过50-100个fasta 文件
点击 Multiple alignment,由于序列大部分是相同的,因此我只需要对不同的进行着色,
根据具体的情况是否需要换行。
就可以看出哪些氨基酸进行突变,那块缺失了。
但是结果不可以导出。
怎么导出结果?
需要关注的2个方面
-- 是否是+61, 人眼观察结果就可以了。
如果序列能对上,并且E的编号是1则是+61
有些E的编号是1则加61
有些E的编号是7【4DI2】则加53
有些没有E是M
-- 保守性的序列的编号,由于不可以导出结果,那就人工记录下存在突变的氨基酸编号
K136A E138A N153Q
T133C Q134K
loop 缺失区域的序列比对往往会存在错位现象。GAP区域的突变往往是错位。
query窗口类似于table窗口,可以进行选择并删除序列。
4. 加入61
借助于脚本step4_renameResi.py 完成
尽管有插入序列,但是PDB库中氨基酸的编号还是统一的,EMVDL 都是从1开始的
5. 检测序列差异,并check 结果是否完整
尽管有插入序列,但是PDB库中氨基酸的编号还是统一的,EMVDL 都是从1开始的
发现里面除了有A chain 还有B chain;
只保留一个chain
6. 保留A 链
发现依旧有问题 ,存在多构象问题
7. 仅仅保留alt A+' '构象
再次check,
发现有4个体系是空的,经过经常发现其没有A链
手动处理 1XN3 1XS7 3LPJ 3LPK
8. 得到result.txt 至此序列分析结束。
EMVDL-GYN446结束,
9. 计算最大alpha C 距离
62到446 最大alphaC距离
10. 选择典型坐标,构建自己的坐标体系,
在每个蛋白上指定2个原子,计算其距离。
距离1: ASP-93 是稳定存在的有354个CA 134GLN [('GLN', 340)] 有340个
距离2: ASN-66是稳定存在的 [('ASN', 354)] 226GLU的距离 [('GLU', 126), ('ILE', 3), ('GLY', 3), ('MET', 1)]
距离3: 329 [('VAL', 345), ('ILE', 2), ('PRO', 2), ('PHE', 1), ('TRP', 1), ('ASN', 1)] THR 375 [('THR', 203), ('SER', 4), ('GLN', 1), ('ALA', 1)]
对应3个loop,
132-138 131-139 loop1
217-231 216-232
370-380 369-381
11. 计算3个距离
pdbid dist1,dist2,dist3
12. 取最小的和最大的距离的2个构象,一共取6个构象
4FRS_nohet_number_A_altA.pdb 5I3X_nohet_number_A_altA.pdb
3VEU_nohet_number_A_altA.pdb 4FM8_nohet_number_A_altA.pdb
4D89_nohet_number_A_altA.pdb 3VG1_nohet_number_A_altA.pdb
# 发现4D85 的编号有问题
check 的时候需要检查最后一个GYN的编号,
总结:
这2个很像 1m4h 1fkn 1XN2 1XS7 1YM2
1XN3
1W51 2b8l 2b8v 2F3E
1tqf
#### 工具
1. https://www.ebi.ac.uk/Tools/msa/clustalo/ 序列比对
2.
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-20 01:53
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社