Sviluppo di un database bioinformatico per la ricerca sui legami disolfuro

Il database bioinformatica della Protein Data Bank (PDB) è il più grande archivio al mondo di strutture di proteine, acidi nucleici e complessi complessi determinate sperimentalmente. Tutti i dati vengono raccolti utilizzando metodi sperimentali come raggi X, spettroscopia, cristallografia, NMR, ecc.

Questo articolo spiega come estrarre, filtrare e pulire i dati dal PDB. Questo, a sua volta, abilita il tipo di analisi spiegato nell'articolo Occurrence of protein disolfuro legami in diversi domini della vita: un confronto di proteine ​​dalla banca dati delle proteine, pubblicato in Protein Engineering, Design and Selection , Volume 27, Issue 3, 1 marzo 2014, pagg. 65–72.

Il PDB ha molte strutture ripetitive con diverse risoluzioni, metodi, mutazioni, ecc. Fare un esperimento con proteine ​​uguali o simili può produrre bias in qualsiasi analisi di gruppo, quindi dovremo scegliere la struttura corretta tra qualsiasi serie di duplicati . A tale scopo, abbiamo bisogno di utilizzare un insieme di proteine ​​non ridondanti (NR).

Ai fini della normalizzazione, consiglio di scaricare il dizionario dei composti chimici per importare i nomi degli atomi in un database che utilizza 3NF o utilizza uno schema a stella e una modellazione dimensionale. (Ho anche usato DSSP per eliminare le strutture problematiche. Non lo tratterò in questo articolo, ma tieni presente che non ho utilizzato altre funzionalità DSSP.)

I dati utilizzati in questa ricerca contengono proteine ​​a unità singola che contengono almeno un legame disolfuro prelevato da specie diverse. Per eseguire un'analisi, tutti i legami disolfuro vengono prima classificati come consecutivi o non consecutivi, per dominio (archea, procariota, virale, eucariota o altro) e anche per lunghezza.

Produzione

Per essere pronto per l'input in R, SPSS o qualche altro strumento, un analista avrà bisogno che i dati si trovino in una tabella di database con questa struttura:

Materiali e metodi

Per raggiungere questo obiettivo, il primo passo è raccogliere dati da rcsb.org. Quel sito contiene strutture PDB scaricabili di esperimenti in vari formati.

Sebbene i dati siano archiviati in più formati, in questo esempio verrà utilizzato solo il formato testuale delimitato da spazio fisso (PDB) formattato. Un'alternativa al formato testuale PDB è la sua versione XML, PDBML, ma a volte contiene voci di denominazione della posizione dell'atomo non corrette, che possono causare problemi per l'analisi dei dati. Potrebbero essere disponibili anche i vecchi formati mmCIF e altri, ma non verranno spiegati in questo articolo.

Il formato PDB

Il formato PDB è un formato testuale frammentato a larghezza fissa che può essere facilmente analizzato da query SQL, plug-in Java o moduli Perl, ad esempio. Ciascun tipo di dati nel contenitore di file è rappresentato da una riga che inizia con il tag appropriato: esamineremo ciascun tipo di tag nelle seguenti sottosezioni. La lunghezza della riga è inferiore o uguale a 80 caratteri, dove un tag richiede sei o meno caratteri più uno o più spazi che insieme occupano otto byte. Esistono anche casi senza spazi tra tag e dati, di solito per i tag CONECT .

TITLE

Il tag TITLE contrassegna una riga come (parte del) titolo dell'esperimento, contenente il nome della molecola e altri dati rilevanti come l'inserimento, la mutazione o l'eliminazione di un amminoacido specifico.

 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

Nel caso in cui ci siano più righe in un record TITLE , allora il titolo deve essere concatenato, ordinandolo per un numero di continuazione, che viene posizionato, allineato a destra, sui byte 9 e 10 (a seconda del numero di queste righe).

ATOM

I dati memorizzati nelle linee ATOM sono dati di coordinate per ciascun atomo in un esperimento. A volte un esperimento ha inserimenti, mutazioni, posizioni alternative o più modelli. Ciò si traduce nella ripetizione dello stesso atomo più volte. La scelta degli atomi giusti sarà spiegata più avanti.

 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

