Développement d'une base de données bioinformatique pour la recherche sur les liaisons disulfure

La base de données bioinformatiques Protein Data Bank (PDB) est le plus grand référentiel au monde de structures de protéines, d'acides nucléiques et d'assemblages complexes déterminés expérimentalement. Toutes les données sont recueillies à l'aide de méthodes expérimentales telles que les rayons X, la spectroscopie, la cristallographie, la RMN, etc.

Cet article explique comment extraire, filtrer et nettoyer les données de la PDB. Ceci, à son tour, permet le type d'analyse expliqué dans l'article Occurrence of protein disulfide bonds in different domains of life: a comparison of protein from the Protein Data Bank, publié dans Protein Engineering, Design and Selection , Volume 27, Issue 3, 1er mars 2014, p. 65–72.

Le PDB a beaucoup de structures répétitives avec différentes résolutions, méthodes, mutations, etc. Faire une expérience avec des protéines identiques ou similaires peut produire un biais dans toute analyse de groupe, nous devrons donc choisir la structure correcte parmi n'importe quel ensemble de doublons . Pour cela, nous devons utiliser un ensemble de protéines non redondant (NR).

Aux fins de normalisation, je recommande de télécharger le dictionnaire des composés chimiques pour importer des noms d'atomes dans une base de données qui utilise 3NF ou utilise un schéma en étoile et une modélisation dimensionnelle. (J'ai également utilisé DSSP pour aider à éliminer les structures problématiques. Je ne couvrirai pas cela dans cet article, mais notez que je n'ai utilisé aucune autre fonctionnalité DSSP.)

Les données utilisées dans cette recherche contiennent des protéines à une seule unité qui contiennent au moins une liaison disulfure provenant de différentes espèces. Pour effectuer une analyse, toutes les liaisons disulfure sont d'abord classées comme consécutives ou non consécutives, par domaine (archaea, procaryote, viral, eucaryote ou autre), et également par longueur.

Sortir

Pour être prêt à être entré dans R, SPSS ou un autre outil, un analyste aura besoin que les données soient dans une table de base de données avec cette structure :

Matériels et méthodes

Afin d'atteindre cet objectif, la première étape consiste à recueillir des données à partir de rcsb.org. Ce site contient des structures PDB téléchargeables d'expériences dans divers formats.

Bien que les données soient stockées dans plusieurs formats, dans cet exemple, seul le format texte délimité par un espace fixe (PDB) formaté sera utilisé. Une alternative au format textuel PDB est sa version XML, PDBML, mais il contient parfois des entrées de nommage de position d'atome malformées, ce qui peut causer des problèmes pour l'analyse des données. Les anciens formats mmCIF et autres peuvent également être disponibles, mais ils ne seront pas expliqués dans cet article.

Le format APB

Le format PDB est un format textuel fragmenté à largeur fixe qui peut être facilement analysé par des requêtes SQL, des plugins Java ou des modules Perl, par exemple. Chaque type de données dans le conteneur de fichiers est représenté sous la forme d'une ligne commençant par la balise appropriée. Nous passerons en revue chaque type de balise dans les sous-sections suivantes. La longueur de la ligne est inférieure ou égale à 80 caractères, où une étiquette prend six caractères ou moins plus un ou plusieurs espaces qui occupent ensemble huit octets. Il existe également des cas sans espaces entre les balises et les données, généralement pour les balises CONECT .

TITLE

La balise TITLE marque une ligne comme étant (partie de) le titre de l'expérience, contenant le nom de la molécule et d'autres données pertinentes telles que l'insertion, la mutation ou la suppression d'un acide aminé spécifique.

 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

Dans le cas où il y a plusieurs lignes dans un enregistrement TITLE , le titre doit être concaténé, classé par un numéro de continuation, qui est placé, aligné à droite, sur les octets 9 et 10 (selon le nombre de ces lignes).

ATOM

Les données stockées dans les lignes ATOM sont des données de coordonnées pour chaque atome dans une expérience. Parfois, une expérience comporte des insertions, des mutations, des emplacements alternatifs ou plusieurs modèles. Cela se traduit par la répétition du même atome plusieurs fois. Le choix des bons atomes sera expliqué plus tard.

 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'exemple ci-dessus est tiré de l'expérience 1BAH . La première colonne marque le type d'enregistrement et la deuxième colonne est le numéro de série de l'atome. Chaque atome d'une structure a son propre numéro de série.

