Импорт библиотек

Испорт необходимых библиотек и функций

In [24]:
import matplotlib.pyplot as plt
import numpy as np
import prody as pr
import seaborn as sns
import pandas as pd
from scipy.stats import pearsonr
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score

Импорт pdf-структуры, получение B-факторов и центров масс

Импорт структуры 7axc (hPXR-LBD в комплексе с ферритином), получение объекта (белка), готового к итерированию.

In [2]:
structure = pr.parsePDB('7axc')
selection = structure.select('protein')
hv = selection.getHierView()
hv

@> PDB file is found in working directory (7axc.pdb.gz).
@> 2242 atoms and 1 coordinate set(s) were parsed in 0.03s.
Out[2]:
<HierView: Selection 'protein' (1 chains, 264 residues)>

Получение pd-dataframe, содержащего средние значения В-факторов для атомов, а также координаты их центров масс.

In [3]:
B_factors = np.zeros((264, 6))
for i, residue in enumerate(hv.iterResidues()):
    B_factors[i][1] = np.mean(residue.getBetas())
    B_factors[i][0] = residue.getResnum()
    string = "resnum "
    string+= str(residue.getResnum())
    group = structure.select(string)
    B_factors[i][2], B_factors[i][3], B_factors[i][4] = pr.calcCenter(group, weights=group.getMasses())
In [4]:
data = pd.DataFrame(B_factors, columns = ["Number", "B-factor", "x_mass", "y_mass", "z_mass", "to_center"])

Поиск остатка, имеющего максимальный средний В-фактор.

In [5]:
data[data["B-factor"] == data["B-factor"].max()]
Out[5]:
Number B-factor x_mass y_mass z_mass to_center
67 230.0 93.3925 26.29246 86.534272 -6.979297 0.0

Поиск остатка, имеющего минимальный средний В-фактор.

In [6]:
data[data["B-factor"] == data["B-factor"].min()]
Out[6]:
Number B-factor x_mass y_mass z_mass to_center
169 340.0 24.345833 3.634301 59.630554 2.529927 0.0

Подсчет центра масс для всей белковой молекулы.

In [7]:
whole_protein_mass = pr.calcCenter(selection, weights=selection.getMasses())
whole_protein_mass
Out[7]:
array([ 4.41446713, 69.47872418,  2.36382885])

Подсчет расстояния

Подсчитаем расстояние от центра масс каждой аминокислоты до центра масс всей белковой молекулы.

In [8]:
for index, row in data.iterrows():
    row["to_center"] = np.sqrt((row["x_mass"] - whole_protein_mass[0])**2 + (row["y_mass"] - whole_protein_mass[1])**2 + (row["z_mass"] - whole_protein_mass[2])**2)
In [9]:
data.head()
Out[9]:
Number B-factor x_mass y_mass z_mass to_center
0 142.0 57.835000 -6.077637 44.922577 12.808048 28.673512
1 143.0 48.076250 -7.249522 47.951080 11.511819 26.137594
2 144.0 51.575714 -10.561891 49.195055 14.039614 27.785653
3 145.0 55.744444 -14.516207 49.532225 10.814803 28.768945
4 146.0 48.913333 -12.882793 53.520639 13.202823 25.910220

Построение графиков, поиск модели распределения

Построим scatter plot, который будет отражать распределение зависимости величины В-фактора от расстояния от центра масс белка.

In [10]:
plt.figure(figsize=(10, 8), dpi=80)
plt.suptitle('График зависимости В-фактора от расстояния до центра масс белка')
ax = []
plt.xlabel("Расстояние, ангстремы")
plt.ylabel("В-фактор")
plt.xlim(0, 35)
plt.ylim(0, 90)
sns.scatterplot(x = data["to_center"], y = data["B-factor"])
Out[10]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f5f9718ce50>

На графике видно, что, конечно, с ростом расстояния от центра увеличивается величина В-фактора, причем, по первому взгляду, зависимость наблюдается практически линейная с присутствием некоторого разброса, который может объясняться иными факторами (например, положение центра масс каждого остатка, а не только расстояние между ним и общим центром масс).
Посмотрим на величину пирсоновского коэффициента корреляции.