L'esempio sopra è tratto dall'esperimento 1BAH . La prima colonna indica il tipo di record e la seconda colonna è il numero di serie dell'atomo. Ogni atomo in una struttura ha il suo numero di serie.

Accanto al numero di serie c'è l'etichetta della posizione dell'atomo, che occupa quattro byte. Da quella posizione dell'atomo è possibile estrarre il simbolo chimico dell'elemento, che non sempre è presente nei dati del record in una propria colonna separata.

Dopo il nome dell'atomo c'è un codice residuo di tre lettere. Nel caso delle proteine, quel residuo corrisponde a un amminoacido. Successivamente, la catena è codificata con una lettera. Per catena si intende una singola catena di amminoacidi, con o senza gap, anche se a volte possono essere assegnati ligandi a una catena; questo caso è rilevabile attraverso lacune molto grandi in una sequenza di amminoacidi, che si trova nella colonna successiva. A volte la stessa struttura può essere scansionata con mutazioni incluse, nel qual caso il codice di inserimento è disponibile in una colonna extra dopo la colonna della sequenza. Il codice di inserimento contiene una lettera che aiuta a distinguere quale residuo è interessato.

Le tre colonne successive sono le coordinate spaziali di ciascun atomo, misurate in Angstrom (Å). Accanto a queste coordinate c'è la colonna di occupazione, che dice qual è la probabilità che l'atomo si trovi in ​​quel punto, sulla solita scala da zero a uno.

La penultima colonna è il fattore di temperatura, che contiene informazioni sul disordine nel cristallo, misurato in Ų. Un valore maggiore di 60Ų indica disordine, mentre uno inferiore a 30Ų indica fiducia. Non è sempre presente nei file PDB perché dipende dal metodo sperimentale.

Le colonne successive, simbolo e carica, di solito mancano. Il simbolo chimico può essere raccolto dalla colonna della posizione dell'atomo, come accennato in precedenza. Quando la carica è presente, al simbolo viene suffisso un numero intero seguito da + o - , ad esempio N1+ .

TER

Questo segna la fine della catena. Anche senza questa linea, è facile distinguere dove finisce una catena. Pertanto, spesso manca la linea TER .

MODEL e ENDMDL

Una riga MODEL indica dove inizia il modello di una struttura e contiene il numero di serie del modello.

Dopo tutte le linee atomiche in quel modello, termina con una linea ENDMDL .

SSBOND

Queste linee contengono legami disolfuro tra gli amminoacidi della cisteina. I legami disolfuro possono essere presenti in altri tipi di residui, ma in questo articolo verranno analizzati solo gli amminoacidi, quindi è inclusa solo la cisteina. L'esempio seguente è tratto dall'esperimento con il codice 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

In questo esempio ci sono quattro legami disolfuro contrassegnati nel file con il loro numero di sequenza nella seconda colonna. Tutti questi legami utilizzano la cisteina (colonne 3 e 6) e tutti sono presenti nella catena A (colonne 4 e 7). Dopo ogni catena c'è un numero di sequenza del residuo che indica la posizione del legame nella catena peptidica. I codici di inserimento sono accanto a ciascuna sequenza di residui, ma in questo esempio non sono presenti perché in quella regione non è stato inserito alcun amminoacido. Le due colonne prima dell'ultima sono riservate alle operazioni di simmetria e l'ultima colonna è la distanza tra gli atomi di zolfo, misurata in Å.

Prendiamoci un momento per dare un contesto a questi dati.

Le immagini sottostanti, scattate utilizzando il visualizzatore NGL di rcsb.org, mostrano la struttura dell'esperimento 132L . In particolare, mostrano una proteina senza ligandi. La prima immagine utilizza una rappresentazione a bastoncino, con la colorazione CPK che mostra gli zolfi ei loro legami in giallo. Le connessioni di zolfo a forma di V rappresentano connessioni di metionina, mentre le connessioni a forma di Z sono legami disolfuro tra le cisteine.

Nella figura successiva, un metodo semplificato di visualizzazione delle proteine ​​chiamato visualizzazione della spina dorsale è mostrato colorato da aminoacidi, dove le cisteine ​​sono gialle. Rappresenta la stessa proteina, con la sua catena laterale esclusa e solo una parte del suo gruppo peptidico incluso, in questo caso, la spina dorsale della proteina. È composto da tre atomi: N-terminale, C-alfa e C-terminale. Questa immagine non mostra i legami disolfuro, ma è semplificata per mostrare la disposizione spaziale della proteina:

I tubi sono creati unendo atomi legati al peptide con un atomo di C-alfa. Il colore della cisteina è lo stesso del colore dello zolfo nel metodo di colorazione CPK. Quando i cistene si avvicinano abbastanza, i loro zolfi creano legami disolfuro e questo rafforza la struttura. Altrimenti questa proteina si legherebbe troppo e la sua struttura sarebbe meno stabile a temperature più elevate.

CONECT

Questi record vengono utilizzati per contrassegnare le connessioni tra atomi. A volte questi tag non sono affatto presenti, mentre altre volte vengono inseriti tutti i dati. Nel caso dell'analisi dei legami disolfuro, questa parte può essere utile, ma non è necessaria. Questo perché, in questo progetto, i legami non contrassegnati vengono aggiunti misurando le distanze, quindi questo sarebbe un sovraccarico e deve anche essere controllato.

SOURCE

Questi record contengono informazioni sull'organismo di origine da cui è stata estratta la molecola. Contengono sottorecord per una più facile localizzazione nella tassonomia e hanno la stessa struttura multilinea che abbiamo visto con i record del titolo.

 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

Formato NR

Questo è un elenco di set PDB a catena non ridondanti (NR). Le sue istantanee possono essere trovate su ftp.ncbi.nih.gov/mmdb/nrtable/. Il suo scopo è evitare inutili pregiudizi causati dalla somiglianza delle proteine. NR ha tre insiemi con diversi livelli di valore p dell'identità creati dal confronto di tutte le strutture PDB. Il risultato viene aggiunto ai file di testo che verranno spiegati in seguito. Non tutte le colonne sono necessarie per questo progetto, quindi verranno spiegate solo quelle importanti.

Le prime due colonne contengono il codice univoco dell'esperimento PDB e l'identificatore di catena come spiegato per i record ATOM sopra. Le colonne 6, 9 e C contengono informazioni sulla rappresentatività del valore p, che è il livello di somiglianza delle sequenze calcolate da BLAST. Se quel valore è zero, allora non è accettato come parte di un insieme; se il valore è 1, allora lo è. Le colonne menzionate rappresentano l'accettazione di insiemi con valori di p di 10e-7, 10e-40 e 10e-80, rispettivamente. Per l'analisi verranno utilizzati solo gli insiemi con un valore di p di 10e-7.

L'ultima colonna contiene informazioni sull'accettabilità di una struttura, dove a è accettabile e n non lo è.

 #--------------------------------------------------------------------------------------------------------------------------- # 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

Costruzione di database e analisi dei dati

Ora che abbiamo un'idea di cosa abbiamo a che fare e di cosa dobbiamo fare, iniziamo.

Download dei dati

Tutti i dati per questa analisi possono essere trovati a questi tre indirizzi:

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

I primi due collegamenti contengono un elenco di archivi. L'ultimo archivio di ciascuno dovrebbe essere utilizzato per evitare problemi derivanti da una mancanza di risoluzione o qualità. Il terzo collegamento contiene direttamente l'archivio della tassonomia più recente.

Analisi dei dati

Di solito l'analisi dei file PDB viene eseguita da plugin o moduli in Java, Perl o Python. Nel caso di questa ricerca, ho scritto un'applicazione Perl personalizzata senza utilizzare un modulo di analisi PDB pre-scritto. Il motivo è che quando si analizza una grande quantità di dati, secondo la mia esperienza, il problema più comune con l'utilizzo di dati sperimentali sono gli errori nei dati. A volte ci sono errori con coordinate, distanze, lunghezze delle linee, commenti in punti in cui non dovrebbero essere, ecc.

Il modo più efficace per affrontare questo problema è archiviare inizialmente tutto nel database come testo non elaborato. I parser comuni vengono scritti per gestire dati ideali che sono completamente conformi alle specifiche. Ma in pratica, i dati non sono l'ideale, e questo sarà spiegato nella sezione dei filtri dove troverai lo script di importazione Perl.

Costruzione di database

