silico.biotoul.fr
 

Atelier Phylogénomique Arbre espèces

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Liens)
m (Liens)
Line 1: Line 1:
==Liens==
==Liens==
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces Arbre espèces]
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces Arbre espèces]
 +
==Introduction==
 +
'''Support de cours :''' [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/2020_supports/ supports]
 +
 +
Nous allons utiliser un sous ensemble de gènes concervés chez ''Prochlorococcus'' et ''Synechococcus'' pour expérimenter les différentes méthodes de reconstruction phylogénomiques. Nous nous initierons à la comparaison d’arbres.
 +
==Extraction des séquences nucléotidiques des gènes orthologues==
 +
Créer un fichier avec toutes les séquences nucléotidiques (cds):
 +
 +
Si le répertoire n'existe pas déjà, le créer.
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
mkdir -p ~/work/ProchlorococcusSynechococcus/DNA
 +
cat ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.ffn ~/work/Synechococcus/prokka/Aaa*/Aaa*.ffn > ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn
 +
 +
grep -c '>' ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn
 +
</pre>
 +
===Extraire les groupes de gènes orthologues===
 +
Exemple de la création d'un script et lancement du job avec sbatch
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA
 +
</pre>
 +
Créer le fichier suivant ''commande.sh''.
 +
<pre>
 +
#!/usr/bin/bash
 +
~/work/scripts/extract_sequences_from_matchtable.pl --fasta ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn --matchtable ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt --outdir ~/work/ProchlorococcusSynechococcus/OG/DNA --quorum 16
 +
</pre>
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">sbatch commande.sh
 +
</pre>
 +
Vérifier que les fichiers contiennent au moins 16 séquences:
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
 +
</pre>
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.1:
 +
A quoi correspond ce quorum. Pourquoi utiliser un seuil de 16?
 +
Combien de fichiers avez-vous obtenus?
 +
Est-ce pertinent dans les situations où vous avez un grand nombre de souches?
 +
Et quand les espèces étudiées sont très éloignées phylogénétiquement parlant ?
 +
</pre>
 +
===Transformation en alignement multiples des séquences ADN des gènes===
 +
Le script ''aa_to_dna_aln.pl'' prend en entrée un fichier fasta contenant des séquences d'ADN. Il traduit des séquences, réalise un alignement multiple avec muscle et retraduit en ADN cet alignement multiple ([https://bioperl.org/ BioPerl]).
 +
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
mkdir ~/work/ProchlorococcusSynechococcus/OG/alignment
 +
</pre>
 +
Attention à faire un script et le lancer en sbatch ou en srun --pty bash
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> ~/work/scripts/aa_to_dna_aln.pl -dna ~/work/ProchlorococcusSynechococcus/OG/DNA/OG_1471.fas --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment
 +
</pre>
 +
 +
La traduction en peptides ajoute un '*' à la fin des séquences. Vous pouvez ignorer le message:
 +
*** WARNING *** Invalid character '*' in FASTA sequence data, ignored
 +
En sortie:
 +
<pre>
 +
pep_OG_1471.fas      peptides
 +
ali_pep_OG_1471.fas  alignement des peptides
 +
ali_dna_OG_1471.fas  alignement des nucléotides
 +
</pre>
 +
Le programme suivant réalise le traitement pour un sous ensemble de ''sample'' fichiers tirés au hasard.
 +
 +
Il soumet sur le cluster donc, merci de vous remettre sur ''genologin'' pour lancer la commande ci-dessous.
 +
 +
Pour monitorer vos jobs : squeue -u $USER
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
~/work/scripts/aa_to_dna_aln_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100
 +
</pre>
 +
Comptez les nombre d'alignements obtenus.
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.2:
 +
Quelle est la raison de passer par un alignement des peptides pour obtenir l'alignement en nucléotides?
 +
Pour quelle raison allons nous travailler sur un sous ensemble d'alignements?
 +
Est-il pertinent de réaliser un échantillonnage des alignements par tirage aléatoire? Quelles autres approches pouvez-vous envisager?
 +
Serait-il pertinent de réaliser plusieurs tirages? Quel usage pourriez-vous en faire?
 +
</pre>
 +
 +
==Construction des arbres des groupes de gènes orthologues==
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.3:
 +
Quels sont les traitements réalisés par le logiciel trimal?
 +
