Atelier Phylogénomique
From silico.biotoul.fr
Matériel pédagogique
Support de cours : Quest for Orthologs
Références
- Kettler et al., PLoS Genet. 2007 Dec;3(12):e231 Patterns and implications of gene gain and loss in the evolution of Prochlorococcus.
- Sun and Blanchard, 2014 Strong Genome-Wide Selection Early in the Evolution of Prochlorococcus Resulted in a Reduced Genome through the Loss of a Large Number of Small Effect Genes
- Yan et al., Appl Environ Microbiol. 2018 Genome rearrangement shapes Prochlorococcus ecological adaptation.
- Biller et al., Nat. Rev. Microbiol. 2015 13(1) 13-27 Prochlorococcus: the structure and function of collective diversity.
- Partensky and Laurence Garczarek Annual Review of Marine Science 2010 Prochlorococcus: Advantages and Limits of Minimalism.
- Prochlorococcus Prochlorococcus.
Logiciels à installer sur vos postes de travail
- seaview : Multiplatform GUI for molecular phylogeny
- mauve : Multiple genome alignment
- Artemis : Genome browser and annotation tool
- Artemis Comparison Tool : Java application for displaying pairwise comparisons between two or more DNA sequences
- splitstree The aim of SplitsTree4 is to provide a framework for evolutionary analysis using both trees and networks.
- FigTree is designed as a graphical viewer of phylogenetic trees and as a program for producing publication-ready figures.
Artemis (art et act) et figtree peuvent être installés avec conda.
Mauve, mauveAligner et progressiveMauve peuvent être installés avec conda mais une erreur peut survenir avec Mauve. En effet, vous pouvez avoir une version trop récente de java (https://edwards.sdsu.edu/research/running-mauve-with-java-10/). Pour y remédier, chercher une "vieille" version de java et remplacer "java" par le chemin de cette version dans le fichier Mauve (ligne JAVA_CMD=java).
Ressources informatiques
Nous allons utiliser les ressources de GenoToul.
Vous avez dans les FAQ, les réponses à toutes vos questions concernant l'utilisation de la ressource.
Sortcuts
Vous allez vous connecter avec un compte fleur:
ssh -Y <login>@genologin.toulouse.inra.fr
Vous pouvez ouvrir deux connections afin de pouvoir lancer qlogin de façon indépendante.
Échange de fichiers:
scp file <login>@genologin.toulouse.inra.fr:~/work
Logiciels disponibles
ou
ou plus directement
ls /usr/local/bioinfo/src/
La documentation est disponible sur le site WEB et dans le répertoire correspondant au logiciel (fichiers How_to_use and/or Readme).
Soumission de jobs avec SLURM
Lancer le job avec sbatch et un script du type "myscript.sh" (les chemins sont à changer).
to submit the job, use the sbatch command line as following
- sbatch (qsub): submit a batch job to slurm (default workq partition()
sbatch myscript.sh
- sarray (qarray): submit a batch job-array to slurm
- scancel (qdel): kill the specified job
INTERACTIVE
- srun [--pty bash] (qrsh): submit an interactive session with a compute node (default workq partition).
- runVisuSession.sh (qlogin): submit a TurboVNC / VirtualGL session with the graphical node (interq partition). Just for graphics jobs.
Pour controler les jobs
- sinfo (qhost): display nodes, partitions, reservations
- squeue (qstat): display jobs and state
- scontrol show : get informations on jobs, nodes, partitions
- sstat (qstat -j): show status of running jobs
- sview (qmon): graphical user interface
- sacct (qacct) : display accounting data
Utiliser R sur le cluster
Tutoriel expliquant comment utiliser R (et compiler des fichiers Rmd) sur le cluster (slurm) de la PF Bioinformatique de Toulouse (Gaëlle Lefort et Nathalie Vialaneix):
Scripts
Aide pour les scripts en bash :https://www.tutorialkart.com/bash-shell-scripting
Des scripts perl ont été écrits pour faciliter certaines étapes du TP. Ils sont disponibles dans le répertoire:
/home/formation/public_html/M2_Phylogenomique/scripts
Vous pouvez les copier dans votre espace de travail et les modifier à votre convenance.
mkdir ~/work/scripts cp /home/formation/public_html/M2_Phylogenomique/scripts/* ~/work/scripts
Copie des fichiers
Rsync (Remote Sync) est une commande couramment utilisée pour copier et synchroniser des fichiers et des répertoires à distance ainsi que localement dans les systèmes Linux/Unix. Avec rsync, vous pouvez copier et synchroniser vos données à distance et localement à travers des répertoires, des disques et des réseaux, effectuer des sauvegardes de données et des mises en miroir entre deux machines Linux.
Nous allons utiliser cette commande pour copier les fichiers de <user>@genologin.toulouse.inra.fr à votre machine.
Exemple:
- Créez un répertoire sur votre machine correspondant à l'atelier.
- Placez vous dans ce répertoire.
rsync --archive --itemize-changes --stats -h -e ssh <user>@genologin.toulouse.inra.fr:/home/<user>/work/Prochlorococcus work
Pour information:
ls -l /home/<user> save -> /save/<user> work -> /work/<user>
Disponibilité des génomes
Caractéristiques des souches étudiées
Table modifiée à partir de la Table 1 de (Kettler et al., 2007).
Cyanobacterium Isolate Light Length(bp) GC% Number Isol. Region Date Accession Code Adap. Genes Depth Prochlorococcus MED4 HL(I) 1,657,990 30.8 1,929 5m Med. Sea Jan. 1989 BX548174 Aaab MIT9515 HL(I) 1,704,176 30.8 1,908 15m Eq. Pacific Jun. 1995 CP000552 Aaag MIT9301 HL(II) 1,642,773 31.4 1,907 90m Sargasso Sea Jul. 1993 CP000576 Aaaj AS9601 HL(II) 1,669,886 31.3 1,926 50m Arabian Sea Nov. 1995 CP000551 Aaaf MIT9215 HL(II) 1,738,790 31.1 1,989 5m Eq. Pacific Oct. 1992 CP000825 Aaak MIT9312 HL(II) 1,709,204 31.2 1,962 135m Gulf Stream Jul. 1993 CP000111 Aaae NATL1A LL(I) 1,864,731 35.1 2,201 30m N. Atlantic Apr. 1990 CP000553 Aaai NATL2A LL(I) 1,842,899 35 2,158 10m N. Atlantic Apr. 1990 CP000095 Aaad CCMP1375/SS120 LL(II) 1,751,080 36.4 1,925 120m Sargasso Sea May 1988 AE017126 Aaaa MIT9211 LL(III)1,688,963 38 1,855 83m Eq. Pacific Apr. 1992 CP000878 Aaah MIT9303 LL(IV) 2,682,807 50.1 3,022 100m Sargasso Sea Jul. 1992 CP000554 Aaal MIT9313 LL(IV) 2,410,873 50.7 2,843 135m Gulf Stream Jul. 1992 BX548175 Aaac Synechococcus CC9311 Syn. 2,606,748 52.5 3017 95m Calif. Current 1993 CP000435 Aaao CC9902 Syn. 2,234,828 54.2 2504 5m Calif. Current 1999 CP000097 Aaam WH8102 Syn. 2,434,428 59.4 2787 Sargasso Sea Mar. 1981 BX548020 Aaap CC9605 Syn. 2,510,659 59.2 2991 51m Calif. Current 1996 CP000110 Aaan
Prochlorococcus
NCBI
Génomes de Prochlorococcus disponibles au NCBI browse
Accessions des génomes utilisés dans Kettler et al., PLoS Genet. 2007:
- accession "BX548174,CP000552,CP000576,CP000551,CP000825,CP000111,CP000553,CP000095,AE017126,CP000878,CP000554,BX548175"
Les fichiers GenBank peuvent être obtenus avec la commande wget en utilisant les chemins enregistrés dans le fichier: accession_file.lst
Pour des raisons de compatibilité avec les programmes, les génomes sont renommés en utilisant un code à quatre lettres. accession_labels_file.lst
GenBank
Les fichiers GenBank renommés sont disponibles dans le répertoire: GenBank
DNA
Les séquences ADN ont été extraites des fichiers GenBank : DNA
Annotation
Prokka
Les réplicons des génomes sont annotés avec le logiciel prokka.
Question 1.1: Pourquoi pensez-vous qu'il soit nécessaire d'annoter les génomes téléchargés du NCBI? Quelles sont les annotations réalisées par Prokka? Quels sont les logiciels utilisés pour réaliser ces annotations ?
Exemple d'utilisation
Nous allons créer un répertoire pour les résultats de prokka et chercher la dernière version disponible de prokka sur le serveur.
mkdir -p ~/work/Prochlorococcus/prokka search_module prokka
srun --pty bash module load bioinfo/prokka-1.12 prokka /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaaa.fas --outdir ~/work/Prochlorococcus/prokka/Aaaa --compliant --addgenes --prefix Aaaa --locustag Aaaa.g --genus Prochlorococcus --species 'Prochlorococcus marinus subsp. marinus' --strain CCMP1375 --kingdom Bacteria --cpus 2
Le programme génère plusieurs fichiers pour chaque réplicon (~/work/Prochlorococcus/prokka/Aaaa), dont:
- annotation en format GenBank AaaaA01.gbk
- annotation en format gff AaaaA01.gff
- annotation en format tabulé AaaaA01.tbl
- les peptides AaaaA01.faa
- les séquences des CDS AaaaA01.ffn
Automatisation des annotations prokka sur l'ensemble des génomes
Les informations sur les génomes sont disponibles dans le fichier : species_strain_names.txt. Ce fichier est lu par le script Perl prokka_loop.pl pour compléter les paramètres de prokka pour chaque génome (--prefix, --locustag, --genus, --species, --strain and --kingdom).
Le script prokka_loop.pl doit être lancé sur le serveur genologin. Il distribue prokka sur les noeuds avec sbatch.
cd ~/work/Prochlorococcus ~/work/scripts/prokka_loop.pl --sample Prochlorococcus squeue -l -u <user>
Une fois les jobs terminés, vérifiez que les fichiers de sortie de prokka existent et ne sont pas vides.
ls -l ~/work/Prochlorococcus/prokka/Aaa*/*.faa
Les fichiers avec le suffixe .err renferment la sortie standard de prokka. Si tout s'est bien passé, vous pouvez supprimer les fichiers .err et .sh.
Question 1.2: Comparez le nombre de gènes obtenus avec ceux reportés dans la publication (Table 1) et commentez les différences observées. Comment faire pour comparer les annotations de ''prokka'' avec celles des fichiers GenBank? Pensez-vous que ''prokka'' soit la meilleure méthode d'annotation? Comment pourriez-vous faire pour évaluer les performances des différentes méthodes d'annotation des génomes?
Visualisation des annotations
Nous pouvons utiliser le logiciel art (Artemis) pour visualiser les annotations des génomes:
Il est fortement recommandé d'utiliser ce logiciel en local sur votre poste de travail.
Conservation de séquence entre souches de Prochlorococcus
Genome pairs
BlastN par pairs
Afin d'estimer les conservations entre les différents génomes, nous allons les comparer par paire de génomes dans l'ordre suivant, à l'aide de blastn:
'Aaab', 'Aaag', 'Aaaj', 'Aaaf', 'Aaak', 'Aaae', 'Aaai', 'Aaad', 'Aaaa', 'Aaah', 'Aaal', 'Aaac'
Les résultats sont dans le repertoire:
mkdir ~/work/Prochlorococcus/BlastN
Nous allons utiliser l'option BLAST-2-Sequences de blastn en précisant -subject <File_In>.
Exemple avec une paire de génomes:
search_module blast srun --pty bash module load bioinfo/ncbi-blast-2.7.1+ blastn -query /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaab.fas -subject /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaag.fas -evalue 1e-05 -outfmt 6 -num_threads 1 -out ~/work/Prochlorococcus/BlastN/Aaab_vs_Aaag.tab
Nous allons exécuter la même commande sur toutes les paires consécutives de génomes dans l'ordre listés ci-dessus. Nous allons utiliser sarray pour soumettre ces commandes en même temps sur le cluster. Vous pouvez vous référer à "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/ pour vous aider dans cette tâche.
Nous allons écrire un script shell pour créer le fichier à soumettre par sarray.
- Dans un premier temps vous devez utiliser une boucle for pour construire les paires de génomes adjacents dans la liste ci-dessus.
- et pour chaque paires reproduire la commande donnée en exemple en changeant les noms des génomes.
MSK
Vérifier le script et lancer le avec sarray
cat blastn_pairs.sh sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blastn_pairs.sh
Vérifier l'exécution des blastn
squeue -l -u <user>
Lister les fichiers obtenus:
ls -l ~/work/Prochlorococcus/BlastN
genoplotR
Nous allons utiliser genoplotR pour visualiser les similarités entre les paires de génomes.
Installation du package genoPlotR
srun --pty bash module load system/R-3.5.1 R install.packages('genoPlotR') ... * installing *source* package ‘genoPlotR’ ... library(genoPlotR)
Sélectionner France (Lyon 2) [https] comme miroir CRAN.
Mise en œuvre
genoplotR nécessite plusieurs objets:
- dna_seg: un objet dna_seg est un ensemble de gènes ou d'éléments le long d'un génome, à représenter sur une carte. Nous allons utiliser les fichiers en format gbk créés par prokka.
- comparison: une comparaison est un ensemble de similitudes, représentant la comparaison entre deux segments d'ADN. Nous allons utiliser les résultats des blastn entre paires de genomes.
- annotation: un objet d'annotation est utilisé pour annoter un segment d'ADN. Nous ne l'utilisons pas ici.
- tree: un arbre au format Newick qui peut être analysé à l'aide du paquetage ade4. Nous l'utiliserons plus tard!
mkdir ~/work/Prochlorococcus/images srun --pty bash module load system/R-3.5.1 Rscript ~/work/scripts/genoplot_blastn_links.R
Pour visualiser les fichiers pdf, il est préférable d'utiliser votre machine en P0. Pensez à faire des rsync avant! Placez-vous dans le répertoire racine de votre TD (au dessus de work).
evince work/Prochlorococcus/images/genoplot_blastn_links.pdf
ACT
Il est également possible d'utiliser le logiciel act (documentation).
Copier les fichiers sur votre porte de travail en P0 et lancer :
act work/Prochlorococcus/prokka/Aaab/Aaab.gbk work/Prochlorococcus/BlastN/Aaab_vs_Aaag.tab work/Prochlorococcus/prokka/Aaag/Aaag.gbk
Vous pouvez aussi utiliser les fichiers en format gff.
act work/Prochlorococcus/prokka/Aaab/Aaab.gff work/Prochlorococcus/BlastN/Aaab_vs_Aaag.tab work/Prochlorococcus/prokka/Aaag/Aaag.gff work/Prochlorococcus/BlastN/Aaag_vs_Aaaj.tab work/Prochlorococcus/prokka/Aaaj/Aaaj.gff work/Prochlorococcus/BlastN/Aaaj_vs_Aaaf.tab work/Prochlorococcus/prokka/Aaaf/Aaaf.gff work/Prochlorococcus/BlastN/Aaaf_vs_Aaak.tab work/Prochlorococcus/prokka/Aaak/Aaak.gff
Question 1.3: Commentez les résultats obtenus avec genoplotR et ACT. Que pensez-vous de la conservation des séquences des génomes?
Groupes de gènes orthologues
Question 1.4: Quelle méthode ont utilisé Kettler et al., 2007 pour identifier les groupes de gènes orthologues?
Préliminaires
Question 1.5: Selon vous qu'est-ce qui guide le choix du type de séquences à utiliser dans les comparaisons (peptides ou nucléotidiques)?
Blast All-All
Nous allons utiliser NCBI_Blast+.
Nous allons copier les fichiers peptides dans un répertoire de travail:
mkdir -p ~/work/Prochlorococcus/peptide cp ~/work/Prochlorococcus/prokka/Aaa*/*.faa ~/work/Prochlorococcus/peptide/. ls -l ~/work/Prochlorococcus/peptide
make blast database
Exemple:
module load bioinfo/ncbi-blast-2.7.1+ srun -n1 -l makeblastdb -in Aaaa.faa -dbtype prot
Vous allez procéder comme précédemment, avec un script donné à sarray, pour réaliser le makeblastdb sur tous les fichiers.
MSK
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 makeblastdb.sh squeue -l -u <user>
Paire de genomes
mkdir -p ~/work/Prochlorococcus/BlastP
Exemple :
module load bioinfo/ncbi-blast-2.7.1+ srun -n1 -l blastp -query ~/work/Prochlorococcus/peptide/Aaaa.faa -db ~/work/Prochlorococcus/peptide/Aaaa.faa -seg yes -dbsize 100000000 -evalue 1e-5 -outfmt 6 -num_threads 1 -out ~/work/Prochlorococcus/BlastP/Aaaa_Aaaa.tab
Question 1.6: Explicitez et justifiez les paramètres de blast utilisés dans votre script.
Boucle sur les genomes
Utilisez sarray pour réaliser les blast en toutes les paires de génomes.
MSK
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh squeue -l -u <user>
vérifiez que les fichiers de sorties de blast sont présents et non vides.
Question 1.7: Combien de fichiers attendez-vous?
PorthoMCL
PorthoMCL: Parallel orthology prediction using MCL for the realm of massive genome availability Ehsan Tabari and Zhengchang Su Big Data Analytics 2017 2:4 DOI: 10.1186/s41044-016-0019-8
Il s'inspire d'orthoMCL
Question 1.8: Rappelez les différentes étapes réalisées par le logiciel. Précisez les paramètres utilisés pour identifier les paires de gènes orthologues.
Le logiciel est disponible ici: /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master. Il est basé sur l'enchainement de programmes python. Ils peuvent être lancés séparément ou automatiquement à l'aide d'un script shell.
mkdir -p ~/work/Prochlorococcus/PorthoMCL_auto12/0.input_faa cd ~/work/Prochlorococcus/PorthoMCL_auto12 cp /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/prokka/A*/*.faa 0.input_faa/. srun --nodes=12 --pty bash /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomcl.sh -t 12 -e 7
Nous allons utiliser une version plus simple:
mkdir ~/work/Prochlorococcus/PorthoMCL cd ~/work/Prochlorococcus/PorthoMCL srun --pty bash ~/work/scripts/short_PorthoMCL.sh
Question 1.9: Editez le script et donnez une description des différentes tâches réalisées. Que proposeriez-vous pour l'améliorer?
Step 4: Parse BLAST results
Comme les différents blastp ont déjà été réalisés, nous passons directement à l'étape 4.
Mais nous devons:
Concaténer les fichiers de résultats blastp par souche:
mkdir -p ~/work/Prochlorococcus/PorthoMCL/3.blastres
MSK
Modifier le nom des protéines dans les fichiers blast et fasta
Pour que les programmes suivants identifient le nom de souche et le nom de gène, il est nécessaire de rajouter un séparateur ('|') dans les noms des protéines: Aaaa.g_00001 devient Aaaa|.g_00001. Pour cela, nous pouvons utiliser la commande sed.
Dans les fichiers ~/work/Prochlorococcus/PorthoMCL/3.blastres/*.tab
MSK
et dans les fichiers ~/work/Prochlorococcus/peptide/*.faa
MSK
Changer le format des fichier blastp
mkdir -p ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq
Exemple pour une souche:
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclBlastParser /home/yquentin/work/Prochlorococcus/PorthoMCL/3.blastres/Aaaa.tab ~/work/Prochlorococcus/PorthoMCL/compliantFasta > ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq/Aaaa.ss.tsv
head ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq/Aaaa.ss.tsv
Format de sortie:
query hit evalueMant evalueExp percentMatch
- evalueMant et evalueExp : mantisse (partie décimale) et exposant de la e-value
- percentMatch = longueur de l'alignement / longueur de la séquence la plus courte (query ou hit) *100
Appliquez à toutes les souches.
MSK
Vérifiez les fichiers obtenus.
Step 5: Finding Best Hits
Les paralogues sont identifiés et un score non normalisé leurs aient assignés. Les scores seront normalisés à l'étape 5.3 pour qu'ils soient comparables entre les génomes.
mkdir -p ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp mkdir -p ~/work/Prochlorococcus/PorthoMCL/5.besthit
Créer un fichier avec la liste des taxon
ls -1 ~/work/Prochlorococcus/peptide/*faa | sed -r 's/.+([A-Z][a-z]{3}).faa/\1/' > ~/work/Prochlorococcus/PorthoMCL/taxon_list
Test avec une souche:
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsBestHit.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -s ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq -b ~/work/Prochlorococcus/PorthoMCL/5.besthit -q ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp -x 1
head ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp/Aaaa.pt.tsv head ~/work/Prochlorococcus/PorthoMCL/5.besthit/Aaaa.bh.tsv
Boucle sur les souches:
MSK
Step 6: Finding Orthologs
A la sortie de cette étape nous obtenons l'ensemble des paires de gènes orthologues. Les fichiers de sortie sont concaténés. Le fichier concaténé est donnée en entrée à MCL.
mkdir -p ~/work/Prochlorococcus/PorthoMCL/6.orthologs
Test avec une souche:
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsOrthologs.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -b ~/work/Prochlorococcus/PorthoMCL/5.besthit -o ~/work/Prochlorococcus/PorthoMCL/6.orthologs -x 1
head ~/work/Prochlorococcus/PorthoMCL/6.orthologs/Aaaa.ort.tsv
Boucle sur les souches:
MSK
Step 7: Finding Paralogs
Les scores calculés à l'étape 5 sont normalisés. Pour la normalisation, PorthoMCL utilise la moyenne des scores des paralogues qui ont des orthologues.
Pour trouver les paralogues, nous devons extraire tous les gènes qui ont une relation orthologique dans chaque génome. Ceci peut être réalisé en utilisant un ensemble de commandes bash comme suit.
mkdir -p ~/work/Prochlorococcus/PorthoMCL/7.ogenes cd ~/work/Prochlorococcus/PorthoMCL
genes in the second column
awk -F'[|\t]' '{print $4 >> ("7.ogenes/"$3".og.tsv")}' 6.orthologs/*.ort.tsv
genes in the first column
awk -F'[|\t]' '{print $2 >> ("7.ogenes/"$1".og.tsv")}' 6.orthologs/*.ort.tsv
Avec ces fichiers (fichiers og) en main, nous pouvons normaliser les relations inparalogiques :
Le résultat de cette étape est l'ensemble des gènes paralogues avec des scores normalisés.
mkdir -p ~/work/Prochlorococcus/PorthoMCL/7.paralogs
Test avec la première souche:
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsInParalogs.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -q ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp -o ~/work/Prochlorococcus/PorthoMCL/7.ogenes -p ~/work/Prochlorococcus/PorthoMCL/7.paralogs -x 1
head ~/work/Prochlorococcus/PorthoMCL/7.paralogs/Aaaa.par.tsv
Boucle sur les souches:
MSK
Step 8: MCL
Question 1.10: Rappelez le principe de MCL et les paramètres utilisés. Quel est l'effet de ces paramètres?
C'est exactement comme la dernière étape d'OrthoMCL. Il suffit d'exécuter MCL sur les sorties du script :
Notez que -t définit le nombre de processeurs que vous devez utiliser pour MCL.
Location: /usr/local/bioinfo/src/MCL
srun --pty bash #Load modules module load bioinfo/mcl-14-137 cat ~/work/Prochlorococcus/PorthoMCL/6.orthologs/*.tsv >> ~/work/Prochlorococcus/PorthoMCL/8.all.ort.tsv mcl ~/work/Prochlorococcus/PorthoMCL/8.all.ort.tsv --abc -I 1.5 -t 4 -o ~/work/Prochlorococcus/PorthoMCL/8.all.ort.group
Taille des OG:
Vous pouvez utiliser awk avec NF qui est le nombre de champs de l'enregistrement courant.
MSK
Que pensez-vous de la distribution des taille des OG? Quelle-est la taille maximum attendue?
Paralogs
cat 7.paralogs/*.tsv >> 8.all.par.tsv mcl 8.all.par.tsv --abc -I 1.5 -t 4 -o 8.all.par.group
Visualisation des OG
srun --pty bash module load system/R-3.5.1
Vision globale avec genoplotR
Rscript ~/work/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_all.pdf evince work/Prochlorococcus/images/8.all.ort.group_all.pdf
Question 1.11: Analysez la figure obtenue. Comparez-là avec celle réalisée avec les blastn.
Restreindre aux OG sans paralogues
Rscript ~/work/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_nopara.pdf --max_paralogs=0 evince work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_nopara.pdf
Question 1.12: Comparez les deux figures, avec et sans les paralogues. Comment se répartissent les régions peu denses en gènes orthologues?
Restreindre à un seul OG
Rscript ~/work/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_OG5.pdf --select_OG=5 evince work/Prochlorococcus/images/8.all.ort.group_OG5.pdf
Liste des membres du groupe 5
Aaaa.g_00266 Aaaa.g_00937 Aaaa.g_01196 Aaab.g_00239 Aaab.g_00820 Aaab.g_00825 Aaac.g_01231 Aaac.g_01233 Aaac.g_02086 Aaad.g_00743 Aaad.g_00747 Aaad.g_00909 Aaae.g_00242 Aaae.g_00786 Aaae.g_00791 Aaaf.g_00249 Aaaf.g_00782 Aaaf.g_00787 Aaag.g_00262 Aaag.g_00862 Aaag.g_00867 Aaah.g_00255 Aaah.g_00856 Aaah.g_01180 Aaai.g_00757 Aaai.g_00761 Aaai.g_00954 Aaaj.g_00252 Aaaj.g_00786 Aaaj.g_00791 Aaak.g_00251 Aaak.g_00834 Aaak.g_00839 Aaal.g_00941 Aaal.g_00943 Aaal.g_02214
Extraction du sous-graphe correspondant à ces gènes:
for i in Aaaa*g_00266 Aaaa\|.g_00937 Aaaa\|.g_01196 Aaab\|.g_00239 Aaab\|.g_00820 Aaab\|.g_00825 Aaac\|.g_01231 Aaac\|.g_01233 Aaac\|.g_02086 Aaad\|.g_00743 Aaad\|.g_00747 Aaad\|.g_00909 Aaae\|.g_00242 Aaae\|.g_00786 Aaae\|.g_00791 Aaaf\|.g_00249 Aaaf\|.g_00782 Aaaf\|.g_00787 Aaag\|.g_00262 Aaag\|.g_00862 Aaag\|.g_00867 Aaah\|.g_00255 Aaah\|.g_00856 Aaah\|.g_01180 Aaai\|.g_00757 Aaai\|.g_00761 Aaai\|.g_00954 Aaaj\|.g_00252 Aaaj\|.g_00786 Aaaj\|.g_00791 Aaak\|.g_00251 Aaak\|.g_00834 Aaak\|.g_00839 Aaal\|.g_00941 Aaal\|.g_00943 Aaal\|.g_02214 do grep $i ~/work/Prochlorococcus/PorthoMCL/8.all.ort.tsv done | sort G5.gr | uniq > ~/work/Prochlorococcus/PorthoMCL/G5.gr
Tracez le graphe de ces gènes avec la librairie igraph (c.f. TP5 de M1).
MSK
Question 1.13: Commentez ces résultats. Pouvez vous formuler des hypothèses sur l'évolution de ces gènes? Comment vérifier ces hypothèses?
Tester l'effet de l'IF sur la taille des OG
Tester l'effet de l'IF sur la paralogie
PanOCT with Prochlorococcus
Pan-genome Ortholog Clustering Tool, is a program written in PERL for pan-genomic analysis of closely related prokaryotic species or strains. Unlike traditional graph-based ortholog detection programs, it uses micro synteny or conserved gene neighborhood (CGN) in addition to homology to accurately place proteins into orthologous clusters.
Vous trouverez le package dans le repertoire suivant:
/home/formation/public_html/M2_Phylogenomique/PanOCT/panoct_v3.23/
Test:
cd /home/formation/public_html/M2_Phylogenomique/PanOCT/panoct_v3.23/example_dir ../bin/panoct.pl -b ../example_dir -t example_blast.txt -.pep -S Y -L 1 -M Y -H Y -V Y -N Y -F 1.33 -G y -c 0,25,50,75,100 -T
Préparation des fichiers d'entrée de PanOCT
srun --pty bash mkdir -p ~/work/Prochlorococcus/panoct/results
gene.att
Créer un fichier avec les coordonnées, noms, fonction et souches des gènes.
cd ~/work/Prochlorococcus for file in peptide/*.faa do prefix=$(basename $file .faa) echo "~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output prokka/$prefix/$prefix.tab" ~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output panoct/results/$prefix.tab done ls panoct/results/*.tab
cat panoct/results/*.tab > panoct/results/combined.att head panoct/results/combined.att
tags.txt
Liste des souches à analyser.
for i in peptide/*.faa do echo $(basename $i .faa) done > panoct/results/genomes.list more panoct/results/genomes.list
peptides
Concaténer les fichier peptides dans un seul fichier.
cat peptide/*.faa > panoct/results/combined.fasta head panoct/results/combined.fasta
blast.txt
Concaténer les résultats des blastp dans un seul fichier.
cat BlastP/*.tab > panoct/results/combined.blast head panoct/results/combined.blast
run panOCT
panoct.pl: avec les paramètres par défaut.
-t: name of btab (wublast-style or ncbi -m8 or -m9) input file [REQUIRED] -f: file containing unique genome identifier tags [REQUIRED] -g: gene attribute file (asmbl_id<tab>protein_identifier<tab>end5<tab>end3<tab>annotation<tab>genome_tag) -P: name of concatinated .pep file [REQUIRED to calc protein lengths]
La commande:
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/panoct.pl -b results -t combined.blast -f genomes.list -g combined.att -P combined.fasta -S yes -L 1 -M Y -H Y -V Y -N Y -F 1.33 -G y -c 0,50,95,100 -T
Lancer avec sbatch et un script du type "panoct_P.csh" (les chemins sont à changer).
sbatch panoct_P.csh squeue -l -u <user>
Analyses pan-génomiques
Les analyses pan-génomiques fournissent un cadre pour déterminer la diversité génomique de l'ensemble des gènomes analysés, mais aussi pour prédire, par extrapolation, combien de séquences génomiques supplémentaires seraient nécessaires pour caractériser l'ensemble du pan-génome ou répertoire génétique.
Inside the Pan-genome - Methods and Software Overview
Définitions
D'après Vernikos et al. 2015
Le pan-génome a été défini pour la première fois par Tettelin et al., 2005.
- le pan-génome englobe l'ensemble du répertoire de gènes accessibles au clade étudié ;
- le génome coeur contient les gènes communs à toutes les souches du clade et comprend généralement des gènes responsables des aspects fondamentaux de la biologie du clade et de ses principaux traits phénotypiques ;
- le génome accessoire est constitué des gènes communs à un sous-ensemble de souches et contribue à la diversité des espèces. Il peut coder des voies biochimiques supplémentaires et des fonctions qui ne sont pas essentielles à la croissance mais qui confèrent des avantages sélectifs, comme l'adaptation à différentes niches, la résistance antibiotique ou la colonisation d'un nouvel hôte (Medini et al.', 2005)
- les gènes spécifiques d’une souche ou singletons désignent des gènes spécifiques à une souche n’ayant par d’orthologues dans les autres souches du clade.
Mise en oeuvre
Ce script permet de prendre en compte les paralogs (Chan et al., 2015).
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/paralog_matchtable.pl -M results/matchtable_0_1.txt -P results/paralogs.txt > results/matchtable_paralog.txt
Celui-ci permet de reformater la table matchtable.
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/matchtable2panoct.result.pl --results_dir results
Ce script calcule un ensemble de statistiques sur les résultats de panOCT.
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/pangenome_statistics.pl --data_file results/combined.att.dat --genomes_list results/genomes.list --method_result results/panoct.result -c results/centroids.fasta -f results/frameshifts.txt --fusion results/fragments_fusions.txt
Un résumé des résultats se trouve dans la fichier: ~/work/Prochlorococcus/panoct/overview_stats.txt
Question 1.14: Commentez les résultats obtenus. Comparez ces résultats à ceux de Kettler et al., 2007.
Tracé de l'histogramme
cat overview_stats.txt | perl -ne 'chomp; if ($p) {print "$_\n";} if (/^Cluster Size Breakdown/) {$p = 1;}' > core_cluster_histogram_data.txt /home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/core_cluster_histogram.R -i core_cluster_histogram_data.txt mv core_cluster_histogram.pdf ~/work/Prochlorococcus/images/.
Question 1.15: Après avoir indiqué sur la figure, le génome cœur, le génome accessoire, les singletons et le pan génome commentez les résultats obtenus. Quels sont les choix méthodologiques qui peuvent influencer la taille du génome cœur (argumentez)?
Taille du pan génome
La taille du pan génome à tendance à augmenter avec le nombre de génomes comparés (pour une revue: Tettelin et al., 2008). Estimer sa taille est un exemple d'une classe générale de mesures où, étant donné une collection d’entités et leurs attributs, le nombre d'attributs distincts observés est fonction du nombre d'entités considérées. C’est par exemple le cas de l’analyse du langage naturel, où les entités sont les textes et les attributs sont les mots. Dans ce contexte, l’augmentation du nombre n d'attributs distincts en fonction du nombre N d'entités considérées suit la loi de Heaps (Heaps 'law). Elle peut être représentée par la formule :
n=k*N-α,
où, dans un contexte génétique, n est le nombre attendu de gènes pour un nombre N de génomes et les paramètres k et α (α=1-γ) sont des paramètres libres qui sont déterminés empiriquement (Tettelin et al., 2008).
Appliqué à la découverte de nouveaux gènes et selon la loi de Heap, lorsque α > 1 (γ < 0), le pan-génome est considéré comme fermé, et l'ajout de nouveaux génomes n'augmentera pas significativement le nombre de nouveaux gènes. Par contre, lorsque α < 1 (0 < γ < 1), le pan-génome est ouvert, et pour chaque nouveau génome ajouté, le nombre de gènes augmente significativement (Tettelin et al., 2008).
Pour déterminer les paramètres k et α nous pouvons calculer toutes les combinaisons de 2 à N génomes et reporter la taille du pan génome pour chaque combinaison. Cependant, comme le nombre de combinaisons augmente très rapidement avec le nombre de génomes (C=N!/(n−1)!⋅(N−n), en pratique un échantillonnage des combinaisons possible est réalisé.
compute_pangenome.R
compute_pangenome.R usage:
-i <input file name for tab delimited 0/1 pan-genome cluster table with column and row headers> -o <output file name for combinatorial counts> -p <percentage of genomes to be considered core - default 100> -q <percentage of genomes to be considered novel - default 0 but clearly at least one genome> -s <maximum number of combinatorial samples to generate>
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/compute_pangenome.R -i results/matchtable_paralog.txt -o results/pangenome_size -p 95 -q 0 -s 250
Ajout de la taille des génomes en première colonne:
for i in results/*.tab; do awk 'END{ printf "1\t"NR"\t0"NR"\t"NR"\t"NR"\t"NR"\t\n"}' $i; done > results/unique.txt cat results/pangenome_size results/unique.txt > results/pangenome_size_comp
Pan génome et génome coeur
Rscript /home/formation/public_html/M2_Phylogenomique/scripts/pan_core_plot.R evince ~/work/Prochlorococcus/images/pangenome.pdf
Question 1.16: Décrire la figure obtenue. Comment sont obtenus les différents points gris? A quoi correspondent les triangles et les carrés de couleurs et les deux courbes associées? Comparez cette figure à la Figure 1 de Kettler et al., 2007.
PanOCT classification des proteins de l'OG 5 de PorthoMCL
Nous avons vu que le OG 5 de PorthoMCL renfermait 36 séquences. Comment se distribuent ces protéines avec panOCT?
MSK
Nous pouvons annoter le graphe de la classe 5 trouvée avec PorthoMCL avec les classes trouvées avec panoct. Comme les noms des séquences ne suivent pas la même règle, nous devons modifier le fichier classes pour ajouter un '|' entre le nom de souche et le nom de gène.
MSK
MSK
Alignement et comparaison de génomes complets
- Jeu de données
Vous pouvez retrouver les informations sur les 12 génomes de Prochlorococcus ici et les données dans le répertoire: /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/
Copiez les 12 génomes de Prochlorococcus dans un répertoire de votre ~/work, par exemple genomes_prochlo:
mkdir -p ~/work/Alignement_genomes/genomes_prochlo cp /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/*.fas ~/work/Alignement_genomes/genomes_prochlo/ ls ~/work/Alignement_genomes/genomes_prochlo/
Exploration de la diversité génomique à partir de l’ANI et des distances Mash
Diversité génomique basée sur l’ANI
- Calculer l’ANI entre toutes les paires de génomes en utilisant la version basée sur Mummer. http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/2018_supports/CoursAligntGenomes2018.pdf
srun --pty bash module load system/Python-3.6.3 module load bioinfo/mummer-4.0.0beta2 average_nucleotide_identity.py -h average_nucleotide_identity.py -v -i ~/work/Alignement_genomes/genomes_prochlo/ -o ~/work/Alignement_genomes/genomes_ANIm_output/ --gformat png,pdf,eps,svg --write_excel -m ANIm
Exemple de script "RunSLURM_ANI.csh" (les chemins sont à changer):
sbatch RunSLURM_ANI.csh squeue -l -u <user>
Question 2.1: Regardez les différents fichiers résultats. Regardez la couverture et le pourcentage d’identité des alignements et commentez les valeurs obtenues. Qu’en concluez-vous sur la diversité des génomes de Prochlorococcus ?
- Construire un arbre de Neighbor-Joining basé sur le ANI (ANIm_percentage_identity et ANIm_alignment_coverage) avec le logiciel de votre choix
Vous pourrez par exemple utiliser la fonction nj du package R ape. Notez que la commande nj prend en entrée une matrice de distance. La fonction heatmap (r-graph-gallery) peut être utile pour visualiser les relations entre les souches.
id_file <- "work/Alignement_genomes/genomes_ANIm_output/ANIm_percentage_identity.tab" id_data <- read.table(file=id_file, header=T, row.names=1) heatmap(as.matrix(id_data), scale="none", symm=T, main="ANIm_percentage_identity") co_file <- "work/Alignement_genomes/genomes_ANIm_output/ANIm_alignment_coverage.tab" co_data <- read.table(file=co_file, header=T, row.names=1) heatmap(as.matrix(co_data), scale="none", symm=T, main="ANIm_alignment_coverage") id_nj <- nj(as.matrix(1-id_data)) plot(id_nj, main="ANIm_percentage_identity") co_nj <- nj(as.matrix(1-co_data)) plot(co_nj, main="ANIm_alignment_coverage")
Question 2.2: Interprétez les résultats.
- Sélectionnez les génomes en ne gardant que le sous-groupe de 6 génomes qui ont au moins 28% de couverture avec tous les autres génomes (pour cela regardez le fichier ANIm_alignement_coverage.tab)
Question 2.3: Citez les.
Distance Mash entre les génomes
Passez à l'étape suivante.
Alignements Mauve et ProgressiveMauve
NB : Commencez par lancer l’alignement ProgressiveMauve (environ 50-60 minutes de temps d’execution) avant de faire la question sur l'alignement Mauve !!!
Alignements Mauve d'un sous-ensemble de 6 génomes
mkdir -p ~/work/Alignement_genomes/cat_genomes_prochlo
Concaténer les 6 génomes sélectionnés à la question précédente dans un fichier multifasta
MSK
Lancement de l’alignement des 6 génomes sur le cluster SLURM
mkdir ~/work/Alignement_genomes/Mauve cd ~/work/Alignement_genomes
Exemple de script "RunSLURM_Mauve_6GProchlo.csh" (les chemins sont à changer):
sbatch RunSLURM_Mauve_6GProchlo.csh squeue -l -u <user>
Analyser et interpréter les résultats en les visualisant via l’interface Mauve (commande Mauve)
Remarques:
- dans le fichier Mauve/6GC_Prochlorococcus_gbk_mauve.xmfa, le chemin du fichier gbk et relatif, penser à lancer Mauve dans le bon répertoire pour avoir les annotations des gènes.
- lien entre le code et le nom de souche: species_strain_names.txt
#FormatVersion Mauve1 #Sequence1File cat_genomes_prochlo/6GC_Prochlorococcus.gbk #Sequence1Entry 1 #Sequence1Format GenBank #Annotation1File cat_genomes_prochlo/6GC_Prochlorococcus.gbk ...
Exploration du contexte génomique
L'outil Sequence Navigator (les jumelles) permet de rechercher un ou plusieurs gènes sur différents critères. Nous allons utiliser cette fonctionnalité pour analyser le contexte génomique des gènes suivants. Les noms des gènes sont accessibles par locus tag. En vous plaçant sur le gène, vous avez ses annotations avec View Genbank.... En quittant les jumelles, vous pouvez analyser la conservation du contexte à différentes échelles.
Aaab.g_00239 Aaab.g_00820 Aaab.g_00825 Aaag.g_00262 Aaag.g_00862 Aaag.g_00867 Aaaj.g_00252 Aaaj.g_00786 Aaaj.g_00791 Aaaf.g_00249 Aaaf.g_00782 Aaaf.g_00787 Aaak.g_00251 Aaak.g_00834 Aaak.g_00839 Aaae.g_00242 Aaae.g_00786 Aaae.g_00791
Question 2.5: Combien y’a-t-il de LCB dans l’alignement ? Quel est leur poids minimal ? Y’a–t-il des réarrangements globaux dans l’alignement et si oui lesquels ? Décrire la structure de l’alignement. Que se passe-t-il si on fait varier le poids des LCB ? Qu'avez-vous appris de l'analyse du contexte génomique des gènes.
Alignement Progressive Mauve de l’ensemble complet des 6 génomes
Afin de comparer les logiciels Mauve et ProgressiveMauve nous allons analyser l'ensemble de 6 génomes avec ProgressiveMauve.
ls mkdir ~/work/Alignement_genomes/ProgressiveMauve cd ~/work/Alignement_genomes
Créer un ficher .csh en prenant pour exemple le fichier "RunSLURM_PMauve_6GProchlo.csh" avec comme ligne de commande:
progressiveMauve --output=ProgressiveMauve/6GC_Prochlorococcus_PMauve.xmfa cat_genomes_prochlo/6GC_Prochlorococcus.gbk
Question 2.6: Comparez et interprétez les résultats obtenus.
Alignement Progressive Mauve de l’ensemble complet des 12 génomes
Concaténer les 12 génomes dans un fichier multifasta
MSK
Lancer l’alignement ProgressiveMauve des 12 génomes sur le cluster SLURM
Exemple de script "RunSLURM_PMauve_12GProchlo.csh" (les chemins sont à changer).
sbatch RunSLURM_PMauve_12GProchlo.csh squeue -l -u <user>
Question 2.7: Analyser et interpréter les résultats en les visualisant via l’interface Mauve Si vous avez l'annotation des gènes dans le résultat de ''Progressive Mauve'', vous pouvez visualiser la conservation du contexte de ces gènes de l'OG 5 de PorthoMCL et proposer une interprétation des liens complexes existant entre ces gènes homologues.
L'article de Yan et al., 2018 Genome rearrangement shapes Prochlorococcus ecological adaptation. Appl Environ Microbiol 84:e01178-18. https://doi.org/10.1128/AEM.01178-18 peut vous aider dans l'interprétation des résultats.
Reconstruction de l’histoire évolutive des réarrangements et des états ancestraux
masqué
Analyses phylogénomiques
Comme dans Kettler et al., 2007, nous allons utiliser quatre génomes de Synechococcus comme groupe externe dans nos analyses.
Mise à jour des scripts!
cp /home/formation/public_html/M2_Phylogenomique/scripts/* ~/work/scripts
Synechococcus
Génomes et annotations
- accession "CP000435, CP000097, BX548020, CP000110"
Les fichiers GenBank peuvent être obtenus avec la commande wget en utilisant les chemins enregistrés dans le fichier: accession_file.lst
Comme précédemment, les génomes sont renommés en utilisant un code à quatre lettres.
Les fichiers GenBank renommés sont disponibles dans le répertoire: GenBank
Les séquences ADN ont été extraites des fichiers GenBank : DNA
Fichier avec les informations sur les souches: species_strain_names.txt
mkdir -p ~/work/Synechococcus/prokka cd ~/work/Synechococcus ~/work/scripts/prokka_loop.pl --sample Synechococcus squeue -l -u <user>
Une fois les jobs terminés, vérifiez que les fichiers de sorties de prokka existent et ne sont pas vides.
ls -l ~/work/Synechococcus/prokka/Aaa*/*.faa
Synechococcus blastp All-All
Nous allons reproduire le même enchaînement de scipts utilisés avec Prochlorococcus (Atelier_Phylogénomique#Blast_All-All) en prenant soin de remplacer Prochlorococcus par Synechococcus dans les noms des répertoires.
Copiez les fichier *.faa de prokka dans ~/work/Synechococcus/peptide/
Créer la blast database.
MSK
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 makeblastdb.sh squeue -l -u <user>
Boucle sur les genomes de Synechococcus.
MSK
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh squeue -l -u <user>
squeue -l -u <user> ls BlastP | wc -l
Prochlorococcus versus Synechococcus
Liens symboliques sur les fichiers peptides
mkdir -p ~/work/ProchlorococcusSynechococcus/peptide ln -s ~/work/Prochlorococcus/peptide/*.faa* ~/work/ProchlorococcusSynechococcus/peptide/. ln -s ~/work/Synechococcus/peptide/*.faa* ~/work/ProchlorococcusSynechococcus/peptide/. ls -l ~/work/ProchlorococcusSynechococcus/peptide/.
Liens symboliques sur les fichiers blastp
mkdir -p ~/work/ProchlorococcusSynechococcus/BlastP ln -s ~/work/Prochlorococcus/BlastP/*.tab ~/work/ProchlorococcusSynechococcus/BlastP/. ln -s ~/work/Synechococcus/BlastP/*.tab ~/work/ProchlorococcusSynechococcus/BlastP/. ls -l ~/work/ProchlorococcusSynechococcus/BlastP/.
Compléter les paires de comparaisons
MSK
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh squeue -l -u <user>
Vérifiez que vous avez bien le nombre de fichiers attendus!
Groupes de gènes orthologues avec PanOCT
Préparation des fichiers combined
srun --pty bash mkdir -p ~/work/Synechococcus/panoct/results mkdir -p ~/work/ProchlorococcusSynechococcus/panoct/results
combined.att Créer un fichier avec les coordonnées, noms, fonction et souches des gènes.
cd ~/work/Synechococcus for file in peptide/*.faa do prefix=$(basename $file .faa) echo "~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output prokka/$prefix/$prefix.tab" ~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output panoct/results/$prefix.tab done squeue -l -u <user> ls panoct/results/*.tab
cat ~/work/Prochlorococcus/panoct/results/*.tab ~/work/Synechococcus/panoct/results/*.tab> ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att head ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att
genomes.list Liste des souches à analyser.
cd ~/work/ProchlorococcusSynechococcus/ for i in peptide/*.faa do echo $(basename $i .faa) done > panoct/results/genomes.list more panoct/results/genomes.list
combined.fasta Concaténer les fichier peptides dans un seul fichier.
cat peptide/*.faa > panoct/results/combined.fasta grep -c '>' panoct/results/combined.fasta
combined.blast Concaténer les résultats des blastp dans un seul fichier.
cat BlastP/*.tab > panoct/results/combined.blast head panoct/results/combined.blast
run panOCT Prochlorococcus vs Synechococcus
Exemple de script "panoct_PS.csh" (les chemins sont à changer).
cd ~/work/ProchlorococcusSynechococcus/panoct mkdir ~/work/ProchlorococcusSynechococcus/images
sbatch panoct_PS.csh squeue -l -u <user>
Question 3.1: Résumez dans une figure les différentes étapes réalisées jusqu'à l’obtention des groupes de gènes orthologues. N'oubliez pas la légende!
Heatmap des groupes de gènes orthologues
MSK
Question 3.2: Décrivez la figure obtenue. Quelles informations vous apporte-t-elle?
Gènes espèces spécifique
Gènes trouvés chez toutes les souches de Prochlorococcus mais dans aucune de Synechococcus:
Vous allez utiliser un script awk pour réaliser ce filtrage.
- Editez le fichier matchtable.txt pour comprendre sa structure.
- Que représente chaque ligne?
- Comment sont organisées les colonnes?
- Identifiez le motif signalant l'absence d'un gène dans un OG.
- Pour chaque OG:
- Compter le nombre de gènes présents chez lesProchlorococcus.
- Compter le nombre de gènes absents chez lesSynechococcus.
- Editer la ligne si elle remplit les conditions
MSK
Quelles-sont les fonctions des gènes retenus?
MSK
Gènes trouvés chez toutes les souches de Synechococcus mais dans aucune de Prochlorococcus:
MSK
Question 3.3: Analysez ces résultats et les confronter à ceux disponibles dans la littérature.
Vous pouvez vous inspirer des articles de Kettler et al. (2007), Garczarek et al. 2007, et Partensky and Garczarek (2010).
Comparaison des résultats de PanOCT avec et sans les génomes de Synechococcus
Ces pairs sont à extraire des fichiers:
- ~/work/Prochlorococcus/panoct/results/matchtable.txt
- ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt
A laide d'un script awk, pour les deux fichiers
- pour chaque OG,
- sélectionner les gènes de Prochlorococcus,
- les assembler par paires et les imprimer.
MSK
La comparaison peut être réalisée graphiquement avec un digramme de Venn.
MSK
Phylogénie basée sur les séquences des ARNr
Question 4.1: Quel-est l’intérêt de réaliser des arbres avec les séquences de l'ARNr? Quels-sont les ARNr présents dans les génomes de procaryotes? A quelle(s) sous-unité(s) ribosomique sont-ils associés?
Annotation des ARNr
Nous utilisons le logiciel rnammer pour annoter les ARNr (lsu, ssu, tsu) dans les génomes.
search_module rnammer srun --pty bash module load bioinfo/rnammer-1.2 rnammer -S bac -m ssu -f ~/work/Prochlorococcus/prokka/Aaaa/Aaaa_ssu.rrna < /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaaa.fas
Vous allez procéder comme précédemment, avec un script donné à sarray, pour réaliser le rnammer sur tous les fichiers et les trois types d'ARNr.
for s in Prochlorococcus Synechococcus do for t in ssu lsu tsu do for i in /home/formation/public_html/M2_Phylogenomique/data/$s/DNA/*.fas do genome=$(basename "$i" .fas) output="~/work/$s/prokka/"$genome"/"$genome"_"$t".rrna" echo "module load bioinfo/rnammer-1.2; rnammer -S bac -m $type -f $output < $i;" done done done > rnammer.sh cat rnammer.sh
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 rnammer.sh squeue -l -u <user>
Vérifiez que les fichiers de sortie ne sont pas vide!
ls -l ~/work/*/prokka/Aaa*/Aaa*su*.rrna
Concaténer les fichiers:
mkdir ~/work/ProchlorococcusSynechococcus/rRNA cat ~/work/*/prokka/Aaa*/Aaa*lsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/lsu.fas cat ~/work/*/prokka/Aaa*/Aaa*ssu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/ssu.fas cat ~/work/*/prokka/Aaa*/Aaa*tsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/tsu.fas
Question 4.2: Combien de gènes codant pour les gènes d'ARNr sont prédits dans les différentes souches? Commentez.
Alignements des ARNr
Mafft comporte deux options, Q-INS-i et X-INS-i, dans lesquelles les informations de structure secondaire de l'ARN sont prises en compte. Ces méthodes sont adaptées à un alignement global de séquences d'ARNc très divergentes. Pour les ARN relativement conservés, tels que les ARNr SSU et LSU, l'avantage de ces méthodes est faible (Katoh et al., 2103). Nous utilisons la version mafft pour des raisons de rapidités.
ssu
srun --pty bash module load bioinfo/mafft-7.313 mafft --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/ssu.fas > ~/work/ProchlorococcusSynechococcus/rRNA/ssu.aln
lsu
srun --pty bash module load bioinfo/mafft-7.313 mafft --globalpair --thread 1 ~/work/ProchlorococcusSynechococcus/rRNA/lsu.fas > ~/work/ProchlorococcusSynechococcus/rRNA/lsu.aln
tsu
srun --pty bash module load bioinfo/mafft-7.313 mafft --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/tsu.fas > ~/work/ProchlorococcusSynechococcus/rRNA/tsu.aln
Question 4.3: Pensez-vous que les alignements auraient été de meilleure qualité avec mafft-qinsi et l'option --maxiterate 1000?
Arbre avec seaview
Utilisez le logiciel seaview pour calculer les arbres avec les trois types ARNr.
Expérimentez plusieurs méthodes avec différents paramètres.
Question 4.4: Comparez les résultats obtenus.
Éditez les fichiers pour ne retenir qu'une seule copie de chaque gènes par souche. Renommer les séquences par le code à quatre lettres.
Concaténez les trois types d'ARNr et calculer l'arbre avec la méthode de votre choix.
Discutez ces résultats.
Code R pour obtenir une illustration des réarrangements présents entre deux arbres (source: phytools blog).
library('phytools') ta <-read.tree(file='all_mod-PhyML_tree.ph') tl <-read.tree(file='lsu_mod-PhyML_tree.ph') ts <-read.tree(file='ssu_mod-PhyML_tree.ph') plot.cophylo(cophylo(ta,tl,rotate=TRUE),fsize=0.7, link.type="curved", link.col="blue") plot.cophylo(cophylo(ta,ts,rotate=TRUE),fsize=0.7, link.type="curved", link.col="blue")
Arbre SSU avec IQ-TREE
IQ-tree doc.
IQ-TREE utilise ModelFinder (Kalyaanamoorthy et al., 2017) pour sélectionner le meilleur modèle adaptés aux données.
Pour seulement trouver le modèle le mieux adapté sans faire de reconstruction d'arbre, utilisez :
module load bioinfo/iqtree-1.6.7 iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy -m MF -redo -AIC
Les résultats sont dans le fichier : ssu_renamed_simplified.phy.iqtree.
grep 'Best-fit model' ssu_renamed_simplified.phy.iqtree
lsu ssu GTR+F+R2 tsu K2P+G4
Évaluation des supports de branches avec approximation bootstrap ultra-rapide (UFBoot):
iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy -pre ssuGTRFR2bb1000bnni -m GTR+F+R2 -bb 1000 -redo -bnni -nt AUTO"
NOTE: les valeurs de support de l'UFBoot ont des interprétations différentes de celles du bootstrap non paramétrique. Suivez le lien UFBoot support values interpretation pour plus d'information.
Évaluer les supports de branche avec des tests de branche simple :
IQ-TREE propose le test du rapport approximatif de vraisemblance de type SH (Guindon et al., 2010).
iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy -pre ssuGTRFR2bbalrt -m GTR+F+R2 -bb 1000 -alrt 1000 -redo -nt AUTO"
Évaluation des supports de branche avec un bootstrap non paramétrique standard :
iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy -pre ssuGTRFR2alrtb -m GTR+F+R2 -alrt 1000 -b 100 -redo -nt AUTO"
Arbre SSU avec FastTree
FastTree doc.
module load bioinfo/FastTree-2.1.10 fasttree -nt -gtr < ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy > ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_fasttree.ph
Comparez et commentez les résultats obtenus avec IQ-TREE et FastTree.
Arbre espèces
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:
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 s'un script et lancement du job avec sbatch
mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA
Créer le fihier 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 "Mon_Login"
~/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?
Edition des alignements avec trimal
Trimer les alignements avec trimal. Nous pouvons utiliser les résultats de trimal pour identifier des alignements de mauvaise qualité.
genologin softwares : trimal
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
Question 5.4: Editez le fichier de résultats et analysez son contenu.
Pour chaque alignement, extraction des % de sites sans indels et de la conservation moyenne des alignements
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
Question 5.5: Vérifiez que vous avez extrait les informations attendues.
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
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()
Question 5.6: Éditez les alignements de plus mauvaises qualités. Commentez.
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}:
~/work/scripts/trimal_selection.pl --alignment ~/work/ProchlorococcusSynechococcus/OG/alignment --gap {0,100} --identity {0, 1}
Question 5.7: Combien d'alignements avez-vous retenu?
Super-alignement
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
~/work/scripts/concat_aligments.pl --alignments ~/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:
~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31.lst
IQ-TREE
genologin softwares : iq-tree
FAQ : http://www.iqtree.org/doc/Frequently-Asked-Questions
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: Lancez le, puis répondez aux questions suivantes car cela va été assez long, vous y reviendrez ensuite. Surtout faites un sbatch --cpus-per-task=4 à partir de genologin. N'oubliez pas de spécifier le modèle dans la ligne de commande vous gagnerez 8 heures de calcul ;-)
module load bioinfo/iqtree-1.6.7 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 reconstruits à partir de la petite sous-unité de l’rRNA. Sur genologin le script é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-1.6.7 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-1.6.7 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) avec la commande cat. Nommez le fichier alltrees.tree.
Si à 17 heure vous n'avez pas obtenu ce fichier, contactez moi par mail à claire.hoede[@]inra.fr, je vous débloquerai.
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("MonPath/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 = "MonPath/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-1.6.7 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 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*.
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-1.6.7
Et à rajouter -nt 1 dans les options.
Question 5.14: Commentez les résultats. Quel est l’arbre le plus différent des autres ? Qu’est-ce qui pourrait l’expliquer ?
Analyse des liens phylogénétiques des gènes du OG5 de PortMCL
Éditez un fichier avec la liste des gènes à étudier.
~/work/scripts/extract_sequences_from_list.pl -fasta ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn --outname ~/work/ProchlorococcusSynechococcus/OG5/OG5_genes.faa --list ~/work/ProchlorococcusSynechococcus/OG5_gene_list.txt ~/work/scripts/aa_to_dna_aln.pl -dna ~/work/ProchlorococcusSynechococcus/OG5/OG5_genes.faa --outdir ~/work/ProchlorococcusSynechococcus/OG5
Réalisez l'analyse phylogénétique avec le logiciel de votre choix. Reportez sur l'arbre obtenu les différents évènements de gain/pert/HT de gènes.
Reconstruction d'états ancestraux
Introduction
Si certaines combinaisons de caractères sont systématiquement associées à plusieurs espèces, cela pourrait suggérer que des forces évolutives, comme la sélection, ont façonné ces associations. Toutefois, les associations non aléatoires de certains traits chez certaines espèces peuvent être dues à l'héritage commun de leurs ancêtres et, par conséquent, les changements concomitants dans le temps ne peuvent être déduits.
Si les caractères ont évolué de façon aléatoire sans association, les espèces plus étroitement apparentées sont plus susceptibles d'être semblables que les autres, créant ainsi des relations apparentes entre les caractères.
Il est donc nécessaire de considérer les relations phylogénétiques entre les espèces lors de l'analyse de leurs caractères.
Deux questions peuvent être abordées lors de l'intégration de la phylogénie dans les données comparatives :
- prise en compte de la non-indépendance inter-espèces dans l'étude des caractères et de leurs relations,
- estimation des paramètres d'évolution des caractères.
Ces deux questions sont étroitement liées. En effet l'impact de la phylogénie sur la distribution des caractères dépend non seulement de la phylogénie mais aussi de la façon dont ces caractères évoluent.
Répertoire
Nous allons travailler dans un nouveau répertoire:
mkdir ~/work/ProchlorococcusSynechococcus/Ancestral_Characters cd ~/work/ProchlorococcusSynechococcus/Ancestral_Characters cp /home/formation/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt ~/work/ProchlorococcusSynechococcus/Ancestral_Characters
Copier l'arbre espèce de référence
Choisissez l'arbre que vous allez utiliser comme référence et copiez le dans le nouveau répertoire.
cp ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31_nuc_GTR+F+R4_rooted.treefile ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph
Composition en GC des génomes
Comme caractère continue, nous allons étudier l'évolution de la composition en G+C des génomes. Le taux de GC est calculé avec gc_freq.pl.
for i in ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.fna do ~/work/scripts/gc_freq.pl --file $i done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt for i in ~/work/Synechococcus/prokka/Aaa*/Aaa*.fna do ~/work/scripts/gc_freq.pl --file $i done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Synechococcus_gc.txt cat ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Synechococcus_gc.txt > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt
Arbre espèces
Nous allons utiliser le logiciel R pour réaliser nos analyses et les librairies ape et phytools.
R library(ape) treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" strains_gc_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt" tr <- read.tree(treefile)
Vérifiez que l'arbre tr est enraciné:
is.rooted(tr) plot(tr) plot(tr <- root(tr, interactive = TRUE)) # pour changer la racine de façon interactive
Informations sur les souches
strains_info_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt" strains_info <- read.table(file=strains_info_file, header=T,sep="\t", row.names=1) strains_info <- strains_info[tr$tip.label,]
Changement des labels de l'arbre
tr$tip.label <- as.character(strains_info$Strain) tipcol <- c() tipcol[1:4] <- 'red' tipcol[5:16] <- 'blue' pdf(file="~/work/ProchlorococcusSynechococcus/images/species_tree.pdf", paper="a4r") plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol) nodelabels(tr$node, bg='white', col='purple', frame='n', adj=c(1,-1), cex=0.5) tiplabels(strains_info$Light, bg='white', col='black', frame='n', adj=c(-0.2, 0.5), cex=1) add.scale.bar(length=0.1) dev.off()
evince ~/work/ProchlorococcusSynechococcus/images/species_tree.pdf
GC et taille des génomes
Lire les fichiers arbre et données:
library(phytools) treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" tr <- read.tree(treefile) strains_gc_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt" data <- read.table(file=strains_gc_file, header=F,sep="\t", row.names=1)
Ajout des noms des colonnes et réordonnement des lignes de strains_gc en fonction des feuilles de l'arbre tr.
pdf(file="~/work/ProchlorococcusSynechococcus/images/length_GC.pdf", paper="a4r") colnames(data) <- c('Length', 'GC') data<- data[tr$tip.label,] plot(data$Length, data$GC, pch=20, ylim=c(0,1), main="Length/GC") dev.off()
Number of genes
Question 6.1: Existe-t-il un lien entre la taille du génome et son contenu en GC? Commentez. Existe-t-il un lien entre la taille du génome et son contenu en gènes?
Cartographie de l'évolution d'un caractère continue sur l'arbre
La fonction trace l'arbre sur lequel est reporté l'évolution d'un caractère continu. La cartographie est réalisée en estimant les états aux nœuds internes à l'aide du ML avec la fonction fastAnc, et l'interpolation des états le long de chaque branche (Felsenstein, 1985).
gc <- data$GC names(gc) <- rownames(data) ln <- data$Length names(ln) <- rownames(data) pdf(file="~/work/ProchlorococcusSynechococcus/images/length_GC_length.pdf", paper="a4r") XX <- contMap(tr, gc, sig=2, outline=F, lwd=5, lims=c(0.3,0.6)) XX <- contMap(tr, ln, sig=2, outline=F, lwd=5) dev.off()
Question 6.2: Commentez les figures obtenues.
Reconstruction par Maximum de vraisemblance
Il existe plusieurs fonctions R utilisant le ML pour estimer les états ancestraux aux nœuds internes pour les caractères continus.: ace(method="ML"), fastAnc, et anc.ML. Les trois méthodes effectuent l'estimation du ML en supposant un modèle brownien pour modéliser l'évolution du caractère. La seule différence est la méthode d'optimisation. fastAnc, par exemple, utilise une méthode de ré-enracinement dans laquelle l'arbre est ré-enraciné à chaque nœud interne et l'algorithme de contraste est utilisé pour obtenir l'état ML pour ce noeud. anc.ML utilise une optimisation numérique, mais initie la recherche avec les valeurs obtenues avec fastAnc.
Comparaison des trois méthodes sur notre exemple (pour chaque méthode '<method>$ace' est un vecteur contenant les états des nœuds internes)
aceML <- ace(gc,tr,type="continuous",method="ML") fAnc <- fastAnc(tr,gc,vars=TRUE,CI=TRUE) ancML <- anc.ML(tr,gc) obj<-cbind(aceML$ace,fAnc$ace,ancML$ace) colnames(obj)<-c("ace(method=\"ML\")","fastAnc","anc.ML") pdf(file="~/work/ProchlorococcusSynechococcus/images/ML_fastAnc_anc.pdf", paper="a4r") pairs(obj,pch=21,bg="grey",cex=1.5) dev.off()
Question 6.3: Que peut on conclure de cette analyse?
Tracé de l'arbre
pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_aceML.pdf", paper="a4r") plotTree(tr) nodelabels(pie=aceML$ace,cex=0.5) dev.off()
Lien entre caractères
Felsenstein a proposé une méthode qui prend en compte la phylogénie dans l'analyse des données comparatives. L'idée qui sous-tend sa méthode des contrastes est que, si nous supposons qu'un trait continu évolue de façon aléatoire dans n'importe quelle direction (modèle de mouvement brownien), le contraste entre deux espèces devrait avoir une distribution normale avec une moyenne nulle et une variance proportionnelle au temps écoulé depuis la divergence. Si les contrastes sont mis à l'échelle avec ces derniers, alors ils ont une variance égale à un.
La fonction pic (phylogenetically independent contrasts)
lns <- ln/10000000 pic.gc <- pic(gc, tr) pic.lns <- pic(lns, tr) pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_contrasts.pdf", paper="a4r") plot(tr) nodelabels(round(gc, 2), adj = c(0, -0.5), frame = "n") nodelabels(round(lns,2), adj = c(0, 1), frame = "n") tiplabels(round(gc, 2), adj = c(-1, -0.5), frame = "n") tiplabels(round(lns,2), adj = c(-1, 1), frame = "n") dev.off()
Lien entre les pairs de contrastes
plot(pic.lns, pic.gc) pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_pairs_contrasts.pdf", paper="a4r") plot(pic.lns, pic.gc, xlim=c(-0.15,0.15), ylim=c(-0.15,0.15)) abline(a = 0, b = 1, lty = 2) # x = y line cor(pic.gc, pic.lns) lm(pic.lns ~ pic.gc) dev.off()
Question 6.4: Existe-t-il un lien entre la taille du génome et son contenu en GC (coefficients de la régression)? Commentez.
Remarque, faire la régression entre les PICs en passant par l'origine est justifiée si les caractères évoluent sous un modèle de mouvement brownien et s'il existe une relation linéaire entre eux.
Autocorrélation phylogénétique
masquée
Habitats
Question 6.6: Rappelez d'ou proviennent les différents isolats analysés et expliquez comment ont été défini les écotypes.
Nous allons maintenant nous intéresser à la profondeur à laquelle les organismes vivent. Cette information est absente pour la souche aaap. Nous allons donc supprimer cette feuille à l'arbre.
tr_d <- drop.tip(tr, "Aaap") strains_info_d <- strains_info[rownames(strains_info) != 'Aaap', ] dp <- strains_info_d$Depth names(dp) <- rownames(strains_info_d) pdf(file="~/work/ProchlorococcusSynechococcus/images/Habitats_profondeur.pdf", paper="a4r") XX <- contMap(tr_d, dp) dev.off()
Question 6.7: Existe-t-il un lien entre profondeur et l'évolution des souches (voir Biller et al., 2015)? Commentez.
Ecotypes
Codage des ecotypes HL et LL
ecotypes <- c() ecotypes[grep('Syn', strains_info$Light)] <- 'Syn' ecotypes[grep('LL', strains_info$Light)] <- 'LL' ecotypes[grep('HL', strains_info$Light)] <- 'HL' names(ecotypes)<-tr$tip.label ecotypes
Reconstruction des états du caractère aux nœuds de l'arbre:
pdf(file="~/work/ProchlorococcusSynechococcus/images/Ecotypes.pdf", paper="a4r") eco_er <- ace(ecotypes, tr, type="discrete", CI=T, model="ER") plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol) tiplabels(strains_info$Light, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1) nodelabels(pie=eco_er$lik.anc,cex=0.5) add.scale.bar(length=0.1) dev.off()
Vous pouvez utiliser toutes l'information des ecotype et refaire l'inférence des états du caractère aux nœuds de l'arbre.
ecotypes <- strains_info$Light
Question 6.8: Sur quelles branches de l'arbre placez vous les transitions d'écotypes? Commentez.
Profils des gènes
R library(ape)
Lecture de l'arbre espèces
treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" tr <- read.tree(treefile)
Matrice de présence/absence
Nous allons utiliser les groupes de gènes orthologues calculés par panOCT.
matchtable_file <- "~/work/ProchlorococcusSynechococcus/panoct/results//matchtable_0_1.txt" genomes.list_file <- "~/work/ProchlorococcusSynechococcus/panoct/results//genomes.list" matchtable <- read.table(matchtable_file, sep="\t", row.names=1) genomes.list <- read.table(genomes.list_file) colnames(matchtable) <- t(genomes.list) head(matchtable)
Transposer la matrice et ordonner les lignes comme les feuilles de l'arbre
t_matchtable <- t(matchtable) t_matchtable <- t_matchtable[tr$tip.label,] nb_strains <- nrow(t_matchtable)
Compiler les différents motifs observés dans une liste
Plusieurs groupes de gènes orthologues peuvent avoir le même profil phylogénétique. Il n'est donc pas nécessaire de faire les calculs sur la totalité.
motifs[[pattern]]$sring : motif associé aux souches motifs[[pattern]]$occurences : occurence du motif
motifs <- list() for ( cluster in 1:ncol(t_matchtable) ) { pattern <- paste(t_matchtable[,cluster], sep='', collapse='') nb <- length(motifs[[pattern]]$occurences) if ( nb == 0 ) { motifs[[pattern]]$sring <- t_matchtable[,cluster] } motifs[[pattern]]$occurences[[nb+1]] <- cluster cat(cluster, pattern, nb, sep="\t", "\n") }
Définir un ensemble de motifs
Pour la mise en main des scripts, il peut être utile de travailler avec un sous ensemble de motifs!
nb_start <- 1 nb_end <- 10 if ( nb_end > length(motifs)) nb_end <- length(motifs)
ou
nb_start <- 1 nb_end <- length(motifs)
Taille des patterns trouvés
nb <- nb_start pat_size <- c() pat_name <- c() i <- 1 while (nb <= nb_end ) { pattern <- names(motifs)[nb] cat(pattern, "\n") size <- length(motifs[[pattern]]$occurences) motifs[[pattern]]$size <- size if ( size > 10 ) { pat_size[i] <- size pat_name[i] <- pattern #cat(pattern, size, "\n") i <- i+1 } nb <- nb+1 } names(pat_size) <- pat_name pat_sort <- sort(pat_size, decreasing=T) pdf(file="~/work/ProchlorococcusSynechococcus/images/taille_des_motifs.pdf", paper="a4r") par(mar=c(5,10,1,1)) barplot(pat_sort, horiz=T, cex.names = 0.6, las=1, cex=1) dev.off()
Question 6.9: Quel sera l'usage de ce calcul?
Inférence des états ancestraux
Les différentes fonctions utilisent des paramètres similaires.
Le nombre d'états est déduit des données.
Le modèle est spécifié avec une matrice d'entiers représentant les indices des paramètres : 1 représente le premier paramètre, 2 le second, et ainsi de suite.
Vous devrez choisir la matrice qui spécifie les probabilités de transition entre les états de votre caractère discret. Ace propose des raccourcis pour
- un modèle à un paramètre (ER) où tous les taux de transitions sont égaux,
- un modèle symétrique (SYM) dans lequel les taux de transitions dans les deux sens sont égaux,
- un modèle où tous les taux sont différents (ARD), toutes les transitions possibles entre les états reçoivent des paramètres distincts.
Vous pouvez également spécifier votre propre matrice. Pour indiquer qu'une transition est impossible, un zéro doit être donné dans la cellule appropriée de la matrice.
La reconstruction la plus parcimonieuse
Cette fonction (Most Parsimonious Reconstruction our MPR) réalise la reconstruction des caractères ancestraux par parcimonie, comme décrit dans Hanazawa et al. (1995) et modifié par Narushima et Hanazawa (1997). Ils ont utilisé l'approche proposée par Farris (1970) et Swofford and Maddison’s (1987) pour reconstruire les états ancestraux. Le caractère est supposé prendre des valeurs entières. L'algorithme recherche les ensembles de valeurs pour chaque nœud sous forme d'intervalles avec des valeurs inférieures et supérieures.
L'arbre doit être non enraciné et entièrement dichotomique.
Exemple
pdf(file="~/work/ProchlorococcusSynechococcus/images/plus_parcimonieuse.pdf", paper="a4r") nb <- 3 pattern <- names(motifs)[nb] mrp <- MPR(motifs[[pattern]]$sring, unroot(tr), "Aaao") plot(unroot(tr), label.offset=0.2, show.node.label=F, cex=1) tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1) nodelabels(paste("[", mrp[, 1], ",", mrp[, 2], "]", sep = "")) dev.off()
Tracé MPR pour chaque motif
pdf(file="~/work/ProchlorococcusSynechococcus/images/plus_parcimonieuse_all.pdf", paper="a4r") nb <- nb_start nb_end < 10 while (nb < nb_end ) { pattern <- names(motifs)[nb] if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) { mrp <- MPR(motifs[[pattern]]$sring, unroot(tr), "Aaao") plot(unroot(tr), label.offset=0.2, show.node.label=F, cex=1, main=paste("MRP", nb)) tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame="n", adj=c(-0.5, 0.5), cex=1) nodelabels(paste("[", mrp[, 1], ",", mrp[, 2], "]", sep = ""), bg="white", col="purple") #cat ("Press [enter] to continue") #line <- readline() } nb <- nb+1 } dev.off() nb_end <- length(motifs)
Marquage stochastique des caractères
La fonction make.simmap de phytools permet de réaliser un marquage stochastique des caractères aux nœuds de l'arbre.
- nsim: number of simulations.
- pi: gives the prior distribution on the root node of the tree, options are equal, estimated, or a vector with the frequencies.
Exemple
motif 3 ("0000001100101111")
pdf(file="~/work/ProchlorococcusSynechococcus/images/marquage_stochastique.pdf", paper="a4r") nb <- 3 pattern <- names(motifs)[nb] ms_er_trees <- make.simmap(tr, motifs[[pattern]]$sring, model="ER", nsim=100, pi="estimated") cols<-setNames(c("blue","red"),c(0,1)) plotSimmap(ms_er_trees[[1]], cols) dev.off()
Calculons maintenant les fréquences d'état à partir des marquages stochastiques pour chaque noeud interne. Ce sont nos probabilités postérieures pour les noeuds.
Fonction pour calculer les états:
map_states <- function(x){ y <- sapply(x$maps,function(x) names(x)[1]) names(y) <- x$edge[,1] y <- y[as.character(length(x$tip)+1:x$Nnode)] return(y) }
Apliquer la fonction à tous les arbres et tracer l'arbre et les probabilités postérieures aux nœuds.
pdf(file="~/work/ProchlorococcusSynechococcus/images/probabilites_posterieures.pdf", paper="a4r") AA <- sapply(ms_er_trees, map_states) piesA <- t(apply(AA,1,function(x,levels,Nsim) summary(factor(x,levels))/Nsim,levels=c(0,1), Nsim=100)) plot(tr, label.offset=0.2, show.node.label=F, cex=1) tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1) nodelabels(pie=piesA, cex=0.5, piecol=cols) dev.off()
Maximum de vraisemblance
La commande ace de la librairie ape peut reconstruire des états ancestraux pour des caractères discrets en utilisant le maximum de vraisemblance (Pagel 1994 Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters. Proceedings of the Royal Society of London. Series B. Biological Sciences, 255, 37–45.).
Par défaut, la probabilité que les différents états ancestraux des caractères discrets sont calculés à l'aide d'une procédure d'estimation conjointe (joint) en utilisant une procédure semblable à celle décrite dans Pupko et al. (2000). Si "marginal = VRAI", une procédure d'estimation marginale est utilisée. Avec cette méthode, les valeurs de vraisemblance à un nœud donné sont calculées en utilisant uniquement les informations provenant des feuilles (et des branches) descendantes de ce nœud. Avec l'estimation conjointe, toutes les informations sont utilisées pour chaque nœud. La mise en œuvre actuelle de l'estimation conjointe utilise un algorithme "en deux passes" qui est beaucoup plus rapide que l'approche stochastique alors que les estimations des deux méthodes sont très proches.
Si l'option CI = TRUE est utilisée, alors la probabilité de chaque état ancestral est retournée pour chaque noeud dans une matrice appelée lik.anc. For discrete characters only maximum likelihood estimation is available
obj<-anc.Bayes(tree,x,ngen=10000) # sample ancestral states print(obj,printlen=20) ## estimates
Remarque: la fonction fitDiscrete de geiger utilise la même syntaxe que ace et donne des résultats similaires, mais les valeurs de log-vraisemblance sont mises à l'échelle différemment cependant le test du rapport de vraisemblance reste le même.
Exemple motif 3 ("0000001100101111")
pdf(file="~/work/ProchlorococcusSynechococcus/images/ML3.pdf", paper="a4r") nb <- 3 pattern <- names(motifs)[nb] er <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ER") sym <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="SYM") ard <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ARD") plot(tr, label.offset=0.2, show.node.label=F, cex=1) tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1) nodelabels(pie=er$lik.anc, cex=0.3, adj=c(0.5, 0.5)) nodelabels(pie=sym$lik.anc, cex=0.3, adj=c(0.53, 0.5)) nodelabels(pie=ard$lik.anc, cex=0.3, adj=c(0.56, 0.5)) dev.off()
La sortie du modèle ER correspond à celle du modèle SYM. Pour un caractère binaire, ces modèles sont identiques, bien qu'ils soient distincts pour les caractères à trois états ou plus. Pour un caractère à trois états, ER est un modèle à un paramètre, SYM un modèle à trois paramètres et ARD un modèle à six paramètres.
Le modèle ARD donne généralement la probabilité la plus élevée. Cependant, il inclut aussi plus de paramètres que le modèle ER. Il est bien connu que l'ajout de paramètres à un modèle augmente généralement sa probabilité. Pour déterminer si l'utilisation du modèle ARD est appropriée, vous devez exécuter un test de vraisemblance.
Test de vraisemblance
La différence de logarithme de probabilité entre deux modèles est distribuée sous forme de chi2, avec des degrés de liberté égaux au nombre de paramètres qui diffèrent entre les modèles. Ainsi, un test de vraisemblance de ces modèles demande si la différence de vraisemblance est suffisamment grande pour se situer dans la queue la plus à droite de la distribution du chi2.
La fonction pchisq(value, df) donne le pourcentage de la fonction de distribution cumulative pour chi2 qui se trouve à gauche de la valeur donnée pour un degré de liberté souhaité(df). Ainsi,
1-pchisq(2*abs(er$loglik - ard$loglik), 1)
donne le pourcentage de la distribution du chi2 qui se trouve à droite de la différence de vraisemblance observée pour les modèles ER et ARD (qui diffèrent d'un paramètre).
Boucle sur tous les motifs
Par curiosité nous allons réaliser le test pour tous les motifs et imprimer les valeurs les plus significatives.
nbnodes <- length(tr$node) nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2, ) nb <- nb_start nbpattern <- 1 #while (nb <= length(motifs) ) { while (nb <= nb_end ) { pattern <- names(motifs)[nb] if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) { er <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ER") motifs[[pattern]]$er <- er ard <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ARD") motifs[[pattern]]$er <- ard ml_test <- 1-pchisq(2*abs(er$loglik - ard$loglik), 1) if ( ml_test < 0.01 ) { cat(nb, pattern, er$loglik, ard$loglik, ml_test, "\n") } } nb <- nb+1 }
Fonction pour annoter les états des feuilles de l'arbre
Les inférences des états du caractère au nœuds sont calculées sur les noeuds interne de l'arbre, nous ajoutons ici les noeuds feuilles.
tip_states <- function(pattern) { leaves <- matrix(c(0), nrow=nb_strains, ncol=2) for ( i in 1:nb_strains ) { leaves[i,1] <- 0 leaves[i,2] <- 0 if (pattern[i] == 0 ) { leaves[i,1] <- 1 } else { leaves[i,2] <- 1 } #cat(pattern[i], leaves[i,1], leaves[i,2], "\n") } return(leaves) }
Exemple d'appel de la fonction
leaves <- tip_states(motifs[[pattern]]$sring)
Concaténation des deux matrices (feuilles et noeuds)
edge_states <- rbind(leaves, er$lik.anc)
Fonction pour annoter les transitions d'états sur les branches de l'arbre
Une transition d'état est identifiée sur une branche (edge) si l'état du caractère présent au nœud parent est différent de l'état du caractère présent au nœud fils.
- En entrée l'arbre tr et les états du caractère inférés aux noeuds (edge_states).
- En sortie une liste avec les listes de transitions (type de tansition, direction et une couleur).
states_transitions <- function(tr, edge_states) { nb_edges <- nrow(tr$edge) direction <- c() transition <- c() transcolor <- c() for ( i in 1:nb_edges ) { child <- tr$edge[i,2] parent <- tr$edge[i,1] if ( edge_states[child, 1] > 0.5 ) { child_states <- 0 } else { child_states <- 1 } if ( edge_states[parent, 1] > 0.5 ) { parent_states <- 0 } else { parent_states <- 1 } if ( parent_states != child_states ) { transition <- c(transition, i) if ( parent_states == 0 ) { direction <- c(direction, '+') transcolor <- c(transcolor, 'blue') } else { direction <- c(direction, '-') transcolor <- c(transcolor, 'red') } cat(i, 'parent', parent , parent_states, '->', child_states, 'child', child, direction, "\n") } } translist <- list(transition=transition, direction=direction, transcolor=transcolor) return(translist) }
Exemple d'appel de la fonction
states_transitions(tr, edge_states)
Tracé pour chaque motif
nbnodes <- length(tr$node) nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2) mycolors <- c('red', 'blue') nb <- nb_start nbpattern <- 1 #while (nb <= length(motifs) ) { while (nb < nb_end ) { pattern <- names(motifs)[nb] if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) { leaves <- tip_states(motifs[[pattern]]$sring) edge_states <- rbind(leaves, motifs[[pattern]]$er$lik.anc) st <- states_transitions(tr, edge_states) plot(tr, label.offset=0.2, show.node.label=F, cex=1) vec <- as.vector(motifs[[pattern]]$sring) bgcol <- unlist(lapply(vec, FUN= function(x) mycolors[x+1])) tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5) nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors) edgelabels(text=st$direction, edge=st$transition, frame="r", bg=st$transcolor, adj=c(0,-0.5),cex=0.2) add.scale.bar(length=0.1) cat ("Press [enter] to continue") line <- readline() } nb <- nb+1 }
Flux de gènes
nbnodes <- length(tr$node) nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2) mycolors <- c('red', 'blue') nb <- nb_start nbpattern <- 1 gain <- vector (mode='numeric',length=nrow(tr$edge)) loss <- vector (mode='numeric',length=nrow(tr$edge)) while (nb < nb_end ) { pattern <- names(motifs)[nb] if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) { leaves <- tip_states(motifs[[pattern]]$sring) edge_states <- rbind(leaves, motifs[[pattern]]$er$lik.anc) st <- states_transitions(tr, edge_states) nb_trans <- length(st$transition) if ( nb_trans > 0 ) { for ( j in 1:nb_trans) { edge <- st$transition[j] cat(nb_trans, j, edge, st$direction[j], "\n") if ( st$direction[j] == "+" ) { cat('gain', nb_trans, j, st$direction[j], "\n") gain[edge] <- gain[edge] + motifs[[pattern]]$size } else if ( st$direction[j] == "-" ) { loss[edge] <- loss[edge] - motifs[[pattern]]$size } } } else { cat("no transition, unexpected!\n") } } nb <- nb+1 }
Gains et pertes de gènes
pdf(file="~/work/ProchlorococcusSynechococcus/images/gains_pertes_genes.pdf", paper="a4r") plot(tr, label.offset=0.2, show.node.label=F, cex=1) tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5) nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors) edgelabels(text=gain, frame="n", col='blue', bg='white', adj=c(0,-0.5),cex=0.8) edgelabels(text=loss, frame="n", col='red', bg='white', adj=c(0,1.5),cex=0.8) add.scale.bar(length=0.1) dev.off()
evince ~/work/ProchlorococcusSynechococcus/images/gains_pertes_genes.pdf
Flux de gènes
pdf(file="~/work/ProchlorococcusSynechococcus/images/flux_genes.pdf", paper="a4r") plot(tr, label.offset=0.2, show.node.label=F, cex=1) flux <- gain + loss plot(tr, label.offset=0.2, show.node.label=F, cex=1) tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5) nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors) edgelabels(text=flux, frame="n", col='black', bg='white', adj=c(0,-0.5),cex=0.8) add.scale.bar(length=0.1) dev.off()
Question 6.10: Commentez les figures obtenues. Comparez ces résultats avec ceux de Sun et Blanchard, 2014)
evince ~/work/ProchlorococcusSynechococcus/images/flux_genes.pdf
Inférence Bayésienne
BayesTraits
BayesTraits is a computer package for performing analyses of trait evolution among groups of species for which a phylogeny or sample of phylogenies is available. This new package incoporates our earlier and separate programes Multistate, Discrete and Continuous. BayesTraits can be applied to the analysis of traits that adopt a finite number of discrete states, or to the analysis of continuously varying traits. Hypotheses can be tested about models of evolution, about ancestral states and about correlations among pairs of traits.
La mise en œuvre de ce logiciel est complexe et ne sera pas traitée dans le cadre de cd TP.
Liens utiles
Mesquite
Mesquite: A modular system for evolutionary analysis