À côté du numéro de série se trouve l'étiquette de position de l'atome, qui prend quatre octets. À partir de cette position d'atome, il est possible d'extraire le symbole chimique de l'élément, qui n'est pas toujours présent dans les données d'enregistrement dans sa propre colonne séparée.

Après le nom de l'atome, il y a un code de résidu à trois lettres. Dans le cas des protéines, ce résidu correspond à un acide aminé. Ensuite, la chaîne est codée avec une lettre. Par chaîne, nous entendons une seule chaîne d'acides aminés, avec ou sans lacunes, bien que parfois des ligands puissent être attribués à une chaîne ; ce cas est détectable à travers de très grandes lacunes dans une séquence d'acides aminés, qui se trouve dans la colonne suivante. Parfois, la même structure peut être scannée avec des mutations incluses, auquel cas le code d'insertion est disponible dans une colonne supplémentaire après la colonne de séquence. Le code d'insertion contient une lettre pour aider à distinguer quel résidu est affecté.

Les trois colonnes suivantes sont les coordonnées spatiales de chaque atome, mesurées en Angströms (Å). À côté de ces coordonnées se trouve la colonne d'occupation, qui indique la probabilité que l'atome se trouve à cet endroit, sur l'échelle habituelle de zéro à un.

L'avant-dernière colonne est le facteur de température, qui contient des informations sur le désordre dans le cristal, mesuré en Ų. Une valeur supérieure à 60Ų signifie un désordre, tandis qu'une valeur inférieure à 30Ų signifie la confiance. Il n'est pas toujours présent dans les fichiers PDB car il dépend de la méthode expérimentale.

Les colonnes suivantes - symbole et charge - sont généralement manquantes. Le symbole chimique peut être recueilli à partir de la colonne de position de l'atome, comme nous l'avons mentionné ci-dessus. Lorsque la charge est présente, elle est suffixée au symbole sous la forme d'un nombre entier suivi de + ou - , par exemple N1+ .

TER

Cela marque la fin de la chaîne. Même sans cette ligne, il est facile de distinguer où se termine une chaîne. Ainsi, souvent la ligne TER est absente.

MODEL et ENDMDL

Une ligne MODEL marque le début du modèle d'une structure et contient le numéro de série du modèle.

Après toutes les lignes atomiques de ce modèle, il se termine par une ligne ENDMDL .

SSBOND

Ces lignées contiennent des liaisons disulfure entre les acides aminés de la cystéine. Des liaisons disulfure peuvent être présentes dans d'autres types de résidus, mais dans cet article, seuls les acides aminés seront analysés, donc seule la cystéine est incluse. L'exemple suivant est issu de l'expérience avec le code 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

Dans cet exemple, il y a quatre liaisons disulfure étiquetées dans le fichier avec leur numéro de séquence dans la deuxième colonne. Toutes ces liaisons utilisent la cystéine (colonnes 3 et 6), et toutes sont présentes dans la chaîne A (colonnes 4 et 7). Après chaque chaîne, il y a un numéro de séquence de résidus indiquant la position de la liaison dans la chaîne peptidique. Les codes d'insertion sont à côté de chaque séquence de résidus, mais dans cet exemple, ils ne sont pas présents car il n'y avait pas d'acide aminé inséré dans cette région. Les deux colonnes avant la dernière sont réservées aux opérations de symétrie, et la dernière colonne est la distance entre les atomes de soufre, mesurée en Å.

Prenons un moment pour donner un peu de contexte à ces données.

Les images ci-dessous, prises à l'aide du visualiseur NGL rcsb.org, montrent la structure de l'expérience 132L . En particulier, ils montrent une protéine sans ligands. La première image utilise une représentation en bâton, avec une coloration CPK montrant les soufres et leurs liaisons en jaune. Les connexions soufre en forme de V représentent les connexions méthionine, tandis que les connexions en forme de Z sont des liaisons disulfure entre les cystéines.

Dans l'image suivante, une méthode simplifiée de visualisation des protéines appelée visualisation du squelette est représentée colorée par les acides aminés, où les cystéines sont jaunes. Il représente la même protéine, avec sa chaîne latérale exclue, et seulement une partie de son groupe peptidique inclus - dans ce cas, le squelette protéique. Il est composé de trois atomes : N-terminal, C-alpha et C-terminal. Cette image ne montre pas de liaisons disulfure, mais elle est simplifiée pour montrer l'arrangement spatial de la protéine :

