E. coli - Gestion de données non structurées - Applications post-génomiques - Mise en oeuvre

#knitr::opts_chunk$set(echo = TRUE)
options(width = 95) # nombre de colonnes pour l'affichage sur la sorties HTML
# knitr::opts_chunk$set(cache=F)
# knitr::opts_chunk$set(fig.width=6, fig.height=6) # for rstudio
# knitr::opts_chunk$set(fig.width=14, fig.height=14) # for html

Environnement conda/mamba (et variable d’environnement ?) :

mamba create -n gdnsap r-tidyverse r-neo4r r-reticulate r-codetools r-dt r-kableextra bioconductor-stringdb r-rmdformats r-factoextra scipy   igraph py2neo monotonic packaging numpy pandas scikit-learn biopython   r-markdown

Librairie utilisée pour lancer quelques commandes bash ou python sous Rstudio pour générer ce RNotebook.

library(reticulate)

Répertoires qui seront utilisés pour les données, les scripts, etc. :

mkdir downloads      # données d'origine téléchargées
mkdir reference.sets # ensembles de référence générés pour la recherche d'enrichissement
mkdir query.sets     # ensembles requête
mkdir neo4j.import   # répertoire pour les fichiers de données à importer dans neo4j
mkdir neo4j.data     # répertoire où neo4j stocke ses données
mkdir scripts        # scripts de recherche d'enrichissement et de prioritization

Illustration sur un organisme. Constitution d’une base dédiée à un organisme et son exploitation.

1 Approche par ensembles : Enrichissement

Il a été identifié un ou quelques ensembles de gènes d’intérêt chez E. coli. A quels processus biologiques, fonctions moléculaires ou localisations sub-cellulaires peut-on les relier ?

Ensembles :

Utiliser l’interface Web http://silico.biotoul.fr/enrichment/ afin de caractériser ces ensembles.

1.1 Ensembles cibles et script de recherche d’enrichissement

Un premier exercice, relativement simple, va être de reproduire ce type d’analyse. Nous allons utiliser les mots-clés associés aux protéines dans UniProt.

Les étapes vont être les suivantes :

  • récupération du protéome : toutes les protéines avec les annotations qui nous intéressent, notamment les keywords pour commencer
  • reformatage : liste des protéines annotées par chaque keyword
  • écriture d’un script python, qui prend en paramètre une liste de protéines et qui cherchent les keywords enrichis en comparant cette liste aux listes associées à chacun des keywords

1.1.1 Génération des ensembles cibles (UniProt keywords, Interpro domains, EcoCyc pathways, EcoCyc transcription units)

Aller sur UniProt et télécharger le protéome correspondant à E. coli K12 MG1655 ; étapes :
→ sélectionner Species Proteomes puis rechercher Escherichia coli K12 MG1655. Identifier le bon résultat (protéome de référence) et cliquer sur son Proteome ID
→ Cliquer ensuite sur le lien pour télécharger toutes les protéines
→ Dans les paramètres pour le téléchargement, spécifier de télécharger toutes les protéines,
→ au format TSV, → sélectionner les colonnes/champs que nous allons utiliser par la suite : Entry name, Gene names (ordered locus), Gene Ontology IDs, Interpro et Keywords → Puis télécharger le fichier (uniprotkb_proteome_UP000000625_2023_10_11tsv.gz) dans un répertoire destiné à garder les fichiers d’origine récupérés (par exemple downloads)

zcat downloads/uniprotkb_proteome_UP000000625_*.tsv.gz | head | column -ts $'\t' 
## Entry   Entry Name  Gene Names                Gene Names (ordered locus)  Gene Ontology IDs                                                                                                                                                       Keywords                                                                                                                                                                                                   InterPro                                                                                    Organism (ID)  Organism
## A5A616  MGTS_ECOLI  mgtS yneM b4599 JW1527.1  b4599 JW1527.1              GO:0005886; GO:0010350                                                                                                                                                  3D-structure;Cell inner membrane;Cell membrane;Membrane;Reference proteome;Stress response;Transmembrane;Transmembrane helix                                                                                                                                                                           83333          Escherichia coli (strain K12)
## O32583  THIS_ECOLI  thiS thiG1 b4407 JW3955   b4407 JW3955                GO:0000166; GO:0009228; GO:0009229; GO:0018117; GO:0052837; GO:0097163; GO:1902503; GO:1990228                                                                          3D-structure;Nucleotide-binding;Phosphoprotein;Reference proteome;Thiamine biosynthesis;Thioester bond                                                                                                     IPR012675;IPR016155;IPR010035;IPR003749;                                                    83333          Escherichia coli (strain K12)
## P00350  6PGD_ECOLI  gnd b2029 JW2011          b2029 JW2011                GO:0004616; GO:0005829; GO:0006098; GO:0009051; GO:0042802; GO:0042803; GO:0046177; GO:0050661; GO:0097216                                                              3D-structure;Gluconate utilization;NADP;Oxidoreductase;Pentose shunt;Reference proteome                                                                                                                    IPR008927;IPR013328;IPR006114;IPR006113;IPR006115;IPR006184;IPR036291;IPR006183;            83333          Escherichia coli (strain K12)
## P00363  FRDA_ECOLI  frdA b4154 JW4115         b4154 JW4115                GO:0000104; GO:0005829; GO:0005886; GO:0006113; GO:0006974; GO:0009055; GO:0009061; GO:0016020; GO:0044780; GO:0045284; GO:0050660; GO:0071949; GO:0102040              3D-structure;Cell inner membrane;Cell membrane;Direct protein sequencing;Electron transport;FAD;Flavoprotein;Membrane;Nucleotide-binding;Oxidoreductase;Reference proteome;Transport                       IPR003953;IPR036188;IPR003952;IPR037099;IPR015939;IPR005884;IPR030664;IPR027477;IPR014006;  83333          Escherichia coli (strain K12)
## P00370  DHE4_ECOLI  gdhA b1761 JW1750         b1761 JW1750                GO:0004354; GO:0005737; GO:0005829; GO:0006536; GO:0006537; GO:0042802; GO:0097216; GO:1990148                                                                          3D-structure;Direct protein sequencing;NADP;Oxidoreductase;Reference proteome                                                                                                                              IPR046346;IPR006095;IPR006096;IPR006097;IPR033524;IPR014362;IPR036291;IPR033922;            83333          Escherichia coli (strain K12)
## P00393  NDH_ECOLI   ndh b1109 JW1095          b1109 JW1095                GO:0005886; GO:0008137; GO:0009060; GO:0009061; GO:0019646; GO:0030964; GO:0050660; GO:0055070                                                                          Cell inner membrane;Cell membrane;Direct protein sequencing;Electron transport;FAD;Flavoprotein;Membrane;NAD;Oxidoreductase;Reference proteome;Transport;Ubiquinone                                        IPR036188;IPR023753;                                                                        83333          Escherichia coli (strain K12)
## P00448  SODM_ECOLI  sodA b3908 JW3879         b3908 JW3879                GO:0003677; GO:0004784; GO:0005737; GO:0005829; GO:0006801; GO:0006979; GO:0009408; GO:0010447; GO:0016209; GO:0019430; GO:0030145; GO:0042803; GO:0046872; GO:0071291  3D-structure;Direct protein sequencing;Manganese;Metal-binding;Oxidoreductase;Reference proteome                                                                                                           IPR001189;IPR019833;IPR019832;IPR019831;IPR036324;IPR036314;                                83333          Escherichia coli (strain K12)
## P00452  RIR1_ECOLI  nrdA dnaF b2234 JW2228    b2234 JW2228                GO:0004748; GO:0005524; GO:0005829; GO:0005971; GO:0006260; GO:0009185; GO:0009263; GO:0009265; GO:0015949; GO:0042802; GO:0044183                                      3D-structure;Acetylation;Allosteric enzyme;Alternative initiation;ATP-binding;Deoxyribonucleotide synthesis;Direct protein sequencing;Disulfide bond;Nucleotide-binding;Oxidoreductase;Reference proteome  IPR005144;IPR013346;IPR000788;IPR013509;IPR008926;IPR039718;                                83333          Escherichia coli (strain K12)
## P00490  PHSM_ECOLI  malP b3417 JW5689         b3417 JW5689                GO:0005737; GO:0005829; GO:0005980; GO:0008184; GO:0030170; GO:0030980; GO:0031220; GO:0042803; GO:0102250; GO:0102499                                                  3D-structure;Acetylation;Allosteric enzyme;Carbohydrate metabolism;Direct protein sequencing;Glycosyltransferase;Pyridoxal phosphate;Reference proteome;Transferase                                        IPR011833;IPR000811;IPR035090;                                                              83333          Escherichia coli (strain K12)

Pour générer un fichier texte au format simple qui donne pour chaque mot-clé, la liste des protéines associées, nous allons utiliser R et les librairies tidyverse et jsonlite.

1.1.1.1 Choix des identifiants de référence

Un des problèmes souvent rencontrés dans ce type d’analyse est d’arriver à identifier : les génomes et les gènes, protéines, etc. En effet, chaque source de données peut utiliser des méthodes de référencement et d’identification qui lui est propre. Pour le protéome d’E. coli, nous allons utiliser un type d’identifiant très souvent utilisé : les bnumbers.

A partir du fichier téléchargé, vous allez donc le charger sous R et le reformater pour avoir un mapping des identifiants UniProt vers les bnumbers.

library(tidyverse)
## ── Attaching core tidyverse packages ─────────────────────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ───────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Chargement du fichier et son contenu :

uniprot <- read_tsv("downloads/uniprotkb_proteome_UP000000625_2023_10_11.tsv.gz", show_col_types = F)
uniprot %>% head

Nombre de lignes attendu : 4403

A partir de ce tibble, il s’agit d’extraire les bnumbers de la colonne Gene names (ordered locus). La fonction str_extract permet d’extraire un pattern à partir d’une expression régulière (ici, ce sera 'b\\d+') d’une chaîne de caractères.

mapping <- uniprot %>% 
  select(Entry, names = `Gene Names (ordered locus)`) %>%
  mutate(bnumber=str_extract(names, 'b\\d+')) %>%
  select(bnumber, uniprotID=Entry) %>%
  filter(!is.na(bnumber)) %>% # 2023 → P0DQD7 and P0A6D5 are lost (no bnumber)
  arrange(bnumber)
mapping %>% head

Nombre de lignes attendu : 4401

1.1.1.2 Reformatage des données

A présent, il s’agit de générer la liste des protéines associées à chaque mot-clé (colonne Keywords). Pour cela, on ne gardera que les colonnes Entry et Keywords, une jointure avec mapping va permettre de faire le mapping avec les bnumbers. La fonction separate_rows va nous servir à découper la colonne Keywords et répartir les mots-clés sur des lignes du tibble :

keywords = uniprot %>% 
  select(uniprotID=Entry, keyword=Keywords) %>%
  right_join(mapping, by="uniprotID") %>% # right join to remove those without bnumber
  separate_rows(keyword, sep=';') %>%
  select(bnumber, keyword) %>%
  arrange(bnumber)
keywords %>% head

Nombre de lignes attendu : 27506

Maintenant, il ne s’agit plus que de regrouper les bnumbers par mot-clé :

ref_sets = keywords %>% 
  group_by(keyword) %>%
  summarise(count=n(), elements = list(bnumber)) %>%
  ungroup %>%
  filter(count>1)  %>%
  select(id=keyword, desc=count, elements)
ref_sets

Pour utiliser “facilement” ces données avec python, nous allons sauvegarder le tibble au format JSON :

library(jsonlite)
## 
## Attaching package: 'jsonlite'
## The following object is masked from 'package:purrr':
## 
##     flatten
ref_sets %>% 
  toJSON %>% 
  write("reference.sets/uniprot.keywords.sets.json")

Apperçu du contenu :