Quando si costruisce il database, tenere presente che questo database è creato per l'elaborazione dei dati. L'analisi successiva verrà eseguita in SPSS o R. Per i nostri scopi qui si consiglia di utilizzare PostgreSQL con almeno la versione 8.4.

La struttura della tabella viene copiata direttamente dai file scaricati con poche piccole modifiche. In questo caso, il numero di record è troppo piccolo perché valga la pena dedicare il nostro tempo alla normalizzazione. Come accennato, questo database è monouso: queste tabelle non sono create per essere pubblicate su un sito Web, ma servono solo per elaborare i dati. Una volta terminato, possono essere eliminati o sottoposti a backup come dati supplementari, magari per ripetere il processo da qualche altro ricercatore.

In questo caso, il risultato finale sarà una tabella che potrà poi essere esportata in un file da utilizzare con alcuni strumenti statistici come SPSS o R.

Tabelle

L'estrazione dei dati dai record ATOM deve essere collegata ai record HEADER o TITLE . La gerarchia dei dati è spiegata nella figura seguente.

Poiché questa immagine è una rappresentazione semplificata di un database nella terza forma normale (3NF), per i nostri scopi contiene un sovraccarico. Il motivo: per calcolare la distanza tra gli atomi per il rilevamento del legame disolfuro, dovremmo fare unioni. In questo caso, avremmo un tavolo unito a se stesso due volte, e anche unito a una struttura secondaria e primaria due volte ciascuna, il che è un processo molto lento. Poiché non tutte le analisi richiedono informazioni sulla struttura secondaria, viene proposto un altro schema nel caso in cui sia necessario riutilizzare i dati o analizzare quantità maggiori di legami disolfuro:

I legami disolfuro non sono così frequenti come altri legami covalenti, quindi non è necessario un modello di magazzino, sebbene possa essere utilizzato. Lo schema a stella e la modellazione dimensionale di seguito richiederanno troppo tempo per essere sviluppati e renderanno le query più complesse:

Nei casi in cui tutti i legami devono essere elaborati, allora raccomando lo schema a stella.

(Altrimenti non è necessario, perché i legami disolfuro non sono comuni come altri legami. Nel caso di questo lavoro, il numero di legami disolfuro è vicino a 30.000, che potrebbe essere abbastanza veloce in 3NF, ma ho deciso di elaborarlo tramite una tabella non normalizzata, quindi non è raffigurata qui.)

Il numero totale previsto di tutti i legami covalenti è almeno il doppio del numero di atomi nella struttura terziaria, e in tal caso 3NF sarebbe molto lento, quindi è necessaria la denormalizzazione utilizzando la forma dello schema a stella. In quello schema, alcune tabelle hanno due controlli di chiave esterna, e questo perché viene creato un legame tra due atomi, quindi ogni atomo deve avere il proprio primary_structure_id , atom_name_id e residue_id .

Ci sono due modi per riempire la tabella delle dimensioni d_atom_name : dai dati e da un'altra fonte, il dizionario dei componenti chimici che ho menzionato prima. Il suo formato è simile al formato PDB: sono utili solo le linee RESIDUE e CONECT . Questo perché la prima colonna di RESIDUE contiene un codice residuo di tre lettere e CONECT contiene il nome dell'atomo e le sue connessioni, che sono anche nomi di atomi. Quindi da questo file possiamo analizzare tutti i nomi degli atomi e includerli nel nostro database, anche se ti consiglio di considerare la possibilità che il database contenga nomi di atomi non elencati.

 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

In questo progetto, la velocità di codifica è più rilevante della velocità di esecuzione e del consumo di memoria. Ho deciso di non normalizzare: dopotutto, il nostro obiettivo è generare una tabella con le colonne menzionate nell'introduzione.

In questa parte verranno spiegate solo le tabelle più importanti.

Le tabelle principali sono:

• proteins ​​: Tabella con nomi e codici degli esperimenti.
• ps : tabella della struttura primaria che conterrà sequence , chain_id e code .
• ts : tabella contenente la struttura terziaria/quaternaria estratta da dati grezzi e trasformata in formato record ATOM . Questo verrà utilizzato come tabella di staging e può essere eliminato dopo l'estrazione. Sono esclusi i ligandi.
• sources : L'elenco degli organismi da cui sono stati derivati ​​​​i dati sperimentali.
• tax_names , taxonomy_path , taxonomy_paths : nomi di tassonomia linneani dal database di tassonomia NCBI, utilizzati per ottenere percorsi di tassonomia da organismi elencati nelle sources .
• nr : Elenco di proteine ​​NCBI non ridondanti estratte dal set NR.
• pdb_ssbond : elenco di legami disolfuro in un determinato file PDB.

