Начало практикума посвящено работе с модулем ProDy. Мы используем его возможности для анализа B-факторов, как в практикуме 3, но уже не «на глаз», по цветовой палитре, а численного анализа.
В первом задании посмотрим на остатки с наибольшим и наименьшим средним B-фактором атомов. Для начала найдём их в белке.
Загрузим модули и файл PDB.
import prody as pd, numpy as np
pdb = pd.parsePDB("7bde")
@> PDB file is found in working directory (7bde.pdb.gz). @> 1140 atoms and 1 coordinate set(s) were parsed in 0.02s.
Пройдёмся по остаткам и вычислим B-факторы.
mean_betas = list()
for residue in pdb.iterResidues():
if "CA" in residue.getNames():
mean_beta = np.mean(residue.getBetas())
mean_betas.append((residue, mean_beta))
min_b_res, max_b_res = sorted(mean_betas, key = lambda elem: elem[1])[0], \
sorted(mean_betas, key = lambda elem: elem[1])[-1]
min_b_res, max_b_res
((<Residue: LEU 22 from Chain A from 7bde (8 atoms)>, 30.37875), (<Residue: ALA 127 from Chain A from 7bde (5 atoms)>, 111.63799999999999))
Таким образом, остаток с наименьшим B-фактором — L22, с наибольшим — A127. Они отображены на рис. .1.
Остаток с минимальным B-фактором входит в альфа-спираль, это лейцин в окружении других гидрофобных разветвлённых остатков, должно быть, его движение стерически затруднено. Остаток с максимальным B-фактором предтерминальный, это небольшой остаток аланина вблизи конца альфа-спирали, погружённой в растворитель. Посмотрим на B-факторы отдельных атомов этого остатка.
max_b_res[0].getNames(), max_b_res[0].getBetas()
(array(['N', 'CA', 'C', 'O', 'CB'], dtype='<U6'), array([103.82, 104.62, 113.84, 136.62, 99.29]))
Можно заметить, что большие B-факторы характерны для атомов остова, особенно для карбонильного кислорода. Возможно, именно поэтому не терминальный серин, а этот аланин имеет наибольший средний B-фактор: доля атомов остова в аланине больше, потому что боковая цепь меньше. Проверим эту гипотезу, посмотрев на B-факторы терминального серина.
pdb["A", 128].getNames(), pdb["A", 128].getBetas()
(array(['N', 'CA', 'C', 'O', 'CB', 'OG'], dtype='<U6'), array([109.41, 95.81, 85.5 , 79.72, 91.31, 84.01]))
Гипотеза не подтвердилась, B-факторы терминального остатка почему-то ниже в целом.
Для наглядности также покажем найденные остатки на одном изображении (рис. 2).
Двадцать второй лейцин повёрнут от своей альфа-спирали в сторону другой, тут можно вспомнить «лейциновые молнии». Видимо, подвижность этого остатка снижена за счёт подобного взаимодействия.
Вычислим положение центра масс белка. Для каждого остатка вычислим расстояние его центра масс от общего центра.
protein_mass_center = pd.calcCenter(pdb, weights=pdb.getMasses())
distances = list()
for residue in pdb.iterResidues():
if "CA" in residue.getNames():
mass_center = pd.calcCenter(residue, weights=residue.getMasses())
distances.append(pd.calcDistance(protein_mass_center, mass_center))
С полученными данными можно визуализировать зависимость между средним B-фактором остатка (эти величины для всех остатков уже вычислены в первом задании) и расстоянием от центра масс белка — см. диаграмму рассеяния ниже.
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
plt.plot(distances, [elem[1] for elem in mean_betas], linestyle="", marker=".")
plt.xlabel("Distance from the protein mass center, Å", fontsize=16)
plt.ylabel("Beta factor", fontsize=16)
plt.show()
По графику видно, что в среднем B-фактор возрастает с удалением от центра масс белка. Кроме того, наблюдается два кластера точек, в которых скорость возрастания отличается. Вероятно, эти кластеры соответствуют C- и N-концевому доменам белка. Наибольший B-фактор наблюдался как раз вблизи C-конца.
Во второй части практикума используются сценарии на языке Python, чтобы смоделировать восстановление положений атомов по данным о гармониках Фурье. Вначале создаём «истинное» распределение электронной плотности, которое в дальнейшем и будем восстанавливать. Вот какие параметры я выбрал:
python3 compile-func.py -g 25,3,3.5+30,3.5,4.5+30,3,9+2,3,10+25,3,15+30,3,16.4+2,3,17.4
Затем вычисляются первые несколько сотен гармоник Фурье для выбранной функции. Оказалось, что для «отличного» восстановления функции достаточно $n_0=37$ первых гармоник. Восстановление по этим гармоникам показано на рис. 3.
Для сравнения приведём восстановление по 36 гармоникам и по избыточному набору из 50 гармоник (рис. 4).
A.
B.
Попробуем восстанавливать функцию по 37 гармоникам, но внесём в данные шум. Попробуем разные варианты: гауссовский шум в $50\,\%$ сигнала, добавленный только к амплитудам, добавленный только к фазам, или $10\,\%$-й шум, добавленный и к фазам, и к амплитудам гармоник. Результаты представлены на рис. 5.
A.
B.
C.
Маленькие пики восстановить не получается ни на одном из графиков. При этом в случаях A и C качество восстановления можно назвать средним, график B же не тянет и на «плохое» — слишком уж велика амплитуда шума.
Вернёмся к незашумлённым данным, однако теперь будем брать не все гармоники подряд, начиная с первой и до какой-то. Попробуем убрать одну гармонику низкого порядка или несколько средних гармоник. Также попытаемся добавить к полному набору из 37 гармоник ещё одну, с номером 47. Результаты представлены на рис. 6.
A.
B.
C.
В случае отсутствия гармоники низкого порядка появилось небольшое низкочастотное искажение. Отсутствие трёх промежуточных гармоник по влиянию похоже на присутствие шума, но, в отличие от шума, появляющееся искажение периодично. Появление одной гармоники высокого порядка не изменило вид графика.
Данные по всем представленным наблюдениям представлены в таблице 1.
| Набор гармоник | Разрешение (Å) | Полнота данных (%) | Шум амплитуды (% от величины $F$) | Шум фазы (% от величины $\varphi$) | Качество восстановления | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| Полный набор гармоник | ||||||||||
| 0–37 | $0{,}81$ | 100 | 0 | 0 | Отличное | |||||
| 0–36 | $0{,}83$ | 100 | 0 | 0 | Хорошее | |||||
| 0–50 | $0{,}60$ | 100 | 0 | 0 | Отличное | |||||
| 0–37 | $0{,}81$ | 100 | 50 | 0 | Среднее | |||||
| 0–37 | $0{,}81$ | 100 | 0 | 50 | Отсутствует | |||||
| 0–37 | $0{,}81$ | 100 | 10 | 10 | Среднее | |||||
| Неполный набор гармоник | ||||||||||
| 0, 2–37 | $0{,}81$ | $97{,}3$ | 0 | 0 | Отличное | |||||
| 0–15, 17–20, 22–25, 27–37 | $0{,}81$ | $91{,}9$ | 0 | 0 | Среднее | |||||
| 0–37, 47 | $0{,}81$ | 100 | 0 | 0 | Отличное |
В разделе с полным набором гармоник разрешение было оценено по гармонике с наибольшим номером. В разделе с неполным набором гармоник в качестве старшей была взята гармоника 37.
Можно заметить, что при полноте около $92\,\%$ качество восстановления структуры всё ещё приемлемое, если при полном наборе оно было отличным. Поэтому я бы предложил задать порог полноты, скажем, $90\,\%$, и оценивать разрешение по самой старшей гармонике, которую удалось измерить и для которой полнота превосходит порог.
В упражнениях мне удалось продемонстрировать, что как зашумлённые данные, так и неполнота измерения гармоник могут вносить погрешности при определении структуры. Хорошо иметь полноту хотя бы $90\,\%$. Кроме того, можно обратить внимание на важность решения фазовой проблемы: при внесении одинакового $50\,\%$-го шума в данные об амплитудах и о фазах искажение фаз делает анализ данных полностью невозможным, в то время как зашумление амплитуд лишь снижает качество до «среднего» уровня.