jq . "reference.sets/uniprot.keywords.sets.json" | head -40
## [
##   {
##     "id": "2Fe-2S",
##     "desc": 26,
##     "elements": [
##       "b0286",
##       "b0724",
##       "b0775",
##       "b0872",
##       "b0873",
##       "b1392",
##       "b1654",
##       "b1684",
##       "b1802",
##       "b1803",
##       "b2236",
##       "b2283",
##       "b2285",
##       "b2525",
##       "b2530",
##       "b2531",
##       "b2538",
##       "b2540",
##       "b2868",
##       "b2881",
##       "b3337",
##       "b3771",
##       "b4063",
##       "b4153",
##       "b4178",
##       "b4367"
##     ]
##   },
##   {
##     "id": "3D-structure",
##     "desc": 1665,
##     "elements": [
##       "b0002",
##       "b0004",
##       "b0006",

1.1.2 Script de recherche d’enrichissement

La recherche d’enrichissement se résume à comparer chacun des ensembles cibles (pour ce premier exemple, constitués des protéines associées à un même mot-clé) au moyen d’un test statistique. Nous utiliserons pour commencer une loi binomiale.

1.1.2.1 Paramètres de la ligne de commande

Ecrire un script python qui prendra en paramètres :

  • la liste ou le nom de fichier contenant les identifiants composant l’ensemble requête query
  • le nom de fichier contenant les ensembles cibles target

On vous fournit le début du script :

#!/usr/bin/env python

import argparse
from os.path import isfile
import json
from scipy.stats import binom, hypergeom

# SCRIPT PARAMETERS
parser = argparse.ArgumentParser(description='Search enriched terms/categories in the provided (gene) set')
parser.add_argument('-q', '--query', required=True, help='Query set.')
parser.add_argument('-t', '--sets', required=True, help='Target sets filename.')
parser.add_argument('-a', '--alpha', required=False, type=float, default=0.05, help='Significance threshold.')
parser.add_argument('-c', '--adjust', required=False, action="store_true", help='Adjust for multiple testing (FDR).')
parser.add_argument('-m', '--measure', required=False, default='binomial', help='Dissimilarity index: binomial (default), hypergeometric, chi2 or coverage. chi2 and coverage are NOT YET IMPLEMENTED')
parser.add_argument('-l', '--limit', required=False, type=int, default=0, help='Maximum number of results to report.')
parser.add_argument('-v', '--verbose', required=False, action="store_true", help='Talk a lot.')
param = parser.parse_args()

1.1.2.2 Ensemble requête

Il s’agit ensuite de charger les identifiants de l’ensemble requête :

# LOAD QUERY
text = param.query
query = set()
if isfile(text):
    with open(text) as f:
        content = ' '.join(f.read().split('\n')).split()
        query |= set(content)
else: # parse string
    query |= set(text.split())

if param.verbose:
  print(f'query set: {query}')

1.1.2.3 Ensembles cibles

Les ensembles cibles peuvent être chargés avec la librairie jsonlite (importée au début):

# LOAD REFERENCE SETS
sets = json.loads(open(param.sets).read())
if param.verbose:
    print('first target sets: ', sets[0:2])

1.1.2.4 Test statistique

Afin d’appliquer un test avec la loi binomiale, nous avons besoin de connaître la taille de la population pour calculer la probabilité de succès ou d’échec à chaque tentative (nombre d’éléments pris au hasard dans la population, avec remise).

Fonction disponible dans scipi.stats → https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom.html

Rappel sur la loi binomiale : \(k\) succès obtenus sur \(n\) tentatives, chaque tentative ayant une probabilité \(p\) de succès. Pour chaque ensemble cible \(T\), ayant \(C\) éléments en commun avec l’ensemble requête \(Q\), il s’agit donc de considérer la probabilité d’avoir au moins \(C\) succès en effectuant \(Q\) tirages aléatoires avec remise, ayant chacun une probabilité de \(T\)/\(G\) de succès, \(G\) étant le nombre de gènes/protéines considérés pour le tirage.

1.1.2.4.1 Taille de la population

Pour appliquer le test, il nous faut déterminer la probabilité de succès et donc le nombre total d’identifiants possibles. Ajouter une partie dans le script qui calcule le nombre d’identifiants “disponibles” parmi les ensembles cibles :

# COMPUTE POPULATION SIZE
population = set()
for s in sets:
  # TO DO
1.1.2.4.2 Comparaison des ensembles cibles à l’ensemble requête

Ensuite, il faut compléter la partie TO DO ci-dessous pour stocker la p-valeur obtenue pour chaque test.

# EVALUATE SETS
results = []
query_size = len(query)
for s in sets:
    elements = set(s['elements' ])
    common_elements = elements.intersection( query )
    if param.measure=='binomial': # par défaut binom.cdf(<=success, attempts, proba). il nous faut calculer p(X>=x)
        pvalue = 100000 # TO DO 
    else:
        print(f'sorry, {param.measure} not (yet) implemented')
        exit(1)
    r = { 'id': s['id'], 'desc': s['desc'], 'common.n':len(common_elements), 'target.n': len(elements), 'p-value': pvalue, 'elements.target': elements, 'elements.common': common_elements }
    results.append( r )
if param.verbose:
  print(results)

1.1.2.5 Affichage des résultats significatifs

Une fois l’ensemble des ensembles cibles comparés et évalués, il s’agit d’afficher uniquement les résultats statistiquement significatifs :

# PRINT SIGNIFICANT RESULTS
results.sort(key=lambda an_item: an_item['p-value'])
for r in results:
    if r['p-value'] > param.alpha: 
      break
    # OUTPUT
    print("{}\t{}\t{}/{}\t{}\t{}".format( r['id'], r['p-value'], r['common.n'], r['target.n'], r['desc'], ', '.join(r['elements.common'])))

Exemple d’utilisation avec set.01.txt :

./scripts/blastsets.json.py -q query.sets/set.01.txt -t reference.sets/uniprot.keywords.sets.json  | column -ts $'\t'
## Amino-acid biosynthesis                 8.294056706277208e-115  102/104  104   b2020, b0273, b0674, b3940, b2414, b0388, b3671, b0073, b1704, b0077, b3941, b4243, b3390, b3359, b3774, b0159, b1261, b3770, b2024, b4024, b4388, b0071, b3281, b0908, b3959, b3672, b3829, b2600, b1263, b2838, b1622, b2019, b2601, b1260, b3773, b2913, b3213, b2818, b3809, b0078, b0242, b2018, b3237, b3772, b1275, b2026, b2599, b4013, b2472, b2413, b0033, b0001, b3960, b0243, b0003, b3607, b3828, b2763, b3769, b2022, b0072, b1264, b3172, b3744, b4019, b2478, b2023, b3938, b3771, b3389, b3958, b3939, b0002, b0075, b0754, b0074, b1693, b4254, b2551, b3766, b1265, b0032, b3212, b0386, b4488, b2021, b0529, b2329, b4054, b3957, b0907, b3433, b3670, b0166, b2421, b3008, b0031, b0261, b2764, b0004, b2025, b1262
## Amino-acid transport                    1.4518272617403864e-85  79/82    82    b2578, b1798, b3653, b0199, b2679, b3709, b3161, b0861, b3454, b3460, b3110, b0874, b0653, b3458, b2677, b2131, b1918, b1533, b0652, b1473, b2308, b0810, b0692, b4115, b3089, b3455, b3456, b3116, b3268, b1492, b0860, b2678, b0401, b1907, b3457, b1296, b0402, b2663, b0655, b1917, b2130, b0864, b3270, b3370, b0862, b0197, b3795, b4141, b1605, b0576, b0654, b4077, b0809, b0198, b2307, b3269, b1920, b2670, b1729, b4132, b2129, b2365, b0486, b0813, b1336, b1453, b0811, b2310, b0112, b2128, b4208, b2014, b0863, b2156, b2306, b2923, b2309, b1015, b3271
## Aromatic amino acid biosynthesis        1.470650209142406e-20   19/19    19    b3281, b1693, b1265, b0908, b2600, b0388, b1264, b4054, b1704, b2329, b1263, b2601, b1260, b3390, b3389, b1261, b2599, b1262, b0754
## Branched-chain amino acid biosynthesis  1.470650209142406e-20   19/19    19    b0074, b0071, b3766, b3672, b3769, b4488, b0072, b3671, b0073, b0077, b4243, b3773, b3774, b3771, b3770, b3670, b0078, b3772, b0075
## Transport                               5.9067056007624344e-18  81/755   755   b2578, b1798, b3653, b0199, b2679, b3709, b3161, b0861, b3454, b3460, b3110, b0874, b0653, b3458, b2677, b2131, b1918, b1533, b0652, b1473, b2308, b0810, b0692, b4115, b3089, b3455, b3456, b3116, b3268, b1492, b0860, b2678, b0401, b1907, b3457, b1296, b0402, b2413, b2663, b0655, b1917, b2130, b0864, b3270, b3370, b0862, b0197, b3795, b4141, b1605, b0576, b0654, b4077, b0809, b0198, b2307, b3269, b1920, b2670, b1729, b4132, b2129, b2365, b0486, b0813, b1336, b1453, b0811, b2310, b0112, b2128, b4208, b2014, b0863, b2156, b2306, b2923, b2309, b2764, b1015, b3271
## Methionine biosynthesis                 1.8411880932039623e-15  14/14    14    b4019, b3829, b3939, b3941, b3008, b0261, b3433, b3938, b3940, b0529, b4013, b0159, b1622, b3828
## Arginine biosynthesis                   1.967682837713794e-13   12/12    12    b3957, b3958, b2818, b0273, b0033, b4254, b3359, b3960, b3237, b3959, b3172, b0032
## Histidine biosynthesis                  2.0855172947473574e-11  10/10    10    b2019, b2020, b2024, b2023, b2022, b2018, b2025, b2021, b2026, b0529
## Lysine biosynthesis                     2.1442313156401837e-10  9/9      9     b2478, b0031, b3809, b3359, b3433, b4024, b2472, b2838, b0166
## Direct protein sequencing               2.4512773987620095e-09  74/940   940   b2663, b0074, b0071, b2020, b0273, b3281, b0033, b1693, b4254, b2551, b1920, b0674, b0243, b0655, b3940, b2679, b3959, b0032, b3212, b2414, b3829, b2763, b1262, b0386, b3769, b0388, b0072, b3089, b0811, b0073, b1264, b0529, b2329, b3172, b1263, b3744, b4019, b0907, b2478, b0077, b2601, b4243, b3460, b1260, b3390, b2023, b2310, b3359, b2913, b3774, b3458, b0860, b3938, b3433, b0863, b3771, b3213, b2156, b1261, b2306, b0166, b2309, b3939, b3008, b0031, b0242, b2764, b0004, b2025, b0002, b1015, b4024, b4013, b0754
## Cysteine biosynthesis                   2.2739095981386226e-08  7/7      7     b2414, b2763, b2421, b2764, b1275, b2413, b3607
## Pyridoxal phosphate                     3.1100476644514014e-08  15/59    59    b1261, b2414, b0907, b3939, b2421, b3008, b2551, b3359, b0004, b3770, b2021, b3772, b4054, b2838, b1622
## Tryptophan biosynthesis                 2.3543720221078866e-07  6/6      6     b1260, b1265, b1264, b1262, b1261, b1263
## Symport                                 2.5393694144588823e-07  14/60    60    b0692, b3653, b1336, b3161, b4132, b1729, b4208, b2014, b3089, b1015, b1296, b4077, b1907, b3116
## Lyase                                   6.433239649528809e-07   23/170   170   b0071, b1693, b3960, b2414, b2022, b0072, b1264, b2329, b2838, b1263, b1622, b2478, b1260, b2023, b3771, b3389, b1261, b3008, b0004, b2025, b3772, b2599, b1262
## Diaminopimelate biosynthesis            2.456742566107881e-06   5/5      5     b2478, b0031, b3433, b2472, b0166
## Leucine biosynthesis                    2.456742566107881e-06   5/5      5     b0074, b0071, b0072, b0075, b0073
## Threonine biosynthesis                  2.456742566107881e-06   5/5      5     b0001, b0004, b0002, b3433, b0003
## NADP                                    5.8304953843610864e-06  15/90    90    b3212, b3958, b2763, b3281, b0386, b0031, b2764, b0243, b3774, b0002, b3433, b3940, b3213, b0529, b2329
## Allosteric enzyme                       4.039027300269103e-05   9/39     39    b2414, b2478, b4254, b0199, b0002, b3772, b4024, b1264, b1263
## Cell inner membrane                     4.568802233070902e-05   64/983   983   b2663, b2578, b3269, b1798, b3653, b0199, b0652, b1473, b2308, b2670, b0810, b1729, b4132, b3709, b1917, b2130, b2365, b0864, b3270, b0486, b0692, b4115, b0813, b1336, b3370, b0862, b3161, b1453, b0861, b3089, b3456, b3795, b3116, b4141, b1605, b0112, b3110, b0874, b0576, b0653, b2128, b4208, b1492, b2014, b2678, b0654, b2677, b4077, b2156, b0401, b1907, b2306, b2923, b0809, b3457, b0198, b2307, b1918, b1015, b3271, b1296, b1533, b0402, b2413
## Leader peptide                          6.032762215227557e-05   6/16     16    b3672, b0001, b3766, b1265, b2018, b0075
## Multifunctional enzyme                  0.00013978810257080063  9/46     46    b2600, b2022, b0002, b3940, b2026, b2599, b0529, b1262, b1263
## Antiport                                0.00014522723523517002  7/27     27    b1798, b4141, b1605, b4115, b0692, b4132, b1492
## Glutamine amidotransferase              0.0001551535184550992   5/12     12    b3212, b2023, b0674, b1263, b0032
## Proline biosynthesis                    0.0002813442221525336   3/3      3     b0386, b0242, b0243
## Serine biosynthesis                     0.0002813442221525336   3/3      3     b0907, b2913, b4388
## Cytoplasm                               0.0004365441569889691   46/686   686   b0074, b0273, b4254, b2551, b0908, b3960, b0243, b0003, b3959, b3607, b3828, b2600, b0386, b2022, b0388, b0073, b4054, b3172, b3744, b3957, b2019, b2478, b0907, b4243, b3773, b3390, b2023, b3359, b3774, b3938, b3389, b0166, b3958, b3939, b2818, b3008, b2024, b0031, b3809, b0242, b3237, b2025, b1275, b2026, b2599, b4013
## Cell membrane                           0.0012588849519370822   65/1123  1123  b2663, b2578, b3269, b1798, b3653, b0199, b0652, b1473, b2308, b2670, b0810, b1729, b4132, b3709, b1917, b2130, b2365, b0864, b3270, b0486, b0692, b4115, b0813, b1336, b3370, b0862, b3161, b0197, b1453, b0861, b3089, b3456, b3795, b3116, b4141, b1605, b0112, b3110, b0874, b0576, b0653, b2128, b4208, b1492, b2014, b2678, b0654, b2677, b4077, b2156, b0401, b1907, b2306, b2923, b0809, b3457, b0198, b2307, b1918, b1015, b3271, b1296, b1533, b0402, b2413
## Aminotransferase                        0.0015335960919950341   5/20     20    b0907, b3359, b2021, b4054, b3770
## Transmembrane helix                     0.0016701525644965663   56/940   940   b2663, b2578, b3269, b1798, b3653, b2308, b1473, b2670, b1729, b0810, b4132, b3709, b2130, b2365, b3270, b0486, b0692, b4115, b0813, b1336, b3370, b0862, b3161, b1453, b0861, b3089, b3456, b3795, b3116, b4141, b1605, b0112, b3110, b0874, b0576, b0653, b2128, b4208, b1492, b2014, b2678, b0654, b2156, b4077, b0401, b1907, b2923, b3457, b0198, b2307, b1918, b1015, b1296, b1533, b0402, b2413
## Transferase                             0.0025021021407739563   38/582   582   b0074, b0273, b4254, b2551, b0908, b0003, b3940, b3959, b3607, b2414, b3829, b3769, b0388, b2021, b3671, b4054, b1704, b1263, b4019, b2019, b0907, b0077, b2601, b3390, b3359, b3770, b3670, b0166, b3939, b2818, b2421, b0078, b0242, b0261, b0002, b4024, b4013, b0754
## Asparagine biosynthesis                 0.003188560885059188    2/2      2     b3744, b0674
## Glutamate biosynthesis                  0.003188560885059188    2/2      2     b3212, b3213
## Isoleucine biosynthesis                 0.003188560885059188    2/2      2     b3772, b4243
## Transmembrane                           0.005765447086765859    56/993   993   b2663, b2578, b3269, b1798, b3653, b2308, b1473, b2670, b1729, b0810, b4132, b3709, b2130, b2365, b3270, b0486, b0692, b4115, b0813, b1336, b3370, b0862, b3161, b1453, b0861, b3089, b3456, b3795, b3116, b4141, b1605, b0112, b3110, b0874, b0576, b0653, b2128, b4208, b1492, b2014, b2678, b0654, b2156, b4077, b0401, b1907, b2923, b3457, b0198, b2307, b1918, b1015, b1296, b1533, b0402, b2413
## Thiamine pyrophosphate                  0.013756919854410646    3/12     12    b3769, b3671, b0077
## Membrane                                0.01796492718671205     65/1253  1253  b2663, b2578, b3269, b1798, b3653, b0199, b0652, b1473, b2308, b2670, b0810, b1729, b4132, b3709, b1917, b2130, b2365, b0864, b3270, b0486, b0692, b4115, b0813, b1336, b3370, b0862, b3161, b0197, b1453, b0861, b3089, b3456, b3795, b3116, b4141, b1605, b0112, b3110, b0874, b0576, b0653, b2128, b4208, b1492, b2014, b2678, b0654, b2677, b4077, b2156, b0401, b1907, b2306, b2923, b0809, b3457, b0198, b2307, b1918, b1015, b3271, b1296, b1533, b0402, b2413
## ATP-binding                             0.02465892142676584     26/422   422   b0199, b0033, b0652, b0674, b0003, b3940, b3959, b1917, b2129, b0864, b0032, b0388, b3455, b3172, b3744, b3454, b2019, b3390, b2677, b2306, b0809, b0242, b0002, b3271, b4024, b2026
## FAD                                     0.026691551197112685    7/70     70    b3212, b0077, b3941, b2764, b4488, b3671, b2329
## One-carbon metabolism                   0.04347117356312915     2/8      8     b2551, b0529