Filtraggio ed elaborazione dei dati

I dati vengono recuperati dagli snapshot dal repository PDB di RCSB.

Ogni file viene importato in una singola tabella raw_pdb nel nostro database PostgreSQL utilizzando uno script Perl. Lo script utilizza transazioni di 10.000 inserimenti per blocco.

La struttura di raw_pdb è questa:

Lo script di importazione è simile al seguente:

 #!/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;

Dopo che le righe sono state importate, vengono analizzate utilizzando le funzioni che definiremo di seguito.

Dai dati raw_pdb , generiamo le tabelle ts , ps , proteins ​​, sources , sources_organela e ss_bond analizzando i record corrispondenti.

La tabella ps ha tre colonne importanti: chain , length e sequence . La sequenza proteica viene generata utilizzando atomi di C-alfa per ciascuna catena e per ogni residuo ordinato per sequenza di residui, prendendo solo il primo inserimento e la prima posizione alternativa. chain è preso dalla colonna TS.chain e la length è semplicemente la lunghezza precalcolata della stringa di sequence . Poiché questo articolo ha lo scopo di analizzare solo catene singole e connessioni intrachain, le proteine ​​a catena multipla vengono saltate nella nostra analisi qui.

All'interno dei record SSBOND , tutti i legami disolfuro sono archiviati nella tabella pdb_ssbond , che eredita dalla tabella pdb_ssbond_extended . pdb_ssbond si presenta così:

Ho aggiunto anche questi indici:

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

Si presume che la distribuzione dei legami disolfuro prima del cutoff sia gaussiana (senza test con, ad esempio, KS-test), quindi le deviazioni standard sono state calcolate su ciascuna distanza tra le cisteine ​​nella stessa proteina per scoprire l'intervallo di lunghezze di legame consentite e confrontarle al taglio. Il cutoff era lo stesso della media calcolata più-meno tre deviazioni standard. Abbiamo esteso il range per introdurre più possibili legami disolfuro che non erano arruolati nel file PDB nelle righe SSBOND ma che abbiamo inserito noi stessi calcolando le distanze tra i record ATOM . L'intervallo scelto per ssbonds è compreso tra 1,6175344456264 e 2,48801947642267 Å, ​​che è la media (2,05) più-meno quattro deviazioni standard:

 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;

La tabella TS contiene le coordinate di tutti gli atomi, ma verranno utilizzate solo le cisteine, con il loro zolfo chiamato " SG " . Quindi viene creata un'altra tabella di staging con solo atomi di zolfo " SG " per accelerare il processo riducendo il numero di record da cercare. Quando si selezionano solo gli zolfi, il numero di combinazioni è molto inferiore che nel caso di tutti gli atomi: 194.574 dei primi contro 122.761.100 dei secondi. All'interno di questa tabella unita a se stessa, le distanze vengono calcolate utilizzando la distanza euclidea ei risultati vengono importati nella tabella pdb_ssbond ma solo dove la distanza è compresa tra le lunghezze definite calcolate in precedenza. Il motivo per fare questo accelerazione è ridurre il tempo necessario per eseguire nuovamente l'intero processo, a scopo di controllo, mantenendolo entro un giorno.

Per pulire i legami disolfuro, utilizziamo il seguente algoritmo:

• Elimina quando entrambi i lati della connessione puntano allo stesso amminoacido
• Elimina le obbligazioni la cui lunghezza non è compresa tra 1.6175344456264 e 2.48801947642267
• Rimuovi inserzioni
• Rimuovi i legami causati da posizioni atomiche alternative, ma lasciando la prima posizione

Il codice per questo (prendendo la tabella pdb_ssbond come primo argomento) è:

 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;

