基于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.