Quelles sont les options possibles pour réaliser ces traitements ? Ne pas lister toutes les options, les regrouper par type de traitement!
 +
Connaissez-vous d'autres logiciels réalisant un traitement comparable des alignements multiples?
 +
</pre>
 +
===Evaluation des alignements avec t_coffee===
 +
<syntaxhighlight lang="bash">
 +
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_*.fas
 +
do
 +
  ip=$(basename $i .aln)   
 +
  outfile="~/work/ProchlorococcusSynechococcus/OG/alignment/"$ip".ali"
 +
  echo "module load  bioinfo/T-COFFEE_11.00.8cbe486; t_coffee -infile $i -output score_ascii, aln, score_html -outfile $outfile;"
 +
done > OG_t_coffee_TCS.sh
 +
</syntaxhighlight>
 +
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 OG_t_coffee_TCS.sh
 +
squeue -l -u $USER
 +
</pre>
 +
 +
Sélection alignements avec un SCORE=1000.
 +
<syntaxhighlight lang="bash">
 +
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/*.ali;
 +
do
 +
ip=$(basename $i .aln)   
 +
grep -H 'SCORE=1000' $i
 +
done
 +
</syntaxhighlight>
 +
 +
===Edition des alignements avec trimal===
 +
MSK
 +
<!--
 +
Trimer les alignements avec trimal. Nous pouvons utiliser les résultats de trimal pour identifier des alignements de mauvaise qualité.
 +
 +
genologin softwares : [http://bioinfo.genotoul.fr/index.php/resources-2/softwares/?searchll=trimal trimal]
 +
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
srun --pty bash
 +
module load bioinfo/trimal-1.4.1
 +
n=1
 +
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_*.fas
 +
do
 +
echo -e "Nb:" $n
 +
echo -e "Alignement:" $i "\n"
 +
dir=$(dirname $i)
 +
prefix=$(basename $i .fas)
 +
trimal -in $dir/$prefix.fas -out $dir/$prefix.trim.aln -sident -sgt
 +
((n++))
 +
done > ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.txt
 +
</pre>
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.4:
 +
Editez le fichier de résultats et analysez son contenu.
 +
</pre>
 +
=====Pour chaque alignement, extraction des % de sites sans indels et de la conservation moyenne des alignements=====
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
awk 'BEGIN {
 +
        i=1;
 +
        col[""]=0;
 +
}
 +
{
 +
        if ( $5==0)
 +
        {
 +
        col[i "," 1]=$2
 +
        }
 +
        else
 +
        if ( $2=="AverageIdentity" && $3!="Average")
 +
        {
 +
        col[i "," 2]=$3
 +
        i = i+1;
 +
        }
 +
}
 +
END {
 +
        max = i
 +
        i=1;
 +
        while (i < max) {
 +
        print i "\t" col[i "," 1] "\t" col[i "," 2];
 +
        i = i+1;
 +
        }
 +
}' < ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.txt > ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.tab
 +
</pre>
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.5:
 +
Vérifiez que vous avez extrait les informations attendues.
 +
</pre>
 +
 +
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
awk 'BEGIN { min = max = "NaN" }
 +
{
 +
    min = (NR==1 || $2<min ? $2 : min)
 +
    max = (NR==1 || $2>max ? $2: max)
 +
}
 +
END { print min, max }' < ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.tab
 +
</pre>
 +
=====Tracé de la distribution des deux paramètres=====
 +
Tracé de la distribution des deux paramètres pour définir des seuils permettant de retenir les meilleurs alignements
 +
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
module load system/R-3.6.1
 +
R
 +
pdf_file <- '~/work/ProchlorococcusSynechococcus/images/alignment_statistics.pdf'
 +
data <- read.table('~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.tab', row.names=1)
 +
pdf(file=pdf_file, paper="a4r")
 +
plot(data, xlim=c(1,100), ylim=c(0,1), pch=20, xlab="% without gap", ylab="AverageIdentity")
 +
text(data, labels=row.names(data), pos=4)
 +
dev.off()
 +
cat(pdf_file, "\n")
 +
quit()
 +
</pre>
 +
 +
Vous pouvez utiliser filezilla (sous Windows) ou un scp (sous linux) pour télécharger le fichier pdf en local.
 +
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.6:
 +
Éditez les alignements de plus mauvaises qualités.
 +
Commentez.
 +
</pre>
 +
 +
=====Sélection des meilleurs alignements=====
 +
Sélection des meilleurs alignements à partir de la distribution du nombre de sites sans indels et de la conservation moyenne des alignements.
 +
 +
Choisir --gap dans {0,100} et --identity dans {0, 1}:
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
~/work/scripts/trimal_selection.pl --alignment ~/work/ProchlorococcusSynechococcus/OG/alignment --gap {0,100} --identity {0, 1}
 +
</pre>
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.7:
 +
Combien d'alignements avez-vous retenu?
 +
</pre>
 +
-->
 +
 +
==Super-alignement (Jeudi)==
 +
Comme dans l’article publié en 2018 de Yan et al., nous utiliserons 31 gènes du core genome tirés au hasard.
 +
Nous partirons des arbres effectués sur le super-alignement protéique de ces gènes et sur celui retranscrits en nucléotides.
 +
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.8:
 +
Proposez une méthode pour obtenir le super-alignement protéique de ces 31 gènes concaténés.
 +
</pre>
 +
<!--
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.9 :
 +
De la même façon comment obtenons nous la « rétro-traduction » en codons ? Dans quels cas cela peut-il être intéressant d’effectuer l’alignement en nucléotides de cette manière ?
 +
</pre>
 +
-->
 +
====Concaténation de 31 alignements====
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
/home/formation/work/scripts/concat_aligments.pl  --alignments /home/formation/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70.lst  --outfile ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 -nb_ali 31
 +
</pre>
 +
Liste des alignements retenus:
 +
/home/formation/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31.lst
 +
Aller regarder le fichier de sortie. Est-il conforme à l'attendu ?
 +
 +
====IQ-TREE====
 +
genologin softwares : [http://bioinfo.genotoul.fr/index.php/resources-2/softwares/?searchll=iq-tree iq-tree]
 +
 +
FAQ : http://www.iqtree.org/doc/Frequently-Asked-Questions
 +
 +
Documentation : http://www.iqtree.org/doc/
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
mkdir ~/work/ProchlorococcusSynechococcus/phyloG/
 +
cd ~/work/ProchlorococcusSynechococcus/phyloG/
 +
cp ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 .
 +
</pre>
 +
Nous allons inférer un arbre à partir du super-alignement en codons:
 +
Créer le fichier condTree.sh contenant les lignes suivantes. Lancez le, puis répondez aux questions suivantes car cela va être assez long, vous y reviendrez ensuite. Surtout faites un sbatch --cpus-per-task=4 à partir de genologin pour le lancer car nous avons demandé 4 slots de calcul.
 +
N'oubliez pas de spécifier le modèle dans la ligne de commande vous gagnerez 8 heures de calcul ;-)
 +
 +
Bonne pratique pour quand vous travaillerez : prenez l'habitude de vérifier la version la plus récente installée sur genologin. Pensez à vérifier si elle doit être mise à jour, si oui, demander la mise à jour. Si non, utilisez la dernière version.
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
#!/bin/bash
 +
module load bioinfo/iqtree-2.0.6
 +
iqtree -s ~/work/ProchlorococcusSynechococcus/phyloG/alignment_0.4_70_31 -redo -bb 1000 -alrt 1000 -st CODON -nt 4 -m KOSI07+F+R6
 +
</pre>
 +
<!--
 +
[http://www.iqtree.org/doc/Substitution-Models Substitution-Models]
 +
<pre>
 +
Akaike Information Criterion:          KOSI07+F+R7
 +
Corrected Akaike Information Criterion: KOSI07+F+R7
 +
Bayesian Information Criterion:        KOSI07+F+R6
 +
Best-fit model: KOSI07+F+R6 chosen according to BIC
 +
</pre>
 +
Empirical codon model ([https://academic.oup.com/mbe/article/24/7/1464/986344 Kosiol et al., 2007]).
 +
-->
 +
Pour vous, nous avons inféré un arbre à partir du super-alignement protéique, la ligne de commande lancée était :
 +
iqtree -s all_pep.fas  -bb 1000 -alrt 1000 -nt 4
 +
 +
Vous trouverez les fichiers générés ici : /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/all_pep.fas.*. Vous pouvez les copier dans votre répertoire ~/work/ProchlorococcusSynechococcus/phyloG/.
 +
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.10:
 +
Pouvez-vous m’expliquer ce qu’on a demandé à IQTREE dans les deux cas ? Quels ont été les modèles choisis ? Pouvez-vous les expliquer ?
 +
Ceci peut vous aider :
 +
http://www.iqtree.org/doc/Substitution-Models
 +
</pre>
 +
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.11:
 +
Comparez les deux arbres (*.treefile) visuellement avec figtree par exemple (ou avec des outils en ligne tels que : phylo.io, PhyD3, iTOL...).
 +
Pour afficher les valeurs de bootstraps, il est nécessaire que figtree les charge comme « labels », ainsi dans la rubrique « node labels » il vous faudra sélectionner « labels » dans le menu déroulant « display ».
 +
Mettez l'arbre généré dans le rapport de TP.
 +
Quelles sont les principales différences ? Que pensez-vous de ces arbres ? Leurs supports ? Leurs congruences ? Les grands clades connus sont-ils retrouvés ? Comparer avec les arbres des articles ou revues suivant(e)s :
 +
La revue de Biller et al : Prochlorococcus  : the structure and function of collective diversity
 +
L’article de Kettler et al. : Patterns and implications of gene gain and loss in the evolution of Prochlorococcus et celui de Yan et al. : Genome rearrangement shapes Prochlorococcus ecological adaptation.
 +
N’oubliez pas de commenter les valeurs de support des nœuds. Pensez à enraciner les arbres en utilisant l’outgroup Synechococcus.
 +
</pre>
 +
 +
NB: Les noms des espèces et les informations des clades correspondants sont disponibles dans le fichier /home/formation/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt ou dans le tableau [http://www.m2p-bioinfo.ups-tlse.fr/p/Atelier_Phylog%C3%A9nomique#Disponibilit.C3.A9_des_g.C3.A9nomes ici].
 +
 +
==Super-arbres==
 +
 +
Nous allons utiliser les arbres individuels protéiques ainsi que celui reconstruit à partir de la petite sous-unité de l’rRNA.
 +
Sur genologin le script ssuTree.sh était celui-là :
 +
#!/bin/bash
 +
iqtree -s ./ssu_renamed_simplified.aln -nt 1 -AIC -bb 1000 -alrt 1000 -redo
 +
 +
Et on a tapé :
 +
module load bioinfo/iqtree-2.0.6
 +
sbatch --mem=20G ssuTree.sh
 +
 +
Voici l’arbre obtenu : /home/formation/work/ProchlorococcusSynechococcus/phyloG/ssu_renamed_simplified.aln.treefile
 +
 +
'''A vous :'''
 +
 +
Loguez-vous sur genologin.
 +
 +
Lançons les arbres protéiques. Pour cela faire un seul script que vous appellerez ind_pep_trees.sh et qui écrira toutes les commandes en bouclant sur chaque fichier d’entrée. Le but étant d’obtenir une ligne par commande iqtree sur un fichier d’alignement protéique. 
 +
Les fichiers d’input sont /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/*_renamed.fas.
 +
 +
Il vous faut les copier dans votre work car sinon iqtree écrira les fichiers d'output sur le répertoire du compte de formation et pas dans le votre...
 +
 +
Pour vous aider inspirez vous de la réponse à la question "How to generate an sarray command file with bash for single (fastq) file ?" sur la page http://bioinfo.genotoul.fr/index.php/faq/bioinfo_tips_faq/.
 +
 +
'''Attention à bien spécifier -nt 1. Si vous souhaitez utiliser plus d’un CPU il faut le réserver avec l’option --cpus-per-task dans la commande sbatch / sarray. Pas besoin d’augmenter la RAM pour les protéines. Pensez aussi à mettre -AIC comme critère de sélection de modèle.'''
 +
 +
Regarder le contenu de ind_pep_trees.sh. Est-il correct ?
 +
 +
Rajouter la première ligne obligatoire sur genologin avec vi par exemple :
 +
#!/bin/bash
 +
 +
Pour le lancer en parallèle sur plusieurs nœuds du cluster :
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
module load bioinfo/iqtree-2.0.6
 +
sarray ind_pep_trees.sh
 +
</pre>
 +
 +
Pour monitorer votre job :
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
squeue -l -u <login>
 +
</pre>
 +
 +
Pour le killer si besoin :
 +
scancel jobid
 +
 +
Quand tous vos jobs sont terminés, vérifiez que les fichiers de sorties ne soient pas vides et que ça s’est bien passé en faisant par exemple :
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
tail *_renamed.fas.log
 +
</pre>
 +
 +
Vous pouvez aussi regarder les modèles sélectionnés :
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
grep 'Best-fit model:' *_renamed.fas.log
 +
</pre>
 +
 +
Concaténer tous les arbres (les 31 arbres protéiques et l’arbre ARNr obtenu hier) avec la commande cat.
 +
Nommez le fichier alltrees.tree.
 +
 +
'''Test de plusieurs méthodes de super-arbres :'''
 +
 +
Commençons par la méthode la plus répandue : le MRP.
 +
 +
Pour aller sur un nœud :
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
srun --pty bash
 +
</pre>
 +
puis taper :
 +
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
module load system/R-3.5.1
 +
R
 +
library(phytools)
 +
trees=read.tree("./alltrees.tree")
 +
supertrees<-mrp.supertree(trees,rearrangements="SPR", start="NJ")
 +
</pre>
 +
Vous avez obtenu les super-arbres les plus parcimonieux. Sauvez-les en utilisant la fonction write.tree de R.
 +
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
write.tree(supertrees, file = "./superTrees.tree")
 +
quit()
 +
</pre>
 +
Nous sommes sortis de R.
 +
 +
Dans notre cas, nous avons un seul arbre le plus parcimonieux, mais nous aurions pu en obtenir plusieurs.
 +
 +
Utiliser iqtree pour obtenir le consensus des 32 arbres avec la règle majoritaire étendue.
 +
 +
Pour cela restez sur le nœud et tapez :
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
module load bioinfo/iqtree-2.0.6
 +
iqtree -con -t alltrees.tree -nt 1
 +
</pre>
 +
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.12:
 +
D’après vous que signifie les valeurs aux nœuds sur cet arbre consensus ?
 +
</pre>
 +
 +
Lancez aussi l’arbre consensus en réseau :
 +
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
iqtree -net -t alltrees.tree -nt 1
 +
</pre>
 +
 +
Et visualisez-le avec splitstree en local.
 +
 +
https://software-ab.informatik.uni-tuebingen.de/download/splitstree5/welcome.html
 +
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.13:
 +
Commentez ce réseau par rapport aux autres arbres obtenus. Qu’en pensez-vous ?
 +
</pre>
 +
 +
==Comparaison des arbres==
 +
 +
Concaténez maintenant les deux arbres de super-matrice (celui sur les codons que vous avez lancé toute à l'heure et celui sur les protéines) ainsi que les deux super-arbres (le consensus et le MRP).
 +
 +
Si vous avez eu des difficultés à obtenir l'arbre à partir de l'alignement en codon il est disponible ici : /home/formation/work/ProchlorococcusSynechococcus/phyloG/alignment_0.4_70_31.fas*.
 +
 +
Attention marquez quelques part l'ordre avec lequel vous les avez concaténés pour créer le fichier d'arbres à comparer afin de vous en souvenir ensuite.
 +
Lancer ensuite le calcul de la distance de Robinson and Foulds sur ce fichier avec iqtree entre les 4 arbres 2 à 2.
 +
 +
Attention à faire le module load bioinfo/iqtree-2.0.6
 +
 +
Et à rajouter -nt 1 dans les options.
 +
Restez sur le nœud ou faites un sbatch, mais ne lancez rien sur le nœud maître.
 +
 +
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 +
Question 5.14:
 +
Commentez les résultats. Quel est l’arbre le plus différent des autres ? Qu’est-ce qui pourrait l’expliquer ?
 +
</pre>
 +
----
----
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces Arbre espèces]
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces Arbre espèces]

Revision as of 15:46, 15 October 2021

Contents

Liens

Introduction

Support de cours : supports

Nous allons utiliser un sous ensemble de gènes concervés chez Prochlorococcus et Synechococcus pour expérimenter les différentes méthodes de reconstruction phylogénomiques. Nous nous initierons à la comparaison d’arbres.

Extraction des séquences nucléotidiques des gènes orthologues

Créer un fichier avec toutes les séquences nucléotidiques (cds):

Si le répertoire n'existe pas déjà, le créer.

mkdir -p ~/work/ProchlorococcusSynechococcus/DNA
cat ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.ffn ~/work/Synechococcus/prokka/Aaa*/Aaa*.ffn > ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn

grep -c '>' ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn

Extraire les groupes de gènes orthologues

Exemple de la création d'un script et lancement du job avec sbatch

mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA

Créer le fichier suivant commande.sh.

#!/usr/bin/bash
~/work/scripts/extract_sequences_from_matchtable.pl --fasta ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn --matchtable ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt --outdir ~/work/ProchlorococcusSynechococcus/OG/DNA --quorum 16
sbatch commande.sh

Vérifier que les fichiers contiennent au moins 16 séquences:

grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
Question 5.1:
A quoi correspond ce quorum. Pourquoi utiliser un seuil de 16? 
Combien de fichiers avez-vous obtenus?
Est-ce pertinent dans les situations où vous avez un grand nombre de souches?
Et quand les espèces étudiées sont très éloignées phylogénétiquement parlant ?

Transformation en alignement multiples des séquences ADN des gènes

Le script aa_to_dna_aln.pl prend en entrée un fichier fasta contenant des séquences d'ADN. Il traduit des séquences, réalise un alignement multiple avec muscle et retraduit en ADN cet alignement multiple (BioPerl).

mkdir ~/work/ProchlorococcusSynechococcus/OG/alignment

Attention à faire un script et le lancer en sbatch ou en srun --pty bash

 ~/work/scripts/aa_to_dna_aln.pl -dna ~/work/ProchlorococcusSynechococcus/OG/DNA/OG_1471.fas --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment

La traduction en peptides ajoute un '*' à la fin des séquences. Vous pouvez ignorer le message:

*** WARNING *** Invalid character '*' in FASTA sequence data, ignored

En sortie:

pep_OG_1471.fas       peptides
ali_pep_OG_1471.fas   alignement des peptides
ali_dna_OG_1471.fas   alignement des nucléotides

Le programme suivant réalise le traitement pour un sous ensemble de sample fichiers tirés au hasard.

Il soumet sur le cluster donc, merci de vous remettre sur genologin pour lancer la commande ci-dessous.

Pour monitorer vos jobs : squeue -u $USER

~/work/scripts/aa_to_dna_aln_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100

Comptez les nombre d'alignements obtenus.

Question 5.2:
Quelle est la raison de passer par un alignement des peptides pour obtenir l'alignement en nucléotides?
Pour quelle raison allons nous travailler sur un sous ensemble d'alignements? 
Est-il pertinent de réaliser un échantillonnage des alignements par tirage aléatoire? Quelles autres approches pouvez-vous envisager?
Serait-il pertinent de réaliser plusieurs tirages? Quel usage pourriez-vous en faire?

Construction des arbres des groupes de gènes orthologues

Question 5.3:
Quels sont les traitements réalisés par le logiciel trimal? 
Quelles sont les options possibles pour réaliser ces traitements ? Ne pas lister toutes les options, les regrouper par type de traitement!
Connaissez-vous d'autres logiciels réalisant un traitement comparable des alignements multiples?