Les tuyaux sont créés en joignant des atomes liés à des peptides avec un atome C-alpha. La couleur de la cystéine est la même que la couleur du soufre dans la méthode de coloration CPK. Lorsque les cystènes se rapprochent suffisamment, leurs soufres créent des liaisons disulfure, ce qui renforce la structure. Sinon, cette protéine se lierait trop et sa structure serait moins stable à des températures plus élevées.

CONECT

Ces enregistrements sont utilisés pour marquer les connexions entre les atomes. Parfois, ces balises ne sont pas présentes du tout, alors que d'autres fois, toutes les données sont renseignées. Dans le cas de l'analyse des liaisons disulfure, cette partie peut être utile, mais elle n'est pas nécessaire. C'est parce que, dans ce projet, des liaisons non étiquetées sont ajoutées en mesurant les distances, donc ce serait une surcharge et doit également être vérifiée.

SOURCE

Ces enregistrements contiennent des informations sur l'organisme source à partir duquel la molécule a été extraite. Ils contiennent des sous-enregistrements pour une localisation plus facile dans la taxonomie et ont la même structure multiligne que celle que nous avons vue avec les enregistrements de titre.

 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

Il s'agit d'une liste d'ensembles de chaînes PDB non redondantes (NR). Ses instantanés peuvent être trouvés sur ftp.ncbi.nih.gov/mmdb/nrtable/. Son but est d'éviter les biais inutiles causés par la similarité des protéines. NR a trois ensembles avec différents niveaux de valeur p d'identité créés par comparaison de toutes les structures PDB. Le résultat est ajouté à des fichiers texte qui seront expliqués plus tard. Toutes les colonnes ne sont pas nécessaires pour ce projet, donc seules les plus importantes seront expliquées.

Les deux premières colonnes contiennent le code d'expérience PDB unique et l'identifiant de la chaîne, comme expliqué pour les enregistrements ATOM ci-dessus. Les colonnes 6, 9 et C contiennent des informations sur la représentativité de la valeur p, qui est le niveau de similarité des séquences calculé par BLAST. Si cette valeur est zéro, alors il n'est pas accepté de faire partie d'un ensemble ; si la valeur est 1, alors c'est le cas. Les colonnes mentionnées représentent l'acceptation d'ensembles avec des seuils de valeurs p de 10e-7, 10e-40 et 10e-80, respectivement. Seuls les ensembles avec un seuil de valeur p de 10e-7 seront utilisés pour l'analyse.

La dernière colonne contient des informations sur l'acceptabilité d'une structure, où a est acceptable et n ne l'est pas.

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

Construction de la base de données et analyse des données

Maintenant que nous avons une idée de ce à quoi nous avons affaire et de ce que nous devons faire, commençons.

Téléchargement de données

Toutes les données pour cette analyse se trouvent à ces trois adresses :

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

Les deux premiers liens contiennent une liste d'archives. La dernière archive de chacun doit être utilisée pour éviter les problèmes résultant d'un manque de résolution ou de qualité. Le troisième lien contient directement l'archive de taxonomie la plus récente.

Analyse des données

Habituellement, l'analyse des fichiers PDB est effectuée par des plugins ou des modules en Java, Perl ou Python. Dans le cas de cette recherche, j'ai écrit une application Perl personnalisée sans utiliser de module d'analyse PDB pré-écrit. La raison en est que lors de l'analyse d'une grande quantité de données, d'après mon expérience, le problème le plus courant avec l'utilisation de données expérimentales est les erreurs dans les données. Parfois, il y a des erreurs avec les coordonnées, les distances, les longueurs de ligne, les commentaires à des endroits où ils ne devraient pas être, etc.

Le moyen le plus efficace de gérer cela est de tout stocker initialement dans la base de données sous forme de texte brut. Les analyseurs communs sont écrits pour traiter des données idéales qui sont entièrement conformes aux spécifications. Mais en pratique, les données ne sont pas idéales, et cela sera expliqué dans la section filtrage où vous trouverez le script d'importation Perl.

Construction de la base de données