Nombre de résultats significatifs obtenus :

./scripts/blastsets.json.py -q query.sets/set.01.txt -t reference.sets/uniprot.keywords.sets.json | wc -l
## 41
1.1.2.5.1 Ajustement dans le cadre de tests multiples → FDR

Il s’agit maintenant d’ajouter une correction du seuil \(\alpha\) dans le cadre de ctests multiples. La FDR consiste à ajuster ce seuil par rapport aux p-valeurs obtenues : Les p-valeurs \(P_i\) sont triées de manière croissante (\(P_1\) à \(P_m\), \(m\) étant le nombre de tests) et on garde les \(k\) plus petites p-valeurs vérifiant \(P_k \le \frac{k}{m}\alpha\)

Nombre de résultats significatifs obtenus :

./scripts/blastsets.json.py --adjust -q query.sets/set.01.txt -t reference.sets/uniprot.keywords.sets.json | wc -l
## 35

A partir de toutes ces informations et du code fourni, écrire le script python permettant de faire la recherche d’ensembles cibles similaires.

1.1.3 Autres ensembles d’ensembles cibles

Même chose à faire pour

Adaptez les traitements effectués pour générer les ensembles cibles de références pour les mots-clés UniProt, afin de générer les ensembles cibles de références pour les domaines InterPro, GOTerms, les pathways et les TUs.

Une autre source qui peut s’avérer intéressante est la littérature biomédicale relative à chacun·e des gènes/protéines. Pour cela Entrez du NCBI va nous permettre de constituer, pour chaque article référençant au moins 2 gènes d’E. coli, un ensemble cible dans la section suivante.

1.1.3.1 Ensembles de gènes cités dans une même publication (PubMed)

1.1.3.1.1 Ficher d’annotation du génome

A partir du site du NCBI, avec Entrez, rechercher E. coli K-12 MG1655 : https://www.ncbi.nlm.nih.gov/genome/?term=escherichia+coli+k-12+MG1655

Dans les NCBI Datasets obtenus, le problème, encore une fois est d’identifier la bonne version des données. Pour cet exemple, nous utiliserons le génome de référence et son annotation RefSeq ASM584v2 → RefSeq:GCF_000005845.2

Télécharger le ficher au format GBFF. Nous aurons besoin ensuite du module biopython pour analyser le contenu et extraire les gènes codants (pour des protéines).

unzip -l downloads/GCF_000005845.2.zip 
## Archive:  downloads/GCF_000005845.2.zip
##   Length      Date    Time    Name
## ---------  ---------- -----   ----
##      1596  10-06-2023 11:56   README.md
##       420  10-06-2023 11:56   ncbi_dataset/data/data_summary.tsv
##      3203  10-06-2023 11:56   ncbi_dataset/data/assembly_data_report.jsonl
##  11870708  10-06-2023 11:56   ncbi_dataset/data/GCF_000005845.2/genomic.gbff
##       390  10-06-2023 11:56   ncbi_dataset/data/dataset_catalog.json
## ---------                     -------
##  11876317                     5 files

Extract gbff to compressed file

unzip -p downloads/GCF_000005845.2.zip ncbi_dataset/data/GCF_000005845.2/genomic.gbff | gzip > downloads/GCF_000005845.2.gbff.gz

Nous allons ici ne conserver que les GeneID (pour chercher les références bibliographiques) et la localisation des gènes sur le chromosome :

Attention : pour ne pas surcharger le serveur du NCBI (et se faire bloquer), il faut passer par leur API et limiter le nombre de requêtes effectuées par seconde.

1.1.3.1.2 Mapping des identifiants

Encore une fois, il faut choisir les “bons” identifiants parmi ceux utiliser. Ici, nous utiliserons les bnumbers.

Pour extraire le mapping à partir du fichier GenBank, créer le script suivant (gbff.to.bnumber.GeneID.py) dans le répertoire data et l’exécuter :

#!/bin/env python
# conda install python3 biopython requests
# ./gbff.to.bnumber.GeneID.py > mapping.bnumber_ncbi.tsv

import gzip
from Bio import SeqIO, Entrez

gbff = 'GCF_000005845.2_ASM584v2_genomic.gbff.gz'

with gzip.open(gbff, "rt") as handle:
    record = SeqIO.read(handle, "genbank")
    #print(record)
    print('bnumber\tdbname\tdbid')
    for f in record.features:
        if f.type=='CDS':
            #print(f)
            f_loc = f.location
            f_xref = f.qualifiers['db_xref']
            f_id = f.qualifiers['locus_tag']
            f_desc = f.qualifiers['product']
            for i in f_id:
                for j in f_xref:
                    (dbsource, dbid) = j.split(':')
                    dbsource = 'UniProtKB' if dbsource.startswith('UniProtKB') else dbsource
                    print(f'{i}\t{dbsource}\t{dbid}')

On obtient à la sortie le mapping entre identifiants :

./scripts/gbff.to.bnumber.GeneID.py downloads/GCF_000005845.2.gbff.gz > generated.data/mapping.bnumber_ncbi.tsv
head generated.data/mapping.bnumber_ncbi.tsv | column -ts $'\t'
## bnumber  dbname     dbid
## b0001    UniProtKB  P0AD86
## b0001    ASAP       ABE-0000006
## b0001    ECOCYC     EG11277
## b0001    GeneID     944742
## b0002    UniProtKB  P00561
## b0002    ASAP       ABE-0000008
## b0002    ECOCYC     EG10998
## b0002    GeneID     945803
## b0003    UniProtKB  P00547
1.1.3.1.3 Génération des ensembles cibles (PubMed ids)

Pour chaque gène/identifiant/GeneID, l’utilitaire elink d’Entrez a été utilisé pour récupérer les publications relatives à chacun des gènes. Ces utilitaires peuvent être utilisés directement en accédant à une certaine URL (comme ici), ou bien sont disponibles dans des librairies R, python ou encore comme commande linux dans un shell.

