Porównywanie struktur 3D

Wprowadzenie

Porównywanie ma na celu ocenę podobieństw i różnic pomiędzy parą lub w obrębie zbioru struktur 3D
Bioinformatyczne metody skupiają się na ocenie ilościowej np. średnia różnica wynosi X lub podobieństwo fragmentu to N procent
Metody porównywania struktur 3D różnią się między sobą:
- Mogą być globalne (oceniają podobieństwo całych struktur) lub lokalne (znajdują różne/podobne fragmenty)
- Mogą mierzyć podobieństwo (im większa wartość, tym bliższy kształt mają struktury) lub odległość (im większa wartość, tym bardziej struktury się różnią)
- Mogą być uniwersalne lub dedykowane RNA lub białkom (lub nawet białkom jednego typu)
Porównywanie struktur przydaje się gdy:
- Chcemy klasyfikować struktury według ich kształtu 3D
- Modelujemy strukturę 3D by (1) najpierw znać podobieństwa między modelami i (2) wyznaczyć jakość modelowania (np. kiedy znamy już strukturę rozwiązaną eksperymentalnie)
- Poszukujemy interesującego fragmentu

Root Mean Square Deviation (RMSD)

RMSD to miara odległości między dwoma zbiorami atomów o tej samej liczności wyrażona w Å:
RMSD = \sqrt{\dfrac{1}{n}\sum_{i=1}^n \delta_{i}^2}
gdzie:
- n to liczba atomów
- \delta_{i} to odległość między i-tą parą atomów
Bardzo istotne jest, żeby przed wyznaczeniem tej miary dokonać superpozycji struktur:
164.9Å 6.4Å
- Traktujemy współrzędne atomów jako macierze N \times 3 (w przykładzie macierze będą się nazywać A i B):
  A = \begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ & \cdots & \\ x_n & y_n & z_n \end{bmatrix}
- Przesuwamy obie struktury do ich centroidów:
  A = \begin{bmatrix} x_1-\bar{x} & y_1-\bar{y} & z_1-\bar{z} \\ x_2-\bar{x} & y_2-\bar{y} & z_2-\bar{z} \\ & \cdots & \\ x_n-\bar{x} & y_n-\bar{y} & z_n-\bar{z} \end{bmatrix} gdzie \bar{x} to średnia ze współrzędnych x (analogicznie \bar{y} i \bar{z})
- Wyznaczamy macierz kowariancji:
  A^TB
- Wyznaczamy rozkład według wartości osobliwych (ang. Singular Value Decomposition (SVD)):
  A^TB = U\Sigma V^T
- Wyznaczamy znak:
  d = \text{sgn}(\det(A^TB))
- Otrzymujemy macierz rotacji:
  R = V \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & d \end{bmatrix} U^T
- Współrzędne macierzy B mnożymy przez macierz rotacji R:
  B' = B R
- Na podstawie macierzy A i B' wyznaczamy RMSD wg wzoru powyżej
RMSD jest najbardziej znaną i powszechnie używaną miarą odległości, ale nie jest pozbawiona wad:
- Struktury wielodomenowe porównywane w całości mogą mieć dowolnie dużą wartość RMSD, nawet gdy domeny analizowane osobno są parami bardzo podobne
- RMSD jest zależne od wielkości cząsteczki, więc przykładowo informacja o tym, że RMSD wynosi 4Å nie może być zinterpretowana bez wiedzy o tym ile nukleotydów lub aminokwasów zostało porównane

Mean of Circular Quantities (MCQ)

MCQ to miara odległości w przestrzeni trygonometrycznej:
MCQ = \arctan \left( \sum_{i=1}^N \sin \Delta_{i}, \sum_{i=1}^N \cos \Delta_{i} \right)
gdzie:
- N to liczba wartości kątowych
- \Delta_{i} to różnica między i-tą parą kątów
Wizualna reprezentacja:
MCQ pozwala uśrednić dowolne wartości kątowe, ale w kontekście struktur 3D wykorzystywane są kąty torsyjne:
W przypadku RNA, każdy nukleotyd opisany jest przez wiele kątów torsyjnych:
Aby wyznaczyć wartość kąta torsyjnego między atomami A, B, C i D:
- Wyznacz kolejne wektory: v_{AB}, v_{BC}, v_{CD}
- Wartość kąta jest równa:
  \arctan \left( |v_{BC}| v_{AB} \cdot (v_{BC} \times v_{CD}), (v_{AB} \times v_{BC}) \cdot (v_{BC} \times v_{CD}) \right) gdzie:
  - \mathbf{x} \times \mathbf{y} to iloczyn wektorowy
  - \mathbf{x} \cdot \mathbf{y} to iloczyn skalarny