Lors de la construction de la base de données, notez que cette base de données est conçue pour le traitement des données. Une analyse ultérieure sera effectuée dans SPSS ou R. Pour nos besoins ici, il est recommandé d'utiliser PostgreSQL avec au moins la version 8.4.

La structure de la table est directement copiée à partir des fichiers téléchargés avec seulement quelques petites modifications. Dans ce cas, le nombre d'enregistrements est bien trop faible pour que cela vaille la peine de consacrer notre temps à la normalisation. Comme mentionné, cette base de données est à usage unique : ces tables ne sont pas conçues pour être servies sur un site Web, elles ne sont là que pour le traitement des données. Une fois cette opération terminée, elles peuvent être supprimées ou sauvegardées en tant que données supplémentaires, peut-être pour répéter le processus par un autre chercheur.

Dans ce cas, le résultat final sera un tableau qui pourra ensuite être exporté vers un fichier pour être utilisé avec un outil statistique comme SPSS ou R.

les tables

L'extraction de données des enregistrements ATOM doit être connectée aux enregistrements HEADER ou TITLE . La hiérarchie des données est expliquée dans l'image ci-dessous.

Étant donné que cette image est une représentation simplifiée d'une base de données sous la troisième forme normale (3NF), pour nos besoins, elle contient trop de temps système. La raison : pour calculer la distance entre les atomes pour la détection des liaisons disulfure, nous aurions besoin de faire des jointures. Dans ce cas, nous aurions une table jointe à elle-même deux fois, et également jointe à une structure secondaire et primaire deux fois chacune, ce qui est un processus très lent. Étant donné que toutes les analyses n'ont pas besoin d'informations sur la structure secondaire, un autre schéma est proposé au cas où vous auriez besoin de réutiliser des données ou d'analyser de plus grandes quantités de liaisons disulfures :

Les liaisons disulfure ne sont pas aussi fréquentes que les autres liaisons covalentes, donc un modèle d'entrepôt n'est pas nécessaire, bien qu'il puisse être utilisé. Le schéma en étoile et la modélisation dimensionnelle ci-dessous prendront trop de temps à développer, et rendront les requêtes plus complexes :

Dans les cas où toutes les obligations doivent être traitées, je recommande le schéma en étoile.