In [11]:
corr, _ = pearsonr(data["to_center"], data["B-factor"])
print('Pearsons correlation: %.3f' % corr)

Pearsons correlation: 0.658

Как и предполагалось, мы наблюдаем достаточно высокий коэффициент корреляции между данными параметрами.

Построение линейных моделей

Высокий коэффициент корреляции, а также оющий вид графика позволяет нам предположить, что линейные модели могут с достаточно высокой точностью оценивать этот параметр. Построим линейную регрессию, которая будет отражать зависимость среднего В-фактора от расстояния до центра белка.

In [12]:
linear_regression = LinearRegression()
X_train, X_test, y_train, y_test = train_test_split(data["to_center"].values, data["B-factor"].values, test_size=0.33)
linear_regression.fit(X_train.reshape(-1, 1), y_train)
print("Coef: ", round(linear_regression.coef_[0], 4), "; Intercept: ", round(linear_regression.intercept_, 4))

Coef:  1.6845 ; Intercept:  13.1496
In [13]:
i = 0
x_reg = [i for i in range(34)]
y_reg = [k*1.6658 + 12.9251 for k in x_reg]
plt.figure(figsize=(10, 8), dpi=80)
plt.suptitle('График зависимости В-фактора от расстояния до центра масс белка с регрессионной прямой')
ax = []
plt.xlabel("Расстояние, ангстремы", fontsize = 13)
plt.ylabel("В-фактор", fontsize = 13)
plt.xlim(0, 35)
plt.ylim(0, 90)
sns.scatterplot(x = data["to_center"], y = data["B-factor"], label = "True points")
sns.lineplot(x = x_reg, y = y_reg, color = 'red', label = "Regression line")
plt.legend(fontsize = 13)
Out[13]:
<matplotlib.legend.Legend at 0x7f5f96e46f10>
In [14]:
linear_regression.score(X_test.reshape(-1, 1), y_test)
Out[14]:
0.3525256232484003

Мы наблюдает достаточно невысокое качество полученной модели, так как входных данных, как мне кажется, не достаточно для более точной модели. Построим линейную модель, которая использует большее количество параметров, а именно, координаты центра масс для каждого остатка.

In [25]:
X_many = data[["x_mass", "y_mass", "z_mass", "to_center"]]
y_many = data["B-factor"]
X_train, X_test, y_train, y_test = train_test_split(X_many, y_many, test_size=0.33)
In [26]:
linear_regression = LinearRegression()
linear_regression.fit(X_train, y_train)
linear_regression.score(X_test, y_test)
Out[26]:
0.6738551783864699
In [27]:
linear_regression.coef_
Out[27]:
array([0.58178854, 0.23341988, 0.12298193, 1.49910339])
In [28]:
y_pred = linear_regression.predict(X_test)
plt.figure(figsize=(10, 8), dpi=80)
plt.suptitle('График зависимости В-фактора от расстояния до центра масс белка с регрессионной прямой')
ax = []
plt.xlabel("Расстояние, ангстремы", fontsize = 13)
plt.ylabel("В-фактор", fontsize = 13)
plt.xlim(0, 35)
plt.ylim(0, 90)
sns.scatterplot(x = X_test["to_center"], y = y_test, label = "True points")
sns.scatterplot(x = X_test["to_center"], y = y_pred, label = "Predicted points", color = 'red')
Out[28]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f5f8a807350>

В данной модели мы видим на порядок более высокое качество по сравнению с обычной линейной регрессией, такую модель, возможно, можно использовать для предсказания B-факторов других структур.

Построение нелинейных моделей

В построенной линейной модели все еще есть недостатки: она немного занижает значения для В-факторов для близких к центру масс остатков и для далеких. Также можно отметить, что модель не до конца способна охватить разброс в значениях. Попробуем применить в данном случае нелинейную модель -- рандомный лес.

