Opracowanie bazy bioinformatycznej do badań wiązań dwusiarczkowych

Bioinformatyczna baza danych Protein Data Bank (PDB) jest największym na świecie repozytorium eksperymentalnie określonych struktur białek, kwasów nukleinowych i złożonych zespołów. Wszystkie dane są gromadzone za pomocą metod eksperymentalnych, takich jak promieniowanie rentgenowskie, spektroskopia, krystalografia, NMR itp.

W tym artykule wyjaśniono, jak wyodrębniać, filtrować i czyścić dane z pliku PDB. To z kolei umożliwia rodzaj analizy wyjaśniony w artykule Występowanie wiązań dwusiarczkowych białek w różnych domenach życia: porównanie białek z Protein Data Bank, opublikowane w Protein Engineering, Design and Selection , tom 27, wydanie 3, 1 marca 2014, s. 65–72.

PDB ma wiele powtarzających się struktur z różnymi rozdzielczościami, metodami, mutacjami itp. Przeprowadzenie eksperymentu z tymi samymi lub podobnymi białkami może spowodować błąd systematyczny w dowolnej analizie grupowej, więc będziemy musieli wybrać prawidłową strukturę spośród dowolnego zestawu duplikatów . W tym celu musimy użyć nienadmiarowego (NR) zestawu białek.

W celu normalizacji zalecam pobranie słownika związków chemicznych do importowania nazw atomów do bazy danych, która używa 3NF lub używa schematu gwiaździstego i modelowania wymiarowego. (Użyłem również DSSP, aby pomóc wyeliminować problematyczne struktury. Nie będę tego omawiał w tym artykule, ale zauważ, że nie korzystałem z żadnych innych funkcji DSSP.)

Dane użyte w tym badaniu zawierają białka jednojednostkowe, które zawierają co najmniej jedno wiązanie dwusiarczkowe pobrane z różnych gatunków. Aby przeprowadzić analizę, wszystkie wiązania dwusiarczkowe są najpierw klasyfikowane jako kolejne lub niekolejne, według domeny (archaea, prokariota, wirus, eukariota lub inne), a także według długości.

Wyjście

Aby być gotowym do wprowadzenia danych do R, SPSS lub innego narzędzia, analityk będzie potrzebował danych znajdujących się w tabeli bazy danych o następującej strukturze:

Materiały i metody

Aby osiągnąć ten cel, pierwszym krokiem jest zebranie danych z rcsb.org. Ta strona zawiera do pobrania struktury PDB eksperymentów w różnych formatach.

Chociaż dane są przechowywane w wielu formatach, w tym przykładzie zostanie użyty tylko sformatowany format tekstowy rozdzielany spacjami (PDB). Alternatywą dla formatu tekstowego PDB jest jego wersja XML, PDBML, ale czasami zawiera on zniekształcone wpisy nazewnictwa pozycji atomów, co może powodować problemy z analizą danych. Starsze formaty mmCIF i inne mogą być również dostępne, ale nie zostaną wyjaśnione w tym artykule.

Format PDB

Format PDB to pofragmentowany format tekstowy o stałej szerokości, który można łatwo analizować na przykład za pomocą zapytań SQL, wtyczek Java lub modułów Perla. Każdy typ danych w kontenerze pliku jest reprezentowany jako wiersz rozpoczynający się odpowiednim znacznikiem — w kolejnych podsekcjach omówimy każdy typ znacznika. Długość linii jest mniejsza lub równa 80 znaków, gdzie znacznik zajmuje sześć lub mniej znaków plus jedną lub więcej spacji, które łącznie zajmują osiem bajtów. Istnieją również przypadki bez spacji między tagami i danymi, zwykle dla tagów CONECT .

TITLE

Znacznik TITLE oznacza linię jako (część) tytułu eksperymentu, zawierającą nazwę cząsteczki i inne istotne dane, takie jak insercja, mutacja lub delecja określonego aminokwasu.

 12345678901234567890123456789012345678901234567890123456789012345678901234567890 TITLE A TWO DISULFIDE DERIVATIVE OF CHARYBDOTOXIN WITH DISULFIDE TITLE 2 13-33 REPLACED BY TWO ALPHA-AMINOBUTYRIC ACIDS, NMR, 30 TITLE 3 STRUCTURES

W przypadku, gdy rekord TITLE ma wiele wierszy, tytuł musi zostać połączony, uporządkowany według numeru kontynuacji, który jest umieszczany wyrównany do prawej w bajtach 9 i 10 (w zależności od liczby tych wierszy).

ATOM

Dane przechowywane w liniach ATOM to współrzędne dla każdego atomu w eksperymencie. Czasami eksperyment zawiera wstawki, mutacje, alternatywne lokalizacje lub wiele modeli. Powoduje to wielokrotne powtarzanie tego samego atomu. Wybór właściwych atomów zostanie wyjaśniony później.

 12345678901234567890123456789012345678901234567890123456789012345678901234567890 ATOM 390 N GLY A 26 -1.120 -2.842 4.624 1.00 0.00 N ATOM 391 CA GLY A 26 -0.334 -2.509 3.469 1.00 0.00 C ATOM 392 C GLY A 26 0.682 -1.548 3.972 1.00 0.00 C ATOM 393 O GLY A 26 0.420 -0.786 4.898 1.00 0.00 O ATOM 394 H GLY A 26 -0.832 -2.438 5.489 1.00 0.00 H ATOM 395 HA2 GLY A 26 0.163 -3.399 3.111 1.00 0.00 H ATOM 396 HA3 GLY A 26 -0.955 -2.006 2.739 1.00 0.00 H

