silico.biotoul.fr
 

Atelier Phylogénomique OrthoFinder

From silico.biotoul.fr

Jump to: navigation, search

Contents

Liens

MSK

Paramètres

Nous allons utiliser le programme OrthoFinder avec les paramètres par défaut (config.json].

  • diamond: sequence search program
  • DendroBLAST: gene tree inference

Mais vous pouvez utiliser n'importe quel programme d'alignement, de reconstruction d'arbre ou de comparaison de séquences que vous préférez. Pour utiliser un autre programme, il suffit de modifier le fichier de configuration appelé config.json dans le répertoire orthofinder. Le mieux est de créer un fichier au même format appelé config_orthofinder_user.json dans votre répertoire utilisateur. 

/usr/local/bioinfo/src/OrthoFinder/OrthoFinder-2.5.2/config.json

Mise en œuvre

How_to_use_SLURM_OrthoFinder

Nous utilisons le script donné en exemple dans /usr/local/bioinfo/src/OrthoFinder/example_on_cluster comme modèle.

Créer un répertoire ~/work/OrthoFinder et copier le script test_OrthoFinder-2.2.1.sh dans ce répertoire en changeant pour la version utilisée. 

cp /usr/local/bioinfo/src/OrthoFinder/example_on_cluster/test_OrthoFinder-2.2.1.sh ~/work/OrthoFinder/prochlo_OrthoFinder-2.5.4.sh

Créer un sous-répertoire Prochlorococcus et copier les fichiers peptides issues de Prokka dans ce repertoire. 

cp ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.faa ~/work/OrthoFinder/Prochlorococcus/.

Le script copié est édité pour changer le répertoire de travail et la version du programme. 

-X: option pour ne pas ajouter les noms des espèces aux noms des séquences.
-n: suffixe à ajouter au répertoire contenant les résultats.

Options à utiliser dans différents runs. 

-S: programme à utiliser pour aligner les séquences.
-M: méthode pour inférer les arbres.

Pro_OrthoFinder-2.5.4.sh 

sbatch Pro_OrthoFinder-2.5.4.sh
 
squeue -l -u $USER

Par défaut, OrthoFinder crée un répertoire OrthoFinder dans le répertoire du protéome d'entrée (Prochlorococcus) et y place les résultats.

Fichiers de sorties

ls -l ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*

OrthoFinder produit un ensemble de fichiers décrivant

  • les orthogroupes,
  • les orthologues,
  • les arbres gènes,
  • les arbres gènes avec les évènements de duplications/délétions,
  • l'arbre des espèces enracinées,
  • les événements de duplication de gènes,
  • les statistiques de génomiques comparatives pour toutes les espèces analysées.

Remarque: les fichiers .tsv d'OrthoFinder peuvent être visualisés dans un tableur comme Excel ou LibreOffice Calc.

Résultats préliminaires

Les résultats de cette première étape sont basés sur les premières versions d'OrthoFinder (Emms and Kelly, 2015).

OrthoFinder est un algorithme qui infère les orthogroupes entre plusieurs espèces. La méthode ne nécessite pas d'informations sur la synténie et peut donc être utilisée avec des séquences de protéines prédites 1) à partir d'ensembles de données sur le génome ou 2) le transcriptome. Un orthogroupe est l'ensemble des gènes dérivés d'un seul gène dans le dernier ancêtre commun de toutes les espèces considérées. C'est l'approche utilisée par OrthoMCL.

Les étapes :

  1. Alignement de toutes les protéines contre toutes les protéines (diamond, BLAST, autre).
  2. Normalisation de la longueur des protéines et de la distance phylogénétique des scores binaires pour obtenir les scores à utiliser pour l'inférence des orthogroupes.
  3. Sélection de paires de protéines apparentées putatifs à partir des scores normalisés.
  4. Construction d'un graphe d'orthogroupe, les protéines sont des sommets dans le graphe et les paires de gènes sont reliées par une arête dont le poids est donné par le score binaire normalisé.
  5. Regroupement des gènes en orthogroupes discrets à l'aide de MCL.

Les résultats de cette première étape sont sauvegardées dans le répertoire Orthogroups.