In [17]:
out = {}
for i in range(1,51):
    regr = RandomForestRegressor(n_estimators=i)
    regr.fit(X_train, y_train)
    pred = regr.predict(X_test)
    out[i] = r2_score(y_test, pred) 
out
Out[17]:
{1: 0.6646441404993266,
 2: 0.7410800106694281,
 3: 0.7600907624047271,
 4: 0.7513012488493311,
 5: 0.7974226125113055,
 6: 0.7791801569799592,
 7: 0.7698337697426869,
 8: 0.7782586574292827,
 9: 0.7936702157407675,
 10: 0.8013017320361556,
 11: 0.7834071746627708,
 12: 0.8249058371548263,
 13: 0.7863346392046178,
 14: 0.7970125609452576,
 15: 0.7971702331623862,
 16: 0.8040783060461811,
 17: 0.777919339539186,
 18: 0.8139437898004227,
 19: 0.779496070404871,
 20: 0.8217979316754725,
 21: 0.7979367168989531,
 22: 0.8300567009501576,
 23: 0.8211913019176172,
 24: 0.79641954378982,
 25: 0.8145109363987233,
 26: 0.8054687332035728,
 27: 0.7895503449656054,
 28: 0.7973407097156635,
 29: 0.8149421067553255,
 30: 0.8206241533095516,
 31: 0.7917443523215173,
 32: 0.8088165532513496,
 33: 0.8141411473010292,
 34: 0.8013516221038114,
 35: 0.8102698284639812,
 36: 0.8272474817979414,
 37: 0.8013562766224192,
 38: 0.8233750437801725,
 39: 0.8182999617581039,
 40: 0.8146867045039261,
 41: 0.7936427453384935,
 42: 0.822389554853144,
 43: 0.8120561137001961,
 44: 0.8182879190015774,
 45: 0.801445506276409,
 46: 0.8027373032152794,
 47: 0.823051489645411,
 48: 0.8151023015231269,
 49: 0.8230986804283724,
 50: 0.8141918718205701}

Можно сказать, что начиная с 6 деревьев качество модели перестает очевидно расти. Остановимся на этом количестве и построим такую модель, оценим качество ее работы численно и графически.

In [41]:
regr = RandomForestRegressor(n_estimators=6)
regr.fit(X_train, y_train)
pred = regr.predict(X_test)
regr.score(X_test, y_test)
Out[41]:
0.8489569066723275
In [42]:
plt.figure(figsize=(10, 8), dpi=80)
plt.suptitle('График зависимости В-фактора от расстояния до центра масс белка с регрессионной прямой')
ax = []
plt.xlabel("Расстояние, ангстремы", fontsize = 13)
plt.ylabel("В-фактор", fontsize = 13)
plt.xlim(0, 35)
plt.ylim(0, 90)
sns.scatterplot(x = X_test["to_center"], y = y_test, label = "True points")
sns.scatterplot(x = X_test["to_center"], y = pred, label = "Predicted points", color = 'red')
Out[42]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f5f89c49b10>

Данная модель работает немного лучше, чем линейная со множественными параметрами; способна не занижать точки на краях и показывает чуть больший разброс. И хотя тест качества проводился на тестовой выборке, это все равно может отражать переобучение под конкретную модель.

Проверка моделей

Проверим данные модели на других белковых структурах -- например, на лизоциме.

In [43]:
structure = pr.parsePDB('7PNH')
selection = structure.select('protein')
hv = selection.getHierView()
hv

