基于人工智能的多肽药物分析问题(五)
2021SC@SDUSC
1. 前言
前面代码已经完成了1T+蛋白质数据的预处理,之后准备构建模型
部分术语:
残基: 在蛋白质的序列中,氨基酸之间的氨基和羧基脱水成键,氨基酸由于其部分基团参与了肽键的形成,剩余的结构部分则称氨基酸残基。
2. 代码分析
2.1 计算原子之间距离图
maps = extract_multi_distance_map(pose)
在将蛋白质数据预处理之后的特征值都存储到pdict变量中,之后,程序中下一步调用了extract_multi_distance_map方法,输入为pose,也就是最开始输入的蛋白质数据,返回的是“不同原子的距离图”
def extract_multi_distance_map(pose):
x1 = get_distmaps(pose, atom1="CB", atom2="CB", default="CA")
x2 = get_distmaps(pose, atom1=dict_3LAA_to_tip, atom2=dict_3LAA_to_tip)
x3 = get_distmaps(pose, atom1="CA", atom2=dict_3LAA_to_tip)
x4 = get_distmaps(pose, atom1=dict_3LAA_to_tip, atom2="CA")
output = np.stack([x1,x2,x3,x4], axis=-1)
return output
如代码所示,x1表示CB到CB的距离map,x2表示的是CA 到CA 的距离图,x3表示的是tip与tip原子之间的距离图,x4表示的是tip与CA的距离图,四种距离图分别保存在x1,x2,x3,x4四个变量中,然后调用np.stack将四个变量合成一个更高维的变量。并返回给最上方的maps变量。
2.1.1 get_distmaps
def get_distmaps(pose, atom1="CA", atom2="CA", default="CA"):
psize = pose.size()
xyz1 = np.zeros((psize, 3))
xyz2 = np.zeros((psize, 3))
for i in range(1, psize+1):
r = pose.residue(i)
if type(atom1) == str:
if r.has(atom1):
xyz1[i-1, :] = np.array(r.xyz(atom1))
else:
xyz1[i-1, :] = np.array(r.xyz(default))
else:
xyz1[i-1, :] = np.array(r.xyz(atom1.get(r.name(), default)))
if type(atom2) == str:
if r.has(atom2):
xyz2[i-1, :] = np.array(r.xyz(atom2))
else:
xyz2[i-1, :] = np.array(r.xyz(default))
else:
xyz2[i-1, :] = np.array(r.xyz(atom2.get(r.name(), default)))
return distance_matrix(xyz1, xyz2)
get_distmaps 的目的是根据传入的pose参数(pose类型在前面的文章有讲过),得到不同原子之间的距离矩阵。
思路解析:
首先获取pose的大小,初始化xyz1,xyz2两个空变量用于保存原子坐标,接着进行一个循环,利用p.residue(i)
方法查看pose实例的第i个残基(返回值是三个大写字母组成的字符串,如ARG),接着判断传入的参数是否为字符串类型,这是因为在extract_multi_distance_map
方法中,传入的参数既有字符串类型,又有字典map类型,其中dict_3LAA_to_tip
是一个转化关系,如下所示:
dict_3LAA_to_tip = {"ALA":"CB", "CYS":"SG", "ASP":"CG", "ASN":"CG", "GLU":"CD",
"GLN":"CD", "PHE":"CZ", "HIS":"NE2", "ILE":"CD1", "GLY":"CA",
"LEU":"CG", "MET":"SD", "ARG":"CZ", "LYS":"NZ", "PRO":"CG",
"VAL":"CB", "TYR":"OH", "TRP":"CH2", "SER":"OG", "THR":"OG1"}
如果为str类型,接着判断残基中是否有方法所输入的原子(也就是输入的字符串),有的话,则将该残基中所有的该原子的坐标数组保存至xyz1,如果atom1参数所表示的原子不存在,则利用默认的"CA"计算。
如果不是str类型,那么计算残基的name根据dict_3LAA_to_tip转化成的原子和default的"CA"之间的坐标
同样,接下来对atom2以上述流程进行计算,最后返回距离矩阵
2.2 计算能量矩阵
_2df, aas = extract_EnergyDistM(pose, energy_terms)
使用extract_EnergyDistM计算能量矩阵以及aa序列
传入的energy是能量
energy_terms = [pyrosetta.rosetta.core.scoring.ScoreType.fa_atr,\
pyrosetta.rosetta.core.scoring.ScoreType.fa_rep,\
pyrosetta.rosetta.core.scoring.ScoreType.fa_sol,\
pyrosetta.rosetta.core.scoring.ScoreType.lk_ball_wtd,\
pyrosetta.rosetta.core.scoring.ScoreType.fa_elec,\
pyrosetta.rosetta.core.scoring.ScoreType.hbond_bb_sc,\
pyrosetta.rosetta.core.scoring.ScoreType.hbond_sc]
2.2.1 extract_EnergyDistM
该方法的思路为,获取残基总数,并利用Pose类内置的api获得能量图
length = int(pose.total_residue())
tensor = np.zeros((1+len(energy_terms)+2, length, length))
energies = pose.energies()
graph = energies.energy_graph()
用graph变量保存
用能量填充矩阵:
aas = []
for i in range(length):
index1 = i + 1
aas.append(pose.residue(index1).name().split(":")[0].split("_")[0])
iru = graph.get_node(index1).const_edge_list_begin()
irue = graph.get_node(index1).const_edge_list_end()
while iru!=irue:
edge = iru.__mul__()
evals = [edge[e] for e in energy_terms]
index2 = edge.get_other_ind(index1)
count = 1
for k in range(len(evals)):
e = evals[k]
t = energy_terms[k]
if t == pyrosetta.rosetta.core.scoring.ScoreType.hbond_bb_sc or t == pyrosetta.rosetta.core.scoring.ScoreType.hbond_sc:
if e != 0.0:
tensor[count, index1-1, index2-1] = 1
else:
tensor[count, index1-1, index2-1] = e
count += 1
iru.plus_plus()
进行残基总数次的循环,每次循环中,将残基相关信息添加至aas变量,即该方法返回的aa序列,之后获取edge迭代器,接着对能量图进行解析,如果该edge边的起点和终点不是同一个节点,则获取该边,然后获取边的能量值保存在evals
变量中,index2是index1节点通过edge这条边对应的另外一个节点的索引
接着,再进行内层的循环,这里只关注hbond_bb_sc和hbond_sc两种能量类型,若t为这两种类型之一,那么对应的tensor三维张量对应的位置存储为1,若不是,则按照原能量e进行存储至相应的位置。最后移除已经计算过的节点。
再往后,对能量项目 energy term进行简单的转化
for i in range(1, 1+len(evals)):
temp = tensor[i]
if i == 1 or i == 2:
tensor[i] = np.arcsinh(np.abs(temp))/3.0
elif i == 3 or i==4 or i==5:
tensor[i] = np.tanh(temp)
利用CB的距离矩阵计算残基之间的distance
for i in range(length):
for j in range(length):
index1 = i + 1
index2 = j + 1
# Calculate distance and store. Use CA if CB does not exist.
if pose.residue(index1).has("CB"):
vector1 = pose.residue(index1).xyz("CB")
else:
vector1 = pose.residue(index1).xyz("CA")
if pose.residue(index2).has("CB"):
vector2 = pose.residue(index2).xyz("CB")
else:
vector2 = pose.residue(index2).xyz("CA")
distance = vector1.distance(vector2)
tensor[0, index1-1, index2-1] = distance
遍历每一条边的两个节点计算残基之间的距离,如果“CB”存在,那么利用"CB"进行计算,如果“CB”不存在,那么就利用“CA”进行距离的计算,保存在tensor的第一个大维度中
2.2.2 获取hbond信息
hbonds = get_hbonds(pose)
for hb in hbonds[0]:
index1 = hb[0]
index2 = hb[1]
tensor[count, index1-1, index2-1] = 1
count +=1
for hb in hbonds[1]:
index1 = hb[0]
index2 = hb[1]
tensor[count, index1-1, index2-1] = 1
这里使用get_hbonds方法得到backbone-to-backbone的氢键,输入参数为pose,输出结果为两个数组,接着,分别对得到的数据进行节点的遍历赋值,并保存至tensor中
下面是get_hbonds()方法的解析:
def get_hbonds(pose):
hb_srbb = []
hb_lrbb = []
hbond_set = pose.energies().data().get(pyrosetta.rosetta.core.scoring.EnergiesCacheableDataType.HBOND_SET)
for i in range(1, hbond_set.nhbonds()):
hb = hbond_set.hbond(i)
if hb:
acceptor = hb.acc_res()
donor = hb.don_res()
wtype = pyrosetta.rosetta.core.scoring.hbonds.get_hbond_weight_type(hb.eval_type())
energy = hb.energy()
is_acc_bb = hb.acc_atm_is_protein_backbone()
is_don_bb = hb.don_hatm_is_protein_backbone()
if is_acc_bb and is_don_bb:
if wtype == pyrosetta.rosetta.core.scoring.hbonds.hbw_SR_BB:
hb_srbb.append((acceptor, donor, energy))
elif wtype == pyrosetta.rosetta.core.scoring.hbonds.hbw_LR_BB:
hb_lrbb.append((acceptor, donor, energy))
return hb_srbb, hb_lrbb
首先初始化变量,之后利用pyrosetta的api获得pose蛋白质的能量种类,其中HBOND_SET是EnergiesCacheableDataType这个类的枚举值,是能量数据类型的枚举值,得到氢键集合hbond_set
随后遍历该氢键集合,获取相关的氢键数据并保存
最后进行判断,判断氢键的weight类型,这里判断了两个类型,分别是hbw_SR_BB和hbw_LR_BB,最后都存储到最开始初始化的两个变量中【但是这里面几个api我查了很久都没有查到有效的解释文档,只查到了输入和输出】
2.3 提取OneBodyTerms
_1df, _ = extractOneBodyTerms(pose)
利用pose信息获取键角和键长,保存至相关变量,在extractOneBodyTerms方法中,首先调用get_feature_matrix方法获取键角和键长:
get_feature_matrix
def get_feature_matrix(mypose, padval=0):
# GIVEN: a pose and a value for when lenghts and angles don't make
# sense at the C and N terminus
# RETURNS: a 2d numpy array
# rows correspond to residues and columns correspond to features
result = []
column_names = ["NcCAc_len", "CAcCc_len", "CcNn_len", "CpNcCAc", "NcCAcCc", "CAcCcNn"]
for res_pos in range(1,len(mypose.sequence())+1):
feature_dict = get_bond_lengths_and_angles(mypose,res_pos)
data_row = []
# "zero padding"
if res_pos == 1:
feature_dict["CpNcCAc"] = padval
if res_pos == len(mypose.sequence()):
feature_dict["CcNn_len"] = padval
feature_dict["CAcCcNn"] = padval
for feature in column_names:
data_row.append(feature_dict[feature])
result.append(data_row)
return np.array(result).T
当长度和角度在C和N端点没有实际的意义的时候的参数pose和padval,返回值是一个二维数组,其中行对应的是残基,列对应的是相应的特征features,初始化了6种features的类型名称
之后遍历pose,调用get_bond_lengths_and_angles
方法,该方法接受的参数为pose对象以及残基的索引,返回值为一个字典类型,该字典类型的key为features,value为对应的值。
之后,对于第一个和最后一个残基,将feature_dict的相应键赋值为0,随后存储到变量中并返回
这里调用的get_bond_lengths_and_angles方法,考虑了三种长度的氢键以及角度,分别是
键长
N(k)-CA(k)
CA(k)-C(k)
Ca(k)-N(K+1) (except for C-term)
键角
C(k-1)-N(k)-CA(k) (except for N-term)
N(k)-CA(k)-C(k)
CA(k)-C(k)-N(k+1) (except for C-term)
该方法具体分析下周再仔细研究一下