Un critère pour évaluer la qualité des groupes de gènes orthologues est de calculer le nombre de paralogues par OG. Cette information peut être calculée à partir du fichier :

  • ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Orthogroups/Orthogroups.txt (Orthogroups.txt).

Taux de paralogie (avec awk mais en vous inspirant du TP5 de M1):

MSK NF is a predefined variable whose value is the number of fields in the current record. 

cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Orthogroups/Orthogroups.txt | awk '{
  j=0;
  SUM = NF-1;
  for (i=2; i<=NF; i++) {
    str=substr($i,1,4);
    strlist[str]++;
    j++;
  }
  if ( j > 12 ) {
    nb=0;
    for (str in strlist ) {nb++}
    par=SUM-nb
    freq=par/SUM
    printf("%s\t%i\t%i\t%i\t%5.3f", $1, SUM, nb, par, freq)
    for (str in strlist ) {
      printf("\t%s:%s", str, strlist[str]);
    }
    printf("\n")
  }
  delete strlist;
}'
Question 1:
Commentez les résultats obtenus.

Vue d'ensemble

Nous avons un résumé de l'analyse réalisée dans le fichier

  • ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Statistics_Overall.tsv (Statistics_Overall.tsv). 
cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Statistics_Overall.tsv

La plupart des termes figurant dans les fichiers 'Statistics_Overall.csv' et 'Statistics_PerSpecies.csv' sont explicites, les autres sont définis ci-dessous.

  • Species-specific orthogroup : Un orthogroupe entièrement constitué de gènes d'une seule espèce.
  • G50 : Le nombre de gènes dans l'orthogroupe tel que 50% des gènes sont dans des orthogroupes de cette taille ou plus.
  • O50 : Le plus petit nombre d'orthogroupes tel que 50% des gènes sont dans des orthogroupes de cette taille ou plus.
  • Single-copy orthogroup : Un orthogroupe comportant exactement un gène (et pas plus) de chaque espèce. Ces orthogroupes sont idéaux pour inférer un arbre des espèces et de nombreuses autres analyses.
  • Unassigned gene : Un gène qui n'a pas été placé dans un orthogroupe avec d'autres gènes.

La première chose à vérifier est de savoir combien de gènes ont été assignés à des orthogroupes. En général, il est bon de voir au moins 80% des gènes assignés à des orthogroupes. 

Question 2:
Quelles remarques pouvez vous faire sur le lien entre la taille des génomes et les différentes statistiques sur les orthologues?

Nous pouvons également obtenir les statistiques pour chaque souche (Statistics_PerSpecies.tsv). 

cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Statistics_PerSpecies.tsv

Liste des OG ayant le plus grand nombre de duplications : 

sort -n -k 2 ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Duplications_per_Orthogroup.tsv

Distribution du nombre de duplications par OG :

MSK 

cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Duplications_per_Orthogroup.tsv | awk '
BEGIN{ bin_width=1; nb=0 }
{
   if ( $1 ~ /OG/ ) { 
      bin=int($2/bin_width );
      if(bin in hist){hist[bin]+=1} else {hist[bin]=1}
      nb=nb+1;
   }
}
END{
    for (h in hist )
        printf " * > %2.2f  ->  %i %f\n", h*bin_width, hist[h], hist[h]/nb
}' | sort -n -k 3

Question 4:
Comparez les taux de paralogies observés avec le nombre de duplications estimées par OG.

Le fichier :

  • ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Orthogroups_SpeciesOverlaps.tsv (Orthogroups_SpeciesOverlaps.tsv])

est un fichier texte séparé par des tabulations qui contient le nombre d'orthogroupes partagés entre chaque paire d'espèces sous forme de matrice carrée. 

mkdir ~/work/OrthoFinder/images
cd ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics

MSK 

srun --pty bash # si nécessaire!
module load system/R-4.1.1_gcc-9.3.0
R
pdf_file <- '~/work/OrthoFinder/images/Orthogroups_SpeciesOverlaps.pdf'
Orthogroups_SpeciesOverlaps <- 'Orthogroups_SpeciesOverlaps.tsv'
 
data <- read.delim(file=Orthogroups_SpeciesOverlaps, sep="\t", header=TRUE, row.names=1)
 
