def detransform_structure(structure_to_align, Q):
'''Применить обратные повороты к моделям в структуре'''
for i,model in enumerate(structure_to_align):
transform_model(model, np.transpose(Q[i]), [0,0,0])
def transform_structure(structure_to_align, Q):
'''Применить повороты к моделям в структуре'''
for i,model in enumerate(structure_to_align):
transform_model(model, Q[i], [0,0,0])
def get_coord_list(model):
'''Получить лист координат атомов модели'''
return [atom.get_coord() for atom in model.get_atoms()]
def transform_model(model, rot, tran):
'''Повернуть/перенести атомы модели по данным матрице поворота rot и вектору пар. переноса tran'''
for atom in model.get_atoms():
atom.transform(rot, tran)
def random_rotate_structure(structure):
'''Случайным образом повернуть все конформации в структуре'''
for model in structure:
Q = get_random_rot_mat_haar(3)
transform_model(model, Q, [0,0,0])
def random_translate_atoms_in_structure(structure, eps):
'''Случайным образом сдвинуть все атомы во всех моделях в структуре'''
#print(len(list(structure.get_atoms())), len(list(structure[0].get_atoms())))
for atom in structure.get_atoms():
atom.set_coord(atom.get_coord()+eps*(np.random.rand(3)-1))
def RMSD_pairwise_avg(structure):
'''Вернуть среднее RMSD по всем парам в структуре'''
#Находим число моделей в структуре
N = len(list(structure))
RMSD = 0.0
for ref_model in structure:
#Получаем координаты атомов первой модели.
ref_atoms = np.asarray(get_coord_list(ref_model))
for alt_model in structure:
#Получаем координаты атомов второй модели.
alt_atoms = np.asarray(get_coord_list(alt_model))
#Находим RMSD для пары моделей с помощью функции Bio.SVDSuperimposer.SVDSuperimposer._rms(...)
sup = Bio.SVDSuperimposer.SVDSuperimposer()
RMSD += sup._rms(ref_atoms, alt_atoms)
#Усредняем RMSD полученное по всем парам моделей
RMSD = RMSD / (N * N)
return RMSD
def get_pair_rot_tran(ref_model, alt_model):
'''Вернуть поворот/трансляцию совмещающие (наилучшим образом) пару моделей (по алгоритму Каббша)'''
#Получаем координаты атомов моделей
ref_atoms = np.asarray(get_coord_list(ref_model))
alt_atoms = np.asarray(get_coord_list(alt_model))
#Создаем суперимпозер
sup = Bio.SVDSuperimposer.SVDSuperimposer()
#Инициализируем наш суперимпозер двумя наборами данных
sup.set(ref_atoms, alt_atoms)
#Находим выравнивание
sup.run()
return (sup.get_rotran()[0], sup.get_rotran()[1])
def translate_to_center(structure):
"""Выполнить трансляции (параллельные переносы) всех моделей в структуре к центру координат"""
#Количество атомов в модели структуры
n = len(list(structure[0].get_atoms()))
for model in structure:
center = np.mean([np.array(atom.get_coord()) for atom in model.get_atoms()], axis = 0)
transform_model(model, np.identity(3), -center)
def get_pair_covar_matr(model1, model2):
'''Получить матрицу ковариаций для двух моделей'''
model1_X = np.asarray(get_coord_list(model1))
model2_X = np.asarray(get_coord_list(model2))
R = (np.transpose(model1_X)).dot(model2_X)
return R
def get_ort_mat_SVD(R):
'''Получить матрицу поворота из матрицы R как результат произведения матриц SVD разложения R'''
U, s, V = np.linalg.svd(R)
E_alter = np.array([1.0, 1.0, np.linalg.det(np.dot(U, V))])
Q_ans = np.dot(U, np.dot(np.diag(E_alter), V))
return np.transpose(Q_ans)
Наша задача ставится так:
$$\min_{Q_1,\ldots,Q_N\in SO(3)}\frac{1}{2}\sum\limits_{i,j=1}^N\sum\limits_{k=1}^n\|Q_ir_i^k-Q_jr_j^k\|^2$$Здесь $Q_i$ - искомые матрицы поворотов, а $r_i^k$ - вектор $k$-го элемента $i$-й структуры.
Считаем что в структурах $n$ элементов, а всего структур - $N$.
Обозначим $X_i\in\mathbb{R}^{3\times n}$ - матрицу координат атомов $i-й$ структуры. Тогда имеющийся функционал перепишется в следующем виде:
$$\min_{Q_1,\ldots,Q_N\in SO(3)}\frac{1}{2}\sum\limits_{i,j=1}^N\|Q_iX_i^k-Q_jX_j^k\|^2_F$$Раскроем теперь эти произведения:
$$\|Q_iX_i^k-Q_jX_j^k\|^2_F=Tr((Q_iX_i^k-Q_jX_j^k,Q_iX_i^k-Q_jX_j^k))=Tr((Q_iX_i^k,Q_iX_i^k))-2Tr((Q_iX_i^k,Q_jX_j^k))+Tr((Q_jX_j^k,Q_jX_j^k))$$Пользуясь тем, что $Tr(Q_iX_i^k,Q_iX_i^k)=Tr((X_i^k,Q_i^TQ_iX_i^k))=Tr((X_i^k,X_i^k))$отбрасываем первое и третье слагаемое как константы, и приходим к следующей задаче максимизации:
$$\max_{Q_1,\ldots,Q_N\in SO(3)}\sum\limits_{i,j=1}^N Tr(X_i^TQ_i^TQ_jX_j) = \max_{Q_1,\ldots,Q_N\in SO(3)}\sum\limits_{i,j=1}^N Tr(Q_i^TQ_jX_jX_i^T)$$Тогда вводя матрицы $R = (Q_1, Q_2, \ldots, Q_N)$ и $C_{ji} = \{X_jX_i^T\}_{ji}$ получим:
$$\max_{Q_1,\ldots,Q_N\in SO(3)} Tr(R^TRC) = \max_{Q_1,\ldots,Q_N\in SO(3)} Tr(GC)$$где $G = R^TR$ - положительно полуопределенная матрица. Данная задача и ее релаксации разобраны в https://web.math.princeton.edu/~amits/publications/sync_ima.pdf Следуя тем же путем, рассматриваем следующую релаксацию:
\begin{align} &\max_{Q_1,\ldots,Q_N\in SO(3)} ~Tr(GC)\\ &\text{subject to} ~G\succeq 0\\ & ~G_{ii} = I_3,~for~i = 1,\ldots,N \end{align}После нахождения решения данной SDP задачи остается один этап - восстановить матрицы поворота из полученной матрицы G. Для этого можно воспользоваться одним из трех способов:
1)Можно положить первую модель неподвижной (т.е. положить $Q_1=I_3$). Тогда, в виду того что первых три столбца матрицы G это $R^TQ_1$ получаем, что искомый вектор повротов $R$ получается транспонированием первых трех столбцов матрицы G.
2)Также мы можем восстановливать матрицы поворота из трех наибольших собственных векторов матрицы G. Из данных трех собственных векторов получается набор матриц $R_i$,из которых $Q_i$ восстанавливаются известным способом через $SVD$ разложение.
3)Еще один способ заключается в домножении матрицы G на вектор случайных матриц поворота $F = (F_1, \ldots, F_N)$, после чего восстановить искомые матрицы поворота из матрицы $FG$ тем же способом что и в пункте 2 (через SVD). Метод генерации данных случайных матриц взят из статьи https://arxiv.org/pdf/math-ph/0609050.pdf
def Wang_Singer_SDP_get_C_matrix(structure, N):
'''Получить матрицу C попарных ковариаций координат атомов моделей'''
C = np.zeros((3*N,3*N))
submatr_list = [0,1,2]
for i,ref_model in enumerate(structure):
for j, alt_model in enumerate(structure):
C[np.ix_([x + 3*i for x in submatr_list],
[x + 3*j for x in submatr_list])] = (get_pair_covar_matr(ref_model,alt_model))
#(rot, tran) = get_pair_rot_tran(ref_model, alt_model)
#C[np.ix_([x + 3*i for x in submatr_list],
# [x + 3*j for x in submatr_list])] = np.transpose(np.asarray(rot))
return C;
from cvxpy import *
def Wang_Singer_SDP_cvxproblem_solve(C, N):
'''Найти решение SDP (с помощью cvxpy)'''
submatr_list = [0,1,2]
#import cvxpy as cvx
#G - положительно полуопределенная матрица, G=R^T*R
G = Semidef(3*N,3*N)
#Условия оптимизации:
sdp_constraints = []
for i in range(0,N):
#Блоки на диагонали - это единичные матрицы:
sdp_constraints.append(G[np.ix_([x+3*i for x in submatr_list],
[x+3*i for x in submatr_list])]==np.identity(3))
#sdp_constraints.append(G>>np.identity(3*N)*(10**(-5)))
#Максимизируем функционал
sdp_objective = Maximize(trace(G*C))
#Постановка оптимизационной задачи
sdp_problem = Problem(sdp_objective, sdp_constraints)
#Находим решение поставленной задачи
#print(G.value)
sdp_problem.solve(solver = MOSEK)
#print(G.value)
return G.value
def change_columns_order_in_matr(R, order):
'''Поменять порядок столбцов в переданном списке матриц, так,
чтобы первым столбцом во всех матрицах списка шел столбец fir, вторым - sec, третьим - thir'''
for (i,R_i) in enumerate(R):
R_new = R_i
R_new[:,[0,1,2]] = R_i[:, order]
R[i]=R_new
return R
def Wang_Singer_SDP_get_rot_matrices_from_eigenvectors(structure_to_align, G, N):
'''Восстановить матрицы поворота из трех наибольших собственных векторов матрицы G'''
import scipy
print("eigenvectors:")
print(scipy.linalg.eigh(G,eigvals=[3*N-6, 3*N-1])[0])
#Берем три наибольших собственных вектора (упорядочены от наименьшего к наибольшему)
G_top3_eigenvectors = scipy.linalg.eigh(G,eigvals=[3*N-3, 3*N-1])
#Получаем набор матриц (еще не поворота)
R = np.split(G_top3_eigenvectors[1], N)
RMSD_min = RMSD_pairwise_avg(structure_to_align)
Q_ans = [np.identity(3) for R_i in R]
#Ищем перестановку столбцов, наилучшим образом совмещающую модели
import itertools
for order in itertools.permutations([0, 1, 2]):
#Создаем временную структуру
structure_temp = structure_to_align.copy()
#Берем собственные вектора в новом порядке
R_temp = change_columns_order_in_matr(R, order)
#Восстанавливаем матрицы поворота из собственных векторов
Q_ans_temp = [np.transpose(get_ort_mat_SVD(R_i)) for R_i in R_temp]
#Применяем полученные повороты к моделям во временной структуре
for i,model in enumerate(structure_temp):
transform_model(model, Q_ans_temp[i], [0,0,0])
#Находим среднее попарное RMSD временной структуры после выравнивания
RMSD_min_new = RMSD_pairwise_avg(structure_temp)
#Если оно меньше - берем это в качестве нового набора матриц поврота
#print(RMSD_min_new, RMSD_min, type(RMSD_min), type(RMSD_min_new))
if RMSD_min_new <= RMSD_min :
RMSD_min = RMSD_min_new
Q_ans = Q_ans_temp
return (Q_ans)
def Wang_Singer_SDP_get_rot_matrices_from_first_columns(G, N):
'''Восстановить матрицы поворота из трех первых столбцов матрицы G'''
print("columns:")
#Берем три первых столбца матрицы G
G_first3_columns = np.asarray(G)[:, [0, 1, 2]]
#Получаем набор матриц (еще не поворота)
R = np.split(G_first3_columns, N)
Q_ans = [np.transpose(get_ort_mat_SVD(R_i)) for R_i in R]
return(Q_ans)
def get_random_rot_mat_haar(n):
from scipy import random, linalg, dot, diagonal, absolute, multiply
'''Получить случайную ортогональную матрицу'''
'''Взято из https://arxiv.org/pdf/math-ph/0609050.pdf раздел 4. Беру только лишь реальную часть.'''
'''A Random matrix distributed with Haar measure'''
z = random.randn(n,n)
q,r = linalg.qr(z)
d = diagonal(r)
ph = d/absolute(d)
q = multiply(q,ph,q)
E_alter = np.array([1.0, 1.0, np.linalg.det(q)])
q = np.dot(np.asarray(q), np.diag(E_alter))
return q
def get_random_unitary_matrices_vector(N):
'''Получить вектор из случайных ортогональных матриц'''
Q_ans = [get_random_rot_mat_haar(3) for i in range(0,N)]
return np.reshape(np.asarray(Q_ans), (3*N, 3))
def Wang_Singer_SDP_get_rot_matrices_stochastically(G, N):
print("stochas:")
'''Восстановить матрицы поворота умножением G на вектор случайных матриц'''
Q = get_random_unitary_matrices_vector(N)
G_min_eig_value = min(list(np.linalg.eig(G)[0]))
'''ЗДЕСЬ ПОДГОН. ПОЧЕМУ? G - ДОЛЖНА(!) БЫТЬ ПОЛ. ОПР. ПОСЛЕ РЕШЕНИЯ SDP.
НО ЭТО ИНОГДА НЕ ТАК (ХОТЯ И С МАЛОЙ ПОГРЕШНОСТЬЮ). А это нужно для разложения Холецкого'''
if(G_min_eig_value <= 10**(-8)):
L = np.linalg.cholesky(np.asarray(G)-2*np.identity(3*N)*G_min_eig_value)
else:
L = np.linalg.cholesky(np.asarray(G))
LQ_list = np.split(np.dot(L, Q), N)
Q_ans = [np.transpose(get_ort_mat_SVD(LQ_i)) for LQ_i in LQ_list]
return(Q_ans)
def Wang_Singer_SDP_get_rot_matrices(structure_to_align, G, N, mode):
'''Восстановить матрицы поворота одним из трех способов:
mode = 0 - из трех наибольших собственных векторов матрицы G
mode = 1 - из трех первых столбцов матрицы G
mode = 2 - умножением вектора случайных матриц на G'''
if(mode == 0):
return Wang_Singer_SDP_get_rot_matrices_from_eigenvectors(structure_to_align, G, N)
elif(mode == 1):
return Wang_Singer_SDP_get_rot_matrices_from_first_columns(G, N)
elif(mode == 2):
return Wang_Singer_SDP_get_rot_matrices_stochastically(G, N)
'''print(mode)
return {
0: Wang_Singer_SDP_get_rot_matrices_from_eigenvectors(structure_to_align, G, N),
1: Wang_Singer_SDP_get_rot_matrices_from_first_columns(G, N),
2: Wang_Singer_SDP_get_rot_matrices_stochastically(G, N)
}[mode]'''
def Wang_Singer_SDP(structure_to_align, rot_matr_recov_mode, G_is_found, G_ext):
'''Получить искомое выравнивание на основе алгоритма из Wang,Singer'''
#Транслируем все структуры к центру координат
translate_to_center(structure_to_align)
#Находим число моделей в структуре
N = len(list(structure_to_align))
if(G_is_found == 0):
#C - матрица попарных ковариаций.
C = Wang_Singer_SDP_get_C_matrix(structure_to_align, N)
#G - положительно полуопределенная матрица, которая теоретически имеет вид R^T*R, где R - столбец из искомых
#матриц поворота.
G = Wang_Singer_SDP_cvxproblem_solve(C, N)
else:
G = G_ext
#Получаем матрицы поворота одним из методов
#mode = 0 - из трех наибольших собственных векторов G
#mode = 1 - из трех первых столбцов матрицы G
#mode = 2 - умножением вектора случайных матриц на G
Q_ans = Wang_Singer_SDP_get_rot_matrices(structure_to_align, G, N, mode = rot_matr_recov_mode)
return (Q_ans, G)
Ниже - алгоритм Каббша попарного выравнивания (2 реализации, встроенная в Bio и моя. Работают одинаково, как и ожидалось)
Встроенная:
#Выравнивание по алгоритму Каббша. Используется встроенный метод bio
def Kabsch(structure):
'''Получить искомое выравнивание на основе алгоритма Каббша через попарное выравнивание'''
#Мы реализуем попарное сравнение, так что просто сравниваем все модели с первой, выбранной статичной.
ref_model = structure[0]
Q_ans = []
for alt_model in structure:
#Получаем попарный поворот/сдвиг
(rot,tran) = get_pair_rot_tran(ref_model,alt_model)
#Применяем их, преобразуя очередную модель
transform_model(alt_model, np.identity(3), tran)
Q_ans.append(rot)
return Q_ans
Моя реализация алгоритма Каббша:
def MyKabsch_get_pair_rot(ref_model, alt_model):
'''Получить матрицу поворота, осуществляющего наложение alt_model на ref_model'''
#Получаем матрицу ковариаций
R = get_pair_covar_matr(ref_model, alt_model)
#Возвращаем матрицу поворта, полученную из матрицы ковариаций
return get_ort_mat_SVD(R)
def MyKabsch(structure):
'''Получить искомое выравнивание на основе алгоритма Каббша через попарное выравнивание'''
#Переносим все модели в центр координат (переносим их центроиды)
translate_to_center(structure)
ref_model = structure[0]
Q_ans = []
for alt_model in structure:
#Получаем попарный поворот
rot = MyKabsch_get_pair_rot(ref_model,alt_model)
#Применяем их, преобразуя очередную структуру
transform_model(alt_model, np.identity(3), [0,0,0])
Q_ans.append(rot)
return Q_ans
# import os, sys
import Bio.PDB.PDBParser
import Bio.SVDSuperimposer
#import Bio.PDB.Polypeptide
import numpy as np
import os
#Количество различных конформаций для данных PDB (беру ровно 20, хотя в некоторых
#имеется и больше, просто чтобы для всех структур было одинаково)
N = 3
#Для того чтобы потом сделать графики.
pdb_codes = []
Start_RMSDs = []
Start_RMSDs_after_translation = []
End_RMSDs_Kabsch = []
End_RMSDs_MyKabsch = []
End_RMSDs_Wang_Singer_0 = []
End_RMSDs_Wang_Singer_1 = []
End_RMSDs_Wang_Singer_2 = []
ns = []
#Вычисление для набора PDB
for file in os.listdir("./Structures/"):
if file.endswith(".pdb"):
#print(os.path.join("/mydir", file))
pdb_code = os.path.splitext(file.upper())[0]
###print(pdb_code)
pdb_filename = "{}.pdb".format(pdb_code)
pdb_out_filename = "%s_aligned.pdb" % pdb_code
print ("Loading PDB file %s" % pdb_filename)
structure = Bio.PDB.PDBParser().get_structure(pdb_code, "./Structures/"+pdb_filename.lower())
structure_first_N_models = Bio.PDB.Structure.Structure(pdb_code)
for model in list(structure)[:N]:
structure_first_N_models.add(model.copy())
Super_Start_RMSD = RMSD_pairwise_avg(structure_first_N_models)
random_rotate_structure(structure_first_N_models)
Start_RMSD = RMSD_pairwise_avg(structure_first_N_models)
eps = (list((structure_first_N_models.get_atoms()))[0]).__sub__(list((structure_first_N_models.get_atoms()))[1])/4
#print(eps)
random_translate_atoms_in_structure(structure_first_N_models, eps)
#io=Bio.PDB.PDBIO()
#io.set_structure(structure)
#io.save(structure.id + '_out.pdb')
#!pymol 1JOY_out.pdb
#Транслируем все модели к центру координат
translate_to_center(structure_first_N_models)
Start_RMSD_after_translation = RMSD_pairwise_avg(structure_first_N_models)
Q_ans = Kabsch(structure_first_N_models)
transform_structure(structure_first_N_models, Q_ans)
End_RMSD_Kabsch = RMSD_pairwise_avg(structure_first_N_models)
detransform_structure(structure_first_N_models, Q_ans)
#Q_ans = MyKabsch(structure_first_N_models)
#transform_structure(structure_first_N_models, Q_ans)
#End_RMSD_MyKabsch = RMSD_pairwise_avg(structure_first_N_models)
#detransform_structure(structure_first_N_models, Q_ans)
G_ext = 0
Q_ans, G_ext = Wang_Singer_SDP(structure_first_N_models, 0, 0, G_ext)
#print(Q_ans)
transform_structure(structure_first_N_models, Q_ans)
End_RMSD_Wang_Singer_0 = RMSD_pairwise_avg(structure_first_N_models)
detransform_structure(structure_first_N_models, Q_ans)
Q_ans, G_ext = Wang_Singer_SDP(structure_first_N_models, 1, 1, G_ext)
transform_structure(structure_first_N_models, Q_ans)
End_RMSD_Wang_Singer_1 = RMSD_pairwise_avg(structure_first_N_models)
detransform_structure(structure_first_N_models, Q_ans)
Q_ans, G_ext = Wang_Singer_SDP(structure_first_N_models, 2, 1, G_ext)
#io=Bio.PDB.PDBIO()
#io.set_structure(structure_first_N_models)
#io.save(structure_first_N_models.id + '_before_out.pdb')
transform_structure(structure_first_N_models, Q_ans)
#io=Bio.PDB.PDBIO()
#io.set_structure(structure_first_N_models)
#io.save(structure_first_N_models.id + '_after_out.pdb')
End_RMSD_Wang_Singer_2 = RMSD_pairwise_avg(structure_first_N_models)
#detransform_structure(structure_first_N_models, Q_ans)
n = sum([len(ref_chain) for ref_chain in structure_first_N_models[0]])
print("pdb=%s tr)%s k)%s ws0)%s ws1)%s ws2)%s"%(pdb_code, Start_RMSD_after_translation,
End_RMSD_Kabsch, End_RMSD_Wang_Singer_0, End_RMSD_Wang_Singer_1, End_RMSD_Wang_Singer_2))
#if(pdb_code == "5KK9"):
# io=Bio.PDB.PDBIO()
# io.set_structure(structure)
# io.save(structure.id + '_out.pdb')
pdb_codes.append(pdb_code)
Start_RMSDs.append(Start_RMSD)
Start_RMSDs_after_translation.append(Start_RMSD_after_translation)
End_RMSDs_Kabsch.append(End_RMSD_Kabsch)
#End_RMSDs_MyKabsch.append(End_RMSD_MyKabsch)
End_RMSDs_Wang_Singer_0.append(End_RMSD_Wang_Singer_0)
End_RMSDs_Wang_Singer_1.append(End_RMSD_Wang_Singer_1)
End_RMSDs_Wang_Singer_2.append(End_RMSD_Wang_Singer_2)
ns.append(n)
Loading PDB file 5N08.pdb eigenvectors: [ -3.42129806e-08 1.70223428e-07 2.72962319e-07 2.99999990e+00 2.99999997e+00 3.00000033e+00] columns: stochas: pdb=5N08 tr)28.5057944661 k)4.98764127368 ws0)4.98401991287 ws1)4.98401991356 ws2)4.98404647741 Loading PDB file 4YWA.pdb eigenvectors: [ 2.57383244e-08 8.62315815e-08 2.72365727e-07 2.99999990e+00 2.99999997e+00 3.00000007e+00] columns: stochas: pdb=4YWA tr)26.5460699467 k)8.53622000849 ws0)8.53173050519 ws1)8.53173050514 ws2)8.53173335475 Loading PDB file 1Z6E.pdb eigenvectors: [ -1.94072340e-08 -1.13893476e-08 4.87726222e-08 3.00000005e+00 3.00000011e+00 3.00000023e+00] columns: stochas: pdb=1Z6E tr)17.9593225437 k)6.99730392219 ws0)6.98289889334 ws1)6.98289889337 ws2)6.98289918975 Loading PDB file 5ICJ.pdb eigenvectors: [ 3.27659777e-08 9.49798448e-08 2.39983530e-07 2.99999990e+00 2.99999995e+00 3.00000023e+00] columns: stochas: pdb=5ICJ tr)18.8135583011 k)7.35535700203 ws0)7.35194114671 ws1)7.35194114667 ws2)7.35195754295 Loading PDB file 5TEC.pdb eigenvectors: [ 2.53746445e-08 4.91699228e-08 1.00686213e-07 2.99999996e+00 2.99999996e+00 3.00000004e+00] columns: stochas: pdb=5TEC tr)20.6959588053 k)6.90141926699 ws0)6.89987652345 ws1)6.89987652296 ws2)6.89989485807 Loading PDB file 5JFE.pdb eigenvectors: [ 4.12240088e-09 6.93582802e-08 2.22790467e-07 2.99999989e+00 2.99999999e+00 3.00000024e+00] columns: stochas: pdb=5JFE tr)22.0417104504 k)6.96044666109 ws0)6.94002108605 ws1)6.94002108615 ws2)6.94002460451 Loading PDB file 5HU4.pdb eigenvectors: [ -1.57237429e-08 -1.20007951e-08 -8.23391459e-09 3.00000003e+00 3.00000004e+00 3.00000004e+00] columns: stochas: pdb=5HU4 tr)12.0051134407 k)5.52613132491 ws0)5.51608041489 ws1)5.51608041497 ws2)5.51607941756 Loading PDB file 5HY9.pdb eigenvectors: [ 5.05825493e-07 9.13321859e-07 1.40859855e-06 2.99999948e+00 3.00000022e+00 3.00000081e+00] columns: stochas: pdb=5HY9 tr)25.3673168458 k)8.13058307775 ws0)8.04898144903 ws1)8.04898143454 ws2)8.04922344314 Loading PDB file 5T4U.pdb eigenvectors: [ -1.20116514e-08 -8.02968865e-10 1.67848332e-08 3.00000001e+00 3.00000002e+00 3.00000005e+00] columns: stochas: pdb=5T4U tr)15.0446426048 k)5.39624673338 ws0)5.3773069702 ws1)5.37730697014 ws2)5.37730849785 Loading PDB file 5T89.pdb eigenvectors: [ -9.01175130e-09 1.40830521e-08 2.27151619e-07 2.99999996e+00 3.00000003e+00 3.00000005e+00] columns: stochas: pdb=5T89 tr)53.2160962024 k)9.94013451856 ws0)9.93605773107 ws1)9.93605773117 ws2)9.93609314488 Loading PDB file 4Z4P.pdb
/usr/lib/python3.6/site-packages/Bio/PDB/Atom.py:98: PDBConstructionWarning: Could not assign element 'Z' for Atom (name=ZN) with given element '' warnings.warn(msg, PDBConstructionWarning)
eigenvectors: [ 3.08927139e-08 1.47178453e-07 3.09170588e-07 2.99999986e+00 2.99999996e+00 3.00000036e+00] columns: stochas: pdb=4Z4P tr)10.4417872192 k)5.6046081278 ws0)5.58476372239 ws1)5.58476372235 ws2)5.58475040861 Loading PDB file 5UGW.pdb eigenvectors: [ 3.87691009e-08 1.36118021e-07 5.18293874e-07 2.99999973e+00 2.99999999e+00 3.00000050e+00] columns: stochas: pdb=5UGW tr)13.640743744 k)6.55559546305 ws0)6.54492806013 ws1)6.54492805515 ws2)6.54501432854 Loading PDB file 4YW6.pdb eigenvectors: [ 1.53047084e-07 2.63920809e-07 5.87674464e-07 2.99999964e+00 2.99999996e+00 3.00000026e+00] columns: stochas: pdb=4YW6 tr)21.074048194 k)9.56716301001 ws0)9.50525789255 ws1)9.50525789012 ws2)9.50526546457 Loading PDB file 5MX3.pdb eigenvectors: [ 6.24937858e-08 1.00085307e-07 2.17756175e-07 2.99999991e+00 2.99999996e+00 3.00000004e+00] columns: stochas: pdb=5MX3 tr)35.9693484792 k)10.145735789 ws0)10.1446591953 ws1)10.1446591953 ws2)10.1446686701 Loading PDB file 5GSV.pdb eigenvectors: [ 1.77632349e-07 3.29619097e-07 4.50271048e-07 2.99999984e+00 3.00000001e+00 3.00000018e+00] columns: stochas: pdb=5GSV tr)21.036648663 k)7.49924747144 ws0)7.48944871679 ws1)7.48944871636 ws2)7.48946062409 Loading PDB file 5HYD.pdb eigenvectors: [ 1.35267486e-07 2.65786150e-07 4.11754262e-07 2.99999958e+00 2.99999996e+00 3.00000032e+00] columns: stochas: pdb=5HYD tr)25.1715904901 k)6.95628913945 ws0)6.93787220631 ws1)6.93787220869 ws2)6.9379069559 Loading PDB file 5MMO.pdb eigenvectors: [ -1.26534981e-09 3.59733124e-08 4.60121486e-08 2.99999999e+00 2.99999999e+00 3.00000002e+00] columns: stochas: pdb=5MMO tr)18.0856957806 k)6.52271528183 ws0)6.49438702824 ws1)6.49438702804 ws2)6.49438693461 Loading PDB file 5ML6.pdb eigenvectors: [ -4.64188766e-08 2.80372789e-07 4.26102676e-07 2.99999984e+00 3.00000004e+00 3.00000031e+00] columns: stochas: pdb=5ML6 tr)14.75831778 k)4.53979446522 ws0)4.53554543338 ws1)4.53554543295 ws2)4.53554216891 Loading PDB file 5GVJ.pdb eigenvectors: [ -7.36328461e-09 9.32391061e-08 1.40302781e-07 2.99999994e+00 3.00000000e+00 3.00000017e+00] columns: stochas: pdb=5GVJ tr)25.7407530131 k)6.71856165134 ws0)6.71217199226 ws1)6.71217200015 ws2)6.71210130342
print(Start_RMSDs_after_translation)
print(End_RMSDs_Kabsch)
print(End_RMSDs_Wang_Singer_0)
print(End_RMSDs_Wang_Singer_1)
print(End_RMSDs_Wang_Singer_2)
[28.505794466139282, 26.546069946701763, 17.959322543716691, 18.813558301135419, 20.695958805312799, 22.041710450426397, 12.005113440736702, 25.367316845834814, 15.044642604847228, 53.216096202400244, 10.441787219160018, 13.640743743987152, 21.074048193972384, 35.969348479166086, 21.036648663020308, 25.171590490115381, 18.085695780621585, 14.758317779968792, 25.740753013120056] [4.9876412736828772, 8.536220008487243, 6.9973039221891362, 7.3553570020265413, 6.9014192669921792, 6.9604466610874418, 5.5261313249113151, 8.1305830777476356, 5.3962467333759339, 9.9401345185645766, 5.6046081278048048, 6.5555954630492703, 9.5671630100072349, 10.145735789022019, 7.4992474714357735, 6.9562891394486757, 6.5227152818251986, 4.539794465216862, 6.7185616513393329] [4.9840199128736788, 8.5317305051911489, 6.9828988933412495, 7.3519411467143421, 6.8998765234474968, 6.9400210860539095, 5.5160804148903013, 8.0489814490286893, 5.3773069701974174, 9.9360577310745928, 5.5847637223914077, 6.5449280601322117, 9.5052578925473981, 10.144659195328995, 7.4894487167893322, 6.9378722063131129, 6.4943870282428202, 4.5355454333847778, 6.7121719922563399] [4.9840199135550263, 8.5317305051412848, 6.9828988933704537, 7.3519411466702378, 6.8998765229579995, 6.9400210861450056, 5.516080414967095, 8.0489814345384687, 5.3773069701415244, 9.9360577311689404, 5.5847637223467057, 6.5449280551486595, 9.5052578901223796, 10.144659195330922, 7.4894487163626353, 6.9378722086871329, 6.4943870280428007, 4.5355454329541622, 6.712172000148497] [4.9840464774101605, 8.5317333547538343, 6.9828991897476236, 7.3519575429537412, 6.8998948580709092, 6.9400246045076743, 5.5160794175622447, 8.0492234431395726, 5.3773084978488477, 9.9360931448784662, 5.5847504086135844, 6.5450143285410292, 9.5052654645681436, 10.144668670075063, 7.4894606240929464, 6.9379069559004813, 6.4943869346102465, 4.5355421689127375, 6.7121013034201749]
import pandas as pd
zipped = (End_RMSDs_Kabsch, Start_RMSDs_after_translation, End_RMSDs_Wang_Singer_0,
End_RMSDs_Wang_Singer_1, End_RMSDs_Wang_Singer_2, pdb_codes)
#zipped.sort()
results = np.column_stack(zipped)
df = pd.DataFrame(results)
print (df.to_latex())
\begin{tabular}{lllllll} \toprule {} & 0 & 1 & 2 & 3 & 4 & 5 \\ \midrule 0 & 4.987641273682877 & 28.50579446613928 & 4.984019912873679 & 4.984019913555026 & 4.9840464774101605 & 5N08 \\ 1 & 8.536220008487243 & 26.546069946701763 & 8.531730505191149 & 8.531730505141285 & 8.531733354753834 & 4YWA \\ 2 & 6.997303922189136 & 17.95932254371669 & 6.9828988933412495 & 6.982898893370454 & 6.982899189747624 & 1Z6E \\ 3 & 7.355357002026541 & 18.81355830113542 & 7.351941146714342 & 7.351941146670238 & 7.351957542953741 & 5ICJ \\ 4 & 6.901419266992179 & 20.6959588053128 & 6.899876523447497 & 6.8998765229579995 & 6.899894858070909 & 5TEC \\ 5 & 6.960446661087442 & 22.041710450426397 & 6.9400210860539095 & 6.940021086145006 & 6.940024604507674 & 5JFE \\ 6 & 5.526131324911315 & 12.005113440736702 & 5.516080414890301 & 5.516080414967095 & 5.516079417562245 & 5HU4 \\ 7 & 8.130583077747636 & 25.367316845834814 & 8.04898144902869 & 8.048981434538469 & 8.049223443139573 & 5HY9 \\ 8 & 5.396246733375934 & 15.044642604847228 & 5.377306970197417 & 5.377306970141524 & 5.377308497848848 & 5T4U \\ 9 & 9.940134518564577 & 53.216096202400244 & 9.936057731074593 & 9.93605773116894 & 9.936093144878466 & 5T89 \\ 10 & 5.604608127804805 & 10.441787219160018 & 5.584763722391408 & 5.584763722346706 & 5.584750408613584 & 4Z4P \\ 11 & 6.55559546304927 & 13.640743743987152 & 6.544928060132212 & 6.5449280551486595 & 6.545014328541029 & 5UGW \\ 12 & 9.567163010007235 & 21.074048193972384 & 9.505257892547398 & 9.50525789012238 & 9.505265464568144 & 4YW6 \\ 13 & 10.14573578902202 & 35.96934847916609 & 10.144659195328995 & 10.144659195330922 & 10.144668670075063 & 5MX3 \\ 14 & 7.4992474714357735 & 21.03664866302031 & 7.489448716789332 & 7.489448716362635 & 7.489460624092946 & 5GSV \\ 15 & 6.956289139448676 & 25.17159049011538 & 6.937872206313113 & 6.937872208687133 & 6.937906955900481 & 5HYD \\ 16 & 6.5227152818251986 & 18.085695780621585 & 6.49438702824282 & 6.494387028042801 & 6.4943869346102465 & 5MMO \\ 17 & 4.539794465216862 & 14.758317779968792 & 4.535545433384778 & 4.535545432954162 & 4.5355421689127375 & 5ML6 \\ 18 & 6.718561651339333 & 25.740753013120056 & 6.71217199225634 & 6.712172000148497 & 6.712101303420175 & 5GVJ \\ \bottomrule \end{tabular}
zipped = zip(End_RMSDs_Kabsch, Start_RMSDs_after_translation, End_RMSDs_Wang_Singer_0,
End_RMSDs_Wang_Singer_1, End_RMSDs_Wang_Singer_2, pdb_codes)
zipped.sort()
End_RMSDs_Kabsch, Start_RMSDs_after_translation, End_RMSDs_Wang_Singer_0, End_RMSDs_Wang_Singer_1, End_RMSDs_Wang_Singer_2, pdb_codes = zip(*zipped)
Start_RMSDs_after_translation = [20.992278940347511, 24.743422863430773, 29.059301765212766, 44.406963477085014, 18.984985814352225, 36.658388466513451, 29.458873049039536, 20.721102462408137, 34.896421088833371, 18.15698054900172, 39.408029099426535, 31.693173368391797, 38.771015519949529, 36.722618495225248, 20.115682305895291, 81.31222477177667, 28.902379594390037, 24.119635878539714, 86.451381076922587, 30.349716857838946]
End_RMSDs_Kabsch = [8.6803676689681168, 8.6581763928494802, 10.914479879555632, 12.90871248456833, 8.6083568383792262, 12.374481146504074, 10.842910144670677, 8.3764871833234764, 12.175707515036379, 7.4839852394140518, 11.755192066926121, 9.4896011381547751, 10.269951741023654, 10.470387305749464, 7.8556185941566943, 13.07394875791489, 11.066477139715442, 8.409734418518342, 13.489417675819135, 11.035510184515156]
End_RMSDs_Wang_Singer_0 = [8.6247685332496218, 8.5624436365943311, 10.845892595869472, 12.850207398806615, 8.5304007765353695, 12.291373572848572, 10.800094744089225, 8.3362052569883858, 12.113353918049711, 7.3809006546917395, 11.701766112678015, 9.4490771711987023, 10.222737265048847, 10.417780392310945, 7.8200220855300024, 13.057382880909028, 11.006972859821948, 8.3708571068165174, 12.773439948599405, 10.983384095694241]
End_RMSDs_Wang_Singer_1 = [8.6247685328724639, 8.5624436363790792, 10.845892595701972, 12.850207398206905, 8.530400776521363, 12.29137357301399, 10.800094743058027, 8.336205256903062, 12.113353917729505, 7.3809006546298681, 11.701766111581168, 9.449077171883534, 10.222737263629325, 10.417780390933466, 7.8200220853674329, 13.057382882177222, 11.00697285981618, 8.370857106543653, 12.773439976684029, 10.983384095662368]
End_RMSDs_Wang_Singer_2 = [8.6247746220855799, 8.5624519109529764, 10.845913579938497, 12.850273104549954, 8.5304026179798758, 12.291454197130331, 10.800177996345488, 8.3362072058476606, 12.113368345124719, 7.3809080407477712, 11.701814134292521, 9.4491940856930103, 10.22289899699104, 10.417866009139972, 7.8200241216908539, 13.057956758621156, 11.006974607217131, 8.3708679716723609, 12.774905857150928, 10.983401974623618]
import matplotlib.pyplot as plt
%matplotlib inline
import numpy
plt.figure(figsize=(12,8))
plt.rcParams.update({'font.size': 20})
##x = numpy.arange(0, 20, 1)
##y = numpy.arange(0, 30000, 1500)
##plt.scatter(x, y)
ax = plt.subplot(111)
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs, color='#000000', label="Start")
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs_after_translation, color='#00bb00', label="Without translations")
plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Kabsch, color='#bb0000', label="Pairwise")
#plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_0, color='#0000aa', label="SDP with eigenvectors")
#plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_1, color='#00aaaa', label="SDP-first3")
plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_2, color='#aa00aa', label="SDP stochastic")
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs, color='#ff0000')
plt.setp(ax.lines, linewidth=2.0);
#fig = plt.figure()
ax = plt.gca()
ax.set_xticks(range(1,len(pdb_codes)+1))
ax.set_xticklabels(pdb_codes, rotation=75, rotation_mode="anchor", ha="right");
plt.legend(bbox_to_anchor=(0.35, 0.95), loc=1, borderaxespad=0.)
plt.grid()
#ax.set_xticks(range(1,len(pdb_codes)+1))
#ax.set_xticklabels(pdb_codes, rotation=75, rotation_mode="anchor", ha="right");
#plt.title('Pairwise/SDP-first3')
plt.title('Pairwise/SDP stochastic')
#plt.title('Pairwise/SDP with eigenvectors')
#plt.title('Pairwise/SDP')
plt.xlabel('PDB codes')
plt.ylabel('APRMSD');
plt.savefig('30_without_tran_3.eps', format='eps', dpi=300)
End_RMSDs_Wang_Singer_1_without = [x/y for (x,y) in zip(End_RMSDs_Wang_Singer_1,End_RMSDs_Wang_Singer_0)]
print(End_RMSDs_Wang_Singer_1_without)
End_RMSDs_Wang_Singer_2_without = [x/y for (x,y) in zip(End_RMSDs_Wang_Singer_2,End_RMSDs_Wang_Singer_0)]
print(End_RMSDs_Wang_Singer_2_without)
End_RMSDs_Wang_Singer_0_without = [x/y for (x,y) in zip(End_RMSDs_Wang_Singer_0,End_RMSDs_Wang_Singer_0)]
print(End_RMSDs_Wang_Singer_0_without)
[0.99999999999161737, 0.99999999997921107, 0.99999999998976463, 0.99999999996740307, 0.99999999999835809, 0.99999999997486089, 0.99999999995627042, 1.000000000072476, 0.99999999986114063, 0.9999999998677761, 0.99999999990451949, 0.99999999998455635, 0.9999999999970981, 0.99999999999947597, 0.99999999990626653, 0.99999999997356581, 1.0000000000134581, 0.99999999995333066, 1.0000000000971248, 1.0000000021986735] [1.0000010006984754, 1.0000002603778901, 1.0000002337825442, 1.0000012979382762, 1.0000002158684631, 1.0000009663548159, 1.000000705970941, 1.0000123731124415, 1.0000158208059151, 1.0000082183369012, 1.0000077084746231, 1.0000019347480016, 1.0000016278160921, 1.0000001587534744, 1.0000041037920295, 1.0000011910058193, 1.0000065594200096, 1.0000051132048924, 1.0000439504391778, 1.0001147622376918] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
print(min([x/y for x,y in zip(End_RMSDs_Kabsch, Start_RMSDs_after_translation)]))
print(min([x/y for x,y in zip(End_RMSDs_Wang_Singer_0, Start_RMSDs_after_translation)]))
print(min([x/y for x,y in zip(End_RMSDs_Wang_Singer_1, Start_RMSDs_after_translation)]))
print(min([x/y for x,y in zip(End_RMSDs_Wang_Singer_2, Start_RMSDs_after_translation)]))
0.169867884985 0.159891343501 0.15989134438 0.159897281275
print(max([(a-b)/a for a,b in zip([x/y for x,y in zip(End_RMSDs_Kabsch, Start_RMSDs_after_translation)],
[x/y for x,y in zip(End_RMSDs_Wang_Singer_0, Start_RMSDs_after_translation)])]))
print(max([(a-b)/a for a,b in zip([x/y for x,y in zip(End_RMSDs_Kabsch, Start_RMSDs_after_translation)],
[x/y for x,y in zip(End_RMSDs_Wang_Singer_1, Start_RMSDs_after_translation)])]))
print(max([(a-b)/a for a,b in zip([x/y for x,y in zip(End_RMSDs_Kabsch, Start_RMSDs_after_translation)],
[x/y for x,y in zip(End_RMSDs_Wang_Singer_2, Start_RMSDs_after_translation)])]))
0.0587311809119 0.0587311757345 0.0586962256652
End_RMSDs_Wang_Singer_1_without = [x/y for (x,y) in zip(End_RMSDs_Wang_Singer_1,End_RMSDs_Wang_Singer_0)]
print(End_RMSDs_Wang_Singer_1_without)
End_RMSDs_Wang_Singer_2_without = [x/y for (x,y) in zip(End_RMSDs_Wang_Singer_2,End_RMSDs_Wang_Singer_0)]
print(End_RMSDs_Wang_Singer_2_without)
End_RMSDs_Wang_Singer_0_without = [x/y for (x,y) in zip(End_RMSDs_Wang_Singer_0,End_RMSDs_Wang_Singer_0)]
print(End_RMSDs_Wang_Singer_0_without)
import matplotlib.pyplot as plt
%matplotlib inline
import numpy
plt.figure(figsize=(12,8))
plt.rcParams.update({'font.size': 20})
##x = numpy.arange(0, 20, 1)
##y = numpy.arange(0, 30000, 1500)
##plt.scatter(x, y)
#zipped = zip(End_RMSDs_Kabsch, Start_RMSDs_after_translation, End_RMSDs_Wang_Singer_0, End_RMSDs_Wang_Singer_1, End_RMSDs_Wang_Singer_2, pdb_codes)
#ax = plt.subplot(111)
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs, color='#000000', label="Start")
#plt.plot(range(1,len(pdb_codes)+1), zipped[1], color='#00bb00', label="Without translations")
#plt.plot(range(1,len(pdb_codes)+1), zipped[0], color='#bb0000', label="Pairwise")
plt.plot(range(1,len(pdb_codes[:7])+1), End_RMSDs_Wang_Singer_0_without[:7], color='#0000aa', label="SDP1")
plt.plot(range(1,len(pdb_codes[:7])+1), End_RMSDs_Wang_Singer_1_without[:7], color='#00aaaa', label="SDP2")
plt.plot(range(1,len(pdb_codes[:7])+1), End_RMSDs_Wang_Singer_2_without[:7], color='#aa00aa', label="SDP3")
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs, color='#ff0000')
plt.setp(ax.lines, linewidth=2.0);
#fig = plt.figure()
ax = plt.gca()
ax.set_xticks(range(1,pdb_codes[:7]+1))
ax.set_xticklabels(pdb_codes, rotation=75, rotation_mode="anchor", ha="right");
plt.legend(bbox_to_anchor=(0.25, 0.95), loc=1, borderaxespad=0.)
plt.grid()
#ax.set_xticks(range(1,len(pdb_codes)+1))
#ax.set_xticklabels(pdb_codes, rotation=75, rotation_mode="anchor", ha="right");
plt.title('APRMSD, diff from SDP1')
plt.xlabel('PDB codes')
plt.ylabel('RMSD');
plt.savefig('30_without_kabsch.eps', format='eps', dpi=300)
[0.99999999999161737, 0.99999999997921107, 0.99999999998976463, 0.99999999996740307, 0.99999999999835809, 0.99999999997486089, 0.99999999995627042, 1.000000000072476, 0.99999999986114063, 0.9999999998677761, 0.99999999990451949, 0.99999999998455635, 0.9999999999970981, 0.99999999999947597, 0.99999999990626653, 0.99999999997356581, 1.0000000000134581, 0.99999999995333066, 1.0000000000971248, 1.0000000021986735] [1.0000010006984754, 1.0000002603778901, 1.0000002337825442, 1.0000012979382762, 1.0000002158684631, 1.0000009663548159, 1.000000705970941, 1.0000123731124415, 1.0000158208059151, 1.0000082183369012, 1.0000077084746231, 1.0000019347480016, 1.0000016278160921, 1.0000001587534744, 1.0000041037920295, 1.0000011910058193, 1.0000065594200096, 1.0000051132048924, 1.0000439504391778, 1.0001147622376918] [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-134-8ca7b35ce390> in <module>() 34 ax = plt.gca() 35 ---> 36 ax.set_xticks(range(1,pdb_codes[:7]+1)) 37 ax.set_xticklabels(pdb_codes, rotation=75, rotation_mode="anchor", ha="right"); 38 TypeError: can only concatenate tuple (not "int") to tuple
'''Просто чтобы посмотреть на тот же график в другом виде'''
import math
for i,a in enumerate(Start_RMSDs):
Start_RMSDs[i]/=Start_RMSDs_after_translation[i]
End_RMSDs_Kabsch[i]/=Start_RMSDs_after_translation[i]
End_RMSDs_MyKabsch[i]/=Start_RMSDs_after_translation[i]
End_RMSDs_Wang_Singer_0[i]/=Start_RMSDs_after_translation[i]
End_RMSDs_Wang_Singer_1[i]/=Start_RMSDs_after_translation[i]
End_RMSDs_Wang_Singer_2[i]/=Start_RMSDs_after_translation[i]
Start_RMSDs_after_translation[i]/=Start_RMSDs_after_translation[i]
import math
for i,a in enumerate(Start_RMSDs):
Start_RMSDs[i] = math.log(Start_RMSDs[i])
Start_RMSDs_after_translation[i] = math.log(Start_RMSDs_after_translation[i])
End_RMSDs_Kabsch[i] = math.log(End_RMSDs_Kabsch[i])
End_RMSDs_MyKabsch[i] = math.log(End_RMSDs_MyKabsch[i])
End_RMSDs_Wang_Singer_0[i] = math.log(End_RMSDs_Wang_Singer_0[i])
End_RMSDs_Wang_Singer_1[i] = math.log(End_RMSDs_Wang_Singer_1[i])
End_RMSDs_Wang_Singer_2[i] = math.log(End_RMSDs_Wang_Singer_2[i])
import math
for i,a in enumerate(Start_RMSDs):
#Start_RMSDs[i] = math.log(Start_RMSDs[i])
#Start_RMSDs_after_translation[i] = math.log(Start_RMSDs_after_translation[i])
End_RMSDs_MyKabsch[i] -= End_RMSDs_Kabsch[i]
End_RMSDs_Wang_Singer_0[i] -= End_RMSDs_Kabsch[i]
End_RMSDs_Wang_Singer_1[i] -= End_RMSDs_Kabsch[i]
End_RMSDs_Wang_Singer_2[i] -= End_RMSDs_Kabsch[i]
End_RMSDs_Kabsch[i] -= End_RMSDs_Kabsch[i]
import matplotlib.pyplot as plt
%matplotlib inline
import numpy
plt.figure(figsize=(12,8))
plt.rcParams.update({'font.size': 20})
##x = numpy.arange(0, 20, 1)
##y = numpy.arange(0, 30000, 1500)
##plt.scatter(x, y)
ax = plt.subplot(111)
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs, color='#000000', label="before")
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs_after_translation, color='#00ff00', label="aft tran")
plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Kabsch, color='#ff0000', label="kabsch")
plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_MyKabsch, color='#0000ff', label="my kabsch")
plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_0, color='#00ffff', label="ws 0")
plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_1, color='#ffff00', label="ws 1")
plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_2, color='#ff00ff', label="ws 2")
#plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs, color='#ff0000')
plt.setp(ax.lines, linewidth=2.0);
#fig = plt.figure()
ax = plt.gca()
ax.set_xticks(range(1,len(pdb_codes)+1))
ax.set_xticklabels(pdb_codes, rotation=75, rotation_mode="anchor", ha="right");
plt.legend(bbox_to_anchor=(0.30, 0.95), loc=2, borderaxespad=0.)
plt.grid()
#ax.set_xticks(range(1,len(pdb_codes)+1))
#ax.set_xticklabels(pdb_codes, rotation=75, rotation_mode="anchor", ha="right");
plt.title('Kabsch + SDP relaxations RMSDs')
plt.xlabel('PDB codes')
plt.ylabel('RMSD');
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-53-7d4b4616b9ed> in <module>() 16 #plt.plot(range(1,len(pdb_codes)+1), Start_RMSDs_after_translation, color='#00ff00', label="aft tran") 17 plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Kabsch, color='#ff0000', label="kabsch") ---> 18 plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_MyKabsch, color='#0000ff', label="my kabsch") 19 plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_0, color='#00ffff', label="ws 0") 20 plt.plot(range(1,len(pdb_codes)+1), End_RMSDs_Wang_Singer_1, color='#ffff00', label="ws 1") /home/kesha/anaconda2/lib/python2.7/site-packages/matplotlib/pyplot.pyc in plot(*args, **kwargs) 3316 mplDeprecation) 3317 try: -> 3318 ret = ax.plot(*args, **kwargs) 3319 finally: 3320 ax._hold = washold /home/kesha/anaconda2/lib/python2.7/site-packages/matplotlib/__init__.pyc in inner(ax, *args, **kwargs) 1890 warnings.warn(msg % (label_namer, func.__name__), 1891 RuntimeWarning, stacklevel=2) -> 1892 return func(ax, *args, **kwargs) 1893 pre_doc = inner.__doc__ 1894 if pre_doc is None: /home/kesha/anaconda2/lib/python2.7/site-packages/matplotlib/axes/_axes.pyc in plot(self, *args, **kwargs) 1404 kwargs = cbook.normalize_kwargs(kwargs, _alias_map) 1405 -> 1406 for line in self._get_lines(*args, **kwargs): 1407 self.add_line(line) 1408 lines.append(line) /home/kesha/anaconda2/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in _grab_next_args(self, *args, **kwargs) 405 return 406 if len(remaining) <= 3: --> 407 for seg in self._plot_args(remaining, kwargs): 408 yield seg 409 return /home/kesha/anaconda2/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in _plot_args(self, tup, kwargs) 383 x, y = index_of(tup[-1]) 384 --> 385 x, y = self._xy_from_xy(x, y) 386 387 if self.command == 'plot': /home/kesha/anaconda2/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in _xy_from_xy(self, x, y) 242 if x.shape[0] != y.shape[0]: 243 raise ValueError("x and y must have same first dimension, but " --> 244 "have shapes {} and {}".format(x.shape, y.shape)) 245 if x.ndim > 2 or y.ndim > 2: 246 raise ValueError("x and y can be no greater than 2-D, but have " ValueError: x and y must have same first dimension, but have shapes (1,) and (0,)
'''
io=Bio.PDB.PDBIO()
io.set_structure(structure)
io.save(structure.id + '_out.pdb')
!pymol 2CH7_100_out_out.pdb
def get_random_rot_mat(dim=3):
#Получить случайную матрицу поворота
random_state = np.random
H = np.eye(dim)
D = np.ones((dim,))
for n in range(1, dim):
x = random_state.normal(size=(dim-n+1,))
D[n-1] = np.sign(x[0])
x[0] -= D[n-1]*np.sqrt((x*x).sum())
# Householder transformation
Hx = (np.eye(dim-n+1) - 2.*np.outer(x, x)/(x*x).sum())
mat = np.eye(dim)
mat[n-1:, n-1:] = Hx
H = np.dot(H, mat)
# Fix the last sign such that the determinant is 1
D[-1] = (-1)**(1-(dim % 2))*D.prod()
# Equivalent to np.dot(np.diag(D), H) but faster, apparently
H = (D*H.T).T
return H'''