Задание 1¶
В этом задании нужно было с помощью модуля Prody найти в структуре 3CUZ остаток с наибольшим и наименьшим B-фактором. Поскольку B-фактор это понятие для атомов, B-фактор для остатка определяется как среднее по B-факторам составляющих его атомов.
Найдя нужные остатки, нужно было проанализировать их расположение в структуре белка и сделать выводы.
import prody as prd
import numpy as np
protein = prd.parsePDB("3CUZ")
betas = []
for residue in protein.iterResidues():
mean_beta = np.mean(residue.getBetas())
betas.append([residue, mean_beta])
Найдем остаток с наименьшим B-фактором
sorted(betas, key = lambda x: x[1])[0]
Наименьший B-фактор имеет остаток изолейцина 261: 6.13625
Остаток ILE 261 показан в виде фиолетовых сфер. Как можно видеть, он находится со стороны альфа-спирали, обращенной внутрь белковой глобулы. Посмотрим подробнее изнутри глобулы.
Наш остаток покрашен фиолетовым. Его окружение показано в виде палочек, атомы остова скрыты. Красными сферами показаны молекулы воды.
Поскольку гидрофобный эффект является главной движущей силой сворачивания белка и поддержания его структуры, такие гидрофобные кластеры устойчивы и мало варьируют от молекулы к молекуле в кристалле. Дополнительно, наш остаток находится в составе вторичной структуры, что структурирует по крайней мере атомы остова.
То же самое, но атомы показаны в виде сфер, видно, радикал остатка плотно запакован и потому стабилен.
Для наглядности уберем часть атомов, с добавленными водородами становится еще понятнее, что места для "телодвижений" у нашего остатка немного в виду его плотной упаковки.
Найдем теперь остаток с наибольшим B-фактором.
sorted(betas, key = lambda x: x[1], reverse = True)[0]
Упс! Нашли воду. Давайте уточним, что ищем по аминокислотам.
res_betas = [res for res in betas if "CA" in res[0].getNames()]
sorted(res_betas, key = lambda x: x[1], reverse = True)[0]
Готово.
Наибольший B-фактор имеет остаток глутамата 497: 23.6580
Как можно видеть, остаток со всех сторон окружен растворителем и настолько подвижен, что не получилось вписать атомы радикала дальше Сb в электронную плотность.
В принципе, нахождение остатка на поверхности белка не обязательно должно приводить к его большой подвижности в случае кристаллических структур, поскольку этот остаток может стабилизироваться кристаллическими контактами. Проверим, есть ли они в случае нашего остатка (хотя по тому, что от радикала нам известен только один атом, уже все понятно).
Сгенерировал соседей по кристаллу, остатки в радиусе 10 ангстрем показал палочками. Glu497 показан в виде сфер. Как можно видеть, вокруг нашего остатка в кристалле никаких контактов нет и смотрит он в растворитель. Между делом можно заметить два немного подозрительных кислорода выше нашего остатка, рядом с радикалом пролина. Что если это атомы нашего глутамата, которые приняли за воду?
Подгрузили плотность, показали на уровне подрезки 0.5 окружение (carve = 8) нашего глутамата. Где у него радикал по плотности правда не понятно. А расстояние между двумя подозирительными кислородами составляет 2.7 ангстрем, что можно интерпретировать как водородную связь между этими кислородами. Загадочно, правда, почему между кислородами так много плотности. Скорее всего, они просто сильно колеблются.
Сравнивая расстояние 2.7 между теми кислородами и расстоянием 2.2 между кислородами карбоксильной группы нормально разрешенного глутамата 452, понимаем, что те две воды это, скорее всего, просто две воды, слишком уж большое расстояние между ними, чтобы быть карбоксильной группой.
Посмотрим как внутри нашего обделенного атомами глутамата распределен B-фактор.
Рядом с атомами подписаны значения B-факторов. Атомы раскрашены по значениям их B-факторов от синего (наименьший) к красному(наибольший). B-факторы не сильно различаются.
Как можно ожидать, Сb атом имеет большую подвижность, чем атомы остова, правда, неожиданно, самым подвижным атомом всего наблюдаемого куска остатка оказался карбонильный кислород. Может, этот кислород не формирует водородной связи, хоть и находится в альфа-спирали?
Фрагмент альфа-спирали с нашим остатком (его раскраска сохранена с прошлой картинки). Отмеченные расстояния измерены между тяжелыми атомами. Протон серина был повернут вручную в сторону соседнего карбонила, так как иначе между кислородом серина и этим карбонильным кислородом было бы слишком малое расстояние.
Видим, что водородная связь у нашего подвижного карбонильного кислорода не оптимальная по углу, и по расстоянию хуже, чем по крайней мере еще одна видимая водородная связь длины 2.9.
Вообще, кажется, что атомы радикала должны быть подвижнее, чем атомы остова, особенно, если остов вовлечен во вторичную структуру. Этот же вывод мы сделали в прошлом практикуме в задании 2. За исключением странного карбонильного кислорода, мы наблюдаем этот тренд, атомы дальше Cb по радикалу вообще настолько подвижны, что их не видно.
Наверное, даже для петли атомы остова будут менее подвижны, чем атомы радикала (если это не фенилаланин какой-нибудь или что-то подобное, где атомы по понятным причинам зафиксированы относительно друг друга), потому что остов представляет из себя систему плоскостей, соответствующих пептидным связям, не все взаимные ориентации которых возможны (карта Рамачандрана), а радикал, который смотрит в раствор в движении не особо ограничен.
Посмотрим еще раз на Ile261 и Glu497 и сделаем выводы:
Ile261 покрашен голубым, Glu497 красным.
Ile261 самый неподвижный, он погружен в глобулу белка, радикал плотно упакован в гидрофобное ядро и неподвижен. Glu497 самый подвижный, он смотрит в растворитель и не контактирует с другими молекулами в кристалле, подвижность его радикала ничего не ограничивает.
Задание 2¶
Целью этого задания являлось построить зависимость B-факторов остатков от расстояния между их центром масс и центром массы белковой глобулы, прокомментировать наблюдения.
protein = prd.parsePDB("3CUZ")
protein_atoms = protein.select("protein") # Only aminoacids
protein_center = prd.calcCenter(protein_atoms, weights=protein_atoms.getMasses())
res_betas_task2 = []
distances = []
for residue in protein.iterResidues():
if "CA" in residue.getNames(): # if residue is an aminoacid
mean_beta = np.mean(residue.getBetas())
res_center = prd.calcCenter(residue, weights=residue.getMasses())
distance = prd.calcDistance(protein_center, res_center)
res_betas_task2.append(mean_beta)
distances.append(distance)
Визуализируем результаты
import seaborn as sns
plot = sns.jointplot(x=distances, y=res_betas_task2, marker = ".")
plot.set_axis_labels('Distances from centers', 'B-factor')
Из графика видим следующее:
- Аминокислот с малым расстоянием от центра масс белка мало и постепенно становится больше. Это связано с тем фактом, что аминокислотные остатки -- не материальные точки, которые могут располагаться как угодно, а имеют вполне определенные размеры. Много их упаковать вокруг заданной точки не получится;
- Явно выделяются остатки, сильно удаленные от центра масс и имеющие большие B-факторы. Возможно, эти остатки смотрят в раствор и, наверное, принадлежат N- и C- концам белка.
Проверим идею, что у концов большие B-факторы. Отсортируем по B-факторам список остатков.
terminality_check_res = []
prot_length = 0
for residue in protein.iterResidues():
if "CA" in residue.getNames(): # if residue is an aminoacid
mean_beta = np.mean(residue.getBetas())
res_center = prd.calcCenter(residue, weights=residue.getMasses())
distance = prd.calcDistance(protein_center, res_center)
terminality_check_res.append([residue, mean_beta, distance])
prot_length += 1 # лень было искать как из prody достать длину белка
print(f"Total protein length is {prot_length}")
terminality_check_res.sort(key = lambda x: x[1], reverse = True)
for elem in terminality_check_res[:15]:
print(elem)
Кажется, мы не правы насчет концов. В белке всего 529 аминокислот, и совсем близко к началу или к концу белка ничего нет. Наверное, мы не правы в случае нашего белка, потому что концы прижаты к белку, а не торчат совсем произвольно. На N-конце паймол рисует маленький кусочек бета-листа, который, вероятно, его стабилизирует, а на С-конце у нас альфа-спираль.
Тогда остается думать что самые подвижные остатки это экспонированные в раствор не обязательно на концах.
Про доступность остатка растворителю у нас будет дальше в курсе лекций, поэтому я пока не могу проверить гипотезу, что самые подвижные = самые экспонированные. Также интересно было бы покрасить точки в зависимости от вторичной структуры остатка, но лекции про вторичную структуру у нас тоже пока не было :D
Сделаем предыдущий скаттерплот гладким и красивым.
plot = sns.jointplot(x=distances, y=res_betas_task2, kind="kde")
plot.set_axis_labels('Distances from centers', 'B-factor')
Видим, что больше всего остатков на расстоянии 20-25 ангстрем от центра белка. Можно попробовать объяснить это тем, что, как я писал ранее, остатки это не точки и не могут произвольно быть разбросанными вокруг заданной материальной точки (центра глобулы). Кажется разумным, что чем дальше мы удаляемся от центра белка (можно представить, что мы находимся при этом на поверхности сферы с центром в центре глобулы), тем больше у нас становится места ($S_{сферы} = 4\pi r^2$) для упаковки центров остатков. То есть, по мере удаления от центра остатков становится больше, потому что их можно больше упаковать, а после какого-то расстояния их становится меньше просто потому что остатки закончились.
По поводу B-факторов (смотрим на распределение справа) можно сказать, что сильно неподвижных остатков мало (белок -- динамичная структура) и в среднем все остатки довольно подвижны. Аномально подвижных остатков тоже немного, но больше, чем малоподвижных (хвост длиннее на соответствующем распределении)
Двумя абзацами выше мы придумали гипотезу:
На небольшом расстоянии от центра глобулы остатки пакуются более-менее плотно и их число пропорционально площади сферы с данным радиусом.
В нашем случае плотная упаковка заканчивается примерно на 22.5 Å (максимум верхнего распределения на графике выше).
Иными словами, $ N_{остатков}^{0..22.5Å} = N_{остатков}^{0..22.5Å}\int_{0}^{22.5}Ar^2 \,dr$, где A - нормировка.
Нормировка определяется тривиально: надо посчитать интеграл и взять обратное число: $ A = (\int_{0}^{22.5}r^2\,dr))^{-1} = \frac{3}{22.5^3 - 0^3}$
def packing_pseudodensity(r, normalization):
return r**2 * normalization
normalization = 3 / (22.5**3 - 0)
# select distances less than or equal to 22.5
close_res = [d for d in distances if d <= 22.5]
N_res = len(close_res)
x_points = np.linspace(0, 22.5, 100)
theoretical_density = np.apply_along_axis(lambda r : packing_pseudodensity(r, normalization),
arr=x_points, axis=0)
theoretical_counts = N_res * theoretical_density
Дальше я встретился с проблемой: если строить гистограмму как отображение числа остатков на данном расстоянии от центра глобулы, то возникает произвол с выбором числа бинов у гистограммы. Варьируя число бинов, можно подогнать, по-моему, любое распределение под кривую.
sns.distplot(close_res, kde=False, bins = 4)
sns.lineplot(x_points, theoretical_counts)
Взяли мало бинов (4).
sns.distplot(close_res, kde=False, bins = 17)
sns.lineplot(x_points, theoretical_counts)
Постарались подогнать бины, чтобы было красиво. Тут 17 бинов.
sns.distplot(close_res, kde=False, bins = 50)
sns.lineplot(x_points, theoretical_counts)
Взяли слишком большое число бинов (50).
Возможно, компромиссным вариантом является построение ядерной оценки плотности, но у нее, опять же, тоже есть параметр (bandwidth), который можно двигать, чтобы подогнать получаемую плотность под желаемую кривую. Однако, тут я параметр не двигал, он определяется независимо от меня, исходя из данных.
sns.distplot(close_res, kde=True)
sns.lineplot(x_points, theoretical_density)
Вроде бы и на плотность похоже и на гистограмму немного. Может, и правда, остатки более-менее плотно пакуются. Понятно, что для фибриллярных или других вытянутых белков это точно неверно, поскольку для них такое сферическое приближение белка неуместно. Может быть их можно фитить эллипсоидом.
Какой глубокий смысл того, что я тут сделал? Не знаю. Было интересно, но получилось ли у меня или я совсем чушь какую-то написал, я не знаю.
Задание 3¶
В это задании мы моделируем одномерный РСА эксперимент: возьмем модельную электронную плотность (ЭП) и разложим ее в ряд Фурье. Далее будем смотреть, как изменение амплитуд и фаз в ряде Фурье влияет на восставление нашей исходной ЭП.
Будем моделировать ЭП молекул $SCN^-$ и $N_2$.
Параметры compile_func.py: -g 16.55,4.5,5+6.1,3,6.6+7.35,3,7.8+7,3,15+7,3,16.1
Нецелые "числа электронов" у атомов в случае $SCN^-$ мотивированы тем, что один электрон делокализован по молекуле. Добавка электронов к каждому атому, а также длины связей были списаны с каких-то картинок в интернете.
import pandas as pd
def parse_file(filename, label):
data = pd.read_table(filename, names = ["x","y"], skiprows=1)
data["label"] = label
return data
true_density = parse_file("true_density.txt", "true density")
sns.lineplot(x = true_density["x"], y = true_density["y"])
Вот они, красавцы. Видно, что азот в тиоцианате имеет чуть большую ЭП, чем азот в молекулярном азоте.
А вот теперь мы над этой плотностью будем издеваться по-всякому.
Полный набор гармоник¶
Посмотрим, как число гармоник влияет на качество восстанавливаемой плотности и сколько гармоник хватает, чтобы все увидеть. Рассмотрим также влияние шума по амплитудам и по фазам.
density_0_10 = parse_file("ft-0-10.txt", "0-10")
sns.lineplot(x =density_0_10["x"], y = density_0_10["y"])
Первые 10 гармоник. Плохое восстановление. Видно, где расположены молекулы, но не более того.
density_0_20 = parse_file("ft-0-20.txt", "0-20")
sns.lineplot(x =density_0_20["x"], y = density_0_20["y"])
Первые 20 гармоник. Что-то среднее между хорошим и отличным восстановлением: в $SCN^-$ видно все атомы, в молекулярном азоте видно очень широкий пик, в который вполне можно положить два атома. Другое дело, что большой пик может принадлежать какому-то сильно подвижному атому с большим числом электронов.
Если подходить к задаче с позиций реального РСА, то там мы априори знаем, из каких аминокислотных остатков состоит белок. Поэтому, если у нас есть знание, что где-то должен быть молекулярный азот, то все атомы мы тут различим.
density_0_30 = parse_file("ft-0-30.txt", "0-30")
sns.lineplot(x =density_0_30["x"], y = density_0_30["y"])
Гармоники 0-30. Отличное восстановление, только по высоте пиков у азотов проврались: все три пика примерно одинаковой высоты. При этом пики в молекулярном азоте разной высоты, что не соответствует реальности.
Теперь про шум. Согласно заданию, работаем с гармониками от 0 до 30, которые у нас отлично восстанавливают исходную плотность. Здесь и далее я буду пользоваться словом "фон", имея в виду высокочастотные колебания в районе базовой линии, которые происходят из-за того, что мы работаем не с всеми гармониками (которых, формально, бесконечность), а только с первыми тридцатью.
def plot_two(f1_df, f2_df):
data = f1_df.append(f2_df)
sns.lineplot(x = data["x"], y = data["y"], hue = data["label"])
F5phi0 = parse_file("F5_phi0.txt", "5% amplitude noise")
plot_two(density_0_30, F5phi0)
5% шум по амплитудам. Отличное восстановление, сложно отличить от исходной функции
F0phi5 = parse_file("F0_phi5.txt", "5% phase noise")
plot_two(density_0_30, F0phi5)
5% шум по фазам. Отличное восстановление, изменения более существенные, чем при 5% шуме по амплитудам.
F25phi0 = parse_file("F25_phi0.txt", "25% amplitude noise")
plot_two(density_0_30, F25phi0)
25% шум по амплитудам. Отличное восстановление, но сильный фоновый сигнал.
F0phi25 = parse_file("F0_phi25.txt", "25% phase noise")
plot_two(density_0_30, F0phi25)
25% шум по фазам. На грани отличного и хорошего восстановления, изменения намного более существенные, чем при 25% шуме по амплитудам. Фон почти сопоставим с сигналом, часть пиков в сигнале просела.
F25phi5 = parse_file("F25_phi5.txt", "25% amplitude noise, 5% phase noise")
plot_two(density_0_30, F25phi5)
25% шум по амплитудам, 5% шум по фазам. Отличное восстановление.
F5phi25 = parse_file("F5_phi25.txt", "5% amplitude noise, 25% phase noise")
plot_two(density_0_30, F5phi25)
5% шум по амплитудам, 25% шум по фазам. Отличное восстановление. (Ситуация чуть лучше просто 25% шума по фазам, потому что для создания шума используются случайные числа. Тут нам больше повезло).
Можно подтвердить тезис с лекций, что фазы нам важнее, чем амплитуды.
Шум по фазам приводит к большему уровню фонового сигнала, чем сопоставимый шум по амплитудам. Если при шуме по амплитудам у нас пики восстановленной плотности остаются над пиками исходной плотности и мы максимум перепутаем какие-то элементы, то при шуме по фазам у нас могут пропасть пики в одном месте или появиться пики в другом месте.
Понятно, что убрать пик или прибавить новый может и шум по амплитуде, но для этого величина шума должна быть сопоставима с величиной сигнала. Кажется, что таких ситуаций немного.
Таким образом, амплитуды определяют, какие по величине мы получаем пики, а фазы определяют то, где эти пики вообще будут находиться. Для восстановления нам нужны и фазы и амплитуды, "важность" фаз или амплитуд, о которой мы говорим здесь, это скорее то, к неточностям в какой части данных более чувствительно восстановление ЭП.
Неполный набор гармоник¶
no_0_harm = parse_file("1-30_harmonics_func.txt", "Without 0 harmonic")
plot_two(density_0_30, no_0_harm)
Убрали нулевую гармонику. Так как нулевая гармоника - это просто константа, график съехал вниз. Восстановление отличное.
no_0_and_1_harm = parse_file("2-30_harmonics_func.txt", "Without 0 and 1 harmonics")
plot_two(density_0_30, no_0_and_1_harm)
Убрали нулевую и первую гармонику. Первая гармоника - это уже синусоида, но такого большого периода, что на детализированность восстановления плотности ее отсутствие никак не повлияло. Отличное восстановление.
no_16_to_18_harm = parse_file("no_16-18_harmonics_func.txt", "Without 16-18 harmonics")
plot_two(density_0_30, no_16_to_18_harm)
Убрали гармоники 16, 17 и 18. Сразу увеличился фон, пострадал пик углерода в $SCN^-$. Несмотря на это, восстановление все еще отличное.
plus_40_harm = parse_file("plus_40_harmonic_func.txt", "Plus harm №40")
plot_two(density_0_30, plus_40_harm)
Добавили гармонику номер 40. Она содержит очень мелкие детали и почти никак не повлияла на общую картину. Отличное восстановление.
Таким образом, отсутствие первых гармоник безусловно, влияет на абсолютные значения восстанавливаемой плотности, но при этом различить детали не составляет трудности.
Отсутствие гармоник в середине набора ведет к увеличению фонового сигнала и может затрагивать высоту пиков сигнала.
Прибавляя новые гармоники сверх уже достаточных для отличного восстановления, мы незначительно улучшаем восстанавливаем картину.
Таблица по заданию 3.¶
| Набор гармоник | Разрешение (Å) | Полнота (%) | Шум амплитуды (% от величины F) | Шум фазы (% от величины phi) | Качество восстановления (отличное, хорошее, среднее, плохое) |
|---|---|---|---|---|---|
| 0-10 | 3 | 100% | 0 | 0 | Плохое |
| 0-20 | 1.5 | 100% | 0 | 0 | Между хорошо и отлично |
| 0-30 | 1 | 100% | 0 | 0 | Отличное |
| 0-30 | 1 | 100% | 5 | 0 | Отличное |
| 0-30 | 1 | 100% | 0 | 5 | Отличное |
| 0-30 | 1 | 100% | 25 | 0 | Отличное |
| 0-30 | 1 | 100% | 0 | 25 | Между хорошо и отлично |
| 0-30 | 1 | 100% | 25 | 5 | Отличное |
| 0-30 | 1 | 100% | 5 | 25 | Отличное (повезло с случайными числами) |
| 1-30 | 1 | 97% | 0 | 0 | Отличное |
| 2-30 | 1 | 94% | 0 | 0 | Отличное |
| 0-15, 19-30 | 1 | 90% | 0 | 0 | Отличное |
| 0-30, 40 | 1* | 100%* | 0 | 0 | Отличное |
* Для диапазона гармоник 0-30 разрешение 1 ангстрем, полнота 100%, для диапазона гармоник 0-40 разрешение 0.75 ангстрем, полнота 78%.