Le script suivant a été utilisé pour récupérer, à partir des GeneID, les identifiants PubMed relatifs à une séquence :

Ceci n’est PAS à réaliser → récupérez plutôt l’ensemble des résultats sur silico.

#!/bin/env python

import sys
import datetime
import time
import requests
import os.path
import json

sleep_for = 1
n=0
last_sleep = time.time()
GeneIDs = set()
mapping = {}
with open('mapping.bnumber_ncbi.tsv', 'r') as f:
    for line in f.readlines():
        (bnumber, dbname, dbid) = line.strip().split('\t')
        if dbname == 'GeneID':
            GeneIDs.add(dbid)
            mapping[dbid] = bnumber
            resfile = 'ncbi/pubmed/GeneID.'+dbid+'.PMIDlinks.json'
            if os.path.isfile(resfile):
                sys.stderr.write('%s file %s exists for %s skipped\n' % (datetime.datetime.now(), resfile, dbid))
            else:
                url = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?retmode=json&dbfrom=gene&db=pubmed&id='+dbid
                sys.stderr.write(f'fetching pmid links of {dbid} from {url}\n')
                r = requests.get(url)
                fh = open(resfile, mode='w')
                fh.write(r.text)
                fh.close()
                # timing
                n += 1
                elapsed = (time.time() - last_sleep)
                rps = n / elapsed
                sys.stderr.write("  %3d requests in %2.2f seconds → %2.1f requests per second" % (n, elapsed, rps))
                if rps > 3 : # max 3 requests/second at NCBI = 9 req in 3 sec
                    sys.stderr.write("sleeping for %s seconds\n" % (sleep_for))
                    time.sleep(sleep_for)
                    last_sleep = time.time()
                    n=0
                elif n>=10: # reset counters after 10 requests
                    last_sleep = time.time()
                    n=0

print('bnumber\tPMID')
for dbid in GeneIDs:
    sys.stderr.write("extracing PMIDs for %s\n" % (dbid))
    # extract PMIDs from XML file
    resfile = 'ncbi/pubmed/GeneID.'+dbid+'.PMIDlinks.json'
    f = open(resfile, 'r')
    d = json.load(f)
    for ls in d['linksets']:    
        for lsd in ls['linksetdbs']:
            if lsd['linkname'] == "gene_pubmed":
                for i in lsd['links']:
                    print(f'{mapping[dbid]}\t{i}')

Nous obtenons donc le fichier associant les identifiants de gène aux identifiants PubMed :

curl http://silico.biotoul.fr/enseignement/m2bbs/gdns-ap/2023/bnumber.PMID.tsv > generated.data/bnumber.PMID.tsv
head generated.data/bnumber.PMID.tsv
##   % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
##                                  Dload  Upload   Total   Spent    Left  Speed
## 
  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0
100 1442k  100 1442k    0     0  32.5M      0 --:--:-- --:--:-- --:--:-- 32.7M
## bnumber  PMID
## b4551    16738553
## b4551    9278503
## b4551    16397293
## b3028    16511004
## b3028    11488591
## b3028    2118386
## b3028    15613473
## b3028    7568050
## b3028    17639604

Utiliser la même technique que pour les mots-clés UniProt afin de générer le fichier au format JSON correspondant qui associe à chaque PMID la liste des gènes (bnumbers) cités dans la publications.

1.2 Données supplémentaires sur la localisation des gènes sur le chromosome

Par la suite, nous pourrons aussi être intéressés par la localisation des gènes sur le chromosome.

Effectuez un travail similaire à partir du fichier d’annotation GenBank afin d’extraire la localisation des gènes sur le chromosome et obtenir le tibble suivant :

1.3 Intégration dans Neo4J

Schéma de la base de données
Schéma de la base de données

1.3.1 Neo4j installation

Sites et documentations :

Récupérer une image docker :

podman pull neo4j:4.4.9-community

Remarques :

Démarrage du conteneur à partir du shell :

mkdir neo4j.data neo4j.import
podman run --rm -it --name neo4j \
   -p 7474:7474 \
   -p 7687:7687 \
   -v ./neo4j.data:/data \
   -v ./neo4j.import:/import \
   -e NEO4J_AUTH=neo4j/omybioinfo \
   --user="$(id -u):$(id -g)" \
   neo4j:4.4.9-community

Le processus est au premier plan donc pour arrêter le serveur il faut faire Ctrl + C dans le terminal.

Remarque : au lancement, le container vérifie à qui appartiennent les fichiers et répertoires (notamment /data et /var/lib/neo4j/import) dont il se sert, et les réattribue si besoin. Pour que vous puissiez créer/ajouter des fichiers dans le répertoire import, il faudra avoir les permissions.

Utilisation depuis le navigateur (vérifier le port renseigné lors de la précédente commande) : http://localhost:7474/

Au premier lancement du container, le mot de passe de l’utilisateur neo4j est redéfini omybioinfo (mais pas les suivantes une fois que la base de données est initialisée).

Aller sur http://localhost:7474 une première fois.

Remarque : le driver pour R neo4r ne semble plus maintenu et ne fonctionne pas avec les dernières versions de neo4j. Cela ne remets pas en question la constitution de la base de données son interrogation mais il faudra adapter les commandes pour un autre langage (python ?) ou une autre interface pour alimenter la base de données.

Ouverture d’un shell cypher sur le conteneur en cours (depuis le shell) :

podman exec -it neo4j ./bin/cypher-shell

1.3.2 Alimentation de la base de données

Librairie neo4r https://github.com/neo4j-rstats/neo4r et https://neo4j-rstats.github.io/user-guide/

if (!require('neo4r')) { # client neo4j 
  install.packages('neo4r')
}
## Loading required package: neo4r
library(neo4r)
neodb = neo4j_api$new(
  url = "http://localhost:7474", 
  user = "neo4j", 
  password = "omybioinfo"
)
cypher = function(query, neo4j=neodb, ...) call_neo4j(query=query, con=neo4j, ...)

1.3.2.1 Nettoyage

Attention, cette partie nettoyage permet de supprimer tous les sommets et relations qu’il y a dans la base de données pour recommencer les commandes. Elle ne supprime pas les étiquettes ou types déjà rencontrés (ex: Gene)

Pour supprimer tous les liens et sommets

'MATCH ()-[r]-() DELETE(r)' %>% cypher
'MATCH ()-[r]-() RETURN count(r)' %>% cypher %>% unlist
'MATCH (n) DELETE(n)' %>% cypher
'MATCH (n) RETURN count(n)' %>% cypher %>% unlist

1.3.2.2 Création du graphe

1.3.2.2.1 Genes

Copier le fichier dans le répertoire d’import de neo4j s’il n’est pas au bon endroit.

# cp data/ncbi.bnumber.location.tsv neo4j.import/
head neo4j.import/ncbi.bnumber.location.tsv
## bnumber  rank    strand  begin   end
## b0001    1   +   189 255
## b0002    2   +   336 2799
## b0003    3   +   2800    3733
## b0004    4   +   3733    5020
## b0005    5   +   5233    5530
## b0006    6   -   5682    6459
## b0007    7   -   6528    7959
## b0008    8   +   8237    9191
## b0009    9   +   9305    9893

Load CSV

# IMPORT
"
LOAD CSV WITH HEADERS FROM 'file:///ncbi.bnumber.location.tsv' AS row FIELDTERMINATOR '\t' 
CREATE (n:Gene)
SET n = row,
 n.id = row.bnumber,
 n.organism = toInteger('511145'),
 n.rank = toInteger(row.rank),
 n.strand = row.strand,
 n.begin = toInteger(row.begin),
 n.end = toInteger(row.end)
"  %>% cypher

Vérification

'MATCH (n:Gene) RETURN count(n)' %>% cypher %>% unlist
## count(n).value 
##           4315

Contenu

"MATCH (n:Gene) RETURN n LIMIT 6" %>% cypher %>% .$n

Index sinon l’ajout de liens et les requêtes peuvent être très longs

'CREATE INDEX ON :Gene(id)' %>% cypher
## $error_code
## [1] "Neo.ClientError.Schema.EquivalentSchemaRuleAlreadyExists"
## 
## $error_message
## [1] "An equivalent index already exists, 'Index( id=3, name='index_817871c2', type='GENERAL BTREE', schema=(:Gene {id}), indexProvider='native-btree-1.0' )'."
1.3.2.2.2 Keywords et liens Keyword → Gene

Load CSV

keywords %>% 
  select(keyword) %>%
  unique %>%
  write_csv("neo4j.import/uniprot.keywords.csv")
'LOAD CSV WITH HEADERS FROM "file:///uniprot.keywords.csv" AS row 
CREATE (n:Keyword)
SET n = row,
 n.id = row.keyword
'  %>% cypher
'CREATE INDEX ON :Keyword(id)' %>% cypher

Vérification

'MATCH (n:Keyword) RETURN count(n)' %>% cypher %>% unlist
## count(n).value 
##            385

Liens Keyword → Gene

keywords %>% write_csv("neo4j.import/uniprot.keywords.genes.csv")

Import du CSV

"MATCH (:Keyword)-[r:describes]->(:Gene) DELETE r" %>% cypher
"
LOAD CSV WITH HEADERS FROM 'file:///uniprot.keywords.genes.csv' AS line 
MATCH (k:Keyword),(g:Gene) 
WHERE k.id=line.keyword AND g.id=line.bnumber
WITH k,g 
MERGE (k)-[:describes]->(g)
" %>% cypher

Vérification

"MATCH (:Keyword)-[r:describes]->(:Gene) RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##          27330
1.3.2.2.3 PubMed
pmids = read_tsv("generated.data/bnumber.PMID.tsv", col_types = "cc")
pmids %>% select(PMID) %>% unique %>% write_csv("neo4j.import/ncbi.pmid.csv")
pmids %>% rename(gene_id = bnumber) %>% write_csv("neo4j.import/ncbi.pmid.genes.csv")

Import du fichier CSV

"MATCH (:PubMed)-[r:cites]->(:Gene) DELETE r" %>% cypher
'MATCH (n:PubMed) DELETE n' %>% cypher
'LOAD CSV WITH HEADERS FROM "file:///ncbi.pmid.csv" AS row 
CREATE (n:PubMed)
SET n = row,
 n.id = row.PMID
'  %>% cypher
'CREATE INDEX ON :PubMed(id)' %>% cypher

Vérification

'MATCH (n:PubMed) RETURN count(n)' %>% cypher %>% unlist
## count(n).value 
##          35320

Import du fichier CSV

"MATCH (:PubMed)-[r:cites]->(:Gene) DELETE r" %>% cypher
"
LOAD CSV WITH HEADERS FROM 'file:///ncbi.pmid.genes.csv' AS line 
MATCH (p:PubMed),(g:Gene) 
WHERE p.id=line.PMID AND g.id=line.gene_id
WITH p,g 
MERGE (p)-[:cites]->(g)
" %>% cypher

Vérification

"MATCH (:PubMed)-[r:cites]->(:Gene) RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##         101841

Complétez cette partie afin de charger dans la base de données :

  • Les domaines InterPro,
  • les pathways,
  • les unités de transcription (TUs),

ainsi que leurs liens avec les gènes.

1.3.2.2.4 Autres sommets et relations : TUs, Pathways
1.3.2.2.4.1 Domaines InterPro

Load CSV

Résultats attendus.

Vérification pour les domaines InterPro.

"MATCH (n:InterPro) RETURN count(n)" %>% cypher %>% unlist
## count(n).value 
##           6679

Liens InterPro → Genes

"MATCH (n:InterPro)-[r:harbored_by]->(:Gene) RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##          15401
1.3.2.2.4.2 Transcription Units d’EcoCyc
Mapping aux identifiants bnumber

Load CSV

Vérification des TUs

"MATCH (n:TU) RETURN count(n)" %>% cypher %>% unlist
## count(n).value 
##           3696

Vérification des liens TU → Gene