Aby wyznaczyć różnicę kątową, spójrzmy najpierw na projekcję Newmana:
W projekcji Newmana wyjątkowo dobrze widać, że porównująć wartości kątowe trzeba wziąć pod uwagę ich okresowość:
Z tego względu, wzór na różnicę kątową to:
dm \left( \angle_{1}, \angle_{2} \right) = \min \left( |\angle_{1} - \angle_{2}|, 360\degree - |\angle_{1} - \angle_{2}| \right)

Przykład

Dane wejściowe:


	1EHZ			1EVV
	x	y	z	x	y	z
P	50.63	49.73	50.57	22.40	4.87	49.78
O5’	50.16	49.14	52.02	22.62	4.85	51.39
C5’	50.22	49.95	53.21	23.91	4.53	51.92
C4’	50.97	49.23	54.31	23.84	4.39	53.42
C3’	52.45	49.03	54.07	25.19	4.02	54.01
O3’	53.20	50.18	54.43	25.76	5.20	54.54

RMSD

Przesunięcie do centroidu:


	1EHZ			1EVV
	x	y	z	x	y	z
P	-0.64	0.19	-2.53	-1.55	0.23	-2.73
O5’	-1.11	-0.40	-1.08	-1.33	0.21	-1.12
C5’	-1.05	0.41	0.11	-0.04	-0.11	-0.59
C4’	-0.30	-0.31	1.21	-0.11	-0.25	0.91
C3’	1.18	-0.51	0.97	1.24	-0.62	1.50
O3’	1.93	0.64	1.33	1.81	0.56	2.03

Macierz kowariancji:
A^TB = \begin{bmatrix} 7.50 & 0.16 & 9.02 \\ 0.78 & 0.67 & -0.06 \\ 8.83 & -0.98 & 13.31 \end{bmatrix}
Dekompozycja macierzy:
U = \begin{bmatrix} -0.59 & -0.63 & -0.50 \\ -0.02 & -0.61 & 0.79 \\ -0.81 & 0.48 & 0.35 \end{bmatrix} \Sigma = \begin{bmatrix} 19.81 & 0 & 0 \\ 0 & 1.56 & 0 \\ 0 & 0 & 0.14 \end{bmatrix} V = \begin{bmatrix} -0.58 & -0.65 & -0.49 \\ 0.03 & -0.63 & 0.78 \\ -0.81 & 0.44 & 0.39 \end{bmatrix}
Znak:
\det(A^TB) = 4.31 d = 1
Macierz rotacji:
R = \begin{bmatrix}0.99988 & 0.0131 & -0.00809 \\ -0.01362 & 0.99829 & -0.05680 \\ 0.00732 & 0.05691 & 0.99835 \end{bmatrix}

Wynik:


	1EHZ			1EVV
	x	y	z	x	y	z	δ²
P	-0.65	0.19	-2.53	-1.58	0.05	-2.73	0.93
O5’	-1.11	-0.41	-1.08	-1.34	0.13	-1.12	0.33
C5’	-1.06	0.41	0.11	-0.05	-0.15	-0.58	1.79
C4’	-0.30	-0.31	1.21	-0.10	-0.20	0.92	0.13
C3’	1.18	-0.51	0.97	1.26	-0.52	1.52	0.31
O3’	1.93	0.63	1.32	1.81	0.70	1.98	0.44
							RMSD = 0.81

MCQ