Powyższy przykład pochodzi z doświadczenia 1BAH . Pierwsza kolumna to typ rekordu, a druga to numer seryjny atomu. Każdy atom w strukturze ma swój numer seryjny.

Obok numeru seryjnego znajduje się etykieta pozycji atom, która zajmuje cztery bajty. Z tej pozycji atomu można wyodrębnić symbol chemiczny pierwiastka, który nie zawsze występuje w danych rekordu w osobnej kolumnie.

Po nazwie atomu znajduje się trzyliterowy kod reszty. W przypadku białek reszta ta odpowiada aminokwasowi. Następnie łańcuch jest kodowany jedną literą. Przez łańcuch rozumiemy pojedynczy łańcuch aminokwasów, z przerwami lub bez, chociaż czasami do łańcucha można przypisać ligandy; ten przypadek jest wykrywalny przez bardzo duże przerwy w sekwencji aminokwasowej, która znajduje się w następnej kolumnie. Czasami tę samą strukturę można przeskanować z uwzględnieniem mutacji, w którym to przypadku kod wstawienia jest dostępny w dodatkowej kolumnie po kolumnie sekwencji. Kod wstawiania zawiera literę, która pomaga odróżnić, której pozostałości dotyczy problem.

Kolejne trzy kolumny to współrzędne przestrzenne każdego atomu mierzone w angstremach (Å). Obok tych współrzędnych znajduje się kolumna zajętości, która mówi, jakie jest prawdopodobieństwo, że atom znajdzie się w tym miejscu, w zwykłej skali od zera do jednego.

Kolumna przedostatnia to współczynnik temperatury, który zawiera informacje o nieuporządkowaniu kryształu, mierzonym w Ų. Wartość większa niż 60 Ų oznacza zaburzenie, a wartość niższa niż 30 Ų oznacza pewność. Nie zawsze występuje w plikach PDB, ponieważ zależy od metody eksperymentalnej.

Zwykle brakuje następnych kolumn — symbolu i ładunku. Symbol chemiczny można pobrać z kolumny pozycji atomu, jak wspomnieliśmy powyżej. Gdy ładunek jest obecny, jest dodawany do symbolu jako liczba całkowita, po której następuje + lub - , np. N1+ .

TER

To oznacza koniec łańcucha. Nawet bez tej linii łatwo odróżnić, gdzie kończy się łańcuch. Dlatego często brakuje linii TER .

MODEL i ENDMDL

Linia MODEL oznacza, gdzie zaczyna się model konstrukcji i zawiera numer seryjny modelu.

Po wszystkich liniach atomowych w tym modelu kończy się linią ENDMDL .

SSBOND

Linie te zawierają wiązania dwusiarczkowe między aminokwasami cysteinowymi. Wiązania disiarczkowe mogą być obecne w innych typach reszt, ale w tym artykule analizowane będą tylko aminokwasy, więc uwzględniona jest tylko cysteina. Poniższy przykład pochodzi z eksperymentu z kodem 132L :

 12345678901234567890123456789012345678901234567890123456789012345678901234567890 SSBOND 1 CYS A 6 CYS A 127 1555 1555 2.01 SSBOND 2 CYS A 30 CYS A 115 1555 1555 2.05 SSBOND 3 CYS A 64 CYS A 80 1555 1555 2.02 SSBOND 4 CYS A 76 CYS A 94 1555 1555 2.02

W tym przykładzie w pliku są oznaczone cztery wiązania dwusiarczkowe, których numer kolejny znajduje się w drugiej kolumnie. Wszystkie te wiązania wykorzystują cysteinę (kolumny 3 i 6) i wszystkie są obecne w łańcuchu A (kolumny 4 i 7). Po każdym łańcuchu znajduje się numer sekwencyjny reszty wskazujący pozycję wiązania w łańcuchu peptydowym. Kody insercji znajdują się obok każdej sekwencji reszt, ale w tym przykładzie nie są one obecne, ponieważ w tym regionie nie wstawiono aminokwasu. Dwie kolumny przed ostatnią są zarezerwowane dla operacji symetrii, a ostatnia kolumna to odległość między atomami siarki, mierzona w Å.

Poświęćmy chwilę, aby nadać kontekst tym danym.

Poniższe zdjęcia, wykonane przy użyciu przeglądarki rcsb.org NGL, pokazują strukturę eksperymentu 132L . W szczególności pokazują białko bez ligandów. Pierwsze zdjęcie wykorzystuje reprezentację sztyftu, w której kolory CPK pokazują siarki i ich wiązania na żółto. Połączenia siarki w kształcie litery V reprezentują połączenia metioniny, podczas gdy połączenia w kształcie litery Z są wiązaniami dwusiarczkowymi między cysteinami.