"MATCH (n:TU)-[r:harbors]->(:Gene) RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##           6841
1.3.2.2.4.3 Pathways d’EcoCyc
pathways = read_tsv('downloads/All-instances-of-Pathways-in-Escherichia-coli-K-12-substr.-MG1655.txt')
pathways %>% 
  select(id=Pathways, name=`Common-Name`, ref_id=`Genes of pathway`) %>%
  unique %>% 
  mutate(organism=511145) %>% 
  write_csv("neo4j.import/ecocyc.pathways.csv")

"MATCH (:Pathway)-[r:requires]->(:Gene) DELETE r" %>% cypher
"MATCH (n:Pathway) DELETE n" %>% cypher
"
LOAD CSV WITH HEADERS FROM 'file:///ecocyc.pathways.csv' AS line 
CREATE (n:Pathway)
SET n = line
"  %>% cypher
"CREATE INDEX ON :Pathway(id)" %>% cypher

pathways %>% 
  select(pathway=Pathways, ref_id=`Genes of pathway`) %>%
  separate_rows(ref_id, sep=' // ') %>%
  left_join(mapping.ecocyc, by = c("ref_id"="ecocyc")) %>% 
  select(pathway, bnumber) %>%
  unique %>% 
  write_csv("neo4j.import/ecocyc.pathways.bnumbers.csv")
"
LOAD CSV WITH HEADERS FROM 'file:///ecocyc.pathways.bnumbers.csv' AS line 
MATCH (n:Pathway),(g:Gene) 
WHERE n.id=line.pathway AND g.id=line.bnumber
WITH n,g 
MERGE (n)-[:requires]->(g)
" %>% cypher

Vérification des pathways

"MATCH (n:Pathway) RETURN count(n)" %>% cypher %>% unlist
## count(n).value 
##            442

Vérification des liens Pathway → Gene

"MATCH (n:Pathway)-[r:requires]->(:Gene) RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##           2745

1.3.3 Adaptation du script pour l’utilisation de neo4j

Modules python nécessaire :

mamba install py2neo monotonic packaging

Créer un nouveau script en adaptant le précédent qui utilisera, non plus les fichiers au format JSON, mais une connexion au serveur Neo4J pour effectuer les recherches d’enrichissement.

Pour pouvoir construire la requête Neo4J qui identifiera les ensembles cibles, nous avons besoin du type de sommets (PubMed, TU, Pathways, …) ainsi que de l’identifiant taxonomique (en prévision de l’ajout d’autres organismes).

Exemple d’utilisation avec set.01.txt :

./scripts/blastsets.neo4j.py -q query.sets/set.01.txt -t Pathway --species 511145 -c
COMPLETE-ARO-PWY    5.38e-20    19/20   superpathway of aromatic amino acid biosynthesis    b1264, b0754, b1704, b1260, b2601, b3390, b1262, b0388, b1693, b4054, b3389, b1263, b2329, b3281, b3770, b0908, b1261, b2600, b2599
P4-PWY  5.55e-19    18/19   superpathway of L-lysine, L-threonine and L-methionine biosynthesis I   b3433, b4013, b0031, b0004, b3939, b0003, b1622, b3359, b4024, b3809, b2472, b2478, b0166, b3829, b4019, b0907, b2838, b3008
PWY0-781    6.26e-17    18/25   aspartate superpathway  b3433, b4013, b0031, b0004, b3939, b0003, b1622, b3359, b4024, b3809, b2472, b2478, b0166, b3829, b4019, b0907, b2838, b3008
BRANCHED-CHAIN-AA-SYN-PWY   6.21e-13    12/13   superpathway of branched chain amino acid biosynthesis  b4054, b0072, b3771, b3770, b3774, b0073, b0078, b0071, b3772, b0074, b3769, b0077
ALL-CHORISMATE-PWY  3.41e-12    19/55   superpathway of chorismate metabolism   b1264, b0754, b1704, b1260, b2601, b3390, b1262, b0388, b1693, b4054, b3389, b1263, b2329, b3281, b3770, b0908, b1261, b2600, b2599
ARGSYN-PWY  6.26e-12    11/12   L-arginine biosynthesis I (via L-ornithine) b2818, b3172, b3359, b4254, b0273, b0033, b3959, b3960, b0032, b3957, b3958
ARO-PWY 2.52e-11    10/10   chorismate biosynthesis I   b3281, b0388, b1693, b0908, b3389, b0754, b2601, b1704, b3390, b2329
DAPLYSINESYN-PWY    2.52e-11    10/10   L-lysine biosynthesis I b0031, b3433, b3359, b4024, b3809, b2472, b0166, b2478, b0907, b2838
ARG+POLYAMINE-SYN   7.63e-10    11/19   superpathway of arginine and polyamine biosynthesis b2818, b3172, b3359, b4254, b0273, b3959, b0033, b3960, b0032, b3957, b3958
METSYN-PWY  2.56e-09    8/8 superpathway of L-homoserine and L-methionine biosynthesis  b3433, b4013, b3939, b1622, b3829, b3940, b4019, b3008

...

1.4 Ajout de Gene Ontology

Récupération de la dernière version au format OBO :

curl http://current.geneontology.org/ontology/go-basic.obo -o downloads/go-basic.obo

Réutilisation du code de la librairie faite en graphes pour générer les fichiers correspondants aux sommets et arcs à charger dans neo4j à partir du script GeneOntology.py

curl http://silico.biotoul.fr/enseignement/m2bbs/gdns-ap/GeneOntology.py > scripts/GeneOntology.py
python scripts/GeneOntology.py downloads/go-basic.obo nodes > neo4j.import/go.nodes.tsv
python scripts/GeneOntology.py downloads/go-basic.obo edges is_a > neo4j.import/go.is_a.edges.tsv
python scripts/GeneOntology.py downloads/go-basic.obo edges part_of > neo4j.import/go.part_of.edges.tsv

Sommets

"MATCH (:GOTerm)-[r]-() DELETE r" %>% cypher
"MATCH (n:GOTerm) DELETE n" %>% cypher
"
LOAD CSV WITH HEADERS FROM 'file:///go.nodes.tsv' AS row FIELDTERMINATOR '\t' 
CREATE (n:GOTerm)
SET n.id = row.id,
n.name = row.desc,
n.desc  = row.def,
n.namespace = row.namespace
"  %>% cypher

Vérification

"MATCH (n:GOTerm) RETURN count(n)" %>% cypher %>% unlist
## count(n).value 
##          42887

Index

"CREATE INDEX ON :GOTerm(id)" %>% cypher
## $error_code
## [1] "Neo.ClientError.Schema.EquivalentSchemaRuleAlreadyExists"
## 
## $error_message
## [1] "An equivalent index already exists, 'Index( id=9, name='index_1a03f871', type='GENERAL BTREE', schema=(:GOTerm {id}), indexProvider='native-btree-1.0' )'."

Extait du contenu :

"MATCH (n:GOTerm) RETURN n LIMIT 6" %>% cypher %>% .$n

Relations is a

"
LOAD CSV WITH HEADERS FROM 'file:///go.is_a.edges.tsv' AS line FIELDTERMINATOR '\t' 
MATCH (t1:GOTerm),(t2:GOTerm) 
WHERE t1.id=line.term1 AND t2.id=line.term2
WITH t1,t2 
MERGE (t2)-[r:is_a]->(t1)
" %>% cypher

Vérification

"MATCH ()-[r:is_a]->() RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##          68521

Relation part of

"
LOAD CSV WITH HEADERS FROM 'file:///go.part_of.edges.tsv' AS line FIELDTERMINATOR '\t' 
MATCH (t1:GOTerm),(t2:GOTerm) 
WHERE t1.id=line.term1 AND t2.id=line.term2
WITH t1,t2 
MERGE (t2)-[r:part_of]->(t1)
" %>% cypher

Vérification

"MATCH ()-[r:part_of]->() RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##           6785

Gene associations

GOTerms = uniprot %>% 
  select(uniprotID=Entry, GOTerm=`Gene Ontology IDs`) %>%
  right_join(mapping, by = join_by(uniprotID)) %>% # right join to remove those without bnumber
  separate_rows(GOTerm, sep='; ') %>%
  select(bnumber, GOTerm) %>%
  arrange(bnumber)
GOTerms

Save CSV file for Neo4J

GOTerms %>% write_csv("neo4j.import/uniprot.GOTerm.bnumber.csv")

Merge into graph

"
LOAD CSV WITH HEADERS FROM 'file:///uniprot.GOTerm.bnumber.csv' AS line
MATCH (t:GOTerm),(g:Gene) 
WHERE t.id=line.GOTerm AND g.id=line.bnumber
WITH t,g
MERGE (t)-[:annotates]->(g)
" %>% cypher

Vérification

"MATCH ()-[r:annotates]->() RETURN count(r)" %>% cypher %>% unlist
## count(r).value 
##          22382

1.4.1 Adaptation du script en conséquence

Adapter le script précédent afin qu’il compare les ensembles cibles correspondants aux GOTerm (en n’oubliant pas les associations GOterm → GeneProduct implicites).

Exemple d’utilisation :

./scripts/blastsets.neo4j.py -q query.sets/set.01.txt -t GOTerm --species 511145 -c
GO:0008652  3.37e-99    101/141 amino acid biosynthetic process b0031, b0166, b2019, b3672, b2599, b1265, b1263, b4254, b3774, b3773, b2601, b3828, b2472, b2020, b0004, b4019, b0908, b0078, b2026, b3771, b3390, b2414, b2024, b0071, b2478, b0388, b3433, b1261, b0003, b4054, b4024, b1262, b3212, b0386, b4013, b0674, b0243, b3829, b2763, b3389, b3213, b3770, b0075, b4243, b0074, b3766, b3744, b3940, b0273, b2818, b0073, b0032, b2838, b2023, b3359, b3008, b2913, b2764, b3938, b3958, b3607, b2600, b1260, b3959, b0077, b0159, b1704, b2022, b3671, b0072, b0001, b0754, b3281, b2329, b0261, b2413, b2018, b3772, b1264, b1275, b0242, b3172, b2551, b1693, b2421, b3237, b3960, b2021, b3809, b0033, b0907, b2025, b0529, b0002, b3769, b3939, b4388, b3670, b1622, b3941, b3957
GO:0046394  2.35e-84    101/200 carboxylic acid biosynthetic process    b0031, b0166, b2019, b3672, b2599, b1265, b1263, b4254, b3774, b3773, b2601, b3828, b0004, b2020, b2472, b4019, b0908, b0078, b2026, b3771, b3390, b2414, b0071, b2024, b2478, b0388, b3433, b0003, b1261, b4054, b4024, b1262, b3212, b0386, b4013, b0674, b0243, b3829, b2763, b3389, b3213, b3770, b0075, b4243, b0074, b3766, b3940, b3744, b0273, b2818, b0073, b0032, b2838, b2023, b3008, b3359, b2764, b2913, b3938, b3958, b3607, b2600, b1260, b3959, b0159, b0077, b1704, b2022, b3671, b0001, b0072, b0754, b3281, b2329, b0261, b2413, b2018, b3772, b0242, b1264, b1275, b3172, b2551, b1693, b2421, b3237, b3960, b2021, b3809, b0033, b0907, b2025, b0529, b0002, b3769, b3939, b3670, b4388, b1622, b3941, b3957
GO:0016053  3.81e-84    101/201 organic acid biosynthetic process   b0031, b0166, b2019, b3672, b2599, b1265, b1263, b4254, b3774, b3773, b2601, b3828, b0004, b2020, b2472, b4019, b0908, b0078, b2026, b3771, b3390, b2414, b0071, b2024, b2478, b0388, b3433, b0003, b1261, b4054, b4024, b1262, b3212, b0386, b4013, b0674, b0243, b3829, b2763, b3389, b3213, b3770, b0075, b4243, b0074, b3766, b3940, b3744, b0273, b2818, b0073, b0032, b2838, b2023, b3008, b3359, b2764, b2913, b3938, b3958, b3607, b2600, b1260, b3959, b0159, b0077, b1704, b2022, b3671, b0001, b0072, b0754, b3281, b2329, b0261, b2413, b2018, b3772, b0242, b1264, b1275, b3172, b2551, b1693, b2421, b3237, b3960, b2021, b3809, b0033, b0907, b2025, b0529, b0002, b3769, b3939, b3670, b4388, b1622, b3941, b3957
GO:0006520  1.32e-78    105/260 amino acid metabolic process    b0031, b0166, b2019, b3672, b2599, b1265, b1263, b4254, b3774, b3773, b2601, b3828, b0004, b2020, b2472, b4019, b0908, b0078, b2663, b2026, b3771, b3390, b2414, b0071, b2024, b2478, b0388, b3433, b0003, b1261, b4054, b4024, b1262, b3212, b0386, b4013, b0674, b0243, b3829, b2763, b3389, b3213, b3770, b0075, b4243, b0074, b3766, b3370, b3940, b3744, b0273, b2818, b0073, b0032, b2838, b2023, b3008, b3359, b2764, b2913, b3938, b3958, b3607, b2600, b1260, b3959, b0159, b0077, b1704, b2022, b3671, b0001, b0072, b0754, b3281, b2329, b0261, b2413, b2018, b3772, b0242, b1264, b1275, b3172, b2551, b1693, b2421, b3110, b3237, b3960, b2021, b3809, b0033, b3709, b0907, b2025, b0529, b0002, b3769, b3939, b3670, b4388, b1622, b3941, b3957
GO:1901607  5.47e-76    81/117  alpha-amino acid biosynthetic process   b0031, b0166, b3672, b2599, b1265, b1263, b4254, b3774, b3828, b2472, b0004, b4019, b0078, b3771, b2414, b2024, b0071, b2478, b3433, b1261, b0003, b4054, b4024, b1262, b3212, b0386, b4013, b0674, b0243, b3829, b2763, b3213, b3770, b0075, b4243, b0074, b3744, b3940, b0273, b2818, b0073, b0032, b2838, b3359, b3008, b2913, b2764, b3938, b3958, b3607, b2600, b1260, b3959, b0077, b0159, b3671, b0072, b0001, b0261, b2413, b3772, b1264, b1275, b0242, b3172, b2551, b2421, b3237, b3960, b3809, b0033, b0907, b0529, b0002, b3769, b3939, b4388, b3670, b1622, b3941, b3957
GO:0006865  1.34e-71    73/93   amino acid transport    b3089, b1473, b0809, b2130, b1729, b0864, b2156, b2678, b1015, b0199, b0813, b2309, b1917, b1336, b0862, b2129, b2014, b0112, b2365, b4115, b2679, b0652, b4132, b2306, b3457, b3110, b0860, b1605, b3795, b0486, b2578, b0811, b0197, b3269, b2670, b2308, b2128, b1453, b3458, b3653, b1533, b1296, b2923, b3456, b0576, b0874, b4208, b3709, b3161, b4141, b1492, b0402, b0654, b3460, b1920, b3270, b0401, b2131, b0653, b2677, b2663, b0861, b1798, b1918, b2310, b3455, b2307, b3116, b0655, b0810, b0198, b0863, b3454
...