Wyznaczymy wartość kąta \beta (P-O5'-C5'-C4') dla przykładu 1EHZ
Wektory wyznaczamy przez różnicę między współrzędnymi:
v_{AB} = \begin{bmatrix} \mathrm{O5'}_{x} - \mathrm{P}_{x} & \mathrm{O5'}_{y} - \mathrm{P}_{y} & \mathrm{O5'}_{z} - \mathrm{P}_{z} \end{bmatrix} v_{AB} = \begin{bmatrix} -0.46 & -0.59 & 1.45 \end{bmatrix} v_{BC} = \begin{bmatrix} 0.05 & 0.81 & 1.19 \end{bmatrix} v_{CD} = \begin{bmatrix} 0.75 & -0.72 & 1.10 \end{bmatrix}
Długość wektora to pierwiastek sumy kwadratów wartości składowych:
|v| = \sqrt{x^2 + y^2 + z^2} |v_{BC}| = 1.44 |v_{BC}| v_{AB} = \begin{bmatrix} -0.67 & -0.85 & 2.09 \end{bmatrix}
Iloczyn wektorowy dla wektorów 3D:
u \times v = \begin{bmatrix} u_{y} v_{z} - u_{z} v_{y} & u_{z} v_{x} - u_{x} v_{z} & u_{x} v_{y} - u_{y} v_{x} \end{bmatrix} v_{AB} \times v_{BC} = \begin{bmatrix} -1.88 & 0.63 & -0.34 \end{bmatrix} v_{BC} \times v_{CD} = \begin{bmatrix} 1.74 & 0.83 & -0.65\end{bmatrix}
Iloczyn skalarny dla wektorów 3D:
u \cdot v = u_{x} v_{x} + u_{y} v_{y} + u_{z} v_{z} |v_{BC}| v_{AB} \cdot (v_{BC} \times v_{CD}) = -3.23 (v_{AB}| \times v_{BC}) \cdot (v_{BC} \times v_{CD}) = -2.53
Wartość kąta:
\beta = \arctan(-3.23, -2.53) = -128.05\degree
Pozostałe wartości kątów i wyliczenia:

1EHZ 1EVV Δ cos(Δ) sin(Δ)

β -128.05 -172.20 44.15 0.72 0.70

γ 67.80 178.70 110.90 -0.36 0.93

δ 82.90 101.10 18.20 0.95 0.31

Σ = 1.31 Σ = 1.94

MCQ = 56.00

	1EHZ	1EVV	Δ	cos(Δ)	sin(Δ)
β	-128.05	-172.20	44.15	0.72	0.70
γ	67.80	178.70	110.90	-0.36	0.93
δ	82.90	101.10	18.20	0.95	0.31
				Σ = 1.31	Σ = 1.94
					MCQ = 56.00

Zadanie

Napisz program, który porówna dwie struktury w formacie PDB przy pomocy miary RMSD oraz MCQ np.
```
$ ./program X.pdb Y.pdb
RMSD = 1.36
MCQ = 12.27
```

Uzyskane punkty będą zależały od tego czy daną funkcjonalność wykorzystano z gotowej biblioteki (B) czy też w postaci własnego kodu (W):


Uzyskane punkty	50%	60%	70%	80%	90%	100%
Wyliczenie MCQ	B	W	W	W	W	W
Wyznaczenie wartości kątów torsyjnych	B	B	B	W	W	W
Wyliczenie RMSD	B	B	W	W	W	W
Superpozycja struktur	B	B	B	B	W	W
Parsowanie PDB	B	B	B	B	B	W

W wersji na 100% dopuszczalne jest użycie wyłącznie biblioteki do wyliczenia dekompozycji macierzy np. scipy.linalg.svg
W wersji na 90% można wykorzystać parser formatu PDB np. Bio.PDB.PDBParser
W wersji na 80% dodatkowo dopuszczalne jest użycie gotowego kodu do superpozycji np. Bio.PDB.Superimposer
W wersji na 70% możesz wyznaczać wartości kątów torsyjnych przy pomocy gotowych funkcji np. z Bio.PDB.vectors
We wszystkich przypadkach, niech MCQ będzie dedykowane strukturom RNA i niech wyznacza tylko wartości kątów \beta (P-O5'-C5'-C4'), \gamma (O5'-C5'-C4'-C3') i \delta (C5'-C4'-C3'-O3')
We wszystkich przypadkach, niech RMSD operuje tylko na atomach fosforu (P)
Jeśli dwie wejściowe struktury mają różną liczbę nukleotydów to należy wypisać komunikat i zakończyć działanie
Dane do testów można pobierać z repozytorium modeli z konkursu RNA-Puzzles np. dla rp08