silico.biotoul.fr
 

Atelier Phylogénomique OrthoFinder

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Digramme de Venn)
m (HOG)
Line 395: Line 395:
cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Phylogenetic_Hierarchical_Orthogroups/N0.tsv | awk -F $'\t' '{
cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Phylogenetic_Hierarchical_Orthogroups/N0.tsv | awk -F $'\t' '{
if ( $1~/N0.HOG/ ) {
if ( $1~/N0.HOG/ ) {
-
if ($2~"OG0000001") {
+
#if ($2~"OG0000001") {
   j=0;
   j=0;
   for (i=4; i<=NF; i++) {
   for (i=4; i<=NF; i++) {
Line 411: Line 411:
  ncol=j;
  ncol=j;
  if ( ncol > 1 ) {
  if ( ncol > 1 ) {
 +
  n=asort(col);
   for (k=1; k<=ncol; k++) {
   for (k=1; k<=ncol; k++) {
   for (m=k+1; m<=ncol; m++) {
   for (m=k+1; m<=ncol; m++) {
Line 419: Line 420:
  delete col
  delete col
}
}
-
}
+
#}
}' > ~/work/OrthoFinder/Figures/HOGN0_pairsOG01.txt
}' > ~/work/OrthoFinder/Figures/HOGN0_pairsOG01.txt
</pre>
</pre>

Revision as of 15:45, 18 October 2021

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 prochlo_OrthoFinder-2.5.2.sh

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

cp ~/work/Prochlorococcus/peptide/*.faa 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.2.sh

sbatch Pro_OrthoFinder-2.5.2.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

ll ~/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 tous les gènes contre tous les gènes (BLAST).
  2. Normalisation de la longueur des gènes et de la distance phylogénétique des scores binaires de BLAST pour obtenir les scores à utiliser pour l'inférence des orthogroupes.
  3. Sélection de paires de gènes apparentés putatifs à partir des scores BLAST normalisés.
  4. Construction d'un graphe d'orthogroupe, les gènes sont des nœuds 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 Orthogroups.txt.

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

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 u résumé de l'analyse réalisée dans le fichier 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 :

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

srun --pty bash # si nécessaire!

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

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 et raciné à l'aide de l'algorithme STRIDE. Il est donc prêt à être interprété (normalement, vous devriez d'abord raciner un arbre vous-même). Les fichiers sont :

ll ~/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 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.

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.

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

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

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

Histogramme du taux de paralogie

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

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

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?
awk '{ if ( $5!~/^Terminal/ && $4 >= 0.5 ) {print $_}}' Gene_Duplication_Events/Duplications.tsv | wc
awk '{ if ( $5~/^Terminal/ && $4 >= 0.5 ) {print $_}}' Gene_Duplication_Events/Duplications.tsv | wc

Comparaison des OG avec les HOG

Diagramme de Venn

OG

cat  ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Orthogroups/Orthogroups.txt | awk '
{
if ($1~"OG0000001") {
 j=0;
 for (i=2; i<=NF; i++) {
  j++;
  col[j]=$i;
 }
 ncol=j;
 if ( ncol > 1 ) {
  for (k=1; k<=ncol; k++) {
   for (m=k+1; m<=ncol; m++) {
    printf("%s%s\n", col[k], col[m]);
   }
  }
 }
}
 delete col
}' > ~/work/OrthoFinder/Figures/Orthogroups_pairsOG01.txt

HOG

cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/Phylogenetic_Hierarchical_Orthogroups/N0.tsv | awk -F $'\t' '{
if ( $1~/N0.HOG/ ) {
#if ($2~"OG0000001") {
  j=0;
  for (i=4; i<=NF; i++) {
    if ( $i~/[A-Z]/ ) {
       split($i, paral, ", ");
       for ( k in paral ) {
          j++;
          col[j]=paral[k];

       }
       delete paral;
    }
  }
  
 ncol=j;
 if ( ncol > 1 ) {
  n=asort(col);
  for (k=1; k<=ncol; k++) {
   for (m=k+1; m<=ncol; m++) {
     printf("%s%s\n", col[k], col[m]);
   }
  }
 }
 delete col
}
#}
}' > ~/work/OrthoFinder/Figures/HOGN0_pairsOG01.txt

Digramme de Venn

srun --pty bash
module load system/R-3.5.1
R
library('gplots')
setwd("~/work/OrthoFinder/Figures")
pdf_file <- 'Venn_Orthogroups_HOGN0OG01.pdf'
Orthogroups <- read.table('Orthogroups_pairsOG01.txt')
HOGN0 <- read.table('HOGN0_pairsOG01.txt')
input <- list(Orthogroups=Orthogroups, HOGN0=HOGN0)
pdf(file=pdf_file, paper="a4r")
venn(input)
dev.off()

MSK

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

MSK

Ajout des Synechococcus

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

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

Quels pourraient être les rôles des génomes de Synechococcus dans l'analyse des génomes de Prochlorococcus ?