(Sinon, ce n'est pas nécessaire, car les liaisons disulfure ne sont pas aussi courantes que les autres liaisons. Dans le cas de ce travail, le nombre de liaisons disulfure est proche de 30 000, ce qui peut être assez rapide dans 3NF, mais j'ai décidé de le traiter via une table non normalisée, elle n'est donc pas illustrée ici.)

Le nombre total attendu de toutes les liaisons covalentes est au moins deux fois le nombre d'atomes dans la structure tertiaire, et dans ce cas, 3NF serait très lent, donc une dénormalisation utilisant la forme de schéma en étoile est nécessaire. Dans ce schéma, certaines tables ont deux vérifications de clé étrangère, et c'est parce qu'une liaison est créée entre deux atomes, donc chaque atome doit avoir ses propres primary_structure_id , atom_name_id et residue_id .

Il existe deux façons de remplir la table de dimension d_atom_name : à partir des données et à partir d'une autre source, le dictionnaire des composants chimiques que j'ai mentionné plus tôt. Son format est similaire au format PDB : Seules les lignes RESIDUE et CONECT sont utiles. En effet, la première colonne de RESIDUE contient un code résiduel à trois lettres et CONECT contient le nom de l'atome et ses connexions, qui sont également des noms d'atomes. Ainsi, à partir de ce fichier, nous pouvons analyser tous les noms d'atomes et les inclure dans notre base de données, bien que je vous recommande de prévoir la possibilité que la base de données contienne des noms d'atomes non répertoriés.

 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

Dans ce projet, la vitesse de codage est plus pertinente que la vitesse d'exécution et la consommation de stockage. J'ai décidé de ne pas normaliser - après tout, notre objectif est de générer un tableau avec les colonnes mentionnées dans l'intro.

Dans cette partie, seuls les tableaux les plus importants seront expliqués.

Les principaux tableaux sont :

• proteins : tableau avec noms et codes des expériences.
• ps : Table de structure primaire qui contiendra sequence , chain_id et code .
• ts : Table contenant la structure tertiaire/quaternaire extraite des données brutes et transformée au format d'enregistrement ATOM . Celle-ci sera utilisée comme table intermédiaire et pourra être supprimée après l'extraction. Les ligands sont exclus.
• sources : La liste des organismes à partir desquels les données expérimentales ont été dérivées.
• tax_names , taxonomy_path , taxonomy_paths : noms de taxonomie linnéenne de la base de données de taxonomie NCBI, utilisés pour obtenir des chemins de taxonomie à partir d'organismes répertoriés dans sources .
• nr : Liste des protéines NCBI non redondantes extraites de l'ensemble NR.
• pdb_ssbond : Liste des liaisons disulfure dans un fichier PDB donné.

Filtrage et traitement des données

Les données sont récupérées à partir d'instantanés du référentiel RCSB PDB.

Chaque fichier est importé dans une seule table raw_pdb dans notre base de données PostgreSQL à l'aide d'un script Perl. Le script utilise des transactions de 10 000 insertions par bloc.

La structure de raw_pdb est la suivante :

Le script d'importation ressemble à ceci :

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

Une fois les lignes importées, elles sont analysées à l'aide de fonctions que nous définirons ci-dessous.

À partir des données raw_pdb , nous générons les tables ts , ps , sources_organela , sources , proteins et ss_bond en parsant les enregistrements correspondants.

La table ps a trois colonnes importantes : chain , length et sequence . La séquence protéique est générée à l'aide d'atomes C-alpha pour chaque chaîne et pour chaque résidu classé par séquence de résidus, en ne prenant que la première insertion et le premier emplacement alternatif. chain est tiré de la colonne TS.chain et length est simplement la longueur précalculée de la chaîne de sequence . Étant donné que cet article est destiné à analyser uniquement les chaînes uniques et les connexions intrachaînes, les protéines à chaînes multiples sont ignorées dans notre analyse ici.

Dans les enregistrements SSBOND , toutes les liaisons disulfure sont stockées dans la table pdb_ssbond , qui hérite de la table pdb_ssbond_extended . pdb_ssbond ressemble à ceci :

J'ai également ajouté ces index:

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

On suppose que la distribution des liaisons disulfure avant la coupure est gaussienne (sans test avec, par exemple, le test KS), de sorte que les écarts-types ont été calculés sur chaque distance entre les cystéines dans la même protéine pour découvrir la gamme de longueurs de liaison autorisées et les comparer. à la coupure. Le seuil était le même que la moyenne calculée plus-moins trois écarts-types. Nous avons étendu la plage pour introduire davantage de liaisons disulfure possibles qui n'étaient pas inscrites dans le fichier PDB dans les lignes SSBOND mais que nous avons insérées nous-mêmes en calculant les distances entre les enregistrements ATOM . La plage choisie pour les ssbonds est comprise entre 1,6175344456264 et 2,48801947642267 Å, ​​qui est la moyenne (2,05) plus-moins quatre écarts-types :

 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 table TS contient les coordonnées de tous les atomes, mais seules les cystéines seront utilisées, avec leur soufre nommé " SG " . Ainsi, une autre table intermédiaire contenant uniquement des atomes de soufre " SG " est créée pour accélérer le processus en réduisant le nombre d'enregistrements à rechercher. Lors de la sélection des soufres uniquement, le nombre de combinaisons est bien inférieur à celui de tous les atomes - 194 574 des premiers contre 122 761 100 des seconds. Dans cette table jointe à elle-même, les distances sont calculées à l'aide de la distance euclidienne et les résultats sont importés dans la table pdb_ssbond mais uniquement lorsque la distance est comprise entre les longueurs définies calculées précédemment. La raison de cette accélération est de réduire le temps nécessaire pour exécuter à nouveau l'ensemble du processus - à des fins de vérification - en le gardant en une journée.

Pour nettoyer les liaisons disulfure, nous utilisons l'algorithme suivant :

• Supprimer lorsque les deux côtés de la connexion pointent vers le même acide aminé
• Supprimer les liaisons dont la longueur n'est pas comprise entre 1,6175344456264 et 2,48801947642267
• Supprimer les insertions
• Supprimer les liaisons causées par des emplacements d'atomes alternatifs, mais en laissant le premier emplacement

Le code pour cela (en prenant la table pdb_ssbond comme premier argument) est :

 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;

Avec cela, l'ensemble non redondant de protéines est importé dans la table nr qui est jointe aux tables ps et proteins , et les ensembles sont marqués ( set7 , set40 et set80 ). À la fin, selon la quantité de protéines, un seul ensemble sera analysé. Les chaînes incompatibles entre PDB et NR sont supprimées de l'analyse.

Les protéines sans liaisons disulfure sont exclues de la recherche, ainsi que les protéines qui n'appartiennent à aucun ensemble. Les données sont traitées avec DSSP, et ces fichiers qui avaient des problèmes de résolution ou trop d'atomes à traiter sont également exclus. Seules les protéines à chaînes simples sont utilisées comme résultat pour l'analyse car les connexions interchaînes n'ont pas été conservées, bien qu'elles soient facilement calculées à partir de la table ssbond en comptant le nombre de connexions qui ont des chaînes différentes.

Pour les protéines restantes, une mise à jour est effectuée pour le nombre total de liaisons et le nombre de liaisons qui se chevauchent, et cela est fait pour chacun des ensembles.

L'organisme source est extrait des enregistrements SOURCE . Dans les cas où il est inconnu, synthétique, conçu, artificiel ou hybride, il est écarté de la recherche. Les protéines à faible résolution sont également exclues uniquement lorsque leur chaîne latérale n'est pas visible.

Les enregistrements SOURCE sont stockés dans la table des sources , qui contient des lignes de taxonomie. Dans certains cas, la taxonomie est manquante ou incorrecte. Dans ces cas, la correction manuelle d'experts est nécessaire.

À partir de la source et de la taxonomie téléchargées depuis NCBI, la classe est attribuée à chaque structure primaire. Dans le cas où une classe ne peut pas être attribuée, la protéine est retirée de la liste d'analyse. Sachant que des bases de données biologiques sont utilisées, une vérification supplémentaire de tous les enregistrements source et des classes de taxonomie NCBI est recommandée par un biologiste, sinon il pourrait y avoir des problèmes de classification entre différents domaines. Pour faire correspondre les emplacements cellulaires sources avec les identifiants de taxonomie, les données de la table source sont extraites dans la table sources_organela où toutes les données sont écrites sous forme de code, de balise et de valeur. Son format est présenté ci-dessous :

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

L'archive de taxonomie que nous souhaitons utiliser est un fichier ZIP contenant sept fichiers de vidage. Parmi ces fichiers, deux des plus importants sont names.dmp et merged.dmp . Les deux fichiers sont des fichiers CSV délimités par des tabulations, comme indiqué dans la documentation :

• Le fichier merged.dmp contient une liste des identifiants de taxonomie précédents et les identifiants de taxonomie actuels dans lesquels chacun a été fusionné.
• names.dmp contient ces colonnes importantes dans cet ordre :
  • tax_id : L'identifiant de la taxonomie
  • name_txt : le nom de l'espèce et, le cas échéant, le nom unique (lorsque les espèces peuvent être trouvées avec plusieurs noms, cela aide à lever l'ambiguïté)