Na następnym rysunku uproszczona metoda wizualizacji białek zwana wizualizacją szkieletu jest pokazana pokolorowana aminokwasami, gdzie cysteiny są żółte. Reprezentuje to samo białko, z wykluczonym łańcuchem bocznym i zawiera tylko część jego grupy peptydowej – w tym przypadku szkielet białkowy. Składa się z trzech atomów: N-końcowego, C-alfa i C-końcowego. Ten obrazek nie pokazuje wiązań dwusiarczkowych, ale jest uproszczony, aby pokazać układ przestrzenny białka:

Rury powstają przez połączenie atomów związanych peptydem z atomem C-alfa. Kolor cysteiny jest taki sam jak kolor siarki w metodzie barwienia CPK. Gdy cysteny zbliżą się wystarczająco blisko, ich siarki tworzą wiązania dwusiarczkowe, co wzmacnia strukturę. W przeciwnym razie białko to wiązałoby się zbyt mocno, a jego struktura byłaby mniej stabilna w wyższych temperaturach.

CONECT

Rekordy te służą do oznaczania połączeń między atomami. Czasami te znaczniki w ogóle nie występują, a innym razem wszystkie dane są wypełnione. W przypadku analizy wiązań dwusiarczkowych ta część może być przydatna – ale nie jest konieczna. Dzieje się tak dlatego, że w tym projekcie nieoznakowane wiązania są dodawane przez pomiar odległości, więc byłoby to nad głową i również musi zostać sprawdzone.

SOURCE

Rekordy te zawierają informacje o organizmie źródłowym, z którego wyodrębniono cząsteczkę. Zawierają podrekordy dla łatwiejszej lokalizacji w taksonomii i mają taką samą wielowierszową strukturę, jaką widzieliśmy w przypadku rekordów tytułowych.

 SOURCE MOL_ID: 1; SOURCE 2 ORGANISM_SCIENTIFIC: ANOPHELES GAMBIAE; SOURCE 3 ORGANISM_COMMON: AFRICAN MALARIA MOSQUITO; SOURCE 4 ORGANISM_TAXID: 7165; SOURCE 5 GENE: GST1-6; SOURCE 6 EXPRESSION_SYSTEM: ESCHERICHIA COLI; SOURCE 7 EXPRESSION_SYSTEM_TAXID: 562; SOURCE 8 EXPRESSION_SYSTEM_STRAIN: BL21(DE3)PLYSS; SOURCE 9 EXPRESSION_SYSTEM_VECTOR_TYPE: PLASMID; SOURCE 10 EXPRESSION_SYSTEM_PLASMID: PXAGGST1-6

Format NR

To jest lista nienadmiarowych (NR) zestawów PDB łańcucha. Jego migawki można znaleźć na ftp.ncbi.nih.gov/mmdb/nrtable/. Jego celem jest uniknięcie niepotrzebnych błędów spowodowanych podobieństwem białek. NR ma trzy zestawy o różnych poziomach wartości p tożsamości utworzonych przez porównanie wszystkich struktur PDB. Wynik jest dodawany do plików tekstowych, które zostaną wyjaśnione później. Nie wszystkie kolumny są potrzebne do tego projektu, więc wyjaśnione zostaną tylko te ważne.

Pierwsze dwie kolumny zawierają unikalny kod eksperymentu PDB i identyfikator łańcucha, jak wyjaśniono powyżej dla rekordów ATOM . Kolumny 6, 9 i C zawierają informacje o reprezentatywności wartości p, która jest poziomem podobieństwa sekwencji obliczonym przez BLAST. Jeśli ta wartość wynosi zero, to nie jest akceptowane jako część zestawu; jeśli wartość wynosi 1, to tak jest. Wspomniane kolumny reprezentują akceptację zestawów z wartościami odcięcia wartości p odpowiednio 10e-7, 10e-40 i 10e-80. Do analizy zostaną użyte tylko zestawy z wartością odcięcia 10e-7.

Ostatnia kolumna zawiera informację o akceptowalności konstrukcji, gdzie a jest akceptowalne, a n nie.

 #--------------------------------------------------------------------------------------------------------------------------- # 1 2 3 4 5 6 7 8 9 ABCDEFGHIJKLMNOPQ #--------------------------------------------------------------------------------------------------------------------------- 3F8V A 69715 1 1 1 1 1 1 1 1 1 9427 1 1 0.00 0.00 0.00 0.00 1.08 1 6 5 164 X a 3DKE X 68132 1 2 0 1 2 0 1 2 0 39139 1 1 0.00 0.00 0.00 0.00 1.25 1 11 7 164 X a 3HH3 A 77317 1 3 0 1 3 0 1 3 0 90 1 0 0.00 0.00 0.00 0.00 1.25 1 5 4 164 X a 3HH5 A 77319 1 4 0 1 4 0 1 4 0 90 2 0 0.00 0.00 0.00 0.00 1.25 1 4 4 164 X a

Budowa bazy danych i analiza danych

