Atelier Phylogénomique
From silico.biotoul.fr
Matériel pédagogique
Support de cours : Quest for Orthologs
Support de cours : Alignements de génomes et Phylogénomique
Support de cours : Reconstruction des états ancestraux des caractères
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.
- 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.
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:
anemone CHERYN arome MARINE aster SOUFIANE chardon EVE camelia Tristan capucine ELISE gerbera CHLOE clematite NICOLAS cobee JULIETTE coquelicot CLARA cosmos GEOFFREY cyclamen MATTHIEU dahlia JULIEN digitale MARION geranium ALEXANDRE
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
1. First write a script (ex: myscript.sh) with the command line as following:
#!/bin/bash #SBATCH -J test #SBATCH -o output.out #SBATCH -e error.out #SBATCH -t 01:00:00 #SBATCH --mem=8G #SBATCH --mail-type=BEGIN,END,FAIL (the email address is automatically LDAP account's one) #Purge any previous modules module purge #Load the application module load bioinfo/ncbi-blast-2.2.29+ # My command lines I want to run on the cluster blastall ...
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):
Installation du package genoPlotR
srun --pty bash module load system/R-3.5.1 R install.packages('genoPlotR') Warning in install.packages("genoPlotR") : 'lib = "/tools/R/R-3.5.1/lib64/R/library"' ne peut être ouvert en écriture Would you like to use a personal library instead? (yes/No/cancel) yes Would you like to create a personal library ‘~/R/x86_64-pc-linux-gnu-library/3.5’ to install packages into? (yes/No/cancel) yes --- SVP sélectionner un miroir CRAN pour cette session --- Secure CRAN mirrors 26: France (Lyon 2) [https] * installing *source* package ‘ade4’ ... * installing *source* package ‘genoPlotR’ ... Les packages source téléchargés sont dans ‘/tmp/RtmpLdNj99/downloaded_packages’ library(genoPlotR)
Scripts
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
srun --pty bash module load bioinfo/prokka-1.12 prokka /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaaa.fas --outdir /tmp/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, 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 programme pour compléter les paramètres de prokka pour chaque génome.
Le programme doit être lancé sur le serveur genologin et prokka est lancé sur les noeuds avec sbatch.
mkdir -p ~/work/Prochlorococcus/prokka 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 sorties de prokka existent et ne sont pas vides.
ls -l ~/work/Prochlorococcus/prokka/Aaa*/*.faa
Les fichiers avec le suffix .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. 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.
Séquence conservation 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.
Exemple avec une paire:
srun --pty bash module load bioinfo/ncbi-blast-2.6.0+ 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 BlastN/Aaab_vs_Aaag.tab
Sur l'ensemble des génomes:
~/work/scripts/blastn_genome_pair.pl squeue -l -u <user> 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
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
Question 1.3: Commentez le résultat obtenu. Que pensez-vous de la conservation des séquences des génomes?
ACT
Il est également possible d'utiliser le logiciel act (documentation).
Copier les fichiers sur votre porte de travail en P0 et lancer :
act prokka/Aaab/Aaab.gbk BlastN/Aaab_vs_Aaag.tab prokka/Aaag/Aaag.gbk
Vous pouvez aussi utiliser les fichiers en format gff.
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 cd ~/work/Prochlorococcus/peptide cp ~/work/Prochlorococcus/prokka/Aaa*/*.faa ~/work/Prochlorococcus/peptide/. ls -l ~/work/Prochlorococcus/peptide
make blast database
cd ~/work/Prochlorococcus/peptide module load bioinfo/ncbi-blast-2.6.0+
Exemple:
srun -n1 -l makeblastdb -in Aaaa.faa -dbtype prot
Boucle sur tous les fichiers
for i in *.faa; do srun -n1 -l makeblastdb -in $i -dbtype prot; done
Intra genomes
Nous vous proposons un script perl pour réaliser les blastp de l'ensemble des génomes avec sbatch:
mkdir -p ~/work/Prochlorococcus/BlastP cd ~/work/Prochlorococcus ~/work/scripts/blastp_intra.pl squeue -l -u <user> ls BlastP | wc -l
Question 1.6: Explicitez les paramètres de blast utilisés dans le script blastp_intra.pl. Combien de fichiers attendez-vous?
Inter genomes
Pour les blastp inter génomes, nous allons utiliser un script similaire au précédent, mais avec une double boucle (sur -query et -db).
cd ~/work/Prochlorococcus ~/work/scripts/blastp_inter.pl squeue -l -u <user> ls BlastP | wc -l
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.
Disponible: /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master.
Le logiciel est basé sur l'enchainement de programmes python. Ils peuvent être lancés séparément ou automatiquement à l'aide d'un script shell.
Version automatique (que nous n'utiliserons pas car nous avons déjà réalisé les blastp):
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?
Running 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 6.orthologs.-5.50/*.tsv >> 8.all.ort.-5.50.tsv mcl 8.all.ort.-5.50.tsv --abc -I 1.5 -t 4 -o 8.all.ort.group.1.5-5.50
Taille des OG:
awk 'BEGIN { n=1; printf("class\tsize\n"); } { printf(n"\t"NF"\n"); n = n+1; }' PorthoMCL/8.all.ort.group.1.5-5.50 > PorthoMCL/8.all.ort.group.1.5-5.50_class_size.tab
Que pensez-vous de la distribution des taille des OG? Quelle-est la taille maximum attendue?
Paralogs
cat 7.paralogs.-5.50/*.tsv >> 8.all.par.tsv mcl 8.all.par.tsv --abc -I 1.5 -t 4 -o 8.all.par.group.I1.-5.50
Visualisation of 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.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_all.pdf evince work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_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.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_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.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_OG5.pdf --select_OG=5 evince work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_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 PorthoMCL/8.all.ort.-5.50.tsv done
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 surla taille des OG
Lancer MCL avec différents IF (respecter la syntaxe des noms!)
Le motif suivant est attendu: 8.all.ort.group.<IF>-5.50.
Rscript ~/work/scripts/MCL_partition_inflate.R --MCL_dir=~/work/Prochlorococcus/PorthoMCL --nb_strains=12 pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_IF.pdf evince ~/work/Prochlorococcus/images/8.all.ort.group_IF.pdf
Code couleurs:
bleu : tous les OG rouge : OG sans paralogues vert : OG constitués uniquement d'orthologues 1:1 orange : OG avec au moins un paralogue(s)
Tester l'effet de l'IF sur la paralogie
Rscript ~/work/scripts/MCL_OG_paralogy.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_paralogs.pdf Rscript ~/work/scripts/MCL_OG_paralogy.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.8.0-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.8.0-5.50_paralogs.pdf
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 squeue -l -u <user> 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
with default paramerters
-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]
cd panoct /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
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
Afficher les courbes
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.
Log x Log plot des nouveaux gènes
Rscript /home/formation/public_html/M2_Phylogenomique/scripts/new_genes_plot.R evince ~/work/Prochlorococcus/images/newgenome.pdf
Question 1.17: Décrire cette figure. Que représente le 'a' (alpha)? Comparez la valeur obtenue avec celle reportée dans Tettelin et al., 2008. Que pourrait être l'effet de l'ajout de nouveaux génomes sur la taille du pan-génome (voir Biller et al., 2015)? Que pouvez-vous en conclure?
PanOCT classification des pro teins 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?
rm ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt 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 awk -v var=$i '{ pattern=var; if (match($0,pattern)) { printf("%s\t%d\n", var,$1); } }' ~/work/Prochlorococcus/panoct/all_clusters_members.txt >> ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt done sort -nk 2 ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt
Question 1.18: Commentez le résultat obtenu. Que pouvez-vous en conclure?
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.
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
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
- Concaténer les 6 génomes sélectionnés à la question précédente dans un fichier multifasta
for i in Aaax Aaay Aaaz do echo $i cat ~/work/Alignement_genomes/genomes_prochlo/$i.fas >> ~/work/Alignement_genomes/genomes_prochlo/6GC_Prochlorococcus.fna cat ~/work/Prochlorococcus/prokka/$i/$i.gbk >> ~/work/Alignement_genomes/genomes_prochlo/6GC_Prochlorococcus.gbk done
- 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_PMauve_6GProchlo.csh" (les chemins sont à changer):
#!/bin/bash #SBATCH --time=02:00:00 #job time limit #SBATCH -J Mauve_6GProclo #SBATCH -o RunSLURM_Mauve_6GProclo.out #SBATCH -e RunSLURM_Mauve_6GProclo.err #SBATCH --mem=8G #SBATCH --cpus-per-task=4 #ncpu on the same node #SBATCH --chdir=/home/<user>/work/Alignement_genomes #SBATCH --mail-type=BEGIN,END,FAIL (email address is LDAP accounts) #Purge any previous modules module purge #Load the application module load bioinfo/mauve_2.4.0 # My command lines I want to run on the cluster mauveAligner --output=Mauve/6GC_Prochlorococcus_gbk.mauve_def --permutation-matrix-output=Mauve/6GC_Prochlorococcus_gbk.permutation_matrix --output-guide-tree=Mauve/6GC_Prochlorococcus_gbk.tree --output-alignment=Mauve/6GC_Prochlorococcus_gbk_mauve.xmfa genomes_prochlo/6GC_Prochlorococcus.gbk
sbatch RunSLURM_PMauve_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 genomes_prochlo/6GC_Prochlorococcus.gbk #Sequence1Entry 1 #Sequence1Format GenBank #Annotation1File genomes_prochlo/6GC_Prochlorococcus.gbk ...
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 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 ?
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_Mauve_6GProchlo.csh avec comme ligne de commande:
progressiveMauve --output=ProgressiveMauve/6GC_Prochlorococcus_PMauve.xmfa genomes_prochlo/6GC_Prochlorococcus.fna
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
for i in Aaaa Aaab Aaac Aaad Aaae Aaaf Aaag Aaah Aaai Aaaj Aaak Aaal do echo $i cat ~/work/Alignement_genomes/genomes_prochlo/$i.fas >> ~/work/Alignement_genomes/genomes_prochlo/12GC_Prochlorococcus.fna done grep -c '>' ~/work/Alignement_genomes/genomes_prochlo/12GC_Prochlorococcus.fna
cat ~/work/Prochlorococcus/prokka/*/*gbk >> 12GC_Prochlorococcus.gbk
Lancer l’alignement ProgressiveMauve des 12 génomes sur le cluster SLURM
Exemple de script "RunSLURM_PMauve_12GProchlo.csh" (les chemins sont à changer):
#!/bin/bash #SBATCH --time=02:00:00 #job time limit #SBATCH -J ProgressiveMauve_12GProclo #SBATCH -o RunSLURM_ProgressiveMauve_12GProclo.out #SBATCH -e RunSLURM_ProgressiveMauve_12GProclo.err #SBATCH --mem=8G #SBATCH --cpus-per-task=4 #ncpu on the same node #SBATCH --chdir=/home/<user>/work/Alignement_genomes #SBATCH --mail-type=BEGIN,END,FAIL (email address is LDAP accounts) #Purge any previous modules module purge #Load the application module load bioinfo/mauve_2.4.0 # My command lines I want to run on the cluster progressiveMauve --output=ProgressiveMauve/12GC_Prochlorococcus_PMauve.xmfa --output-guide-tree=ProgressiveMauve/12GC_Prochlorococcus_PMauve.ph genomes_prochlo/12GC_Prochlorococcus.fna
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
Mauve 12GC_Prochlorococcus_pmauve.xmfa
Si vous avez l'annotation des gènes dans le résultat de ''Progressive Mauve'', vous pouvez tenter de visualiser la conservation du contexte de ces gènes de l'OG 5 de PorthoMCL et proposer une interprétation des liens complexes existants entre ces gènes homologues.
Reconstruction de l’histoire évolutive des réarrangements et des états ancestraux
- Exporter le fichier de permutation des LCBs de l’alignement généré à la question précédente avec ProgressiveMauve en cliquant dans l’interface Mauve sur Tools => Export => Export Permutation. Ainsi vous pouvez créer les fichiers permutations.
- 6GC_Prochlorococcus_pmauve.permutation (pour se familiariser sur les résultats!)
- 12GC_Prochlorococcus_pmauve.permutation
NB1: Il n’est pas possible de lancer ProgressiveMauve directement avec l’option –permutation-matrix (ne fonctionne qu’avec Mauve).
NB2: Il est nécessaire d'éditer le fichier 12GC_Prochlorococcus_pmauve.permutation pour le mettre au format MLGO : ajouter une ligne « >id genome » dans le même ordre que le fichier d’entrée 12GC_Prochlorococcus.fna et remplacer toutes les virgules par des espaces. Ajouter un espace avant le $ de chaque fin de ligne
grep ">" 12GC_Prochlorococcus.fna
Le fichier doit avoir le format suivant
>Aaaa 1 56 55 54 26 5 52 53 27 $ >Aaab 1 56 55 54 26 5 52 53 27 $ Etc...
- Utiliser le logiciel MLGO pour inférer l’histoire évolutive des réarrangements à partir de la matrice de permutation modifiée (n'oubliez pas de calculer les bootstraps ).
Lien vers MLGO
Question 2.6: Visualiser les arbres produits par MLGO en utilisant l’outil de votre choix (par exemple FigTree ou iTOL) et interpréter les résultats en fonction des écotypes.
NB1: Le fichier gene_order.tree.PMAG contient l’arbre avec les nœuds ancestraux (A1,A2,…etc) et le fichier gene_order.out contient l’ordre des gènes (ici LCBs) dans ces génomes ancestraux.
NB2: Le fichier de sortie MLGO gene_order_bs.tree qui contient les valeurs de bootstrap entre [] est non standard et n’est pas lisible par FigTree. Pour pouvoir le lire effectuer la transformation suivante:
more gene_order_bs.tree (Aaab:0.00038261697365918421,(((Aaal:0.00110779153484151733,Aaac:0.00045702512028280982):0.01128794608013189253[100],((Aaah:0.00195679558538791641,Aaaa:0.00249248732070895488):0.00170192911132310224[98],(Aaad:0.00000001574704787426,Aaai:0.00000001574704787426):0.00384727219137918474[100]):0.00019414433073579879[65]):0.00328248641555709160[100],(Aaaj:0.00000001574704787426,(Aaaf:0.00000001574704787426,(Aaae:0.00000001574704787426,Aaak:0.00000001574704787426):0.00000001574704787426[34]):0.00000001574704787426[33]):0.00056868823569016999[98]):0.00018460659686999882[50],Aaag:0.00040691005233194506); Est à transformer en: (Aaab:0.00038261697365918421,(((Aaal:0.00110779153484151733,Aaac:0.00045702512028280982)100:0.01128794608013189253,((Aaah:0.00195679558538791641,Aaaa:0.00249248732070895488)98:0.00170192911132310224,(Aaad:0.00000001574704787426,Aaai:0.00000001574704787426)100:0.00384727219137918474)65:0.00019414433073579879)100:0.00328248641555709160,(Aaaj:0.00000001574704787426,(Aaaf:0.00000001574704787426,(Aaae:0.00000001574704787426,Aaak:0.00000001574704787426)34:0.00000001574704787426)33:0.00000001574704787426)98:0.00056868823569016999)50:0.00018460659686999882,Aaag:0.00040691005233194506);
Ou plus simplement, vous pouvez utiliser iTOL pous visualiser l'arbre directement avec les valeurs de bootstraps et éventuellement l'annoter.
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:
mkdir -p ~/work/Synechococcus/peptide cd ~/work/Synechococcus/peptide cp ~/work/Synechococcus/prokka/Aaa*/*.faa ~/work/Synechococcus/peptide/. ls -l ~/work/Synechococcus/peptide
module load bioinfo/ncbi-blast-2.6.0+ for i in *.faa; do srun -n1 -l makeblastdb -in $i -dbtype prot; done
mkdir -p ~/work/Synechococcus/BlastP cd ~/work/Synechococcus ~/work/scripts/blastp_intra.pl
~/work/scripts/blastp_inter.pl
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
cd ~/work/ProchlorococcusSynechococcus ~/work/scripts/blastp_inter.pl squeue -l -u <user> ls BlastP | wc -l
Groupes de gènes orthologues
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
cd ~/work/ProchlorococcusSynechococcus/panoct mkdir ~/work/ProchlorococcusSynechococcus/images /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
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
srun --pty bash # si nécessaire! module load system/R-3.5.1; R pdf_file <- '~/work/ProchlorococcusSynechococcus/images/matchtable_heatmap.pdf' strain_file <- '~/work/ProchlorococcusSynechococcus/panoct/results/genomes.list' matchtable_0_1 <- '~/work/ProchlorococcusSynechococcus/panoct/results/matchtable_0_1.txt' strains <- read.delim(file=strain_file, header=FALSE) data <- read.delim(file=matchtable_0_1, , sep="\t", header=FALSE, row.names=1) colnames(data) <- t(strains) pdf(file=pdf_file, paper="a4r") heatmap(t(as.matrix(data)), scale='none', col=c('white', 'darkblue'), xlab="Strains", labCol=NA) dev.off() cat(pdf_file, "\n")
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:
awk 'BEGIN { h=0; } { mp=0 for (i=2; i<=13; i++) { if( $i !~ /-/ ) mp++ } sp=0 for (i=14; i<=17; i++) { if( $i !~ /-/ ) sp++ } if ( mp==12 && sp ==0 ) { h++; printf("%d\t%d", h, $1); for (i=2; i<=17; i++) { if( $i !~ /-/ ) printf("\t%s", $i); } printf("\n"); } }' < ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt > ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt
Expliquez les différentes étapes réalisées par ce script. Quelles-sont les fonctions des gènes retenus?
awk '{ system("grep "$3" ~/work/Prochlorococcus/prokka/Aaaa/Aaaa.gff >> ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt")}' < ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt
Question 3.3: Analysez ces résultats et les confronter à ceux disponibles dans la littérature.
Gènes trouvés chez toutes les souches de Synechococcus mais dans aucune de Prochlorococcus:
A vous de jouer!
Comparaison des résultats de PanOCT avec et sans les génomes de Synechococcus
Pour cela, nous allons identifier les pairs de gènes orthologues obtenues par chaque approche et les comparer.
Ces pairs sont à extraire des fichiers:
- ~/work/Prochlorococcus/panoct/results/matchtable.txt
- ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt
La comparaison peut être réalisée graphiquement avec un digramme de Venn.
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?
Nous utilisons le logiciel rnammer pour annoter les ARNr (lsu, ssu, tsu) dans les génomes.
/home/formation/public_html/M2_Phylogenomique/scripts/rnammer_loop.pl --prokka_dir ~/work/Synechococcus/prokka --model ssu /home/formation/public_html/M2_Phylogenomique/scripts/rnammer_loop.pl --prokka_dir ~/work/Prochlorococcus/prokka --model ssu
Vérifiez que les fichiers de sortie ne sont pas vide!
ls -l ~/work/*/prokka/Aaa*/Aaa*ssu*.rrna
Changer le motif pour obtenir les deux autres ARNr:
ls -l ~/work/*/prokka/Aaa*/Aaa*lsu*.rrna ls -l ~/work/*/prokka/Aaa*/Aaa*tsu*.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 ~/work/ProchlorococcusSynechococcus mkdir ~/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 vi commande.sh i #!/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 :wq! 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
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()
Vous pouvez utiliser filezilla (sous Windows) ou un scp (sous linux) pour télécharger le fichier pdf en local.
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 l’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.
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