阈值扫描矩阵(Csm):通过蛋白内部残基距离来进行结构分类和功能预测

文献:Cutoff Scanning Matrix (CSM): structural classification and function prediction by protein inter-residue distance patterns

1. 文章的主要实验流程

首先进行数据的清洗，之后进行CSM的生成，第三步使用SVD进行数据的降维。第四步使用机器学习的方法进行分类，最后进行模型的评价。

2. CSM计算方法

首先对于每个蛋白的数据集，作者生成一个特征向量(feature vector). 1). 计算所有的Cα原子对的欧几里得距离，并且定义了要考虑的距离范围和距离步长。作者对距离进行扫描，并且根据这些距离的阈值，计算残基的频率。伪代码如下:

1
2
3
4
5
6
7
8


function GENERATECSM(ProteinSet,CSM,Distance_min,Distance_max,Distance_step):
	for all protein i ∈(ProteinSet)do
		j=0
		Calculate the idstances between all pairs of Ca
		for dist <- Distance_min;to Distance_max;step Distance_step do
			CSM[i][j]<-Get frequency of pairs of Ca within a distance dist
			j++
	return CSM

矩阵的每一行代表了一个蛋白质，每一列代表一定距离内的频率。或者，该频率可以看作是蛋白质中的触点数量，用于某个截止距离或者使用该距离阈值定义的接触曲线的边缘计数。

用python实现一下，首先查看一个蛋白的情况,首先导入蛋白质

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14



from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Selection import unfold_entities

parser = PDBParser(PERMISSIVE=1)

structure_id = "2b3p"
filename = "pdb/2b3p.pdb"

structure = parser.get_structure(structure_id, filename)

model = structure[0]
chain = model["A"]
residue_list = unfold_entities(chain, "R")

去除异质残基:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12



residuedatas = []
#剔除异质残基
for residue in residue_list:
    residue_id = residue.get_id()
    # ('H_ CD', 801, ' ')
    hetfield = residue_id[0]
    # 确认其是否是异质残基
    if hetfield[0] != "H" and hetfield[0] != "W":
        if residue.has_id("CA"):
            # 获取残基编号
            residuedatas.append(residue['CA'])

进行原子配对:

1
2
3
4
5
6
7


import itertools

#参考 https://blog.csdn.net/sinat_28371057/article/details/114274832
atom_pairs = list(itertools.combinations(residuedatas, 2))
distances = []
for i in atom_pairs:
    distances.append(i[0] - i[1])

进行阈值的筛选，在网上查找了有两种方法一种是numpy，一种是原生方法，进行一下计算的速度比较:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


import time

time1 = time.time()
result = sum(i > 5 for i in distances)
print(result)
time2 = time.time()
# print(len(distances))

import numpy as np

j = np.array(distances)
print(sum(j > 5))
time3 = time.time()

print(time2-time1)
print(time3-time2)

返回的结果如下:

1
2
3
4


25563
25563
0.09600973129272461
0.0540010929107666

做了一下封装:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95



# 输入残基列表，返回原子实例
def getAtoms(residue_list, atom_name="CA"):
    '''
    :param residue_list: list;残基列表
    :param atom_name: string;获取的原子名，如CA，CB，O
    :return: list;原子实例
    '''
    residuedatas = []
    # 剔除异质残基
    for residue in residue_list:
        residue_id = residue.get_id()
        # ('H_ CD', 801, ' ')
        hetfield = residue_id[0]
        # 确认其是否是异质残基
        if hetfield[0] != "H" and hetfield[0] != "W":
            if residue.has_id(atom_name):
                # 获取残基编号
                residuedatas.append(residue[atom_name])
    return residuedatas


# 输入原子对,返回距离数组
def getDistance(atom_pairs):
    '''
    获得距离的数组
    :param atom_pairs:list;原子对
    :return: list;距离数组
    '''
    import itertools

    atom_pairs = list(itertools.combinations(atom_pairs, 2))
    distances = []
    for i in atom_pairs:
        distances.append(i[0] - i[1])
    return distances


# 统计距离
def sumCutoffDistance(distances, distance):
    '''
    获取距离对
    :param distances:list;距离数组
    :param distance:float;阈值距离
    :return:int;返回数目
    '''
    import numpy as np
    j = np.array(distances)
    return sum(j < distance)


def frequency(pdb):
    from Bio.PDB.PDBParser import PDBParser
    from Bio.PDB.Selection import unfold_entities

    parser = PDBParser(PERMISSIVE=1)

    structure_id = pdb
    filename = "pdb/" + pdb + ".pdb"
    structure = parser.get_structure(structure_id, filename)
    model = structure[0]
    chain = model["A"]
    residue_list = unfold_entities(chain, "R")
    atom_pairs = getAtoms(residue_list, "CA")
    distances = getDistance(atom_pairs)
    return sumCutoffDistance(distances, 50) / len(distances)


def frequences(pdb, distance_min, distance_max, distance_step):
    '''
    获取频率
    :param pdb:PDB_id
    :param distance_min:最小值(埃)
    :param distance_max: 最大值(埃)
    :param distance_step: 步长(埃)
    :return:频率数组
    '''
    from Bio.PDB.PDBParser import PDBParser
    from Bio.PDB.Selection import unfold_entities
    import numpy as np
    parser = PDBParser(PERMISSIVE=1)

    structure_id = pdb
    filename = "pdb/" + pdb + ".pdb"
    structure = parser.get_structure(structure_id, filename)
    model = structure[0]
    chain = model["A"]
    residue_list = unfold_entities(chain, "R")
    atom_pairs = getAtoms(residue_list, "CA")
    distances = getDistance(atom_pairs)
    totals = len(distances)
    frequences = []
    for i in np.arange(distance_min, distance_max + distance_step, distance_step):
        frequences.append(sumCutoffDistance(distances, i) / totals)
    return frequences

可以执行测试一下:

1

frequences("2b3p", 0.0, 30.0, 0.2)