Teraz, gdy mamy już pomysł na to, z czym mamy do czynienia i co musimy zrobić, zacznijmy.

Pobieranie danych

Wszystkie dane do tej analizy można znaleźć pod tymi trzema adresami:

• ftp.wwpdb.org/pub/pdb/data/structures/divided/pdb/
• ftp.ncbi.nih.gov/mmdb/nrtable/
• ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zip

Pierwsze dwa linki zawierają listę archiwów. Należy używać najnowszego archiwum z każdego z nich, aby uniknąć problemów wynikających z braku rozdzielczości lub jakości. Trzeci link zawiera bezpośrednio najnowsze archiwum taksonomii.

Analiza danych

Zwykle parsowanie plików PDB odbywa się za pomocą wtyczek lub modułów w Javie, Perlu lub Pythonie. W przypadku tego badania napisałem niestandardową aplikację Perl bez użycia wcześniej napisanego modułu parsowania PDB. Powodem tego jest to, że przy parsowaniu dużej ilości danych, z mojego doświadczenia wynika, że ​​najczęstszym problemem przy korzystaniu z danych eksperymentalnych są błędy w danych. Czasami pojawiają się błędy ze współrzędnymi, odległościami, długościami linii, komentarzami w miejscach, w których nie powinno ich być itp.

Najskuteczniejszym sposobem radzenia sobie z tym jest początkowo przechowywanie wszystkiego w bazie danych jako nieprzetworzonego tekstu. Popularne parsery są napisane tak, aby radzić sobie z idealnymi danymi, które są całkowicie zgodne ze specyfikacjami. Ale w praktyce dane nie są idealne i zostanie to wyjaśnione w sekcji filtrowania, gdzie znajdziesz skrypt importu Perla.

Budowa bazy danych

Podczas konstruowania bazy danych należy pamiętać, że ta baza danych jest przeznaczona do przetwarzania danych. Późniejsza analiza zostanie przeprowadzona w SPSS lub R. Dla naszych celów tutaj zaleca się użycie PostgreSQL w wersji co najmniej 8.4.

Struktura tabeli jest bezpośrednio kopiowana z pobranych plików z kilkoma niewielkimi zmianami. W tym przypadku liczba rekordów jest zdecydowanie za mała, aby warto było poświęcać czas na normalizację. Jak wspomniano, ta baza danych jest przeznaczona tylko do jednorazowego użytku: te tabele nie są przeznaczone do obsługi w witrynie internetowej — służą tylko do przetwarzania danych. Po zakończeniu można je usunąć lub utworzyć kopię zapasową jako dane uzupełniające, być może w celu powtórzenia procesu przez innego badacza.

W takim przypadku końcowym wynikiem będzie jedna tabela, którą można następnie wyeksportować do pliku w celu użycia z jakimś narzędziem statystycznym, takim jak SPSS lub R.

Stoły

Ekstrakcja danych z rekordów ATOM musi być połączona z rekordami HEADER lub TITLE . Hierarchię danych wyjaśniono na poniższym obrazku.

Ponieważ ten obraz jest uproszczoną reprezentacją bazy danych w trzeciej postaci normalnej (3NF), dla naszych celów zawiera zbyt dużo narzutu. Powód: Aby obliczyć odległość między atomami w celu wykrycia wiązania dwusiarczkowego, musielibyśmy wykonać złączenia. W tym przypadku mielibyśmy stół połączony ze sobą dwukrotnie, a także połączony ze strukturą drugorzędową i pierwotną dwa razy, co jest bardzo powolnym procesem. Ponieważ nie każda analiza wymaga informacji o strukturze drugorzędowej, proponowany jest inny schemat na wypadek konieczności ponownego wykorzystania danych lub analizy większych ilości wiązań dwusiarczkowych:

Wiązania dwusiarczkowe nie są tak częste, jak inne wiązania kowalencyjne, więc model magazynowy nie jest potrzebny, chociaż można go zastosować. Opracowanie poniższego schematu gwiaździstego i modelowania wymiarowego zajmie zbyt dużo czasu i sprawi, że zapytania będą bardziej złożone:

W przypadkach, w których wszystkie obligacje muszą zostać przetworzone, polecam schemat gwiazdy.

(W przeciwnym razie nie jest to potrzebne, ponieważ wiązania dwusiarczkowe nie są tak powszechne, jak inne wiązania. W przypadku tej pracy liczba wiązań dwusiarczkowych wynosi blisko 30 000, co może być wystarczająco szybkie w 3NF, ale zdecydowałem się przetworzyć za pomocą tabela nieznormalizowana, więc nie jest tutaj przedstawiona.)

Oczekiwana całkowita liczba wszystkich wiązań kowalencyjnych jest co najmniej dwukrotnością liczby atomów w strukturze trzeciorzędowej, a w takim przypadku 3NF byłby bardzo powolny, więc potrzebna jest denormalizacja przy użyciu postaci schematu gwiazdy. W tym schemacie niektóre tabele mają dwa sprawdzenia kluczy obcych, a to dlatego, że wiązanie jest tworzone między dwoma atomami, więc każdy atom musi mieć swój własny primary_structure_id , atom_name_id i residue_id .