• division.dmp contient les noms des domaines de premier niveau que nous utiliserons comme classes.
• nodes.dmp est la table qui contient une structure hiérarchique d'organismes utilisant des identifiants de taxonomie.
  • Il contient un identifiant de taxonomie parent, une clé étrangère qui se trouve dans names.dmp .
  • Il contient également un ID de division qui est important pour nous afin de stocker correctement les données de domaine supérieur pertinentes.

Avec toutes ces données et les corrections manuelles (définir les domaines de vie corrects), nous avons pu créer la structure de la table taxonomy_path . Un échantillon de données ressemble à ceci :

 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 en tant qu'intermédiaire de traitement

Dans ce travail, nous avons montré comment traiter les données, de la récupération à l'analyse. Lorsque vous travaillez avec des données scientifiques, une normalisation est parfois nécessaire, et parfois non. Lorsque l'on travaille avec de petites quantités de données qui ne seront pas réutilisées pour une autre analyse, il suffit alors de les laisser dénormalisées là où le traitement est assez rapide.

La raison pour laquelle tout cela a été fait dans une seule base de données bioinformatique est que PostgreSQL est capable d'intégrer de nombreux langages. Cela inclut R, qui peut effectuer des analyses statistiques directement dans la base de données - le sujet d'un futur article sur les outils bioinformatiques.

Un merci spécial aux collègues de Toptal Stefan Fuchs et Aldo Zelen pour leurs précieuses consultations.