pdf(file=pdf_file, paper="a4r")
heatmap(t(as.matrix(data)), scale='none', xlab="Strains", labCol=NA)
dev.off()
cat(pdf_file, "\n")
Question 5:
Quel type de lien entre souches ces données peuvent-elles vous suggérer?

Rappel sur OrthoFinder2

Les nouvelles version d'OrthoFinder infère les HOG, les orthogroupes à chaque niveau hiérarchique (c'est-à-dire à chaque nœud de l'arbre des espèces) en analysant les arbres des gènes enracinés. Cette méthode d'inférence des orthogroupes est beaucoup plus précise que l'approche basée sur la similarité des gènes/graphes utilisée par toutes les autres méthodes (Emms and Kelly, 2019).

Il repose sur les étapes suivantes :

  1. Inférence d'orthogroupes à l'aide de l'algorithme original OrthoFinder.
  2. Inférence des arbres gènes pour chaque orthogroupe.
  3. Analyse de ces arbres de gènes pour en déduire l'arbre des espèces enraciné.
  4. Enracinement des arbres de gènes à l'aide de l'arbre des espèces enraciné.
  5. Inférence des orthologues et des événements de duplication de gènes à l'aide des arbres gènes.
  6. Cartographie des événements de duplication de gènes sur les arbres des espèces et des gènes).

Arbre des espèces

Cet arbre a été inféré par OrthoFinder à l'aide de l'algorithme STAG (Species Tree inference from All Genes) et raciné à l'aide de l'algorithme STRIDE (Species Tree Root Inference from Gene Duplication Events). Il est donc prêt à être interprété (normalement, vous devriez d'abord raciner un arbre vous-même). Les fichiers sont : 

ls -l ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Species_Tree
  • SpeciesTree_rooted.txt : Un arbre des espèces STAG déduit de tous les orthogroupes, contenant des valeurs de support STAG aux nœuds internes et enraciné à l'aide de STRIDE.
  • SpeciesTree_rooted_node_labels.txt Le même arbre mais avec des étiquettes aux noeuds pour permettre de faire des références croisées entre les branches/noeuds de l'arbre des espèces (par exemple, l'emplacement des événements de duplication de gènes).

Pour visualiser ces arbres, vous pouvez utiliser iTOL, figtree, seaview ou la librairie ape de R. 

library(ape)
treefile <-"SpeciesTree_rooted.txt" 
tr <- read.tree(treefile)
plot(tr)

MSK

Arbres des gènes

Nous allons nous intéresser à l'orthogroup OG0000017 qui renferme 20 gènes dont 8 paralogues.

  • ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Gene_Trees/OG0000017_tree.txt 
Question 6:
Quel est le nombre minimum de duplication qui permet d'expliquer la composition de cet OG?

Tracez l'arbre des espèces en vis-à-vis de l'arbre des gènes du OG OG0000017 et positionnez les évènements de duplication et de perte de gènes sur l'arbre espèce.

MSK

Le fichier suivant (Duplications_per_Orthogroup.tsv)vous donne une estimation du nombre de duplications pour chaque OG. 

grep OG0000017 ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Duplications_per_Orthogroup.tsv

et le fichier Duplications.tsv, la prédiction de la localisation des duplications sur l'arbre espèces et les arbres gènes : 

grep OG0000017 ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_*/Comparative_Genomics_Statistics/Gene_Duplication_Events/Duplications.tsv
Commentez ces différents résultats.

Orthogroupes à chaque niveau hiérarchique

Le répertoire Phylogenetic_Hierarchical_Orthogroups renferme un fichier par nœud de l'arbre espèce. Le fichier N0.tsv correspond à la racine de l'arbre. Chaque ligne renferme les gènes qui appartiennent à un unique orthogroup identifié par un HOG ID. Ce fichier remplace avantageusement les orthogroups identifiés par la première étape (MCL)!

Le nombre de HOG du fichier N0.tsv est à comparer aux données du fichier Statistics_Overall.tsv. 

awk '{ if ($2~"OG0000017") print $0}' ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Phylogenetic_Hierarchical_Orthogroups/N*.tsv

nous donne la décomposition de OG0000017 en HOG sur les différentes feuilles de l'arbre espèces.

MSK 

