silico.biotoul.fr
 

M1 BBS Graphes TP Librairies R

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (ggtree)
m (Taxonomie du NCBI)
(22 intermediate revisions not shown)
Line 3: Line 3:
Il existe plusieurs librairies disponibles : <tt>ape</tt>, <tt>genoPlotR</tt>, <tt>ggtree</tt>, ... Celle choisie pour ce TP est <tt>ggtree</tt>.
Il existe plusieurs librairies disponibles : <tt>ape</tt>, <tt>genoPlotR</tt>, <tt>ggtree</tt>, ... Celle choisie pour ce TP est <tt>ggtree</tt>.
-
=== ggtree ===
 
-
[https://bioconductor.org/packages/release/bioc/html/ggtree.html ggtree] est une librairie R/Bioconductor pour la '''visualisation''' et l''''annotation''' d'arbres phylogénétiques. Elle est basée sur <tt>ggplot2</tt>, une librairie pour faire des graphiques plus élaborés qu'avec la fonction <tt>plot</tt> de R, qu'il vous est conseillé d'expérimenter par vous-même.
+
=== Taxonomie du NCBI ===
 +
Nous allons utiliser la [https://www.ncbi.nlm.nih.gov/guide/taxonomy/ taxonomy du NCBI]. Pour cela, dans l'onglet <tt>Downloads</tt> on peut aller sur le [ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/ FTP] pour télécharger la banque. Le fichier README (taxdump_readme.txt) contient les information suivantes (avec en gras, celles qui vont nous intéresser tout particulièrement) :
-
Installation de la librairie (nécessaire le 12/12/2017 et ne fonctionne que sous windows cette année)
+
This directory contains the following NCBI Taxonomy database dump files:
-
<source lang='rsplus'>
+
-
source("https://bioconductor.org/biocLite.R")
+
  taxdmp.zip
-
biocLite("ggtree")
+
  taxdump.tar.Z
-
</source>
+
  '''taxdump.tar.gz'''
 +
 +
All these files containes exactly the same information and are arranged so
 +
for the convenience of unpacking them on various operating environments.
 +
In addition there are files:
 +
 +
  taxdmp.zip.md5
 +
  taxdump.tar.Z.md5
 +
  taxdump.tar.gz.md5
 +
 +
which contain MD5 sums for the corresponding archive files. These files
 +
might be used to check correctness of the download of corresponding
 +
archive file.
 +
 +
taxdmp.zip
 +
----------
 +
 +
Is intended for zip-capable utilities such as pkunzip, unzip, and WinZip.
 +
These utilities are widely available in almost all operating environments.
 +
To unpack it command-line pkunzip and unzip:
 +
 +
        pkunzip taxdmp.zip
 +
or
 +
        unzip taxdmp.zip
 +
 +
Note: pkunzip and/or unzip executables must be in the executable search path
 +
and taxdmp.zip must be in the current directory. Files will be unzipped into
 +
current directory. For desired dump files placement and more please refer to
 +
the manual and/or option descriptions of pkunzip and unzip utilities.
 +
 +
taxdump.tar.Z
 +
-------------
 +
 +
This file is to be unpacked by uncompress utility and subsequent tar
 +
archiver. These utilities are usually used in UNIX-like environment.
 +
Unpacking instructions follows:
 +
 +
          uncompress -c taxdump.tar.Z | tar xf -
 +
 +
'''taxdump.tar.gz'''
 +
--------------
 +
 +
This file is to be unpacked by GNU unzip utility and subsequent tar
 +
archiver. These utilities are usually used in UNIX-like environment.
 +
Unpacking instructions follows:
 +
 +
          gunzip -c taxdump.tar.gz | tar xf -
 +
 +
The content of the archive
 +
--------------------------
 +
 +
It may look like this:
 +
 +
citations.dmp
 +
delnodes.dmp
 +
division.dmp
 +
gencode.dmp
 +
merged.dmp
 +
'''names.dmp'''
 +
'''nodes.dmp'''
 +
readme.txt
 +
 +
The readme.txt file gives a brief description of *.dmp files. These files
 +
contain taxonomic information and are briefly described below. Each of the
 +
files store one record in the single line that are delimited by "\t|\n"
 +
(tab, vertical bar, and newline) characters. Each record consists of one
 +
or more fields delimited by "\t|\t" (tab, vertical bar, and tab) characters.
 +
The brief description of field position and meaning for each file follows.
 +
 +
'''nodes.dmp'''
 +
---------
 +
 +
This file represents taxonomy nodes. The description for each node includes
 +
the following fields:
 +
 +
<b>tax_id -- node id in GenBank taxonomy database
 +
parent tax_id -- parent node id in GenBank taxonomy database</b>
 +
rank -- rank of this node (superkingdom, kingdom, ...)
 +
embl code -- locus-name prefix; not unique
 +
division id -- see division.dmp file
 +
inherited div flag  (1 or 0) -- 1 if node inherits division from parent
 +
genetic code id -- see gencode.dmp file
 +
inherited GC  flag  (1 or 0) -- 1 if node inherits genetic code from parent
 +
mitochondrial genetic code id -- see gencode.dmp file
 +
inherited MGC flag  (1 or 0) -- 1 if node inherits mitochondrial gencode from parent
 +
GenBank hidden flag (1 or 0)            -- 1 if name is suppressed in GenBank entry lineage
 +
hidden subtree root flag (1 or 0)      -- 1 if this subtree has no sequence data yet
 +
comments -- free-text comments and citations
 +
 
 +
  <b>names.dmp
 +
  ---------
 +
  Taxonomy names file has these fields:
 +
 
 +
tax_id -- the id of node associated with this name
 +
name_txt -- name itself
 +
unique name -- the unique variant of this name if name not unique
 +
name class -- (synonym, common name, ...)</b> 
 +
 +
division.dmp
 +
------------
 +
Divisions file has these fields:
 +
division id -- taxonomy database division id
 +
division cde -- GenBank division code (three characters)
 +
division name -- e.g. BCT, PLN, VRT, MAM, PRI...
 +
comments
 +
 
 +
gencode.dmp
 +
-----------
 +
Genetic codes file:
 +
 +
genetic code id -- GenBank genetic code id
 +
abbreviation -- genetic code name abbreviation
 +
name -- genetic code name
 +
cde -- translation table for this genetic code
 +
starts -- start codons for this genetic code
 +
 
 +
delnodes.dmp
 +
------------
 +
Deleted nodes (nodes that existed but were deleted) file field:
 +
 +
  tax_id -- deleted node id
 +
 
 +
merged.dmp
 +
----------
 +
Merged nodes file fields:
 +
 +
old_tax_id                              -- id of nodes which has been merged
 +
new_tax_id                              -- id of nodes which is result of merging
 +
 +
citations.dmp
 +
-------------
 +
Citations file fields:
 +
 +
cit_id -- the unique id of citation
 +
cit_key -- citation key
 +
        medline_id                              -- unique id in MedLine database (0 if not in MedLine)
 +
pubmed_id -- unique id in PubMed database (0 if not in PubMed)
 +
url -- URL associated with citation
 +
text -- any text (usually article name and authors)
 +
-- The following characters are escaped in this text by a backslash:
 +
-- newline (appear as "\n"),
 +
-- tab character ("\t"),
 +
-- double quotes ('\"'),
 +
-- backslash character ("\\").
 +
taxid_list -- list of node ids separated by a single space
 +
Extrait de <tt>nodes.dmp</tt> (rappel des colonnes : tax_id, parent tax_id, rank, ...) :
 +
<pre>
 +
1      |      1      |      no rank |              |      8      |      0      |      1      |      0      |      0      |      0      |      0      |      0      |              |
 +
2      |      131567  |      superkingdom    |              |      0      |      0      |      11      |      0      |      0      |      0      |      0      |      0      |              |
 +
6      |      335928  |      genus  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
7      |      6      |      species |      AC      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
9      |      32199  |      species |      BA      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
10      |      1706371 |      genus  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
11      |      1707    |      species |      CG      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
13      |      203488  |      genus  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
14      |      13      |      species |      DT      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
16      |      32011  |      genus  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
17      |      16      |      species |      MM      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
18      |      213421  |      genus  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
19      |      18      |      species |      PC      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
20      |      76892  |      genus  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
21      |      20      |      species |      PI      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
22      |      267890  |      genus  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
23      |      22      |      species |      SC      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
24      |      22      |      species |      SP      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
25      |      22      |      species |      SH      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
27      |      49928  |      species |      HE      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
28      |      49928  |      species |      HE      |      0      |      1      |      11      |      1      |      0      |      1      |      1      |      0      |              |
 +
29      |      28221  |      order  |              |      0      |      1      |      11      |      1      |      0      |      1      |      0      |      0      |              |
 +
</pre>
 +
 +
Ce fichier nous donne donc la structure de l'arbre taxonomique avec les couples <tt>''tax_id''</tt> et <tt>''parent tax_id''</tt>.
 +
 +
Extrait de <tt>names.dmp</tt> (rappel des colonnes : tax_id, name_txt, unique name, name class) :
 +
<pre>
 +
1      |      root    |              |      scientific name |
 +
2      |      Bacteria        |      Bacteria <bacteria>    |      scientific name |
 +
2      |      bacteria        |              |      blast name      |
 +
2      |      eubacteria      |              |      genbank common name    |
 +
2      |      Monera  |      Monera <bacteria>      |      in-part |
 +
2      |      Procaryotae    |      Procaryotae <bacteria>  |      in-part |
 +
2      |      Prokaryotae    |      Prokaryotae <bacteria>  |      in-part |
 +
2      |      Prokaryota      |      Prokaryota <bacteria>  |      in-part |
 +
2      |      prokaryote      |      prokaryote <bacteria>  |      in-part |
 +
2      |      prokaryotes    |      prokaryotes <bacteria>  |      in-part |
 +
6      |      Azorhizobium Dreyfus et al. 1988 emend. Lang et al. 2013        |              |      authority      |
 +
6      |      Azorhizobium    |              |      scientific name |
 +
7      |      ATCC 43989      |      ATCC 43989 <type strain>        |      type material  |
 +
</pre>
 +
 +
Ce fichier correspond donc aux attributs des sommets de l'arbre avec notamment le nom scientifique d'un taxon.
 +
 +
Récupérer le fichier <tt>taxdump.tar.gz</tt> (copie [[silico:enseignement/m1/graph/taxdump.tar.gz|ici]] du 2019-12-03 contenant seulement les fichiers nécessaires).
 +
 +
'''Combien de sommets a l'arbre ?'''
 +
 +
'''Utiliser les commandes shell (notamment par exemple <tt>sed</tt>) pour se débarasser des espaces et remplacer les | par des tabulations dans le fichier <tt>nodes.dmp</tt>'''.
 +
 +
Récupérer le fichier <tt>Tree.py</tt> sur [https://gitlab.com/rbarriot/graph gitlab] pour continuer le TP. A utiliser pour charger le fichier <tt>nodes</tt> reformaté.
 +
 +
'''Ajouter une fonction <tt>height</tt> qui détermine la hauteur/profindeur de l'arbre à partir d'un sommet donné ou à partir de la racine sinon.'''
 +
 +
Lors d'analyses, comme ici, il arrive souvent de ne s'intéresser qu'à une partie de cet arbre. Pour ce TP, nous allons nous intéresser aux firmicutes. Pour cela, il faut retrouver l'identifiant de ce taxon (tax_id) dans le fichier <tt>names.dmp</tt> (par exemple avec la commande <tt>grep</tt>), afin de ne conserver que le sous-arbre de racine firmicutes.
 +
 +
'''Combien y a-t-il de firmicutes dans l'arbre ? Pour répondre à cette question, jouter une fonction <tt>nb_leaves</tt> qui détermine le nombre de feuilles à partir d'un sommet donné ou à partir de la racine sinon.'''
 +
 +
Dans l'étude qui nous concerne, toutes les bactéries firmicutes n'ont pas été considérées. Les nombres de gènes codant pour ''hscC'' et ''djlB/C'' ont été analysés. Vous trouverez ces données dans le fichier [[silico:enseignement/m1/graph/firmicutes.selection.tsv|firmicutes.selection.tsv]]
 +
 +
'''Charger ce fichier (code à la fin du fichier Tree.py) pour ne considérer que ces firmicutes, et utiliser la fonction <tt>prune</tt> pour élaguer l'arbre et ne conserver que le sous-arbre des firmicutes spécifié dans ce fichier.'''
 +
 +
'''Analyser le code de la fonction <tt>as_newick</tt> et l'utiliser la fonction pour reformater l'arbre au format newick i) en utilisant le tax_id pour les sommets et ii) en utilisant <tt>scientific name</tt> pour les sommets.'''
 +
 +
===iTOL===
 +
 +
Dans cette partie, nous allons principalement visualiser l'arbre précédemment généré.
 +
 +
Aller sur le site https://itol.embl.de
 +
 +
Cliquer sur le bouton <tt>Upload</tt> et téléverser le fichier newick précédent.
 +
 +
Il est possible d'ajouter des annotations (bouton + en bas à droite) sur l'arbre pour
 +
* colorer certains clades
 +
* indiquer la présence/absence de certaines propriétés sur les feuilles
 +
* ...
 +
 +
'''Remarque :''' il est possible d'utiliser l'arbre non élaguer (toute la taxonomie du NCBI) mais la visualisation est peu "réactive" et peu exploitable.
 +
 +
 +
(à faire : transformer le fichier firmicutes.selection.tsv en fichier d'annotation pour iTOL)
 +
 +
Cliquer sur <tt>Tree of Life</tt> et '''Télécharger (Controls → Export) le fichier newick correspondant.'''
 +
 +
=== ggtree ===
 +
 +
==== Prise en main ====
 +
 +
[https://bioconductor.org/packages/release/bioc/html/ggtree.html ggtree] est une librairie R/Bioconductor pour la '''visualisation''' et l''''annotation''' d'arbres phylogénétiques. Elle est basée sur <tt>ggplot2</tt>, une librairie pour faire des graphiques plus élaborés qu'avec la fonction <tt>plot</tt> de R, qu'il vous est conseillé d'expérimenter par vous-même.
[[Image:Tree_toy_example.png|thumb|Toy example. Télécharger au format [[Media:Tree_toy_example.svg|SVG]] ou [[Media:Tree_toy_example.nw|Newick]].]]
[[Image:Tree_toy_example.png|thumb|Toy example. Télécharger au format [[Media:Tree_toy_example.svg|SVG]] ou [[Media:Tree_toy_example.nw|Newick]].]]
Line 19: Line 255:
  ((B:0.2,J:0.4,(C:0.2,D:0.4,H:0.2,I:0.4)E:0.2)F:0.4,(K:0.2,L:0.4)G:0.2)A:0.2;
  ((B:0.2,J:0.4,(C:0.2,D:0.4,H:0.2,I:0.4)E:0.2)F:0.4,(K:0.2,L:0.4)G:0.2)A:0.2;
-
'''Remarque :''' <tt>read.tree</tt> au lien de <tt>read.newick</tt> (ape plus ancienne)
+
Chargement de l'arbre
-
 
+
<source lang="rsplus">
 +
t = read.tree('http://silico.biotoul.fr/site/images/f/f0/Tree_toy_example.nw')
 +
t
 +
</source>
Premier plot avec les paramètres par défaut
Premier plot avec les paramètres par défaut
-
```{r}
+
<source lang="rsplus">
ggtree(t)
ggtree(t)
-
```
+
</source>
N'affiche que l'arbre, sans les étiquettes ni autre.
N'affiche que l'arbre, sans les étiquettes ni autre.
Line 38: Line 277:
ggtree(t)+ geom_tiplab(size=3) + geom_text2(aes(subset=!isTip, label=nodelabels), hjust=-.5, size=3)  
ggtree(t)+ geom_tiplab(size=3) + geom_text2(aes(subset=!isTip, label=nodelabels), hjust=-.5, size=3)  
</source>
</source>
 +
'''Remarque :''' Il doit bien y avoir une méthode plus indiquée mais je ne l'ai pas trouvée en si peu de temps.
-
Affichage sans prendre en compte les longueurs de branche
+
Affichage sans prendre en compte les longueurs de branche (cladogramme)
<source lang="rsplus">
<source lang="rsplus">
ggtree(t, branch.length="none")
ggtree(t, branch.length="none")
Line 68: Line 308:
</source>
</source>
-
'''Ajout d'annotations'''
+
==== Tree of Life ====
-
Pour cette partie, allez sur https://itol.embl.de/itol.cgi pour télécharger l'arbre au format Newick.
+
Charger l'arbre iTOL précédemment exporté :
<source lang="rsplus">
<source lang="rsplus">
-
t = read.newick('iTOL.export.newick.txt')
+
t = read.tree('iTOL.export.newick.txt')
ggtree(t, layout="circular") + ggtitle("iTOL") + geom_tiplab(aes(angle=angle))+ geom_treescale()
ggtree(t, layout="circular") + ggtitle("iTOL") + geom_tiplab(aes(angle=angle))+ geom_treescale()
</source>
</source>
Line 102: Line 342:
   geom_hilight(node=archaea, fill="green", alpha=.1)  
   geom_hilight(node=archaea, fill="green", alpha=.1)  
</source>
</source>
 +
 +
Vous devriez obtenir la figure suivante :
 +
 +
[[Image:iTOL.annotated.png|800px]]
 +
 +
 +
==== Firmicutes ====
 +
 +
Librairies R nécessaires
 +
<source lang='rsplus'>
 +
library(ggtree)
 +
library(tidyverse)
 +
library(treeio)
 +
</source>
 +
 +
'''Charger et afficher les arbres des firmicutes sélectionnés dans la première partie''' avec les ''tax_id'' et puis avec les ''scientific name''.
 +
 +
Chargement des données associées dans un ''tibble'' ;
 +
<source lang='rsplus'>
 +
tree = read.tree("firmicutes.tax_id.nw")
 +
tb = read_tsv("firmicutes.selection.tsv") %>% select(label=taxonID, proteomeID, species:djlC)
 +
tb
 +
</source>
 +
 +
Pour afficher un label différent (nom des souches) sur les feuilles :
 +
<source lang='rsplus'>
 +
ggtree(tree, layout="circular") %<+% tb + geom_tiplab(aes(label=species), size=2)
 +
</source>
 +
 +
Pour le proteomeID :
 +
<source lang='rsplus'>
 +
ggtree(tree, layout="circular") %<+% tb + geom_tiplab(aes(label=proteomeID), size=2)
 +
</source>
 +
 +
Ajout d'annotations (nombre de ''hscC'' et 'djlB/C'' dans chaque génome) :
 +
<source lang='rsplus'>
 +
df = firm %>% select(hscC,djlC) %>% mutate(hscC=hscC>0, djlC=djlC>0) %>% as.data.frame
 +
df
 +
 +
p.rec = ggtree(tree) %<+% tb + geom_tiplab(aes(label=proteomeID), size=2)
 +
p.circ = ggtree(tree, layout="circular") %<+% tb + geom_tiplab(aes(label=proteomeID), size=2)
 +
 +
gheatmap(p.rec, df, offset=3, width=.25, font.size=3, colnames_angle=90, hjust=-0.2)
 +
gheatmap(p.circ, df, offset=3, width=.25, font.size=3, colnames_angle=90, hjust=-0.2)
 +
 +
</source>
 +
 +
 +
Pour aller plus loin : https://yulab-smu.github.io/treedata-book/
 +
 +
 +
==== Annexe ====
 +
 +
Installation de la librairie (nécessaire le 12/12/2017 et ne fonctionne que sous windows cette année)
 +
<source lang='rsplus'>
 +
source("https://bioconductor.org/biocLite.R")
 +
biocLite("ggtree")
 +
</pre>
== Graphes ==
== Graphes ==
Line 171: Line 469:
* Arbre couvrant de poids minimum (<tt>minimum.spanning.tree</tt>)
* Arbre couvrant de poids minimum (<tt>minimum.spanning.tree</tt>)
-
Fichier dressing au format ncol : [[Media:M1BBS_Graphe_dressing.ncol|dressing.ncol]] et celui pour Floyd-Warshall [[Media:M1BBS_Graphe_Graphe.Floyd-Warshall.ncol|Floyd-Warshall.ncol]]
+
Fichier dressing au format ncol : [[Media:M1BBS_Graphe_dressing.ncol|dressing.ncol]], [[Media:M1BBS_Graphe_Bellman-Ford.ncol|Bellman-Ford.ncol]] et celui pour Floyd-Warshall [[Media:M1BBS_Graphe_Graphe.Floyd-Warshall.ncol|Floyd-Warshall.ncol]]
-
* Parcours en largeur et en profondeur, et Belllman-Ford (<tt>graph.bfs</tt>, <tt>graph.dfs</tt>, <tt>distances</tt>)
+
[[Image:igraph.dressing.png|thumb|dressing]]
 +
[[Image:igraph.Bellman-Ford.png|thumb|Graphe utilisé pour l'algorithme Bellman-Ford]]
 +
 +
[[Image:igraph.Floyd-Warshall.png|thumb|Graphe utilisé pour l'algorithme Floyd-Warshall]]
 +
 +
* Effectuer le tri topologique du graphe dressing.
 +
 +
* Parcours en largeur et en profondeur, et Bellman-Ford (<tt>graph.bfs</tt>, <tt>graph.dfs</tt>, <tt>distances</tt>)
 +
 +
* la betweenness d'un sommet ou d'une arête est le nombre de plus courts chemins passant par le sommet ou l'arête. Utiliser les fonctions <tt>betweenness</tt> et  <tt>edge.betweenness</tt> pour calculer cette valeur et l'ajouter au dessin.
 +
 +
Pour la centralité des sommets : [[Image:igraph.node.betweenness.png|400px]]  Pour la centralité des arêtes : [[Image:igraph.edge.betweenness.png|400px]]
-
* la betweenness d'une arête est le nombre de plus courts chemins passant par cette arête. Utiliser la fonction <tt>edge.betweenness</tt> pour calculer cette valeur et l'ajouter au dessin.
 
==== Partitionnement de graphe ====
==== Partitionnement de graphe ====
Line 193: Line 501:
plot(com, g, vertex.size=6, vertex.label=NA)
plot(com, g, vertex.size=6, vertex.label=NA)
</source>
</source>
 +
 +
[[Image:igraph.edge.betweenness.communities.png|600px]]
décomposition en recherchant le partitionnement ayant meilleure modularité
décomposition en recherchant le partitionnement ayant meilleure modularité
Line 240: Line 550:
</source>
</source>
 +
Affichage interactif
<source lang='rsplus'>
<source lang='rsplus'>
visNetwork(nodes, links, width="1024",height="768", main="graph with visNetwork")
visNetwork(nodes, links, width="1024",height="768", main="graph with visNetwork")
Line 249: Line 560:
<!--[[File:R.libs.TP.Rmd]]-->
<!--[[File:R.libs.TP.Rmd]]-->
 +
 +
commandes
 +
cat taxdump/nodes.dmp  | sed 's/\s//g' | sed 's/|/\t/g' > taxdump/nodes.tsv

Revision as of 08:35, 16 January 2020

Contents

Arbres

Il existe plusieurs librairies disponibles : ape, genoPlotR, ggtree, ... Celle choisie pour ce TP est ggtree.


Taxonomie du NCBI

Nous allons utiliser la taxonomy du NCBI. Pour cela, dans l'onglet Downloads on peut aller sur le FTP pour télécharger la banque. Le fichier README (taxdump_readme.txt) contient les information suivantes (avec en gras, celles qui vont nous intéresser tout particulièrement) :

This directory contains the following NCBI Taxonomy database dump files:

  taxdmp.zip
  taxdump.tar.Z
  taxdump.tar.gz

All these files containes exactly the same information and are arranged so
for the convenience of unpacking them on various operating environments.
In addition there are files:

  taxdmp.zip.md5
  taxdump.tar.Z.md5
  taxdump.tar.gz.md5

which contain MD5 sums for the corresponding archive files. These files
might be used to check correctness of the download of corresponding 
archive file.

taxdmp.zip
----------

Is intended for zip-capable utilities such as pkunzip, unzip, and WinZip.
These utilities are widely available in almost all operating environments.
To unpack it command-line pkunzip and unzip:

       pkunzip taxdmp.zip
or
       unzip taxdmp.zip

Note: pkunzip and/or unzip executables must be in the executable search path
and taxdmp.zip must be in the current directory. Files will be unzipped into
current directory. For desired dump files placement and more please refer to
the manual and/or option descriptions of pkunzip and unzip utilities.

taxdump.tar.Z
-------------

This file is to be unpacked by uncompress utility and subsequent tar 
archiver. These utilities are usually used in UNIX-like environment. 
Unpacking instructions follows:

          uncompress -c taxdump.tar.Z | tar xf - 

taxdump.tar.gz
--------------

This file is to be unpacked by GNU unzip utility and subsequent tar 
archiver. These utilities are usually used in UNIX-like environment. 
Unpacking instructions follows:

          gunzip -c taxdump.tar.gz | tar xf - 

The content of the archive
--------------------------

It may look like this:

citations.dmp
delnodes.dmp
division.dmp
gencode.dmp
merged.dmp
names.dmp
nodes.dmp
readme.txt

The readme.txt file gives a brief description of *.dmp files. These files
contain taxonomic information and are briefly described below. Each of the
files store one record in the single line that are delimited by "\t|\n"
(tab, vertical bar, and newline) characters. Each record consists of one 
or more fields delimited by "\t|\t" (tab, vertical bar, and tab) characters.
The brief description of field position and meaning for each file follows.

nodes.dmp
---------

This file represents taxonomy nodes. The description for each node includes 
the following fields:

	tax_id					-- node id in GenBank taxonomy database
	parent tax_id				-- parent node id in GenBank taxonomy database
	rank					-- rank of this node (superkingdom, kingdom, ...) 
	embl code				-- locus-name prefix; not unique
	division id				-- see division.dmp file
	inherited div flag  (1 or 0)		-- 1 if node inherits division from parent
	genetic code id				-- see gencode.dmp file
	inherited GC  flag  (1 or 0)		-- 1 if node inherits genetic code from parent
	mitochondrial genetic code id		-- see gencode.dmp file
	inherited MGC flag  (1 or 0)		-- 1 if node inherits mitochondrial gencode from parent
	GenBank hidden flag (1 or 0)            -- 1 if name is suppressed in GenBank entry lineage
	hidden subtree root flag (1 or 0)       -- 1 if this subtree has no sequence data yet
	comments				-- free-text comments and citations
 
 names.dmp
 ---------
 Taxonomy names file has these fields:
 
	tax_id					-- the id of node associated with this name
	name_txt				-- name itself
	unique name				-- the unique variant of this name if name not unique
	name class				-- (synonym, common name, ...)  

division.dmp
------------
Divisions file has these fields:
	division id				-- taxonomy database division id
	division cde				-- GenBank division code (three characters)
	division name				-- e.g. BCT, PLN, VRT, MAM, PRI...
	comments
 
gencode.dmp
-----------
Genetic codes file:

	genetic code id				-- GenBank genetic code id
	abbreviation				-- genetic code name abbreviation
	name					-- genetic code name
	cde					-- translation table for this genetic code
	starts					-- start codons for this genetic code
 
delnodes.dmp
------------
Deleted nodes (nodes that existed but were deleted) file field:

 	tax_id					-- deleted node id
 
merged.dmp
----------
Merged nodes file fields:

	old_tax_id                              -- id of nodes which has been merged
	new_tax_id                              -- id of nodes which is result of merging

citations.dmp
-------------
Citations file fields:

	cit_id					-- the unique id of citation
	cit_key					-- citation key
       medline_id                              -- unique id in MedLine database (0 if not in MedLine)
	pubmed_id				-- unique id in PubMed database (0 if not in PubMed)
	url					-- URL associated with citation
	text					-- any text (usually article name and authors)
						-- The following characters are escaped in this text by a backslash:
						-- newline (appear as "\n"),
						-- tab character ("\t"),
						-- double quotes ('\"'),
						-- backslash character ("\\").
	taxid_list				-- list of node ids separated by a single space

Extrait de nodes.dmp (rappel des colonnes : tax_id, parent tax_id, rank, ...) :

1       |       1       |       no rank |               |       8       |       0       |       1       |       0       |       0       |       0       |       0       |       0       |               |
2       |       131567  |       superkingdom    |               |       0       |       0       |       11      |       0       |       0       |       0       |       0       |       0       |               |
6       |       335928  |       genus   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |
7       |       6       |       species |       AC      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
9       |       32199   |       species |       BA      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
10      |       1706371 |       genus   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |
11      |       1707    |       species |       CG      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
13      |       203488  |       genus   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |
14      |       13      |       species |       DT      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
16      |       32011   |       genus   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |
17      |       16      |       species |       MM      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
18      |       213421  |       genus   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |
19      |       18      |       species |       PC      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
20      |       76892   |       genus   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |
21      |       20      |       species |       PI      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
22      |       267890  |       genus   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |
23      |       22      |       species |       SC      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
24      |       22      |       species |       SP      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
25      |       22      |       species |       SH      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
27      |       49928   |       species |       HE      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
28      |       49928   |       species |       HE      |       0       |       1       |       11      |       1       |       0       |       1       |       1       |       0       |               |
29      |       28221   |       order   |               |       0       |       1       |       11      |       1       |       0       |       1       |       0       |       0       |               |

Ce fichier nous donne donc la structure de l'arbre taxonomique avec les couples tax_id et parent tax_id.

Extrait de names.dmp (rappel des colonnes : tax_id, name_txt, unique name, name class) :

1       |       root    |               |       scientific name |
2       |       Bacteria        |       Bacteria <bacteria>     |       scientific name |
2       |       bacteria        |               |       blast name      |
2       |       eubacteria      |               |       genbank common name     |
2       |       Monera  |       Monera <bacteria>       |       in-part |
2       |       Procaryotae     |       Procaryotae <bacteria>  |       in-part |
2       |       Prokaryotae     |       Prokaryotae <bacteria>  |       in-part |
2       |       Prokaryota      |       Prokaryota <bacteria>   |       in-part |
2       |       prokaryote      |       prokaryote <bacteria>   |       in-part |
2       |       prokaryotes     |       prokaryotes <bacteria>  |       in-part |
6       |       Azorhizobium Dreyfus et al. 1988 emend. Lang et al. 2013        |               |       authority       |
6       |       Azorhizobium    |               |       scientific name |
7       |       ATCC 43989      |       ATCC 43989 <type strain>        |       type material   |

Ce fichier correspond donc aux attributs des sommets de l'arbre avec notamment le nom scientifique d'un taxon.

Récupérer le fichier taxdump.tar.gz (copie ici du 2019-12-03 contenant seulement les fichiers nécessaires).

Combien de sommets a l'arbre ?

Utiliser les commandes shell (notamment par exemple sed) pour se débarasser des espaces et remplacer les | par des tabulations dans le fichier nodes.dmp.

Récupérer le fichier Tree.py sur gitlab pour continuer le TP. A utiliser pour charger le fichier nodes reformaté.

Ajouter une fonction height qui détermine la hauteur/profindeur de l'arbre à partir d'un sommet donné ou à partir de la racine sinon.

Lors d'analyses, comme ici, il arrive souvent de ne s'intéresser qu'à une partie de cet arbre. Pour ce TP, nous allons nous intéresser aux firmicutes. Pour cela, il faut retrouver l'identifiant de ce taxon (tax_id) dans le fichier names.dmp (par exemple avec la commande grep), afin de ne conserver que le sous-arbre de racine firmicutes.

Combien y a-t-il de firmicutes dans l'arbre ? Pour répondre à cette question, jouter une fonction nb_leaves qui détermine le nombre de feuilles à partir d'un sommet donné ou à partir de la racine sinon.

Dans l'étude qui nous concerne, toutes les bactéries firmicutes n'ont pas été considérées. Les nombres de gènes codant pour hscC et djlB/C ont été analysés. Vous trouverez ces données dans le fichier firmicutes.selection.tsv

Charger ce fichier (code à la fin du fichier Tree.py) pour ne considérer que ces firmicutes, et utiliser la fonction prune pour élaguer l'arbre et ne conserver que le sous-arbre des firmicutes spécifié dans ce fichier.

Analyser le code de la fonction as_newick et l'utiliser la fonction pour reformater l'arbre au format newick i) en utilisant le tax_id pour les sommets et ii) en utilisant scientific name pour les sommets.

iTOL

Dans cette partie, nous allons principalement visualiser l'arbre précédemment généré.

Aller sur le site https://itol.embl.de

Cliquer sur le bouton Upload et téléverser le fichier newick précédent.

Il est possible d'ajouter des annotations (bouton + en bas à droite) sur l'arbre pour

  • colorer certains clades
  • indiquer la présence/absence de certaines propriétés sur les feuilles
  • ...

Remarque : il est possible d'utiliser l'arbre non élaguer (toute la taxonomie du NCBI) mais la visualisation est peu "réactive" et peu exploitable.


(à faire : transformer le fichier firmicutes.selection.tsv en fichier d'annotation pour iTOL)

Cliquer sur Tree of Life et Télécharger (Controls → Export) le fichier newick correspondant.

ggtree

Prise en main

ggtree est une librairie R/Bioconductor pour la visualisation et l'annotation d'arbres phylogénétiques. Elle est basée sur ggplot2, une librairie pour faire des graphiques plus élaborés qu'avec la fonction plot de R, qu'il vous est conseillé d'expérimenter par vous-même.

Toy example. Télécharger au format SVG ou Newick.

Au format Newick :

((B:0.2,J:0.4,(C:0.2,D:0.4,H:0.2,I:0.4)E:0.2)F:0.4,(K:0.2,L:0.4)G:0.2)A:0.2;

Chargement de l'arbre

t = read.tree('http://silico.biotoul.fr/site/images/f/f0/Tree_toy_example.nw')
t

Premier plot avec les paramètres par défaut

ggtree(t)

N'affiche que l'arbre, sans les étiquettes ni autre.

ggtree(t)+ geom_tiplab(size=3)

Et les étiquettes des noeuds internes :

nodelabels = c(t$tip.label, t$node.label)
ggtree(t)+ geom_tiplab(size=3) + geom_text2(aes(subset=!isTip, label=nodelabels), hjust=-.5, size=3)

Remarque : Il doit bien y avoir une méthode plus indiquée mais je ne l'ai pas trouvée en si peu de temps.

Affichage sans prendre en compte les longueurs de branche (cladogramme)

ggtree(t, branch.length="none")

Différents algorithmes de dessin

ggtree(t) + ggtitle("default: rectangular")
ggtree(t, layout="slanted") + ggtitle("slanted")
ggtree(t, layout="circular") + ggtitle("circular")

Avec ou sans racine

ggtree(t, layout="unrooted") + ggtitle("unrooted layout")

Avec l'échelle pour les longueurs de branche

ggtree(t) + geom_treescale()

La même chose avec une thème différent et en stockant le plot dans une variable

p = ggtree(t) + theme_tree2() + ggtitle("rectangular tree with branch lengths and scale")
p
p + geom_tiplab(size=3, color="orange")

Tree of Life

Charger l'arbre iTOL précédemment exporté :

t = read.tree('iTOL.export.newick.txt')
ggtree(t, layout="circular") + ggtitle("iTOL") + geom_tiplab(aes(angle=angle))+ geom_treescale()

annotations :

t$tip.label[1:10]
t$node.label[1:10]

Recherche des n° de sommets correspondants aux ancêtres des bactéries, eucaryotes, et archées :

all=c(t$tip.label, t$node.label)
which(all=='Bacteria')
bacteria = which(all=='Bacteria')[1]
which(all=='Eukaryota')
eukaryota = which(all=='Eukaryota')[1]
which(all=='Archaea')
archaea = which(all=='Archaea')[1]

plot avec l'annotation

ggtree(t, layout="circular") + 
  ggtitle("iTOL") + 
  geom_tiplab(aes(angle=angle))  + 
  geom_hilight(node=bacteria, fill="steelblue", alpha=.1) + 
  geom_hilight(node=eukaryota, fill="pink", alpha=.1) + 
  geom_hilight(node=archaea, fill="green", alpha=.1)

Vous devriez obtenir la figure suivante :


Firmicutes

Librairies R nécessaires

library(ggtree)
library(tidyverse)
library(treeio)

Charger et afficher les arbres des firmicutes sélectionnés dans la première partie avec les tax_id et puis avec les scientific name.

Chargement des données associées dans un tibble ;

tree = read.tree("firmicutes.tax_id.nw")
tb = read_tsv("firmicutes.selection.tsv") %>% select(label=taxonID, proteomeID, species:djlC)
tb

Pour afficher un label différent (nom des souches) sur les feuilles :

ggtree(tree, layout="circular") %<+% tb + geom_tiplab(aes(label=species), size=2)

Pour le proteomeID :

ggtree(tree, layout="circular") %<+% tb + geom_tiplab(aes(label=proteomeID), size=2)

Ajout d'annotations (nombre de hscC et 'djlB/C dans chaque génome) :

df = firm %>% select(hscC,djlC) %>% mutate(hscC=hscC>0, djlC=djlC>0) %>% as.data.frame
df
 
p.rec = ggtree(tree) %<+% tb + geom_tiplab(aes(label=proteomeID), size=2) 
p.circ = ggtree(tree, layout="circular") %<+% tb + geom_tiplab(aes(label=proteomeID), size=2) 
 
gheatmap(p.rec, df, offset=3, width=.25, font.size=3, colnames_angle=90, hjust=-0.2)
gheatmap(p.circ, df, offset=3, width=.25, font.size=3, colnames_angle=90, hjust=-0.2)


Pour aller plus loin : https://yulab-smu.github.io/treedata-book/


Annexe

Installation de la librairie (nécessaire le 12/12/2017 et ne fonctionne que sous windows cette année)

source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")
</pre>
 
== Graphes ==
 
=== igraph ===
 
La librairie [http://igraph.sourceforge.net iGraph] met à disposition tout un ensemble de fonctions pour le traitement et la visualisation de graphes. Nous allons utiliser aujourd'hui son interfaçage avec R. Pour la charger :
<source lang='rsplus'>
library(igraph)

Pour charger un graphe (différents format possibles : pajek, newick, ...) :

g = read.graph("http://silico.biotoul.fr/site/images/9/9f/Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.tgr", directed=FALSE)

Consulter l'aide de la fonction (?read.graph) pour voir les autres formats supportés.


Pour l'afficher, il faut au préalable en effectuer le dessin (layout) :

# soit en une ligne en passant la fonction de dessin :
plot(g, layout=layout.fruchterman.reingold)
 
# soit en sauvegardant ce layout dans une variable :
g.FR = layout.fruchterman.reingold(g)
plot(g, layout=g.FR, vertex.size=3, vertex.label=NA)

Consulter l'aide des fonction plot.igraph et layout.fructhterman.reingold pour voir les options ainsi que les autres algorithmes de dessin disponibles.

Vous trouverez la documentation de la librairie sur le site dédié. Pour celle de l'interface R en ligne : http://igraph.org/r/

  • Pour obtenir la liste des sommets : V(g)
  • la liste des arêtes : E(g)


  • Quel est l'ordre du graphe ? Combien a-t-il d'arêtes ?


On peut assigner des étiquettes aux sommets : V(g)$name = vector_of_labels

Charger les étiquettes des sommets :

V(g)$name=as.character(read.table("http://silico.biotoul.fr/site/images/6/61/Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.cod")[,2])
plot(g, layout=g.FR, vertex.size=3, vertex.label=V(g)$name)

Les paramètres de la fonction plot de igraph sont décrits dans la documentation.

Il est possible de stocker de l'information sur le graphe, les sommets et/ou les arêtes avec les fonctions dédiées :

vertex_attr(g, name="name")
vertex_attr(g, name="name", index=10)


  • Quel est le diamètre du graphe ?
  • Lister les points d'articulation (articulation.points)
  • Longueur moyenne des plus courts chemins (sans valuation) (average.path.length)
  • Afficher sa représentation canonique (canonical.permutation) et la matrice d'adjacence correspondante
  • Lister les cliques maximales (maximal.cliques)
  • Lister les composantes connexes (clusters) avant et après suppression du point d'articulation
  • Obtenir le line graph (line.graph)
  • Arbre couvrant de poids minimum (minimum.spanning.tree)

Fichier dressing au format ncol : dressing.ncol, Bellman-Ford.ncol et celui pour Floyd-Warshall Floyd-Warshall.ncol

dressing
Graphe utilisé pour l'algorithme Bellman-Ford
Graphe utilisé pour l'algorithme Floyd-Warshall
  • Effectuer le tri topologique du graphe dressing.
  • Parcours en largeur et en profondeur, et Bellman-Ford (graph.bfs, graph.dfs, distances)
  • la betweenness d'un sommet ou d'une arête est le nombre de plus courts chemins passant par le sommet ou l'arête. Utiliser les fonctions betweenness et edge.betweenness pour calculer cette valeur et l'ajouter au dessin.

Pour la centralité des sommets : Pour la centralité des arêtes :


Partitionnement de graphe

Pour partitionner le graphe en communautés, différentes méthodes sont disponibles. Vous allez utilisez la betweenness des arêtes pour effectuer le partitionnement du graphe. Cette méthode sélectionne les arêtes dont la betweenness est la plus importantes afin de former des communautés (clusters).

# community detection with edge-betweenness
com = edge.betweenness.community(g)
par(cex=.5)
plot(as.dendrogram(com))

Affichage du meilleur partitionnement

plot(com, g, vertex.size=6, vertex.label=NA)

décomposition en recherchant le partitionnement ayant meilleure modularité

com$best_step=0
com$best_modularity = -Inf
for (i in 1:nrow(com$merges)) {
  #igraph old version: ctm = community.to.membership(g, com$merges, steps=i)
  #igraph old version: mod = modularity(g, ctm$membership+1) # faire ctm$membership+1 si crash il y a (dépend de la version igraph, membership commence à 1 ou 0)
  ctm = cut_at(com, steps=i)
  mod = modularity(g, ctm)
  print(paste("step: ",i,", modularity:",mod))
  if (mod>com$best_modularity) {
    com$best_step = i
    com$best_modularity = mod
    #igraph old version: com$membership=ctm$membership
    com$membership=ctm
  }
}
com$best_step ; com$best_modularity

Autre : marquage de certains sommets (au hasard la 1ère communauté)

gr1 = groups(com)[1]
plot(g, mark.groups=gr1, vertex.label=NA)

librairie visNetwork

library(visNetwork)


conversion du graphe précédent

nodes = data.frame(id = as.vector(V(g)$name))
links = as.data.frame(get.edgelist(g))
colnames(links) = c('from','to')

paramètres pour l'affichage

links$smooth = F
nodes$color.background <- c("red", "green", "blue", "orange", "yellow")[membership(com)]

Affichage interactif

visNetwork(nodes, links, width="1024",height="768", main="graph with visNetwork")


Pour aller plus loin : http://kateto.net/network-visualization


commandes

cat taxdump/nodes.dmp  | sed 's/\s//g' | sed 's/|/\t/g' > taxdump/nodes.tsv