Con questo, l'insieme di proteine ​​non ridondante viene importato nella tabella nr che viene unita alle tabelle ps e proteins , e gli insiemi vengono contrassegnati ( set7 , set40 e set80 ). Alla fine, in base alla quantità di proteine ​​verrà analizzato solo un set. Le catene non corrispondenti tra PDB e NR vengono rimosse dall'analisi.

Le proteine ​​prive di legami disolfuro sono escluse dalla ricerca, insieme alle proteine ​​che non appartengono a nessun insieme. I dati vengono elaborati con DSSP e vengono esclusi anche questi file che hanno avuto problemi di risoluzione o troppi atomi da elaborare. Solo le proteine ​​con catene singole vengono utilizzate come risultato per l'analisi perché le connessioni interchain non sono state mantenute, sebbene siano facilmente calcolabili dalla tabella ssbond contando il numero di connessioni che hanno catene diverse.

Per le restanti proteine, viene eseguito un aggiornamento per il numero totale di legami e il numero di legami sovrapposti, e questo viene fatto per ciascuno degli insiemi.

L'organismo di origine viene estratto dai record SOURCE . Nei casi in cui è sconosciuto, sintetico, progettato, artificiale o ibrido, viene scartato dalla ricerca. Anche le proteine ​​a bassa risoluzione sono escluse solo quando la loro catena laterale non è visibile.

I record SOURCE vengono archiviati nella tabella delle sources , che contiene le righe della tassonomia. In alcuni casi, la tassonomia è mancante o errata. In questi casi è necessaria la correzione manuale degli esperti.

Dalla fonte e dalla tassonomia scaricate da NCBI, la classe viene assegnata a ciascuna struttura primaria. Nel caso in cui non sia possibile assegnare una classe, la proteina viene rimossa dall'elenco di analisi. Sapendo che vengono utilizzati database biologici, si consiglia di eseguire un controllo aggiuntivo di tutti i record di origine e delle classi di tassonomia NCBI da un biologo, altrimenti potrebbero esserci problemi con le classificazioni tra domini diversi. Per abbinare le posizioni cellulari di origine con gli ID tassonomia, i dati dalla tabella di origine vengono estratti nella tabella sources_organela in cui tutti i dati vengono scritti come codice, tag e valore. Il suo formato è mostrato di seguito:

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

L'archivio della tassonomia che vogliamo utilizzare è un file ZIP contenente sette file dump. Tra questi file, due dei più importanti sono names.dmp e merged.dmp . Entrambi i file sono file CSV delimitati da tab-pipe come dettagliato nella documentazione:

• Il file merged.dmp contiene un elenco di ID tassonomia precedenti e gli ID tassonomia correnti in cui ciascuno è stato unito.
• names.dmp contiene queste colonne importanti in questo ordine:
  • tax_id : l'ID tassonomia
  • name_txt : il nome della specie e, se applicabile, il nome univoco (se le specie possono essere trovate con più nomi, questo aiuta a chiarire le ambiguità)
• division.dmp contiene i nomi dei domini di primo livello che useremo come nostre classi.
• nodes.dmp è la tabella che contiene una struttura gerarchica di organismi che utilizzano ID tassonomia.
  • Contiene un ID tassonomia padre, una chiave esterna che può essere trovata in names.dmp .
  • Contiene anche un ID divisione che è importante per noi per archiviare correttamente i dati del dominio principale pertinenti.

Con tutti questi dati e le correzioni manuali (impostando domini di vita corretti) siamo stati in grado di creare la struttura della tabella taxonomy_path . Un campione di dati si presenta così:

 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 come intermediario di elaborazione

In questo lavoro, abbiamo mostrato come elaborare i dati, dal recupero all'analisi. Quando si lavora con dati scientifici, a volte è necessaria la normalizzazione, a volte no. Quando si lavora con piccole quantità di dati che non verranno riutilizzate per un'altra analisi, è sufficiente lasciarli denormalizzati dove l'elaborazione è abbastanza veloce.

Il motivo per cui tutto ciò è stato fatto in un database bioinformatico è che PostgreSQL è in grado di integrare molti linguaggi. Ciò include R, che può eseguire analisi statistiche direttamente nel database, argomento di un futuro articolo sugli strumenti di bioinformatica.

Un ringraziamento speciale ai colleghi di Toptal Stefan Fuchs e Aldo Zelen per le loro preziose consulenze.