Istnieją dwa sposoby wypełnienia tabeli wymiarów d_atom_name : z danych iz innego źródła, ze słownika składników chemicznych, o którym wspomniałem wcześniej. Jego format jest podobny do formatu PDB: przydatne są tylko wiersze RESIDUE i CONECT . Dzieje się tak, ponieważ pierwsza kolumna RESIDUE zawiera trzyliterowy kod reszty, a CONECT zawiera nazwę atomu i jego połączeń, które są również nazwami atomów. Tak więc z tego pliku możemy przeanalizować wszystkie nazwy atomów i uwzględnić je w naszej bazie danych, chociaż zalecam dopuszczenie możliwości bazy danych zawierającej niewymienione nazwy atomów.

 RESIDUE PRO 17 CONECT N 3 CA CD H CONECT CA 4 NC CB HA CONECT C 3 CA O OXT CONECT O 1 C CONECT CB 4 CA CG HB2 HB3 CONECT CG 4 CB CD HG2 HG3 CONECT CD 4 N CG HD2 HD3 CONECT OXT 2 C HXT CONECT H 1 N CONECT HA 1 CA CONECT HB2 1 CB CONECT HB3 1 CB CONECT HG2 1 CG CONECT HG3 1 CG CONECT HD2 1 CD CONECT HD3 1 CD CONECT HXT 1 OXT END HET PRO 17 HETNAM PRO PROLINE FORMUL PRO C5 H9 N1 O2

W tym projekcie szybkość kodowania jest bardziej istotna niż szybkość wykonywania i zużycie pamięci. Postanowiłem nie normalizować – w końcu naszym celem jest wygenerowanie tabeli z kolumnami wymienionymi we wstępie.

W tej części wyjaśnione zostaną tylko najważniejsze tabele.

Główne tabele to:

• proteins : Tabela z nazwami i kodami eksperymentów.
• ps : Podstawowa tabela struktury, która będzie zawierać sequence , chain_id i code .
• ts : Tabela zawierająca trzeciorzędową/czwartorzędową strukturę wyodrębnioną z surowych danych i przekształconą w format rekordu ATOM . Będzie ona używana jako tabela pomostowa i może zostać usunięta po wyodrębnieniu. Ligandy są wykluczone.
• sources : Lista organizmów, z których uzyskano dane eksperymentalne.
• tax_names , taxonomy_path , taxonomy_paths : Linneowskie nazwy taksonomii z bazy taksonomii NCBI, używane do pobierania ścieżek taksonomii z organizmów wymienionych w sources .
• nr : Lista nieredundantnych białek NCBI wyekstrahowanych z zestawu NR.
• pdb_ssbond : Lista wiązań dwusiarczkowych w danym pliku PDB.

Filtrowanie i przetwarzanie danych

Dane są pobierane z migawek z repozytorium RCSB PDB.

Każdy plik jest importowany do pojedynczej tabeli raw_pdb w naszej bazie danych PostgreSQL za pomocą skryptu Perla. Skrypt wykorzystuje transakcje 10 000 wstawek na porcję.

Struktura raw_pdb jest następująca:

Skrypt importu wygląda tak:

 #!/usr/bin/perl -w use FindBin '$Bin'; use DBI; $dbName = 'bioinf'; $dbLogin = 'ezop'; $dbPass = 'XYZ'; $conn = DBI->connect("DBI:Pg:database=$dbName;host=localhost", "$dbLogin", "$dbPass", {'RaiseError' => 1, 'AutoCommit' => 0}); die "./pdb_lines_unos.pl <table> <file>" if not defined($ARGV[0]); $recordCount = 0; $table = $ARGV[0]; $fName = $ARGV[1]; open F, "zcat $fName|"; while (<F>) { chomp; $linija = $_; $recordCount += 1; $sql = "insert into $table (code, line_num, line_cont) values (?, ?, ?)"; $conn->do($sql, undef, $fName, $recordCount, $linija); $conn->commit() if ($recordCount%10000 == 0); } close F; $conn->commit(); 1;

Po zaimportowaniu wierszy są one analizowane za pomocą funkcji, które zdefiniujemy poniżej.

Z danych raw_pdb generujemy tabele ts , ps , proteins , sources , sources_organela i ss_bond , analizując odpowiednie rekordy.

Tabela ps ma trzy ważne kolumny: chain , length i sequence . Sekwencja białka jest generowana przy użyciu atomów C-alfa dla każdego łańcucha i dla każdej reszty uporządkowanej według sekwencji reszt, biorąc tylko pierwszą insercję i pierwszą alternatywną lokalizację. chain jest pobierany z kolumny TS.chain , a length jest po prostu wstępnie obliczoną długością ciągu sequence . Ponieważ ten artykuł ma na celu analizę tylko pojedynczych łańcuchów i połączeń wewnątrzłańcuchowych, w naszej analizie pominięto białka wielołańcuchowe.

W rekordach SSBOND wszystkie wiązania dwusiarczkowe są przechowywane w tabeli pdb_ssbond , która dziedziczy z tabeli pdb_ssbond_extended . pdb_ssbond wygląda tak:

Dodałem również te indeksy:

 "pdb_ssbond_pkey" PRIMARY KEY, btree (id) "ndxcode1" btree (code, chain_id1, res_seq1) "ndxcode2" btree (code, chain_id2, res_seq2)

Zakłada się, że rozkład wiązań dwusiarczkowych przed wartością graniczną jest gaussowski (bez testowania np. testem KS), więc odchylenia standardowe obliczono dla każdej odległości między cysteinami w tym samym białku, aby odkryć zakres dozwolonych długości wiązań i je porównać do granicy. Wartość graniczna była taka sama jak obliczona średnia plus minus trzy odchylenia standardowe. Rozszerzyliśmy zakres, aby wprowadzić więcej możliwych wiązań dwusiarczkowych, które nie zostały wymienione w pliku PDB w wierszach SSBOND , ale które wprowadziliśmy samodzielnie, obliczając odległości między rekordami ATOM . Wybrany zakres dla obligacji ssbonds wynosi od 1,6175344456264 do 2,48801947642267 Å, ​​co stanowi średnią (2,05) plus minus cztery odchylenia standardowe:

 select count(1) as cnt , stddev(dist) as std_dev , avg(dist) as avg_val , -stddev(dist) + avg(dist) as left1 , stddev(dist) + avg(dist) as right1 , -2 * stddev(dist) + avg(dist) as left2 , 2 * stddev(dist) + avg(dist) as right2 , -3 * stddev(dist) + avg(dist) as left3 , 3 * stddev(dist) + avg(dist) as right3 , -4 * stddev(dist) + avg(dist) as left4 , 4 * stddev(dist) + avg(dist) as right4 , min(dist) , max(dist) from pdb_ssbond_tmp where dist > 0 and dist < 20;

Tabela TS zawiera współrzędne wszystkich atomów, ale używane będą tylko cysteiny, których siarka będzie nazywana " SG " . Tak więc utworzona została kolejna tabela stopniowania zawierająca tylko atomy siarki " SG " w celu przyspieszenia procesu poprzez zmniejszenie liczby rekordów do przeszukania. Przy selekcji samych siarki liczba kombinacji jest znacznie mniejsza niż w przypadku wszystkich atomów – 194 574 pierwszego w porównaniu z 122 761 100 drugich. W tej tabeli połączonej ze sobą odległości są obliczane przy użyciu odległości euklidesowej, a wyniki są importowane do tabeli pdb_ssbond , ale tylko wtedy, gdy odległość mieści się między zdefiniowanymi długościami obliczonymi wcześniej. Powodem tego przyspieszenia jest skrócenie czasu ponownego uruchamiania całego procesu — w celach kontrolnych — utrzymując go w ciągu jednego dnia.

Do czyszczenia wiązań dwusiarczkowych stosujemy następujący algorytm:

• Usuń, gdy obie strony połączenia wskazują ten sam aminokwas
• Usuń obligacje, których długość nie mieści się w przedziale od 1.6175344456264 do 2.48801947642267
• Usuń wstawki
• Usuń wiązania spowodowane przez alternatywne lokalizacje atomów, ale pozostawiając pierwszą lokalizację