cat  ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Phylogenetic_Hierarchical_Orthogroups/N0.tsv | awk -F $'\t' '{
if ( $1~/N0.HOG/ ) {
  nb=0;
  HOG=$1
  for (i=4; i<=NF; i++) {
    if ( $i~/[A-Z]/ ) {
       n=split($i, paral, ",");
       str=substr(paral[1],1,4);
       strlist[str]=n;
	   #print $i, str, strlist[str]
       delete paral;
       nb=nb+1;
	}
  }
 
  if ( nb > 0 ) {
     sum=0;
     for (str in strlist ) {sum=sum+strlist[str]}
     par=sum-nb
     freq=par/sum
     printf("%s\t%i\t%i\t%i\t%5.3f", HOG, sum, nb, par, freq)
     for (str in strlist ) {
       printf("\t%s:%s", str, strlist[str]);
     }
     printf("\n")
  }
  delete strlist;
}
}' | sort -g -k 5 > tmp


Histogramme du taux de paralogie

MSK 

cat tmp | awk '
BEGIN{ bin_width=0.1; nb=0 }
{
   if ( $1 ~ /OG/ ) { 
      bin=int($5/bin_width );
      if(bin in hist){hist[bin]+=1} else {hist[bin]=1}
      nb=nb+1;
   }
}
END{
    for (h in hist )
        printf " * > %2.2f  ->  %i %f\n", h*bin_width, hist[h], hist[h]/nb
}' | sort -n -k 3

--> MSK 

 * > 0.00  ->  2896 0.899379
 * > 0.10  ->  63 0.019565
 * > 0.20  ->  61 0.018944
 * > 0.30  ->  68 0.021118
 * > 0.40  ->  13 0.004037
 * > 0.50  ->  80 0.024845
 * > 0.60  ->  22 0.006832
 * > 0.70  ->  5 0.001553
 * > 0.80  ->  10 0.003106
 * > 0.90  ->  2 0.000621

-->

Somme des évènements de duplication

Le fichier SpeciesTree_Gene_Duplications_0.5_Support.txt indique à chaque noeud de l'arbre le nombre d'évènements de duplication (numero du noeud suivit du nombre). Les duplications sont considérées comme bien supportées si au moins 50% des espèces présentes chez les deux descendants ont retenues les deux copies du gène dupliqué.

MSK 

library(ape)
treefile <-"/home/quentin/Atelier_phylo/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Gene_Duplication_Events/SpeciesTree_Gene_Duplications_0.5_Support.txt" 
tr <- read.tree(treefile)
pdf_file <- "/home/quentin/Atelier_phylo/work/OrthoFinder/Figures/SpeciesTree_Gene_Duplications_0.5_Support.pdf"
pdf(file=pdf_file, paper="a4r")
plot(tr, show.node.label=T)
dev.off()

Comparer le nombre de duplication estimée au niveau des noeuds de l'arbre par rapport à celles estimée au niveau des feuilles. 

Que pouvez-vous en conclure?

cd ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Gene_Duplication_Events
awk '{ if ( $5!~/^Terminal/ && $4 >= 0.5 ) {print $_}}' Duplications.tsv | head -n -1 | wc -l 
awk '{ if ( $5~/^Terminal/ && $4 >= 0.5 ) {print $_}}' Duplications.tsv | wc -l

Comparaison des OG avec les HOG

MSK

Comparaison des résultats obtenus avec différents paramètres

You can use any alignment or tree inference program you like the best! For example, to you muscle and iqtree, the command like arguments you need to add are: "-M msa -A muscle -T iqtree" Paramètres: 

-M msa : trees are inferred using multiple sequence alignments (MSA) by using the option. By default MAFFT is used to generate the MSAs and FastTree to generate the gene trees.


MSK

Ajout des Synechococcus

OrthoFinder permet l'ajout de nouveaux génomes sans avoir à recalculer les comparaisons all vs all déjà réalisées.

ProSyn_OrthoFinder-2.5.4.sh 

cd ~/work/OrthoFinder
mkdir Synechococcus
cp ~/work/Synechococcus/peptide/*.faa Synechococcus/.

MSK

MSK 

Quels pourraient être les rôles des génomes de Synechococcus dans l'analyse des génomes de Prochlorococcus ?
Retrieved from "http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_OrthoFinder"