@> Connecting wwPDB FTP server RCSB PDB (USA).
@> pdb7pnh.ent.gz download failed. 7pnh does not exist on ftp.wwpdb.org.
@> PDB download via FTP completed (0 downloaded, 1 failed).
@> Downloading PDB files via FTP failed, trying HTTP.
@> WARNING 7pnh download failed ('https://files.rcsb.org/download/7PNH.pdb.gz' could not be opened for reading, invalid URL or no internet connection).
@> PDB download via HTTP completed (0 downloaded, 1 failed).
@> Trying to parse mmCIF file instead
@> Connecting wwPDB FTP server RCSB PDB (USA).
@> 7pnh downloaded (7pnh.cif)
@> PDB download via FTP completed (1 downloaded, 0 failed).
@> 1227 atoms and 1 coordinate set(s) were parsed in 0.08s.
Out[43]:
<HierView: Selection 'protein' (1 segments, 1 chains, 129 residues)>
In [44]:
B_factors = np.zeros((129, 6))
for i, residue in enumerate(hv.iterResidues()):
    B_factors[i][1] = np.mean(residue.getBetas())
    B_factors[i][0] = residue.getResnum()
    string = "resnum "
    string+= str(residue.getResnum())
    group = structure.select(string)
    B_factors[i][2], B_factors[i][3], B_factors[i][4] = pr.calcCenter(group, weights=group.getMasses())
data = pd.DataFrame(B_factors, columns = ["Number", "B-factor", "x_mass", "y_mass", "z_mass", "to_center"])
whole_protein_mass = pr.calcCenter(selection, weights=selection.getMasses())
for index, row in data.iterrows():
    row["to_center"] = np.sqrt((row["x_mass"] - whole_protein_mass[0])**2 + (row["y_mass"] - whole_protein_mass[1])**2 + (row["z_mass"] - whole_protein_mass[2])**2)
data.head()
Out[44]:
Number B-factor x_mass y_mass z_mass to_center
0 1.0 20.115333 0.582051 -9.995531 -9.606926 14.127372
1 2.0 17.589429 -1.611204 -13.124220 -12.057850 14.509819
2 3.0 16.756000 2.923835 -13.908956 -9.981033 12.181549
3 4.0 17.194750 3.925419 -16.574243 -13.228098 14.389623
4 5.0 18.762727 5.908467 -21.807352 -14.091008 15.480201
In [45]:
plt.figure(figsize=(10, 8), dpi=80)
plt.suptitle('График зависимости В-фактора от расстояния до центра масс белка')
ax = []
plt.xlabel("Расстояние, ангстремы")
plt.ylabel("В-фактор")
plt.xlim(0, 30)
plt.ylim(0, 60)
sns.scatterplot(x = data["to_center"], y = data["B-factor"])
Out[45]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f5f89beb650>

Запустим линейную модель.

In [46]:
X = data[["x_mass", "y_mass", "z_mass", "to_center"]]
y = data["B-factor"]
In [47]:
y_pred = linear_regression.predict(X)
plt.figure(figsize=(10, 8), dpi=80)
plt.suptitle('График зависимости В-фактора от расстояния до центра масс белка с регрессионной прямой')
ax = []
plt.xlabel("Расстояние, ангстремы", fontsize = 13)
plt.ylabel("В-фактор", fontsize = 13)
plt.xlim(0, 35)
plt.ylim(-5, 90)
sns.scatterplot(x = X["to_center"], y = y, label = "True points")
sns.scatterplot(x = X["to_center"], y = y_pred, label = "Predicted points", color = 'red')
Out[47]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f5f89b7e850>

Запустим случайный лес.

In [48]:
y_pred = regr.predict(X)
plt.figure(figsize=(10, 8), dpi=80)
plt.suptitle('График зависимости В-фактора от расстояния до центра масс белка с регрессионной прямой')
ax = []
plt.xlabel("Расстояние, ангстремы", fontsize = 13)
plt.ylabel("В-фактор", fontsize = 13)
plt.xlim(0, 35)
plt.ylim(-5, 90)
sns.scatterplot(x = X["to_center"], y = y, label = "True points")
sns.scatterplot(x = X["to_center"], y = y_pred, label = "Predicted points", color = 'red')
Out[48]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f5f89aa7d50>

Как мы видим, в данном случае и линейная модель, и случайный лес, работают достаточно плохо, хотя множественная линейная модель имеет визуально более высокое качество. Наверно, для описания такой взаимосвязи можно использовать линейные модели со множеством параметров, однако, обучать их стоит на совокупности данных большого количества белковых структур. Более сложные модели также могут лучше справиться с описанием данной зависимости.