1.4.2 Distance sémantique entre GOTerms sur-représentés

Analysez les résultats du script (GOTerm et p-valeurs) avec REVIGO : http://revigo.irb.hr/

Requête sur ReViGO BP scatterplot

BP treemap
BP treemap
CC scatterplot
CC scatterplot
CC treemap
CC treemap
MF scatterplot
MF scatterplot
MF treemap
MF treemap

2 Approches par graphe ou matrice de distances/dissemblances/similarités/probabilités

2.1 Génération de matrices de similarité

A partir des liens Keyword, Pathway, PubMed et GOTerm, il s’agit de générer des matrices pivots individus(bnumber)–variables(Keyword ou pahtways ou …) pour ensuite calculer la distance/similarité entre paires de gènes.

Ajout d’une colonne avec des 1 pour les association présentes dans les données de départ :

keyword.mat = keywords %>% 
  unique %>% 
  mutate(asso=1)
keyword.mat

Génération d’une table pivot (les valeurs de la colonne keyword deviennent des colonnes) :

keyword.tab = keyword.mat %>% 
  pivot_wider(names_from = keyword, values_from = asso, values_fill = 0) %>% 
  as.data.frame
keyword.tab %>% head

Transfert des identifiants de gènes de la colonne bnumber aux noms de lignes (nécessaires pour qu’ils soient reportés sur les noms des lignes/colonnes) de la matrice de distance :

rownames(keyword.tab) = keyword.tab$bnumber
keyword.tab %>% head

Calcul de la matrice de distance (entre paires de lignes/gènes) avec la méthode binary (cf. ?dist) :

keyword.tab = keyword.tab %>% dplyr::select(-bnumber)
keyword.dist = keyword.tab %>% dist(method='binary')
keyword.dist[1:10]
##  [1] 0.8000000 0.6666667 0.6250000 0.8000000 0.8333333 0.9090909 0.9000000 0.8888889 0.9230769
## [10] 0.7500000

Sous forme de scores à la StringDB

keyword.scores = round( (1-keyword.dist) * 1000) %>% as.matrix
keyword.scores[1:10,1:10]
##       b0001 b0002 b0003 b0004 b0005 b0006 b0007 b0008 b0009 b0010
## b0001     0   200   333   375   200   167    91   100   111    77
## b0002   200     0   467   312    67   133    48   235   333    91
## b0003   333   467     0   250   111   100    67   250   400    59
## b0004   375   312   250     0   125   250    71   273   182   133
## b0005   200    67   111   125     0   250   111   125   143    91
## b0006   167   133   100   250   250     0   100   250   286   182
## b0007    91    48    67    71   111   100     0    71    77   800
## b0008   100   235   250   273   125   250    71     0   300   133
## b0009   111   333   400   182   143   286    77   300     0   143
## b0010    77    91    59   133    91   182   800   133   143     0

2.2 Client R à STRINGdb

Client R Bioconductor :

library(STRINGdb)

Autres librairies utiles pour la gestion et l’affichage de tables :

# mise en forme des tableaux:
library(kableExtra) # nicer tables
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
kabex = function(o) o %>% kable(format="html", escape=F) %>% kable_styling(bootstrap_options = c("striped"))
library(DT) 
dted = function(o) o %>% datatable(rownames=F, filter="top", options=list(pageLength = 10), escape=F)

