---
title: "*E. coli* knowledge base"
author: "RB"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    code_folding: show
    toc: TRUE
    toc_float: TRUE
    toc_depth: 4
    theme: paper
    highlight: tango
editor_options: 
  chunk_output_type: console
---


```{r setup, include=TRUE}
#knitr::opts_chunk$set(echo = TRUE)
options(width = 120)
# 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
``` 

# Prérequis

Installation éventuelle des librairies R dans un environnement conda.

shell:
```{bash eval=F}
conda create -n idh
conda deactivate
conda activate idh
conda install r-base=4.0.2 r-reticulate r-tidyverse r-kableextra r-dt r-rmysql bioconductor-stringdb openjdk=11.0.8 r-knitr r-igraph scipy numpy openjdk
```

# Librairies et fonctions utilitaires

```{r chargement librairie}
library(reticulate) # to execute bash or python chunks (only in notebooks)
library(tidyverse)
library(magrittr)

library(kableExtra) # nicer tables
kabex = function(o) o %>% kable(format="html", escape=F) %>% kable_styling(bootstrap_options = c("striped"))
library(DT) # other nice tables
dted = function(o) o %>% datatable(rownames=F, filter="top", options=list(pageLength = 10), escape=F)

library(igraph)
igraph.options(vertex.color=NA)
igraph.options(vertex.label.cex=.6) # font size
igraph.options(vertex.label.family='sans')
igraph.options(vertex.size=2)
igraph.options(edge.label.cex=.6)
igraph.options(edge.label.family='sans')

library(STRINGdb) # client STRING
```


# Génome/Protéome

La première étape consiste à récupérer l'ensemble du génome/protéome qui nous intéresse. 

## STRINGdb Proteins

Création du répertoire pour les données de STRING (shell) :
```{bash eval=F}
mkdir string.data
```