Kod do tego (przyjmując tabelę pdb_ssbond jako pierwszy argument) to:

 CREATE OR REPLACE FUNCTION pdb_ssbond_clean2( clean_icodes boolean, clean_altloc boolean, mark_reactive boolean, mark_consecutive boolean) RETURNS void AS $$ declare _res integer; BEGIN delete from pdb_ssbond b where exists ( select a.id from pdb_ssbond a where a.code = b.code and a.id > b.id and ( ( a.chain_id1 = b.chain_id1 and a.res_seq1 = b.res_seq1 and a.chain_id2 = b.chain_id2 and a.res_seq2 = b.res_seq2) or ( a.chain_id1 = b.chain_id2 and a.res_seq1 = b.res_seq2 and a.chain_id2 = b.chain_id1 and a.res_seq2 = b.res_seq1) ) ) ; delete from pdb_ssbond where chain_id1 = chain_id2 and res_seq1 = res_seq2 and i_code1 = i_code2; update pdb_ssbond set is_consecutive = true , is_reactive = false; delete from pdb_ssbond where not dist between 1.6175344456264 and 2.48801947642267; if clean_icodes then delete from pdb_ssbond where i_code1 not in ('', ' ', 'A') or i_code2 not in ('', ' ', 'A') ; end if; if clean_altloc then delete from pdb_ssbond a where exists ( select 1 from pdb_ssbond b where b.code = a.code and b.chain_id1 = a.chain_id1 and b.res_seq1 = a.res_seq1 and b.i_code1 = a.i_code1 and b.ser_num < a.ser_num ) ; delete from pdb_ssbond a where exists ( select 1 from pdb_ssbond b where b.code = a.code and b.chain_id2 = a.chain_id2 and b.res_seq2 = a.res_seq2 and b.i_code2 = a.i_code2 and b.ser_num < a.ser_num ) ; end if; if mark_reactive then update pdb_ssbond set is_reactive = false ; update pdb_ssbond set is_reactive = true where exists ( select 1 from pdb_ssbond b where b.code = pdb_ssbond.code and b.chain_id1 = pdb_ssbond.chain_id1 and b.res_seq1 = pdb_ssbond.res_seq1 and b.i_code1 = pdb_ssbond.i_code1 and b.ser_num < pdb_ssbond.ser_num ) ; update pdb_ssbond set is_reactive = true where exists ( select 1 from pdb_ssbond b where b.code = pdb_ssbond.code and b.chain_id2 = pdb_ssbond.chain_id2 and b.res_seq2 = pdb_ssbond.res_seq2 and b.i_code2 = pdb_ssbond.i_code2 and b.ser_num < pdb_ssbond.ser_num ) ; update pdb_ssbond set is_reactive = true where (code, chain_id1, res_seq1, i_code1) in ( select code, chain_id1 as c, res_seq1 as r, i_code1 as i from pdb_ssbond intersect select code, chain_id2 as c, res_seq2 as r, i_code2 as i from pdb_ssbond ) ; update pdb_ssbond set is_reactive = true where (code, chain_id2, res_seq2, i_code2) in ( select code, chain_id1 as c, res_seq1 as r, i_code1 as i from pdb_ssbond intersect select code, chain_id2 as c, res_seq2 as r, i_code2 as i from pdb_ssbond ); end if; if mark_consecutive then update pdb_ssbond set is_consecutive = false; update pdb_ssbond set is_consecutive = true where not exists ( select 1 from pdb_ssbond a where a.code = pdb_ssbond.code and ( (a.chain_id1 = pdb_ssbond.chain_id1 and a.chain_id2 = pdb_ssbond.chain_id2) or (a.chain_id1 = pdb_ssbond.chain_id2 and a.chain_id2 = pdb_ssbond.chain_id1) ) and ( (a.res_seq1 between pdb_ssbond.res_seq1 and pdb_ssbond.res_seq2) or (a.res_seq2 between pdb_ssbond.res_seq1 and pdb_ssbond.res_seq2) or (pdb_ssbond.res_seq1 between a.res_seq1 and a.res_seq2) or (pdb_ssbond.res_seq2 between a.res_seq1 and a.res_seq2) ) and not ( a.code = pdb_ssbond.code and a.chain_id1 = pdb_ssbond.chain_id1 and a.chain_id2 = pdb_ssbond.chain_id2 and a.res_seq1 = pdb_ssbond.res_seq1 and a.res_seq2 = pdb_ssbond.res_seq2 ) ); end if; END; $$ LANGUAGE plpgsql;

Dzięki temu nieredundantny zestaw białek jest importowany do tablicy nr , która jest dołączona do tablic ps i proteins , a zestawy są oznaczane ( set7 , set40 i set80 ). Na koniec, w zależności od ilości białka, tylko jeden zestaw zostanie przeanalizowany. Niedopasowane łańcuchy między PDB i NR są usuwane z analizy.

Białka bez wiązań dwusiarczkowych są wyłączone z badań, podobnie jak białka, które nie należą do żadnego zestawu. Dane są przetwarzane za pomocą DSSP, a te pliki, które miały problemy z rozdzielczością lub zbyt wiele atomów do przetworzenia, również są wykluczone. Jako wynik analizy wykorzystuje się tylko białka z pojedynczymi łańcuchami, ponieważ połączenia międzyłańcuchowe nie zostały zachowane, chociaż można je łatwo obliczyć z tabeli ssbond , licząc liczbę połączeń, które mają różne łańcuchy.

W przypadku pozostałych białek aktualizacja jest wykonywana dla całkowitej liczby wiązań i liczby nakładających się wiązań i jest to wykonywane dla każdego z zestawów.

Organizm źródłowy jest wyodrębniany z rekordów SOURCE . W przypadkach, gdy jest nieznany, syntetyczny, zaprojektowany, sztuczny lub hybrydowy, jest odrzucany z badań. Białka o niskiej rozdzielczości są również wykluczone tylko wtedy, gdy ich łańcuch boczny nie jest widoczny.

Rekordy SOURCE są przechowywane w tabeli sources , która zawiera wiersze taksonomii. W niektórych przypadkach taksonomii brakuje lub jest ona nieprawidłowa. W takich przypadkach konieczna jest ręczna korekta ekspertów.

Na podstawie źródła i taksonomii pobranej z NCBI, klasa jest przypisywana do każdej struktury podstawowej. W przypadku, gdy nie można przypisać klasy, białko jest usuwane z listy analiz. Wiedząc, że wykorzystywane są biologiczne bazy danych, zaleca się, aby biolog przeprowadził dodatkową kontrolę wszystkich rekordów źródłowych i klas taksonomii NCBI, w przeciwnym razie mogą wystąpić problemy z klasyfikacją między różnymi dziedzinami. Aby dopasować źródłowe lokalizacje komórkowe do identyfikatorów taksonomii, dane z tabeli źródłowej są wyodrębniane do tabeli sources_organela , w której wszystkie dane są zapisywane jako kod, znacznik i wartość. Jego format pokazano poniżej:

 select * from sources_organela where code = '1rav';