Connexion au serveur de STRING (NCBI taxon id de Ecoli est 511145 : https://string-db.org/cgi/input.pl?input_page_active_form=organism )

confidence_threshold = 333
stringdb = STRINGdb$new(version='12', species=511145, score_threshold=confidence_threshold, input_directory='downloads')

Téléchargement du génome/protéome :

proteins = stringdb$get_proteins() %>% as_tibble
proteins %>% dted

Fichier récupéré

zcat downloads/511145.protein.info.v12.0.txt.gz | head | column -ts $'\t'
## #string_protein_id  preferred_name  protein_size  annotation
## 511145.b0001        thrL            21            Thr operon leader peptide; This protein is involved in control of the biosynthesis of threonine.
## 511145.b0002        thrA            820           Bifunctional: aspartokinase I (N-terminal); homoserine dehydrogenase I (C-terminal); Protein involved in threonine biosynthetic process, methionine biosynthetic process and homoserine biosynthetic process.
## 511145.b0003        thrB            310           Homoserine kinase; Catalyzes the ATP-dependent phosphorylation of L-homoserine to L-homoserine phosphate. Is also able to phosphorylate the hydroxy group on gamma-carbon of L-homoserine analogs when the functional group at the alpha-position is a carboxyl, an ester, or even a hydroxymethyl group. Neither L-threonine nor L-serine are substrates of the enzyme.
## 511145.b0004        thrC            428           L-threonine synthase; Catalyzes the gamma-elimination of phosphate from L- phosphohomoserine and the beta-addition of water to produce L- threonine. To a lesser extent, is able to slowly catalyze the deamination of L-threonine into alpha-ketobutyrate and that of L-serine and 3-chloroalanine into pyruvate. Is also able to rapidly convert vinylglycine to threonine, which proves that the pyridoxal p-quinonoid of vinylglycine is an intermediate in the TS reaction.
## 511145.b0005        yaaX            98            DUF2502 family putative periplasmic protein; To E.coli YpeC.
## 511145.b0006        yaaA            258           Peroxide resistance protein, lowers intracellular iron; Involved in the cellular response to hydrogen peroxide (H(2)O(2)) stress. During H(2)O(2) stress, prevents oxidative damage to both DNA and proteins by diminishing the amount of unincorporated iron within the cell.
## 511145.b0007        yaaJ            476           Putative transporter; Inner membrane transport protein.
## 511145.b0008        talB            317           Transaldolase B; Transaldolase is important for the balance of metabolites in the pentose-phosphate pathway.
## 511145.b0009        mog             195           Molybdochelatase incorporating molybdenum into molybdopterin; Catalyzes the adenylation of molybdopterin as part of the biosynthesis of the molybdenum-cofactor; Belongs to the MoaB/Mog family.

Chargement de set.01.txt pour analyse.

s1 = scan('query.sets/set.01.txt', character()) 
s1
##   [1] "b0001" "b0002" "b0003" "b0004" "b0031" "b0032" "b0033" "b0071" "b0072" "b0073" "b0074"
##  [12] "b0075" "b0077" "b0078" "b0159" "b0166" "b0242" "b0243" "b0261" "b0273" "b0386" "b0388"
##  [23] "b0529" "b0674" "b0754" "b0907" "b0908" "b1260" "b1261" "b1262" "b1263" "b1264" "b1265"
##  [34] "b1275" "b1622" "b1693" "b1704" "b2018" "b2019" "b2020" "b2021" "b2022" "b2023" "b2024"
##  [45] "b2025" "b2026" "b2329" "b2413" "b2414" "b2421" "b2472" "b2478" "b2551" "b2599" "b2600"
##  [56] "b2601" "b2763" "b2764" "b2818" "b2838" "b2913" "b3008" "b3172" "b3212" "b3213" "b3237"
##  [67] "b3281" "b3359" "b3389" "b3390" "b3433" "b3607" "b3670" "b3671" "b3672" "b3744" "b3766"
##  [78] "b3769" "b3770" "b3771" "b3772" "b3773" "b3774" "b3809" "b3828" "b3829" "b3938" "b3939"
##  [89] "b3940" "b3941" "b3957" "b3958" "b3959" "b3960" "b4013" "b4019" "b4024" "b4054" "b4243"
## [100] "b4254" "b4388" "b4488" "b0112" "b0197" "b0198" "b0199" "b0401" "b0402" "b0486" "b0576"
## [111] "b0652" "b0653" "b0654" "b0655" "b0692" "b0809" "b0810" "b0811" "b0813" "b0860" "b0861"
## [122] "b0862" "b0863" "b0864" "b0874" "b1015" "b1296" "b1336" "b1453" "b1473" "b1492" "b1533"
## [133] "b1605" "b1729" "b1798" "b1907" "b1917" "b1918" "b1920" "b2014" "b2128" "b2129" "b2130"
## [144] "b2131" "b2156" "b2306" "b2307" "b2308" "b2309" "b2310" "b2365" "b2578" "b2663" "b2670"
## [155] "b2677" "b2678" "b2679" "b2923" "b3089" "b3110" "b3116" "b3161" "b3268" "b3269" "b3270"
## [166] "b3271" "b3370" "b3454" "b3455" "b3456" "b3457" "b3458" "b3460" "b3653" "b3709" "b3795"
## [177] "b4077" "b4115" "b4132" "b4141" "b4208"

Mapping des identifiants vers ceux de StringDB

s1.mapped = stringdb$mp(s1)
s1.mapped
##   [1] "511145.b0001" "511145.b0002" "511145.b0003" "511145.b0004" "511145.b0031" "511145.b0032"
##   [7] "511145.b0033" "511145.b0071" "511145.b0072" "511145.b0073" "511145.b0074" "511145.b0075"
##  [13] "511145.b0077" "511145.b0078" "511145.b0159" "511145.b0166" "511145.b0242" "511145.b0243"
##  [19] "511145.b0261" "511145.b0273" "511145.b0386" "511145.b0388" "511145.b0529" "511145.b0674"
##  [25] "511145.b0754" "511145.b0907" "511145.b0908" "511145.b1260" "511145.b1261" "511145.b1262"
##  [31] "511145.b1263" "511145.b1264" "511145.b1265" "511145.b1275" "511145.b1622" "511145.b1693"
##  [37] "511145.b1704" "511145.b2018" "511145.b2019" "511145.b2020" "511145.b2021" "511145.b2022"
##  [43] "511145.b2023" "511145.b2024" "511145.b2025" "511145.b2026" "511145.b2329" "511145.b2413"
##  [49] "511145.b2414" "511145.b2421" "511145.b2472" "511145.b2478" "511145.b2551" "511145.b2599"
##  [55] "511145.b2600" "511145.b2601" "511145.b2763" "511145.b2764" "511145.b2818" "511145.b2838"
##  [61] "511145.b2913" "511145.b3008" "511145.b3172" "511145.b3212" "511145.b3213" "511145.b3237"
##  [67] "511145.b3281" "511145.b3359" "511145.b3389" "511145.b3390" "511145.b3433" "511145.b3607"
##  [73] "511145.b3670" "511145.b3671" "511145.b3672" "511145.b3744" "511145.b3766" "511145.b3769"
##  [79] "511145.b3770" "511145.b3771" "511145.b3772" "511145.b3773" "511145.b3774" "511145.b3809"
##  [85] "511145.b3828" "511145.b3829" "511145.b3938" "511145.b3939" "511145.b3940" "511145.b3941"
##  [91] "511145.b3957" "511145.b3958" "511145.b3959" "511145.b3960" "511145.b4013" "511145.b4019"
##  [97] "511145.b4024" "511145.b4054" "511145.b4243" "511145.b4254" "511145.b4388" "511145.b0112"
## [103] "511145.b0197" "511145.b0198" "511145.b0199" "511145.b0401" "511145.b0402" "511145.b0486"
## [109] "511145.b0576" "511145.b0652" "511145.b0653" "511145.b0654" "511145.b0655" "511145.b0692"
## [115] "511145.b0809" "511145.b0810" "511145.b0811" "511145.b0813" "511145.b0860" "511145.b0861"
## [121] "511145.b0862" "511145.b0863" "511145.b0864" "511145.b0874" "511145.b1015" "511145.b1296"
## [127] "511145.b1336" "511145.b1453" "511145.b1473" "511145.b1492" "511145.b1533" "511145.b1605"
## [133] "511145.b1729" "511145.b1798" "511145.b1907" "511145.b1917" "511145.b1918" "511145.b1920"
## [139] "511145.b2014" "511145.b2128" "511145.b2129" "511145.b2130" "511145.b2131" "511145.b2156"
## [145] "511145.b2306" "511145.b2307" "511145.b2308" "511145.b2309" "511145.b2310" "511145.b2365"
## [151] "511145.b2578" "511145.b2663" "511145.b2670" "511145.b2677" "511145.b2678" "511145.b2679"
## [157] "511145.b2923" "511145.b3089" "511145.b3110" "511145.b3116" "511145.b3161" "511145.b3269"
## [163] "511145.b3270" "511145.b3271" "511145.b3370" "511145.b3454" "511145.b3455" "511145.b3456"
## [169] "511145.b3457" "511145.b3458" "511145.b3460" "511145.b3653" "511145.b3709" "511145.b3795"
## [175] "511145.b4077" "511145.b4115" "511145.b4132" "511145.b4141" "511145.b4208"

Plot du graphe correspondant

stringdb$plot_network(s1.mapped)

Enrichment

enrichment = stringdb$get_enrichment(s1.mapped)
enrichment %>% dted

Sommets voisins d’un ou plusieurs sommets donnés

stringdb$get_neighbors(s1.mapped[1:3])
##   [1] "511145.b0002" "511145.b0003" "511145.b0004" "511145.b0005" "511145.b0007" "511145.b0075"
##   [7] "511145.b0379" "511145.b1265" "511145.b1715" "511145.b2018" "511145.b2598" "511145.b2979"
##  [13] "511145.b3672" "511145.b3766" "511145.b4246" "511145.b4403" "511145.b4637" "511145.b0001"
##  [19] "511145.b0003" "511145.b0004" "511145.b0005" "511145.b0026" "511145.b0031" "511145.b0032"
##  [25] "511145.b0033" "511145.b0047" "511145.b0048" "511145.b0071" "511145.b0072" "511145.b0073"
##  [31] "511145.b0074" "511145.b0077" "511145.b0078" "511145.b0082" "511145.b0085" "511145.b0086"
##  [37] "511145.b0088" "511145.b0095" "511145.b0120" "511145.b0121" "511145.b0154" "511145.b0159"
##  [43] "511145.b0166" "511145.b0242" "511145.b0243" "511145.b0261" "511145.b0268" "511145.b0269"
##  [49] "511145.b0273" "511145.b0311" "511145.b0333" "511145.b0352" "511145.b0386" "511145.b0388"
##  [55] "511145.b0408" "511145.b0414" "511145.b0436" "511145.b0438" "511145.b0471" "511145.b0474"
##  [61] "511145.b0523" "511145.b0529" "511145.b0600" "511145.b0674" "511145.b0720" "511145.b0726"
##  [67] "511145.b0728" "511145.b0774" "511145.b0784" "511145.b0800" "511145.b0813" "511145.b0821"
##  [73] "511145.b0828" "511145.b0870" "511145.b0871" "511145.b0888" "511145.b0907" "511145.b0908"
##  [79] "511145.b0911" "511145.b0924" "511145.b0928" "511145.b0930" "511145.b1033" "511145.b1091"
##  [85] "511145.b1096" "511145.b1131" "511145.b1136" "511145.b1190" "511145.b1211" "511145.b1232"
##  [91] "511145.b1241" "511145.b1260" "511145.b1261" "511145.b1262" "511145.b1263" "511145.b1264"
##  [97] "511145.b1281" "511145.b1302" "511145.b1378" "511145.b1380" "511145.b1387" "511145.b1479"
## [103] "511145.b1622" "511145.b1692" "511145.b1713" "511145.b1719" "511145.b1748" "511145.b1761"
## [109] "511145.b1766" "511145.b1767" "511145.b1778" "511145.b1779" "511145.b1800" "511145.b1812"
## [115] "511145.b1814" "511145.b1831" "511145.b1851" "511145.b1852" "511145.b1864" "511145.b1945"
## [121] "511145.b1983" "511145.b2019" "511145.b2020" "511145.b2021" "511145.b2022" "511145.b2024"
## [127] "511145.b2025" "511145.b2026" "511145.b2041" "511145.b2114" "511145.b2160" "511145.b2171"
## [133] "511145.b2290" "511145.b2295" "511145.b2312" "511145.b2319" "511145.b2320" "511145.b2329"
## [139] "511145.b2341" "511145.b2370" "511145.b2373" "511145.b2379" "511145.b2414" "511145.b2421"
## [145] "511145.b2463" "511145.b2472" "511145.b2474" "511145.b2478" "511145.b2496" "511145.b2499"
## [151] "511145.b2507" "511145.b2508" "511145.b2514" "511145.b2530" "511145.b2551" "511145.b2552"
## [157] "511145.b2557" "511145.b2560" "511145.b2567" "511145.b2568" "511145.b2574" "511145.b2599"
## [163] "511145.b2600" "511145.b2606" "511145.b2609" "511145.b2662" "511145.b2687" "511145.b2696"
## [169] "511145.b2699" "511145.b2708" "511145.b2741" "511145.b2763" "511145.b2784" "511145.b2786"
## [175] "511145.b2797" "511145.b2818" "511145.b2838" "511145.b2871" "511145.b2913" "511145.b2926"
## [181] "511145.b2927" "511145.b2935" "511145.b2937" "511145.b2942" "511145.b2957" "511145.b3008"
## [187] "511145.b3018" "511145.b3030" "511145.b3067" "511145.b3073" "511145.b3117" "511145.b3164"
## [193] "511145.b3168" "511145.b3170" "511145.b3172" "511145.b3177" "511145.b3189" "511145.b3199"
## [199] "511145.b3210" "511145.b3212" "511145.b3225" "511145.b3236" "511145.b3237" "511145.b3258"
## [205] "511145.b3281" "511145.b3287" "511145.b3290" "511145.b3340" "511145.b3350" "511145.b3359"
## [211] "511145.b3368" "511145.b3389" "511145.b3390" "511145.b3414" "511145.b3433" "511145.b3461"
## [217] "511145.b3549" "511145.b3607" "511145.b3616" "511145.b3650" "511145.b3670" "511145.b3671"
## [223] "511145.b3699" "511145.b3701" "511145.b3702" "511145.b3729" "511145.b3738" "511145.b3744"
## [229] "511145.b3769" "511145.b3770" "511145.b3771" "511145.b3772" "511145.b3774" "511145.b3783"
## [235] "511145.b3804" "511145.b3809" "511145.b3823" "511145.b3829" "511145.b3846" "511145.b3847"
## [241] "511145.b3859" "511145.b3866" "511145.b3870" "511145.b3888" "511145.b3916" "511145.b3919"
## [247] "511145.b3927" "511145.b3938" "511145.b3939" "511145.b3940" "511145.b3941" "511145.b3956"
## [253] "511145.b3958" "511145.b3959" "511145.b3960" "511145.b3962" "511145.b3972" "511145.b3973"
## [259] "511145.b3987" "511145.b3988" "511145.b4013" "511145.b4018" "511145.b4019" "511145.b4024"
## [265] "511145.b4025" "511145.b4053" "511145.b4054" "511145.b4139" "511145.b4141" "511145.b4146"
## [271] "511145.b4147" "511145.b4177" "511145.b4219" "511145.b4221" "511145.b4244" "511145.b4245"
## [277] "511145.b4255" "511145.b4258" "511145.b4297" "511145.b4298" "511145.b4388" "511145.b4471"
## [283] "511145.b0001" "511145.b0002" "511145.b0004" "511145.b0005" "511145.b0026" "511145.b0031"
## [289] "511145.b0033" "511145.b0038" "511145.b0060" "511145.b0067" "511145.b0071" "511145.b0072"
## [295] "511145.b0073" "511145.b0074" "511145.b0077" "511145.b0078" "511145.b0088" "511145.b0134"
## [301] "511145.b0149" "511145.b0166" "511145.b0242" "511145.b0243" "511145.b0261" "511145.b0386"
## [307] "511145.b0452" "511145.b0641" "511145.b0674" "511145.b0757" "511145.b0813" "511145.b0870"
## [313] "511145.b0907" "511145.b0908" "511145.b1131" "511145.b1261" "511145.b1262" "511145.b1263"
## [319] "511145.b1333" "511145.b1602" "511145.b1622" "511145.b1761" "511145.b1813" "511145.b1814"
## [325] "511145.b1831" "511145.b1856" "511145.b2019" "511145.b2020" "511145.b2021" "511145.b2022"
## [331] "511145.b2023" "511145.b2024" "511145.b2025" "511145.b2026" "511145.b2295" "511145.b2319"
## [337] "511145.b2320" "511145.b2329" "511145.b2472" "511145.b2478" "511145.b2507" "511145.b2551"
## [343] "511145.b2557" "511145.b2599" "511145.b2600" "511145.b2682" "511145.b2687" "511145.b2797"
## [349] "511145.b2818" "511145.b2838" "511145.b2913" "511145.b2942" "511145.b3008" "511145.b3117"
## [355] "511145.b3199" "511145.b3212" "511145.b3237" "511145.b3389" "511145.b3390" "511145.b3433"
## [361] "511145.b3572" "511145.b3607" "511145.b3670" "511145.b3671" "511145.b3769" "511145.b3770"
## [367] "511145.b3771" "511145.b3772" "511145.b3774" "511145.b3823" "511145.b3829" "511145.b3938"
## [373] "511145.b3939" "511145.b3940" "511145.b3941" "511145.b3956" "511145.b3958" "511145.b3959"
## [379] "511145.b4013" "511145.b4019" "511145.b4024" "511145.b4255" "511145.b4388" "511145.b4471"

Interactions entre un ensemble de sommets donnés

stringdb$get_interactions(s1.mapped)

2.3 Données de coexpression à partir des données complètes

Nous allons utiliser les données d’expression déjà intégrées dans STRINGdb. Il est bien sûr possible de calculer un score de co-expression entre les paires de gènes à partir de données de microarray et/ou de RNA-Seq.

2.4 Autres scores StringDB à intégrer

Faire de même avec les liens correspondant aux

  • interactions protéines–protéines → colonne experimental,
  • les associations basées sur la conservation des paires de gènes dans le voisinage sur le chromosome dans d’autres génomes → neighborhood,
  • les associations basées sur la littérature bio-médicale → textmining, (que l’on pourra comparer aux identifiants PubMed→Gene )
  • les associations basées sur les annotations → database, (que l’on pourra comparer aux liens Keyword→Gene, Pathway→Gene et GOTerm→Gene)
  • le score combiné → combined_score

clean up STRINGdb

"MATCH ()-[r:STRINGdb]-() DELETE(r)" %>% cypher
"MATCH ()-[r:STRINGdb]-() RETURN count(r)" %>% cypher

as function

ds='combined_score'
string.to.neo = function(ds) {
  # links import file
  links.detailed %>%
  filter(.[, ds]>0 & protein1 < protein2) %>%
  dplyr::select(protein1, protein2, all_of(ds)) %>%
  mutate(organism=str_extract(protein1, '^\\d+'), id1=str_extract(protein1, 'b\\d+'), id2=str_extract(protein2, 'b\\d+')) %>%
  dplyr::select(organism:id2, all_of(ds)) %>%
  write_csv(paste0("neo4j.import/string.",ds,".csv"))
  # clean up ?
  
  # import
  paste0('
  LOAD CSV WITH HEADERS FROM "file:///string.',ds,'.csv" AS line
  MATCH (g1:Gene),(g2:Gene) WHERE g1.id=line.id1 AND g2.id=line.id2
  WITH g1, g2, toInteger(line.',ds,') AS value 
  MERGE (g1)-[r:STRINGdb]-(g2) SET r.',ds,'=value
  ')  %>% cypher
  # check
  paste0('MATCH ()-[r:STRINGdb]-() WHERE r.',ds,'>0 RETURN count(r) as ',ds) %>% cypher %>% unlist
}

experimental

string.to.neo('experimental')

neighborhood

string.to.neo('neighborhood')

textmining

string.to.neo('textmining')

database

string.to.neo('database')

combined_score

string.to.neo('combined_score')

2.5 Combinaison de matrices de similarité

Calcul du score combiné par STRINGdb (source: http://string-db.org/help/faq/#how-are-the-scores-computed).

En utilisant la librairie igraph et après avoir chargé le graphe \(g\) avec les liens de coexpression, neighborhood et experiment à partir de neo4j :

prior=0.041
no_prior = function(x, prior = 0.041) (ifelse(is.na(x), 0, x) / 1000 - prior) / (1-prior)
s_coexp_nop = no_prior(E(g)$coexpression)
s_ppi_nop = no_prior(E(g)$ppi)
s_neighborhood_nop = no_prior(E(g)$neighborhood)
s_tot_nop = 1 - (1 - s_coexp_nop) * (1 - s_ppi_nop) * (1 - s_neighborhood_nop)
E(g)$combined_score = round(1000 * (s_tot_nop + prior *(1 - s_tot_nop)))

2.6 Prioritisation de gènes

Gènes “training”. Choix du cycle des citrates (pathway TCA). Récupération des gènes impliqués dans le pathway

training.genes = "MATCH (Pathway {id: 'TCA'})-[:requires]->(g: Gene {organism: 511145}) return g.bnumber " %>% cypher  %>% .$g %>% unlist
training.genes
##  value1  value2  value3  value4  value5  value6  value7  value8  value9 value10 value11 
## "b0116" "b0118" "b0720" "b0721" "b0722" "b0723" "b0724" "b0726" "b0727" "b0728" "b0729" 
## value12 value13 value14 value15 value16 value17 value18 value19 value20 value21 
## "b1136" "b1276" "b1479" "b1611" "b1612" "b1675" "b2210" "b2929" "b3236" "b4122"

Gènes “candidats” (ensemble du génome) :

candidates = "MATCH (g:Gene {organism: 511145}) RETURN g.bnumber" %>% cypher %>% .$g 
candidates

Distance d’un gène à l’ensemble de référence

score=function(gene.ids, ref.genes, datasource) {
  ref_str = paste("'", ref.genes,"'", sep = '', collapse = ',')
  sapply(gene.ids, function(gene.id) {
    query = paste0("MATCH (g1:Gene {id: '", gene.id, "', organism: 511145})-[r:STRINGdb]-(g2:Gene {organism: 511145}) WHERE g2.id IN [",ref_str,"] RETURN r.",datasource)
    res = query %>% cypher %>% .$r
    ifelse(is.null(res), 0, replace_na(res, list(value=0)) %>% unlist %>% mean)
  })
}
# test
gene.ids='b0001'
gene.ids='b0116'
ref.genes = training.genes
datasource='coexpression'
score(gene.ids, training.genes, datasource)
##    b0116 
## 233.3333

Application à l’ensemble des gènes

scores = candidates %>%
  mutate(
    coexpression = score(gene.id=value, ref.genes=training.genes, datasource='coexpression'),
    experimental  = score(gene.id=value, ref.genes=training.genes, datasource='experimental'),
    neighborhood  = score(gene.id=value, ref.genes=training.genes, datasource='neighborhood'),
    textmining  = score(gene.id=value, ref.genes=training.genes, datasource='textmining'),
    database  = score(gene.id=value, ref.genes=training.genes, datasource='database')
         )
scores = scores  %>%  replace(is.na(.), 0)
scores

Plots

scores %>%
  ggplot(aes(x=coexpression, y=neighborhood)) +
  geom_point(alpha=0.2, aes(color=value %in% training.genes, shape=value %in% training.genes))

ACP pour l’affichage

pca = scores %>%
  dplyr::select(-value) %>%
  prcomp
pca %>% summary
## Importance of components:
##                             PC1     PC2     PC3      PC4      PC5
## Standard deviation     111.4271 61.9468 56.1815 47.05207 30.54795
## Proportion of Variance   0.5504  0.1701  0.1399  0.09815  0.04137
## Cumulative Proportion    0.5504  0.7206  0.8605  0.95863  1.00000

Visualisation

library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
pca %>% fviz_pca_ind(col.ind = scores$value %in% training.genes, addEllipses = T, alpha.ind = .2, geom='point')

On peut déjà observer une petite distinction entre les gènes de référence et les autres.

Contribution des sources de données

pca %>% fviz_pca_biplot(col.ind = scores$value %in% training.genes)

Analyse discriminante linéaire

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
mat = scores %>% dplyr::select(-value) %>% as.matrix
model = lda(x=mat, grouping = scores$value %in% training.genes)
model
## Call:
## lda(mat, grouping = scores$value %in% training.genes)
## 
## Prior probabilities of groups:
##       FALSE        TRUE 
## 0.995133256 0.004866744 
## 
## Group means:
##       coexpression experimental neighborhood textmining   database
## FALSE     28.57813     8.591391     16.71926   51.94859   7.220236
## TRUE     389.95476    72.936772    129.07804  591.82725 392.427249
## 
## Coefficients of linear discriminants:
##                      LD1
## coexpression 0.006007943
## experimental 0.000204320
## neighborhood 0.002933981
## textmining   0.003821577
## database     0.012523662

Remarque : attention lors de l’utilisation des librairies MASS et tidyverse, notamment à l’ordre de leur chargement car les 2 fournissent la fonction select. Par exemple ici, nous avons d’abord chargé tidyverse puis MASS. Les bouts de code précédents qui utilise select ne vont plus fonctionner sans préciser qu’il s’agit de dplyr::select que l’on souhaite appeler et non celle de MASS.

Distribution des scores sous forme de boîtes à moustache pour les 2 classes :

mat %>% 
  scale( center=T, scale=F ) %*% model$scaling %>% 
  as_tibble %>%
  mutate(gene.id=scores$value) %>%
  ggplot(aes(x=LD1, y=gene.id %in% training.genes, color=gene.id %in% training.genes, shape=gene.id %in% training.genes)) +
  geom_violin() +
  geom_boxplot(varwidth = T) +
  geom_jitter(height = 0.2, alpha=0.1, color='grey') +
  theme_light()

Density

mat %>% 
  scale( center=T, scale=F ) %*% model$scaling %>% 
  as_tibble %>%
  mutate(gene.id=scores$value) %>%
  ggplot(aes(LD1)) + #, color=gene.id %in% training.genes, shape=gene.id %in% training.genes)) +
  geom_density()

On remarquera que faire >4k requête (1 par gène candidat) prend du temps afin d’obtenir la moyenne (ou autre) entre le candidat et les gènes de référence.

Tentative d’optimisation.

Récupération des “scores” entre les candidats et les gènes de référence en une requête

train.str = paste("'", training.genes,"'", sep = '', collapse = ',')
datasources = c('coexpression', 'experimental', 'neighborhood', 'textmining', 'database' )
ds.str=paste(" r.", datasources, sep='', collapse=',')
pri = paste0("MATCH (candidate:Gene {organism: 511145})-[r:STRINGdb]-(training:Gene {organism: 511145}) WHERE training.id IN [", train.str,"]  RETURN DISTINCT candidate.id, training.id, ", ds.str) %>% cypher
tb = tibble(candidate=pri$candidate.id$value, ref.gene=pri$training.id$value)
scores = sapply(datasources, function(ds) {
  dsr=paste0('r.', ds)
  tb[ds] =  pri[[dsr]]$value
})
tb = tb %>% bind_cols(scores)
tb

Replace NAs with 0s and compute mean (or other) by candidate

scores = tb %>% 
  dplyr::select(-ref.gene) %>%
  replace(is.na(.), 0) %>%
  pivot_longer(cols = datasources, names_to='datasource') %>%
  group_by(candidate, datasource) %>%
  summarise(score=mean(value)) %>%
  pivot_wider(names_from = datasource, values_from = score) %>% 
  ungroup
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(datasources)
## 
##   # Now:
##   data %>% select(all_of(datasources))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
## `summarise()` has grouped output by 'candidate'. You can override using the `.groups`
## argument.
scores

Analyse discriminante linéaire

mat = scores %>% dplyr::select(-candidate) %>% as.matrix
model = lda(x=mat, grouping = scores$candidate %in% training.genes)
model
## Call:
## lda(mat, grouping = scores$candidate %in% training.genes)
## 
## Prior probabilities of groups:
##      FALSE       TRUE 
## 0.98952096 0.01047904 
## 
## Group means:
##       coexpression  database experimental neighborhood textmining
## FALSE     61.78146  15.62323     18.49974      35.9322   111.8704
## TRUE     389.95476 392.42725     72.93677     129.0780   591.8272
## 
## Coefficients of linear discriminants:
##                      LD1
## coexpression 0.005540553
## database     0.007714013
## experimental 0.002397328
## neighborhood 0.004865056
## textmining   0.004661039

Distribution des scores sous forme de boîtes à moustache pour les 2 classes :

mat %>% 
  scale( center=T, scale=F ) %*% model$scaling %>% 
  as_tibble %>%
  mutate(gene.id=scores$candidate, training=gene.id %in% training.genes) %>%
  ggplot(aes(x=LD1, y=training,  color=training, shape=training)) +
  geom_violin() +
  geom_jitter(alpha=.5, height = .2)

3 Pour aller plus loin … types d’analyse possibles qui exploitent Neo4J

3.1 Comparaison de la co-expression entre les gènes de mêmes opérons et d’opérons différents

Ici, on peut considérer les TUs comme un partitionnement des gènes. On pourra donc considérer une mesure de qualité de clustering telle que le coefficient de silhouette.

3.2 Quelles sont les unités de transcription co-exprimées ?

A compléter par soi-même.

3.3 Correspondent-elles à certains pathways ? partagent-elles des annotations GO ?

A compléter.