Connexion au serveur de STRING (NCBI taxon id de Ecoli est 511145: https://string-db.org/cgi/input.pl?input_page_active_form=organism )
```{r echo=TRUE}
confidence_threshold = 700
stringdb = STRINGdb$new(version='11', species=511145, score_threshold=confidence_threshold, input_directory='string.data')
```

Téléchargement du génome/protéome :
```{r}
proteins = stringdb$get_proteins() %>% as_tibble
proteins %>% dted
```

## Neo4j installation

Sites et documentations :

  * https://neo4j.com/
  * https://neo4j.com/developer/get-started/ 

Récupérer la dernière version : Téléchargement linux https://neo4j.com/download-center/#releases (onglet *Community server*)

Création des répertoires et désarchivage (shell):
```{bash eval=F}
mkdir -p neo4j/download
cd neo4j/download
version=4.1.3
tar xf neo4j-community-$version-unix.tar.gz
ln -s neo4j-community-$version neo4j-community
cd neo4j-community/
```

Démarrage et arrêt du serveur (shell):
```{bash eval=F}
./bin/neo4j console
```

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

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

A la première connexion, le mot de passe est `neo4j`, le système demande ensuite de changer le mot de passe. Explorez l'interface Web de Neo4j browser, notamment le côté gauche avec les paramètres, et les informations sur la base de données. 

Utilisation depuis le shell :
```{bash eval=F}
./bin/cypher-shell
```

Nous allons maintenant créer un fichier CSV pour le protéome d'Ecoli pour l'importer dans Neo4j (créer les sommets correspondant aux protéines) :
```{r}
proteins %>%
  mutate(taxon_id = str_extract(protein_external_id, "([0-9]+)"), protein_id= str_match(protein_external_id, "\\.(.+)$")[,2]) %>%
  select(protein_external_id, taxon_id, protein_id, name=preferred_name, length=protein_size, annotation) %>%
  write_csv("neo4j/download/neo4j-community/import/ecoli.proteins.csv")
```

<!-- plus nécessaire depuis la 4.1.3
Pour pouvoir importer les données depuis un fichier local, il faut modifier la configuration de neo4j : dans le fichier `conf/neo4j.conf`, ajouter vers la ligne 227 :
```
dbms.security.allow_csv_import_from_file_urls=true
```
-->

Librairie neo4r
```{r}
if (!require('neo4r')) { # client neo4j (pas dispo avec conda)
  install.packages('neo4r')
}
library(neo4r)
neodb = neo4j_api$new(
  url = "http://localhost:7474", 
  user = "neo4j", 
  password = "bioinfo"
)
cypher = function(query, neo4j=neodb, ...) call_neo4j(query=query, con=neo4j, ...)
# 'MATCH (p:Protein {name:"dnaJ"}) RETURN p.protein_id, p.annotation' %>% cypher
```


Import dans neo4j
```{r eval=F}
# DELETE PREVIOUS ATTEMPTS ON Protein
'MATCH ()-[r:STRINGdb]-() DELETE r' %>% cypher 
'MATCH ()-[r:GO_annotation]-() DELETE r' %>% cyper
'MATCH (n:Protein) DELETE n' %>% cypher
# IMPORT
'LOAD CSV WITH HEADERS FROM "file:///ecoli.proteins.csv" AS row 
CREATE (n:Protein)
SET n = row,
 n.protein_external_id = row.protein_external_id,
 n.taxon_id = toInteger(row.taxon_id),
 n.protein_id = row.protein_id,
 n.name = row.name,
 n.length = toInteger(row.length),
 n.annotation = row.annotation
'  %>% cypher
```

Vérification *via* des requêtess cyper:
```{r}
'MATCH (p:Protein {name:"dnaJ"}) RETURN p.protein_id, p.annotation'  %>% cypher(type="row") 
'MATCH (p:Protein) WHERE p.length <= 40 RETURN p'  %>% cypher %>% extract2('p') %>% dted
'MATCH (p:Protein) RETURN count(p)'  %>% cypher 
```

Ajout d'un index (pour accélérer la création des liens par la suite):
```{r eval=F}
'CREATE INDEX ON :Protein(protein_external_id)' %>% cypher
```

# Données d'expression

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 RNAseq.

## STRINGdb detailed links → neo4j

Il faut passer par la page de téléchargement sur https://string-db.org/cgi/download.pl en restreignant l'organisme d'intérêt. 

Le fichier complet faisant 67.6Go (au 15/09/2020), il est conseillé de d'abord sélectionner un organisme avant le téléchargement. Une copie est disponible sur silico à récupérer sur http://silico.biotoul.fr/enseignement/m2bbs/idh/

Qu'y a-t-il dedans ? shell:
```{bash}
zcat string.data/511145.protein.links.detailed.v11.0.txt.gz | head
```

Chargement sous forme de tibble :
```{r}
links.detailed = read_delim("string.data/511145.protein.links.detailed.v11.0.txt.gz", delim=" ", col_types = "ccnnnnnnnn")
links.detailed
```

C'est la **colonne coexpression** qui nous intéresse.

Génération des fichiers CSV pour l'import
```{r}
links.detailed %>%
  filter(coexpression>0 & protein1 < protein2) %>%
  select(protein1, protein2, coexpression) %>%
  write_csv("neo4j/download/neo4j-community/import/ecoli.string.coexpression.csv")
```

Import dans neo4j
```{r eval=F}
'
LOAD CSV WITH HEADERS FROM "file:///ecoli.string.coexpression.csv" AS line
MATCH (p1:Protein),(p2:Protein)
WHERE p1.protein_external_id=line.protein1 AND p2.protein_external_id=line.protein2
WITH p1,p2, toInteger(line.coexpression) AS value 
MERGE (p1)-[:STRINGdb {coexpression: value}]-(p2)
'  %>% cypher
```


Vérification
```{r}
'MATCH ()-[r:STRINGdb]-() RETURN count(r)' %>% cypher
```

Extraction des liens dont la probabilité est supérieure à 0.7 sous forme de tibble
```{r}
foo = 'MATCH (p1)-[r:STRINGdb]-(p2) WHERE r.coexpression>=700 RETURN p1.name, p2.name, r.coexpression' %>% cypher 
tibble(protein1=foo$p1.name$value, protein2=foo$p2.name$value, coexpression=foo$r.coexpression$value) %>% dted
```

Extraction d'un graphe
```{r}
g.coexpr = 'MATCH p = ()-[r:STRINGdb]-() WHERE r.coexpression>=700 RETURN p'  %>% cypher(type="graph")
```

Manipulation pour `igraph` pour récupérer les propriétés des sommets (https://github.com/neo4j-rstats/neo4r)
```{r}
g.coexpr$nodes = g.coexpr$nodes %>% unnest_nodes(what="properties")
g.coexpr$nodes %>% dted
```

Besoin de réordonner les liens
```{r}
g.coexpr$relationships = g.coexpr$relationships %>% 
  unnest_relationships %>% 
  select(startNode, endNode, type, everything()) %>% 
  mutate(coexpression = unlist(coexpression))
g.coexpr$relationships %>% dted
```

Création dans igraph
```{r}
g.coexpr = graph_from_data_frame(d=g.coexpr$relationships, directed=F, vertices = g.coexpr$nodes)
g.coexpr
```

Plot
```{r}
plot(g.coexpr, vertex.label=NA, main='coexpression')
```

# Données d'interaction protéine-protéine
## STRINGdb protein-protein interactions → neo4j

Pour le réseau d'interaction ppi, la colonne `experimental` contient des interactions physiques et fonctionnelles. Un autre fichier est donc nécessaire : protein.actions (binding)

```{r}
actions = read_delim("string.data/511145.protein.actions.v11.0.txt.gz", delim="\t")
actions
```

Seulement les *binding*
```{r}
links.ppi = actions %>% 
  filter(mode=='binding' & item_id_a<item_id_b) %>% 
  dplyr::select(protein1 = item_id_a, protein2 = item_id_b, ppi = score) %>%
  filter(ppi >= confidence_threshold)
```

`experimetal:` BIND, DIP, GRID, HPRD, IntAct, MINT, and PID. (src: https://string-db.org/help/faq/#which-databases-does-string-extract-experimental-data-from)

Comme précédemment, création du fichier CSV puis import dans neo4j
```{r eval=F}
links.ppi %>% write_csv("neo4j/download/neo4j-community/import/ecoli.string.ppi.csv")
'LOAD CSV WITH HEADERS FROM "file:///ecoli.string.ppi.csv" AS line
MATCH (p1:Protein)-[r:STRINGdb]-(p2:Protein)
WHERE p1.protein_external_id=line.protein1 AND p2.protein_external_id=line.protein2
WITH p1,p2, toInteger(line.ppi) AS value
MERGE (p1)-[r:STRINGdb]-(p2)
SET r.ppi=value
'  %>% cypher
```


```{r}
'MATCH p = (p1:Protein {name:"dnaJ"})-[r]-(p2:Protein {name:"dnaK"}) return p' %>% cypher
```

# Graphe de coexpression & ppi

Extraction du graphe et visualisation
```{r}
g.ppi.coexpr = 'MATCH p = ()-[r:STRINGdb]-() WHERE r.coexpression>=700 AND r.ppi>=700 RETURN p'  %>% cypher(type="graph")
g.ppi.coexpr$nodes = g.ppi.coexpr$nodes %>% unnest_nodes(what="properties")
g.ppi.coexpr$relationships = g.ppi.coexpr$relationships %>% 
  unnest_relationships %>% 
  select(startNode, endNode, type, everything()) %>% 
  mutate(coexpression = unlist(coexpression), ppi = unlist(ppi))
g.ppi.coexpr = graph_from_data_frame(d=g.ppi.coexpr$relationships, directed=F, vertices = g.ppi.coexpr$nodes)
g.ppi.coexpr
plot(g.ppi.coexpr, vertex.label=NA, main='coexpr & ppi')
```

Les contraintes de coexpression > 0.7 et de ppi > 0.7 semblent un peu trop restrictives. Nous allons prendre un OR à la place pour avoir un graphe plus dense :
```{r}
g.ppi.coexpr = 'MATCH p = ()-[r:STRINGdb]-() WHERE r.coexpression>=700 OR r.ppi>=700 RETURN p'  %>% cypher(type="graph")
g.ppi.coexpr$nodes = g.ppi.coexpr$nodes %>% unnest_nodes(what="properties")
g.ppi.coexpr$relationships = g.ppi.coexpr$relationships %>% 
  unnest_relationships %>% 
  select(startNode, endNode, type, everything()) %>% 
  mutate(coexpression = unlist(coexpression), ppi = unlist(ppi))
g.ppi.coexpr = graph_from_data_frame(d=g.ppi.coexpr$relationships, directed=F, vertices = g.ppi.coexpr$nodes)
g.ppi.coexpr
plot(g.ppi.coexpr, vertex.label=NA, vertex.size=1, main="coexpr | ppi")
```



# Données phylogénomiques

Dans STRINGdb, ces données sont disponibles pour les bactéries. Elles résultent de l'analyse de la conservation des gènes otrhologues dans le voisinage sur les chromosomes d'autres organismes.

## Ajout des liens *conserved genomic neighborhood*

Génération des fichiers CSV pour l'import comme pour *coexpression*
```{r echo=F}
links.detailed %>%
  filter(neighborhood>0 & protein1 < protein2) %>%
  select(protein1, protein2, neighborhood) %>%
  write_csv("neo4j/download/neo4j-community/import/ecoli.string.neighborhood.csv")
```

Ajout des lien dans neo4j
```{r eval=F}
'LOAD CSV WITH HEADERS FROM "file:///ecoli.string.neighborhood.csv" AS line
MATCH (p1:Protein)-[r:STRINGdb]-(p2:Protein)
WHERE p1.protein_external_id=line.protein1 AND p2.protein_external_id=line.protein2
WITH p1,p2, toInteger(line.neighborhood) AS value
MERGE (p1)-[r:STRINGdb]-(p2)
SET r.neighborhood=value
'  %>% cypher
```

Extraction du graphe et visualisation
```{r}
g.string = 'MATCH p = ()-[r:STRINGdb]-() WHERE r.coexpression>=700 OR r.ppi>=700 OR r.neighborhood>=700 RETURN p'  %>% cypher(type="graph")
g.string$nodes = g.string$nodes %>% unnest_nodes(what="properties")
g.string$relationships = g.string$relationships %>% 
  unnest_relationships %>% 
  select(startNode, endNode, type, everything()) %>% 
  mutate(coexpression = unlist(coexpression), ppi = unlist(ppi), neighborhood = unlist(neighborhood))
g.string = graph_from_data_frame(d=g.string$relationships, directed=F, vertices = g.string$nodes)
g.string
plot(g.string, vertex.label=NA, vertex.size=1, main='coexpr | pppi | neighborhood')
```

Tailles des composantes connexes
```{r}
clusters(g.string)$csize
```

Extraction des sous-graphes correspondant aux CC de taille minimum 10
```{r}
CCs.num = which(clusters(g.string)$csize >=10)
CCs.size = clusters(g.string)$csize[CCs.num]
V.ids = V(g.string)[ clusters(g.string)$membership %in% CCs.num ]
g = induced_subgraph(g.string, V.ids)
plot(g, vertex.label=NA)
```


Calcul du score combiné avant détection de communautés (source: http://version11.string-db.org/help/faq/#how-are-the-scores-computed)
```{r}
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)))
```

Détection de communautés prenant en compte le score combiné
```{r}
g.communities = cluster_louvain(g, weights=E(g)$combined_score)
plot(g.communities, g, vertex.label=NA)
```

Nombre de communautés
```{r}
length(g.communities)
```

Tailles des communautés
```{r}
sizes(g.communities)
```

Modularité
```{r}
modularity(g, g.communities$membership, weights=E(g)$combined)
```

# Génération des fichiers par clusters/communautés pour la recherche d'enrichissement en termes GO
```{r}
tibble(protein_id = V(g)$protein_id, name=V(g)$name, cluster = g.communities$membership) %>% dted
```

```{r}
groups(g.communities) %>% head(2)
```

Remplacement des noms par les identifiants "stringdb"
```{r}
g.clusters = lapply(groups(g.communities), function(x) V(g)[ name %in% x]$protein_id)
head(g.clusters, 2)
```

Génération du fichier
```{r}
lines = sapply(g.clusters, function(x) paste(x, collapse = '\t')) %>% simplify2array 
write_lines(lines, "enrichment.coli/g.communities.txt")
```

Split by lines (shell)
```{bash eval=F}
split -l 1 --numeric-suffix=1 --additional-suffix=.txt enrichment.coli/g.communities.txt g.com.
mv g.com.*.txt enrichment.coli/
```


# Données d'annotation
## Gene Ontology → neo4j


## Restauration de la base fournie par GO

* télécharger le dump sur go.org http://archive.geneontology.org/latest-termdb/go_daily-termdb-data.gz
* restaurer la base
* extraire le contenu -> go_daily-termdb-tables

Librairie MariaDb/MySQL
```{r}
library(RMySQL) # client MariaDb
dbh=dbConnect(MySQL(), user="guest", password="bioinfo", dbname="go_idh_2020", host="127.0.0.1")
dbget = function(query) as_tibble(dbGetQuery(dbh, query))
```


Variables d'environnement (shell) qui facilitent les changements de nom de DB *etc.* (shell)
```{bash eval=F}
dbname=go_idh_2020
mysql="mysql -uguest -pbioinfo -h 127.0.0.1 $dbname"
mysql -uguest -pbioinfo -e "CREATE DATABASE $dbname"
```

Restauration de GO db (shell)
```{bash eval=F}
mkdir go.data
cd go.data
wget http://archive.geneontology.org/latest-termdb/go_daily-termdb-data.gz
gunzip go_daily-termdb-data.gz
$mysql < go_daily-termdb-data
```

## Extraction de `go.terms.csv` pour import dans neo4j
```{r}
"SELECT id, acc, term_type, name FROM term 
WHERE acc LIKE 'GO:%' AND is_obsolete=0 
LIMIT 10;
" %>% dbget %>% head %>% kabex
```

Export dans un shell
```{bash eval=F}
$mysql -e "SELECT id, acc, term_type, name FROM term WHERE acc LIKE 'GO:%' AND is_obsolete=0;" > neo4j/download/neo4j-community/import/go.terms.tsv
```

## Ajout des GOTerm dans Neo4j

* put data into neo4j subdir "import" qui correspond dans neoj browser à file://
* import des termes GO dans Neo4j 

```{r eval=F}
"MATCH ()-[r:is_a]->() DELETE r" %>% cypher
"MATCH ()-[r:part_of]->() DELETE r" %>% cypher
"MATCH (n:GOTerm) DELETE n" %>% cypher
"
LOAD CSV WITH HEADERS FROM 'file:///go.terms.tsv' AS row FIELDTERMINATOR '\t'
CREATE (n:GOTerm)
SET n = row,
  n.id = toInteger(row.id),
  n.acc = row.acc,
  n.term_type = row.term_type,
  n.name = row.name  
" %>% cypher
```

```{r}
"MATCH (n:GOTerm) RETURN count(n)" %>% cypher
```

Création des index pour accélérer la suite
```{r eval=F}
'CREATE INDEX ON :GOTerm(id)' %>% cypher

'CREATE INDEX ON :GOTerm(acc)' %>% cypher
```

Exploration
```{r}
"MATCH (n:GOTerm { name: 'reproduction'} ) RETURN n" %>% cypher
```

## Ajout des relations entre GOTerms

Il faut identifier les types is_a (1) et part_of (20) dans la table term
```{r}
"SELECT * FROM term WHERE name = 'is_a' OR name = 'part_of';" %>% dbget %>% kabex
```

```{bash eval=F}
is_a=1
part_of=20
$mysql -B -e "SELECT term1_id, term2_id FROM term2term WHERE relationship_type_id = $is_a;" > neo4j/download/neo4j-community/import/go.rel.is_a.dump.tsv
```

```{r eval=F}
"MATCH ()-[r:is_a]->() DELETE r" %>% cypher
"
LOAD CSV WITH HEADERS FROM 'file:///go.rel.is_a.dump.tsv' AS line FIELDTERMINATOR '\t' 
MATCH (t1:GOTerm),(t2:GOTerm) 
WHERE t1.id=toInteger(line.term1_id) AND t2.id=toInteger(line.term2_id) 
WITH t1,t2 
MERGE (t2)-[:is_a]->(t1)
" %>% cypher
```


```{r}
"MATCH ()-[r:is_a]->() RETURN count(r)" %>% cypher
```

```{bash eval=F}
$mysql -B -e "SELECT term1_id, term2_id FROM term2term WHERE relationship_type_id = $part_of;" > neo4j/download/neo4j-community/import/go.rel.part_of.dump.tsv
```

```{r eval=F}
"MATCH ()-[r:part_of]->() DELETE r" %>% cypher
"
LOAD CSV WITH HEADERS FROM 'file:///go.rel.part_of.dump.tsv' AS line FIELDTERMINATOR '\t'
MATCH (t1:GOTerm),(t2:GOTerm)
WHERE t1.id=toInteger(line.term1_id) AND t2.id=toInteger(line.term2_id)
WITH t1,t2
MERGE (t2)-[:part_of]->(t1)
" %>% cypher
```


```{r}
"MATCH ()-[r:part_of]->() RETURN count(r)" %>% cypher
```

## Ajout de l'annotation d'un génome/protéome


### GO Annotations hébergées à l'EBI

Fichier GOA de l'EBI : ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/ 

Celui d'Ecoli 
```{bash eval=F}
mkdir go.data
curl ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/18.E_coli_MG1655.goa -o go.data/18.E_coli_MG1655.goa
```

```{bash}
head go.data/18.E_coli_MG1655.goa
```

Le fichier est au format GAF2.1

Chargement (noms des colonnes de http://geneontology.org/docs/go-annotation-file-gaf-format-2.1/)
```{r}
gaf.cnames = c('DB', 'DB Object ID','DB Object Symbol', 'Qualifier', 'GO ID', 'DB:Reference', 'Evidence Code', 'With (or) From', 'Aspect', 'DB Object Name', 'DB Object Synonym', 'DB Object Type', 'Taxon', 'Date', 'Assigned By', 'Annotation Extension', 'Gene Product Form ID')

go.ecoli = read_tsv('go.data/18.E_coli_MG1655.goa', skip = 2, col_names = gaf.cnames)
go.ecoli %>% head %>% kabex
```

Les identifiants utilisés dans le graphe de STRING sont sur la colonne `DB Object Synonym`. mais nécessitent d'être parsés :
```{r}
go.ecoli %>% select(toparse=`DB Object Synonym`) %>% mutate(bnumber=str_extract(toparse, 'b\\d+'))
```

Ok, on garde seulement les colonnes qui nous intéressent
```{r}
go.ecoli = go.ecoli %>% 
  mutate(bnumber=str_extract(`DB Object Synonym`, 'b\\d+')) %>%
  mutate(protein_external_id=paste0('511145.',bnumber)) %>%
  select(GOTerm=`GO ID`, evidence_code=`Evidence Code`,protein_external_id)
go.ecoli
```

Génération du fichier pour l'import dans neo4j
```{r}
go.ecoli %>% write_csv("neo4j/download/neo4j-community/import/ecoli.goannotation.csv")
```


Ajout des lien dans neo4j
```{r eval=F}
'MATCH ()-[r:GO_annotation]-() DELETE r' %>% cypher
'
LOAD CSV WITH HEADERS FROM "file:///ecoli.goannotation.csv" AS line
MATCH (p:Protein),(t:GOTerm)
WHERE p.protein_external_id=line.protein_external_id AND t.acc=line.GOTerm
WITH p,t,line
MERGE (p)-[r:GO_annotation]-(t)
SET r.evidence_code=line.evidence_code
'  %>% cypher
```


```{r}
'MATCH ()-[r:GO_annotation]-() RETURN count(r)' %>% cypher
```

## Exploration

Annotations directes
```{r}
"
MATCH (:Protein {name:'dnaJ'})-[:GO_annotation]-(t:GOTerm) RETURN t ORDER BY t.name
" %>% cypher
```

Annotations implicites
```{r}
"
MATCH (:Protein {name:'dnaJ'})-[:GO_annotation]-(:GOTerm)-[*]->(t:GOTerm) 
RETURN DISTINCT t 
ORDER BY t.name
" %>% cypher  %>% 
  extract2('t') %>% 
  select(-id) %>%
  dted
```

Annotations directes et implicites
```{r}
terms = "
MATCH (:Protein {name:'dnaJ'})-[:GO_annotation]-(t_direct:GOTerm)-[*]->(t_implicit:GOTerm) 
RETURN t_direct, t_implicit 
" %>% cypher
terms$t_direct %>% unique %>% arrange(name) %>% select(-id) %>% kabex
```

```{r}
terms$t_implicit %>% unique %>% arrange(name) %>% select(-id) %>% dted
```

Protéines annotées directement par un terme
```{r}
"
MATCH (p:Protein)-[:GO_annotation]-(:GOTerm {name: 'chaperone binding'}) 
RETURN DISTINCT p 
" %>% cypher %>% extract2('p') %>% select(-annotation) %>% kabex
```


Terme *chaperone binding* trop spécifique, prenons-en un plus générique
```{r}
"
MATCH (p:Protein)-[:GO_annotation]-(:GOTerm)-[:is_a|part_of*]->(:GOTerm {name: 'protein folding'}) 
RETURN DISTINCT p
ORDER BY p.name
" %>% cypher %>% extract2('p') %>% select(-annotation) %>% dted
```


Protéines annotées directement ou implicitement
```{r}
"
MATCH (p:Protein)-[:GO_annotation|is_a|part_of*]->(:GOTerm {name: 'protein folding'}) 
RETURN DISTINCT p
ORDER BY p.name
" %>% cypher %>% extract2('p') %>% select(-annotation) %>% dted
```

# Génération des fichiers pour la recherche de termes GO sur-représentés

Fonction pour récupérer les termes directs
```{r}
get_GeneProducts = function(goterm, all=F) {
  query = paste0("MATCH (p:Protein)-[:GO_annotation]-(:GOTerm {acc: '", goterm, "'}) RETURN DISTINCT p.protein_id ORDER BY p.protein_id")
  if (all) {
  query = paste0("MATCH (p:Protein)-[:GO_annotation|is_a|part_of*]->(:GOTerm {acc: '", goterm,"'}) RETURN DISTINCT p.protein_id ORDER BY p.protein_id")
  }
  res = query %>% cypher 
  if (length(res) > 0) {
    res = res %>% extract2(1) %>% unlist
  }
  else {
    res = c()
  }
  res
}
```


```{r}
get_GeneProducts("GO:0006457", all=F)
```


```{r}
get_GeneProducts("GO:0006457", all=T)
```

Identification des termes GO concernant *Ecoli*
```{r}
#go.ecoli %>% select(GOTerm) %>% unique %>% arrange(GOTerm)
ecoli.terms = "
MATCH (:Protein {taxon_id:511145})-[:GO_annotation|is_a|part_of*]->(t:GOTerm) 
RETURN DISTINCT t.acc
" %>% cypher %>% unlist
```

Direct GeneProducts
```{r warning=F, message=F, eval=F}
go.sets.direct = lapply(ecoli.terms, get_GeneProducts)
names(go.sets.direct)=ecoli.terms
go.sets.direct[['GO:0006457']]
```

Fichiers *sets*
```{r, eval=F}
header = paste0("# format: sets
# version: 1.3
# strain: Escherichia coli K-12 MG1655
# date: ", format(Sys.time(), '%d %B, %Y'), "
# comment: Gene Ontology terms")
sets = lapply(ecoli.terms, function(x) paste(go.sets.direct[[x]], collapse = '|')) %>% unlist
tb = tibble(GOTerm=ecoli.terms, GeneProducts = sets)
go.terms = read_tsv("neo4j/download/neo4j-community/import/go.terms.tsv")
lines = tb %>% inner_join(go.terms, by=c('GOTerm'='acc')) %>% mutate(desc=paste0(term_type,': ', name)) %>% mutate(txt=paste(GOTerm, desc, GeneProducts, sep='\t')) %>% filter(str_detect(GeneProducts,"\\|"))
write_lines(header,"enrichment.coli/ecoli.go.direct.sets")
write_lines(lines$txt,"enrichment.coli/ecoli.go.direct.sets", append=T)
```



Direct and implicit GeneProducts
```{r warning=F, message=F, eval=F}
go.sets.implicit = lapply(ecoli.terms, function(x) get_GeneProducts(x, all=T))
names(go.sets.implicit)=ecoli.terms
go.sets.implicit[['GO:0006457']]
```

Fichiers *sets*
```{r, eval=F}
header = paste0("# format: sets
# version: 1.3
# strain: Escherichia coli K-12 MG1655
# date: ", format(Sys.time(), '%d %B, %Y'), "
# comment: Gene Ontology terms")
sets = lapply(ecoli.terms, function(x) paste(go.sets.implicit[[x]], collapse = '|')) %>% unlist
tb = tibble(GOTerm=ecoli.terms, GeneProducts = sets)
lines = tb %>% inner_join(go.terms, by=c('GOTerm'='acc')) %>% mutate(desc=paste0(term_type,': ', name)) %>% mutate(txt=paste(GOTerm, desc, GeneProducts, sep='\t')) %>% filter(str_detect(GeneProducts,"\\|"))
write_lines(header,"enrichment.coli/ecoli.go.implicit.sets")
write_lines(lines$txt,"enrichment.coli/ecoli.go.implicit.sets", append=T)
```


# Gene Set Enrichment Analysis

```{bash eval=F}
git clone git@gitlab.com:rbarriot/enrichment.git
```

```{bash}
for i in $(seq 1 9); do 
  ./enrichment/blastset.py -q enrichment.coli/g.com.0$i.txt -t enrichment.coli/ecoli.go.direct.sets -c > enrichment.coli/g.com.0$i.enriched.tsv
done
for i in $(seq 10 42); do 
  ./enrichment/blastset.py -q enrichment.coli/g.com.$i.txt -t enrichment.coli/ecoli.go.direct.sets -c > enrichment.coli/g.com.$i.enriched.tsv
done
```

Résultats
```{r}
read_tsv("enrichment.coli/g.com.01.enriched.tsv") %>% kabex
```


```{r warning=F, message=F,  results='asis'}
cat(sapply(2:9, function(i) kabex(read_tsv(paste0("enrichment.coli/g.com.0",i,".enriched.tsv")))))
cat(sapply(10:42, function(i) kabex(read_tsv(paste0("enrichment.coli/g.com.",i,".enriched.tsv")))))
```


# DELETE ALL nodes & edges
```{r eval=F}
MATCH ()-[r]-() DELETE r; # DELETE ALL RELATIONSHIPS
MATCH (n) DELETE n;       # DELETE ALL NODES
```