Evaluation des alignements avec t_coffee

for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_*.fas 
do 
  ip=$(basename $i .aln)     
  outfile="~/work/ProchlorococcusSynechococcus/OG/alignment/"$ip".ali"
  echo "module load  bioinfo/T-COFFEE_11.00.8cbe486; t_coffee -infile $i -output score_ascii, aln, score_html -outfile $outfile;" 
done > OG_t_coffee_TCS.sh
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 OG_t_coffee_TCS.sh
squeue -l -u $USER

Sélection alignements avec un SCORE=1000.

for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/*.ali; 
do 
 ip=$(basename $i .aln)     
 grep -H 'SCORE=1000' $i
done

Edition des alignements avec trimal

MSK

Super-alignement (Jeudi)

Comme dans l’article publié en 2018 de Yan et al., nous utiliserons 31 gènes du core genome tirés au hasard. Nous partirons des arbres effectués sur le super-alignement protéique de ces gènes et sur celui retranscrits en nucléotides.

Question 5.8: 
Proposez une méthode pour obtenir le super-alignement protéique de ces 31 gènes concaténés.

Concaténation de 31 alignements

/home/formation/work/scripts/concat_aligments.pl  --alignments /home/formation/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70.lst  --outfile ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 -nb_ali 31

Liste des alignements retenus:

/home/formation/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31.lst

Aller regarder le fichier de sortie. Est-il conforme à l'attendu ?

IQ-TREE

genologin softwares : iq-tree

FAQ : http://www.iqtree.org/doc/Frequently-Asked-Questions

Documentation : http://www.iqtree.org/doc/

mkdir ~/work/ProchlorococcusSynechococcus/phyloG/
cd ~/work/ProchlorococcusSynechococcus/phyloG/
cp ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 .

Nous allons inférer un arbre à partir du super-alignement en codons: Créer le fichier condTree.sh contenant les lignes suivantes. Lancez le, puis répondez aux questions suivantes car cela va être assez long, vous y reviendrez ensuite. Surtout faites un sbatch --cpus-per-task=4 à partir de genologin pour le lancer car nous avons demandé 4 slots de calcul. N'oubliez pas de spécifier le modèle dans la ligne de commande vous gagnerez 8 heures de calcul ;-)

Bonne pratique pour quand vous travaillerez : prenez l'habitude de vérifier la version la plus récente installée sur genologin. Pensez à vérifier si elle doit être mise à jour, si oui, demander la mise à jour. Si non, utilisez la dernière version.

#!/bin/bash
module load bioinfo/iqtree-2.0.6
iqtree -s ~/work/ProchlorococcusSynechococcus/phyloG/alignment_0.4_70_31 -redo -bb 1000 -alrt 1000 -st CODON -nt 4 -m KOSI07+F+R6

Pour vous, nous avons inféré un arbre à partir du super-alignement protéique, la ligne de commande lancée était :

iqtree -s all_pep.fas  -bb 1000 -alrt 1000 -nt 4

Vous trouverez les fichiers générés ici : /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/all_pep.fas.*. Vous pouvez les copier dans votre répertoire ~/work/ProchlorococcusSynechococcus/phyloG/.

Question 5.10: 
Pouvez-vous m’expliquer ce qu’on a demandé à IQTREE dans les deux cas ? Quels ont été les modèles choisis ? Pouvez-vous les expliquer ?
Ceci peut vous aider : 
http://www.iqtree.org/doc/Substitution-Models
Question 5.11:
Comparez les deux arbres (*.treefile) visuellement avec figtree par exemple (ou avec des outils en ligne tels que : phylo.io, PhyD3, iTOL...).
Pour afficher les valeurs de bootstraps, il est nécessaire que figtree les charge comme « labels », ainsi dans la rubrique « node labels » il vous faudra sélectionner « labels » dans le menu déroulant « display ».
Mettez l'arbre généré dans le rapport de TP.
Quelles sont les principales différences ? Que pensez-vous de ces arbres ? Leurs supports ? Leurs congruences ? Les grands clades connus sont-ils retrouvés ? Comparer avec les arbres des articles ou revues suivant(e)s :
La revue de Biller et al : Prochlorococcus  : the structure and function of collective diversity
L’article de Kettler et al. : Patterns and implications of gene gain and loss in the evolution of Prochlorococcus et celui de Yan et al. : Genome rearrangement shapes Prochlorococcus ecological adaptation.
N’oubliez pas de commenter les valeurs de support des nœuds. Pensez à enraciner les arbres en utilisant l’outgroup Synechococcus.

NB: Les noms des espèces et les informations des clades correspondants sont disponibles dans le fichier /home/formation/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt ou dans le tableau ici.

Super-arbres

Nous allons utiliser les arbres individuels protéiques ainsi que celui reconstruit à partir de la petite sous-unité de l’rRNA. Sur genologin le script ssuTree.sh était celui-là :

#!/bin/bash
iqtree -s ./ssu_renamed_simplified.aln -nt 1 -AIC -bb 1000 -alrt 1000 -redo

Et on a tapé :

module load bioinfo/iqtree-2.0.6
sbatch --mem=20G ssuTree.sh

Voici l’arbre obtenu : /home/formation/work/ProchlorococcusSynechococcus/phyloG/ssu_renamed_simplified.aln.treefile

A vous :

Loguez-vous sur genologin.

Lançons les arbres protéiques. Pour cela faire un seul script que vous appellerez ind_pep_trees.sh et qui écrira toutes les commandes en bouclant sur chaque fichier d’entrée. Le but étant d’obtenir une ligne par commande iqtree sur un fichier d’alignement protéique. Les fichiers d’input sont /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/*_renamed.fas.

Il vous faut les copier dans votre work car sinon iqtree écrira les fichiers d'output sur le répertoire du compte de formation et pas dans le votre...

Pour vous aider inspirez vous de la réponse à la question "How to generate an sarray command file with bash for single (fastq) file ?" sur la page http://bioinfo.genotoul.fr/index.php/faq/bioinfo_tips_faq/.

Attention à bien spécifier -nt 1. Si vous souhaitez utiliser plus d’un CPU il faut le réserver avec l’option --cpus-per-task dans la commande sbatch / sarray. Pas besoin d’augmenter la RAM pour les protéines. Pensez aussi à mettre -AIC comme critère de sélection de modèle.

Regarder le contenu de ind_pep_trees.sh. Est-il correct ?

Rajouter la première ligne obligatoire sur genologin avec vi par exemple :

#!/bin/bash 

Pour le lancer en parallèle sur plusieurs nœuds du cluster :

module load bioinfo/iqtree-2.0.6
sarray ind_pep_trees.sh

Pour monitorer votre job :

squeue -l -u <login>

Pour le killer si besoin :

scancel jobid

Quand tous vos jobs sont terminés, vérifiez que les fichiers de sorties ne soient pas vides et que ça s’est bien passé en faisant par exemple :

tail *_renamed.fas.log

Vous pouvez aussi regarder les modèles sélectionnés :

grep 'Best-fit model:' *_renamed.fas.log

Concaténer tous les arbres (les 31 arbres protéiques et l’arbre ARNr obtenu hier) avec la commande cat. Nommez le fichier alltrees.tree.

Test de plusieurs méthodes de super-arbres :

Commençons par la méthode la plus répandue : le MRP.

Pour aller sur un nœud :

srun --pty bash

puis taper :

module load system/R-3.5.1
R
library(phytools)
trees=read.tree("./alltrees.tree")
supertrees<-mrp.supertree(trees,rearrangements="SPR", start="NJ")

Vous avez obtenu les super-arbres les plus parcimonieux. Sauvez-les en utilisant la fonction write.tree de R.

write.tree(supertrees, file = "./superTrees.tree")
quit()

Nous sommes sortis de R.

Dans notre cas, nous avons un seul arbre le plus parcimonieux, mais nous aurions pu en obtenir plusieurs.

Utiliser iqtree pour obtenir le consensus des 32 arbres avec la règle majoritaire étendue.

Pour cela restez sur le nœud et tapez :

module load bioinfo/iqtree-2.0.6
iqtree -con -t alltrees.tree -nt 1
Question 5.12:
D’après vous que signifie les valeurs aux nœuds sur cet arbre consensus ?

Lancez aussi l’arbre consensus en réseau :

iqtree -net -t alltrees.tree -nt 1

Et visualisez-le avec splitstree en local.

https://software-ab.informatik.uni-tuebingen.de/download/splitstree5/welcome.html

Question 5.13:
Commentez ce réseau par rapport aux autres arbres obtenus. Qu’en pensez-vous ?

Comparaison des arbres

Concaténez maintenant les deux arbres de super-matrice (celui sur les codons que vous avez lancé toute à l'heure et celui sur les protéines) ainsi que les deux super-arbres (le consensus et le MRP).

Si vous avez eu des difficultés à obtenir l'arbre à partir de l'alignement en codon il est disponible ici : /home/formation/work/ProchlorococcusSynechococcus/phyloG/alignment_0.4_70_31.fas*.

Attention marquez quelques part l'ordre avec lequel vous les avez concaténés pour créer le fichier d'arbres à comparer afin de vous en souvenir ensuite. Lancer ensuite le calcul de la distance de Robinson and Foulds sur ce fichier avec iqtree entre les 4 arbres 2 à 2.

Attention à faire le module load bioinfo/iqtree-2.0.6

Et à rajouter -nt 1 dans les options. Restez sur le nœud ou faites un sbatch, mais ne lancez rien sur le nœud maître.

Question 5.14: 
Commentez les résultats. Quel est l’arbre le plus différent des autres ? Qu’est-ce qui pourrait l’expliquer ?