Archiwum taksonomii, którego chcemy użyć, to plik ZIP zawierający siedem plików zrzutów. Wśród tych plików dwa najważniejsze to names.dmp i merged.dmp . Oba pliki są plikami CSV rozdzielanymi tabulatorami, zgodnie z opisem w dokumentacji:

• Plik merged.dmp zawiera listę identyfikatorów poprzednich taksonomii oraz identyfikatorów obecnych taksonomii, z którymi każdy z nich został scalony.
• names.dmp zawiera te ważne kolumny w następującej kolejności:
  • tax_id : identyfikator taksonomii
  • name_txt : nazwa gatunku i, jeśli ma to zastosowanie, unikalna nazwa (gdzie gatunki można znaleźć z wieloma nazwami, to pomaga ujednoznacznieć)
• division.dmp zawiera nazwy domen najwyższego poziomu, których będziemy używać jako naszych klas.
• nodes.dmp to tabela zawierająca hierarchiczną strukturę organizmów przy użyciu identyfikatorów taksonomii.
  • Zawiera nadrzędny identyfikator taksonomii, klucz obcy, który można znaleźć w names.dmp .
  • Zawiera również identyfikator oddziału, który jest dla nas ważny, aby poprawnie przechowywać odpowiednie dane domeny górnej.

Dzięki tym wszystkim danym i ręcznym poprawkom (ustawieniu właściwych dziedzin życia) byliśmy w stanie stworzyć strukturę tabeli taxonomy_path . Próbka danych wygląda tak:

 select * from taxonomy_path order by length(path) limit 20 offset 2000;

Before any analysis, to avoid biases, sequences have to be checked for their level of identity. Although the NR set contains sequences which are already compared between each other, an extra check is always recommended.

For each disulfide bond's prior statistical analysis, data is marked if it is reactive or overlapping. After marking overlaps, that information automatically reveals how many consecutive and non-consecutive bonds are inside each protein, and that data is stored in the proteins table from which all protein complexes are excluded in final result.

Each disulfide bond is marked also for its association to sets, by checking both bond sides to see if they belong to the same NR set. Where that is not the case, the disulfide bond is skipped for that set analysis.

To analyze the quantity of bonds by their variance, each length has to be put in a specific class. In this case, only five classes are chosen as written in the function below.

 CREATE OR REPLACE FUNCTION len2class(len integer) RETURNS integer AS $BODY$ BEGIN return case when len <= 100 then 1 when len > 100 and len <= 200 then 2 when len > 200 and len <= 300 then 3 when len > 300 and len <= 400 then 4 else 5 end; END; $BODY$ LANGUAGE plpgsql VOLATILE COST 100;

Most of the protein sizes are less than 400 amino acids, so length classification is done by splitting lengths into five ranges, each taking 100 amino acids, except the last one which takes the rest. Below you can see how to use this function to extract data for analysis:

 SELECT p.code, p.title, p.ss_bonds, p.ssbonds_overlap, p.intra_count, p.sci_name_src, p.len, p.tax_path, p.pfam_families, p.src_class, ( SELECT s.id FROM src_classes s WHERE s.src_class::text = p.src_class::text) AS src_class_id, p.len_class7, ( SELECT s.val FROM sources_organela s WHERE s.code = p.code::bpchar AND s.tag::text = 'EXPRESSION_SYSTEM_CELLULAR_LOCATION'::text) AS expression_system_cellular_location, ( SELECT s.val FROM sources_organela s WHERE s.code = p.code::bpchar AND s.tag::text = 'CELLULAR_LOCATION'::text) AS cellular_location, ps.sequence, ps.uniprot_code, ps.accession_code, ps.cc, ps.seq_uniprot, ps.chain_id FROM proteins p JOIN nr n ON n.code::text = p.code::text AND n.rep7 = true JOIN ps ps ON ps.code::text = n.code::text AND ps.chain_id = n.chain_id AND ps.het = false WHERE p.is_excluded = false AND p.chain_cnt = 1 AND p.is_set7 = true AND p.reactive_cnt7 = 0 ORDER BY p.code;

An example result with corrected titles and some manually added columns is shown below.

PostgreSQL jako pośrednik przetwarzania

W tej pracy pokazaliśmy, jak przetwarzać dane, od pobierania po analizę. Podczas pracy z danymi naukowymi czasami potrzebna jest normalizacja, a czasami nie. W przypadku pracy z małymi ilościami danych, które nie zostaną ponownie wykorzystane do kolejnej analizy, wystarczy pozostawić je zdenormalizowane, gdy przetwarzanie jest wystarczająco szybkie.

Powodem, dla którego zrobiono to wszystko w jednej bioinformatycznej bazie danych, jest fakt, że PostgreSQL jest w stanie zintegrować wiele języków. Obejmuje to R, który może przeprowadzać analizy statystyczne bezpośrednio w bazie danych — temat przyszłego artykułu na temat narzędzi bioinformatycznych.

Specjalne podziękowania dla kolegów z Toptal Stefana Fuchsa i Aldo Zelen za ich nieocenione konsultacje.