# Поставить свежий проди локально:
# pip install biopython
# pip install prody
import prody as pd
import numpy as np # Модуль для векторных и матричных операций и всякой линейной алгебры, которая нам тут не пригодится
Я скачал только ту из двух рибосом, в которой есть амикумацин (на сайте PDB можно скачать Biological Assembly, это то, в каком виде эта молекула по мнению ученых должна выглядеть и работать в оригинале, а не в условиях кристаллографического эксперимента, которые зачастую очень отличаются от физиологических).
R = pd.parseCIF('../../Downloads/5i4l-assembly2.cif')
206к атомов это просто запредельно для нашей задачи. В идеале всего атомов, вместе с водой, которую мы зальем вокруг, должно быть около 100к, а еще лучше, если меньше. А тут только одного белка 200к. Надо что-то делать!
ami = R.select('resname UAM') # выделил только атомы амикумацина
Читать про выделения тут: http://prody.csb.pitt.edu/manual/reference/atomic/select.html
?ami # это свой тип данных (Selection), по сути это итератор по атомам. Его можно преобразовать в список: list(ami)
Чтобы посмотреть список методов объекта, надо вбить его название и нажать Tab.
ami_c = ami.getCoords()
print(ami_c) # список координат всех атомов амикумацина
ami_center = np.mean(ami_c, axis = 0) # возьмем среднее по строкам
ami_center # и вот наш центер амикумацина (не центр масс, конечно, но тоже сгодится)
within_30 = R.select('within 30 of center', center = ami_center)
len(within_30) # на расстоянии 30 ангстрем от центра амикумацина находится 3345 атома
def list_chain_breaks(atom_list):
resindices = set(atom_list.getResindices())
# сохранил сет индексов остатков, к которым принадлежит исходный набор атомов
# важно: Resindex это не то же, что и Resnum. Resnum может повторяться,
# Resindex уникален и вообще говоря не совпадает с Resnum.
breaks = []
br_start, br_end = 0, 0
inside_break = False
for i in range(min(resindices), max(resindices) + 1):
if i in resindices:
if inside_break:
# это значит, что мы встретили первый остаток после разрыва
inside_break = False
br_end = i-1
breaks.append((br_start, br_end))
else:
continue
else:
if inside_break:
# это значит, что мы продолжаем идти по разрыву
continue
else:
# мы встретили первый остаток разрыва
br_start = i
inside_break = True
return breaks
breaks_chain_5 = list_chain_breaks(within_30.select('chain 5'))
print(breaks_chain_5)
Voila, индексы остатков разрывов. На наше исходное выделение приходится 11 разрывов только по цепи "5". Теперь узнаем, сколько в них атомов:
for br in breaks_chain_5:
start, end = br
loop = R.select('resindex %d to %d' % (start, end))
print(start, end, len(loop))
within_30_chains = set(within_30.getChids()) # Какие цепи попали в выделение?
within_30_chains
within_30_breaks = {}
for chain in within_30_chains:
atoms = within_30.select('chain %s' % chain)
within_30_breaks[chain] = list_chain_breaks(atoms)
for chain,br in within_30_breaks.items():
print(chain, len(br))
print('Всего возможных разрывов %d' % (sum([26,11,4,8,3,1])))
Всего вариантов попарного их сочетания 2**53, что есть
2**53
Пожалуй, сразу надо бы ввести какие-то ограничения на возможную длину заполняемого разрыва. Допустим, оставим только такие разрывы, длина которых не превосходит 50.
w30br_short50 = {} # horrible variable name, I know
for chain,breaks in within_30_breaks.items():
filtered_breaks = []
for br in breaks:
start, end = br
if end - start + 1 <= 50:
filtered_breaks.append(br)
w30br_short50[chain] = filtered_breaks
N_breaks = 0
for chain,br in w30br_short50.items():
N_breaks += len(br)
print(chain, len(br))
print(2**N_breaks)
Все еще дохрена вариантов. Что ж, поиграемся с порогами длины.
breaks_short = {'5': {},
'10': {},
'15': {},
'20': {},
'30': {},
'40': {}}
for br_maxlen in breaks_short:
N_breaks = 0
for chain,breaks in within_30_breaks.items():
filtered_breaks = []
for br in breaks:
start, end = br
if end - start + 1 <= int(br_maxlen):
filtered_breaks.append(br)
breaks_short[br_maxlen][chain] = filtered_breaks
N_breaks += len(filtered_breaks)
print(br_maxlen, N_breaks, 2**N_breaks)
Варианты 5-20 выглядят подъемными. Посмотрим, как распределятся значения при отсечке в 10.
chains = sorted(list(breaks_short['10'].keys()))
chains
lengths = [len(breaks_short['10'][x]) for x in chains]
lengths
chains = ['5','6','c','d','q']
# зафиксировал определенный порядок цепей и убрал 's', так как в ней ноль разрывов
# config = {'chain': [break_number_1, break_number_3,...]}
def configuration_to_values(config, breaks):
n_atoms = 0
n_breaks = 25 # потому что при отсечке в 10 у нас всего 25 разрывов
for chain, c_indices in config.items():
for i in c_indices:
br = breaks[chain][i]
start, end = br
n_atoms += end - start + 1 # спойлер: это не число атомов, а число остатков, но осознал я это только сильно позже
n_breaks -= 1
return n_atoms, n_breaks
# binary_code = '0110110001...01' - строчка длиной в 25, кодирующая конфигурацию разрывов
bool('0'), bool('1')
bool(int('0')), bool(int('1'))
def code_to_config(code, chains, breaks):
i = -1
config = {}
for c in chains:
coded_breaks = []
for bi,b in enumerate(breaks[c]):
i += 1
if int(code[i]):
coded_breaks.append(bi)
config[c] = coded_breaks
return config
code_to_config('0100011001001001001101011', chains, breaks_short['10'])
Все такие коды можно сгенерировать с помощью itertools
from itertools import product
all_25_codes = product(['0','1'], repeat=25)
test = product(['a','b'], repeat = 3)
list(test)
codes = []
N_atoms = []
N_breaks = []
int('10011', 2) # чтобы не забивать память, буду сохранять коды как числа
i = 0
for code in all_25_codes:
codestring = ''.join(code) # перевел код в строку
config = code_to_config(codestring, chains, breaks_short['10'])
na,nb = configuration_to_values(config, breaks_short['10'])
codes.append(int(codestring,2))
N_atoms.append(na)
N_breaks.append(nb)
i += 1
if i % 1000000 == 0:
print("Processed %d combinations" % i)
len(codes)
Уии, теперь есть статистика по 33 миллионам вариантов рибовырезки. (2 Гб оперативки и около 10 минут счета)
Есть смысл сохранить данные, чтобы не пересчитывать. Для этого у Питона есть возможность записать компактный бинарный файл со своими объектами внутри.
import pickle as pkl
pkl.dump(codes, open('./ribo_30_10_codes.pkl', 'wb'))
pkl.dump(N_atoms, open('./ribo_30_10_nA.pkl', 'wb'))
pkl.dump(N_breaks, open('./ribo_30_10_nB.pkl', 'wb'))
Суммарно файлы заняли около 300 Мб
import matplotlib.pyplot as plt # Время графиков!
%matplotlib notebook
fig,ax = plt.subplots()
ax.hist(N_atoms)
fig.show()
На этом месте я осознал, что накодил не число атомов, а число остатков, что конечно не одно и то же, и нам интереснее было бы все же число атомов. Но пока можно проанализировать данные и так.
fig,ax = plt.subplots()
ax.hist(N_breaks)
fig.show()
Я боюсь, что полный поточечный график на 33 миллиона точек - перебор, поэтому сначала я построю 2д-гистограмму и определюсь, какие значения меня больше всего интереуют
fig,ax = plt.subplots()
ax.hist2d(N_atoms, N_breaks, bins=[max(N_atoms)+1,max(N_breaks)+1])
fig.show()
Зависимость линейная. Отнормируем все и посчитаем например произведение двух параметров.
N_atoms_np = np.array(N_atoms)
N_atoms_np = N_atoms_np / np.max(N_atoms_np)
N_breaks_np = np.array(N_breaks)
N_breaks_np = N_breaks_np / np.max(N_breaks_np)
N_breaks_np
fig,ax = plt.subplots()
ax.hist2d(N_atoms_np, N_breaks_np, bins=[100,100])
fig.show()
t1 = np.array([0,1,2,3,4])
t1 * t1
Просто произведение нам не даст того, чего хотим, так как случай с 0 разрывов при умножении на явно большое число остатков даст 0, равно как и случай с 25 разрывами и 0 атомами.
Изолируем эти два случая
N_breaks_np[N_breaks_np == 0.] = 1.0
N_atoms_np[N_atoms_np == 0.] = 1.0
nA_nB = N_breaks_np * N_atoms_np
np.argmin(nA_nB)
nA_nB[2]
codes[2]
N_atoms[2]
N_breaks[2]
Такое себе, метрика с умножением делает лучшим вариант, когда явно много разрывов.
N_atoms_np = np.array(N_atoms)
N_atoms_np = N_atoms_np / np.max(N_atoms_np)
N_breaks_np = np.array(N_breaks)
N_breaks_np = N_breaks_np / np.max(N_breaks_np)
nA_plus_nB = N_atoms_np + N_breaks_np
np.argmin(nA_plus_nB)
codes[9344503]
N_atoms[9344503],N_breaks[9344503]
Вот это больше похоже на то, что хочется. 10 разрывов и 38 новых остатков.