Atelier Phylogénomique
From silico.biotoul.fr
Matériel pédagogique
- Présentation 'Quest for Orthologs'.
- Références
- Présentation 'Reconstruction of ancestral traits'.
- Ancestral State Reconstruction with phytools
- Ancestral State Reconstruction with phytools2
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.
- 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.
- Chisholm Current Biology 27, R431–R510 Prochlorococcus.
Ressources informatiques
Nous allons utiliser les ressources de GenoToul.
Vous avec dans les FAQ, les réponses à toutes vos questions concernant l'utilisation de la ressource.
Ask for Install
trimal Gblocks PhyML IQ-tree
Sortcuts
Vous allez vous connecter avec un compte fleur:
ssh -Y <login>@genologin.toulouse.inra.fr
Vous pouvez ouvrir deux connections afin de pouvoir lancer qlogin de façon indépendante.
Échange de fichiers:
scp file login@genologin.toulouse.inra.fr:~/work
Logiciels disponibles
ou 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).
SLURM cluster softwares
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 29: 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)
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
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.
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:
rsync --archive --itemize-changes --stats -h -e ssh <user>@genologin.toulouse.inra.fr:/home/<user>/work/ProchlorococcusSynechococcus work
Disponibilité des génomes
Caractéristiques des souches étudiées
Table 1. General Characteristics of the Prochlorococcus and Synechococcus Isolates Used in ([https://www.ncbi.nlm.nih.gov/pubmed/18159947 Kettler et al., 2007).
Cyanobacterium Isolate Light Length(bp) GC% Number Isol. Region Date Accession Adap. Genes Depth Prochlorococcus MED4 HL(I) 1,657,990 30.8 1,929 5m Med. Sea Jan. 1989 BX548174 MIT9515 HL(I) 1,704,176 30.8 1,908 15m Eq. Pacific Jun. 1995 CP000552 MIT9301 HL(II) 1,642,773 31.4 1,907 90m Sargasso Sea Jul. 1993 CP000576 AS9601 HL(II) 1,669,886 31.3 1,926 50m Arabian Sea Nov. 1995 CP000551 MIT9215 HL(II) 1,738,790 31.1 1,989 5m Eq. Pacific Oct. 1992 CP000825 MIT9312 HL(II) 1,709,204 31.2 1,962 135m Gulf Stream Jul. 1993 CP000111 NATL1A LL(I) 1,864,731 35.1 2,201 30m N. Atlantic Apr. 1990 CP000553 NATL2A LL(I) 1,842,899 35 2,158 10m N. Atlantic Apr. 1990 CP000095 CCMP1375/SS120 LL(II) 1,751,080 36.4 1,925 120m Sargasso Sea May 1988 AE017126 MIT9211 LL(III)1,688,963 38 1,855 83m Eq. Pacific Apr. 1992 CP000878 MIT9303 LL(IV) 2,682,807 50.1 3,022 100m Sargasso Sea Jul. 1992 CP000554 MIT9313 LL(IV) 2,410,873 50.7 2,843 135m Gulf Stream Jul. 1992 BX548175 Synechococcus CC9311 Syn. 2,606,748 52.5 3017 95m Calif. Current 1993 CP000435 CC9902 Syn. 2,234,828 54.2 2504 5m Calif. Current 1999 CP000097 WH8102 Syn. 2,434,428 59.4 2787 Sargasso Sea Mar. 1981 BX548020 CC9605 Syn. 2,510,659 59.2 2991 51m Calif. Current 1996 CP000110
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.
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 /home/formation/public_html/M2_Phylogenomique/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 c'est bien passé, vous pouvez supprimer les fichiers .err et .sh.
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.
runVisuSession.sh art ~/work/Prochlorococcus/prokka/Aaaa/Aaaa.gff
Séquence conservation entre souches de Prochlorococcus
Genome 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')
Nous allons utiliser l'option BLAST-2-Sequences de blastn. Exemple avec une pair
srun -n1 -l 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
mkdir ~/work/Prochlorococcus/BlastN /home/formation/public_html/M2_Phylogenomique/scripts/blastn_genome_pair.pl ls -l ~/work/Prochlorococcus/BlastN
genoplotR
Nous allons utiliser genoplotR pour visualiser les similarités entre les pairs 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 fichier en format gbk crées 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 pairs de genomes.
- annotation: Un objet d'annotation est utilisé pour annoter un segment d'ADN. Nous ne l'utilisons ici.
- tree: un arbre au format Newick qui peut être analysée à l'aide du paquetage ade4. Nous ne l'utilisons plus tard!
mkdir ~/work/Prochlorococcus/images srun --pty bash module load system/R-3.5.1 Rscript /home/formation/public_html/M2_Phylogenomique/scripts/genoplot_blastn_links.R evince ~/work/Prochlorococcus/images/genoplot_blastn_links.pdf
ACT
Il est également possible d'utiliser le logiciel act (documentation).
Copier les fichiers sur votre porte de travail en P0 et lancer :
act 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
Préliminaires
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+ 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
blastp -query {input} -db {input} -seg {SEG} -dbsize 100000000 -evalue {EVALUE} -outfmt 6 -num_threads 1 -out {output}
Un script perl peut réaliser les blastp pour l'ensemble des génomes avec sbatch:
mkdir -p ~/work/Prochlorococcus/BlastP cd ~/work/Prochlorococcus /home/formation/public_html/M2_Phylogenomique/scripts/blastp_intra.pl squeue -l -u <user> ls BlastP | wc -l
Combien de fichiers attendez-vous?
Inter genomes
Nous allons utiliser un script similaire au précédant, mais avec une double boucle (sur -query et -db).
cd ~/work/Prochlorococcus /home/formation/public_html/M2_Phylogenomique/scripts/blastp_inter.pl squeue -l -u <user> ls BlastP | wc -l
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
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 /home/formation/public_html/M2_Phylogenomique/scripts/short_PorthoMCL.sh
old stuff
Running MCL
C'est exactement comme la dernière étape d'OrthoMCL. Il suffit d'exécuter MCL sur les sorties des étapes 6 et 7 :
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
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
Global view
Rscript /home/formation/public_html/M2_Phylogenomique/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
Restrict view to OG without paralogs
Rscript /home/formation/public_html/M2_Phylogenomique/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
Restrict view to one OG
Rscript /home/formation/public_html/M2_Phylogenomique/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
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 /home/formation/public_html/M2_Phylogenomique/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
Tester l'effet de l'IF sur la paralogie
Rscript /home/formation/public_html/M2_Phylogenomique/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 /home/formation/public_html/M2_Phylogenomique/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
old stuff2
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 "/home/formation/public_html/M2_Phylogenomique/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output prokka/$prefix/$prefix.tab" /home/formation/public_html/M2_Phylogenomique/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output panoct/results/$prefix.tab done ls panoct/results/*.tab
cat panoct/results/*.tab > panoct/results/combined.att head panoct/results/combined.att
tags.txt
Liste des souches à analyser.
for i in peptide/*.faa do echo $(basename $i .faa) done > panoct/results/genomes.list more panoct/results/genomes.list
peptides
Concaténer les fichier peptides dans un seul fichier.
cat peptide/*.faa > panoct/results/combined.fasta head panoct/results/combined.fasta
blast.txt
Concaténer les résultats des blastp dans un seul fichier.
cat BlastP/*.tab > panoct/results/combined.blast head panoct/results/combined.blast
run panOCT
panoct.pl
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’orthologs 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
Commentez les résultats obtenus. Comparez ces résultats à ceux de [https://www.ncbi.nlm.nih.gov/pubmed/18159947 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/.
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. 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 (https://en.wikipedia.org/wiki/Heaps%27_law 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
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
run_panoct.pl
--genome_list_file,-g : File containing list of genome names --att_dir, -A : Directory containing genome .att files --gene_att_file. -a : Use this gene_att_file instead of creating one from att_dir --blast_file,-b : Path to a file of all-vs-all BLAST results. Must be -m 8 or -m 9 tab delimited output. --working_dir, -w : Path to be consedered the working directory for the pipeline. Defaults to cwd ('./'). --no_grid : Run all subtasks locally instead of in parallel on an SGE/Univa grid.
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/run_panoct.pl --genome_list_file tags.txt --gene_att_file gene_att --blast_file combined.blast --no_grid
results/R.plots/core_cluster_histogram.pdf
compute_pangenome.R
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/compute_pangenome.R -i /share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/matchtable_paralog.txt -o /share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/pangenome_size -p 95 -q 0 -s 250
Afficher les courbes:
R pangenome_size <- '/share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/pangenome_size' pdf_file <- '/share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/pangenome.pdf' data <- read.table(file=pangenome_size, header=T) data ymax <- max(data$Genome, data$Pan_genome) ymin <- 0 pdf(file=pdf_file, paper="a4r") plot(data$Genome, data$Pan_genome, col='blue',xlim=c(2, 12), xlab='Number of genomes', ylab='Number of genes', ylim=c(ymin, ymax), pch=20) points(data$Genome, data$Core, col='red', pch=20) dev.off() cat(pdf_file, "\n")
Analyses phylogénomiques
Comme dans Kettler et al., 2007, nous allons utiliser quatre génomes de Synechococcus comme groupe externe dans nos analyses.
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 /home/formation/public_html/M2_Phylogenomique/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 /home/formation/public_html/M2_Phylogenomique/scripts/blastp_intra.pl
/home/formation/public_html/M2_Phylogenomique/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 pairs de comparisons
cd ~/work/ProchlorococcusSynechococcus /home/formation/public_html/M2_Phylogenomique/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 "/home/formation/public_html/M2_Phylogenomique/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output prokka/$prefix/$prefix.tab" /home/formation/public_html/M2_Phylogenomique/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output panoct/results/$prefix.tab done 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
Heatmap des groupes de gènes orthologues
srun --pty bash 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")
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"); } }' < /home/formation/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt > /home/formation/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt
Quelles sont les fonctions de ces gènes?
awk '{ system("grep "$3" ../Prochlorococcus/prokka/Aaaa/Aaaa.gff >> Prochlorococcus_specific_annotations.txt")}' < /home/formation/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt
Gènes trouvés chez toutes les souches de Synechococcus mais dans aucune de Prochlorococcus:
A vous de jouer!
Extraction des séquences nucléotidiques des gènes orthologues
Créer un fichier avec toutes les séquences nucléotidiques:
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
mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA /home/formation/public_html/M2_Phylogenomique/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
Vérifier que les fichiers contiennent au moins 16 séquences:
grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
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 /home/formation/public_html/M2_Phylogenomique/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.
/home/formation/public_html/M2_Phylogenomique/scripts/aa_to_dna_aln_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100
Construction des arbres des groupes de gènes orthologues
Edition des alignements avec trimal
Trimer les alignements avec trimal (pas sur genologin!). Nous pouvons utiliser les résultats de trimal pour identifier des alignements de mauvaise qualité.
ssh -Y <user>@genotoul.toulouse.inra.fr
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_*.fas do dir=$(dirname $i) prefix=$(basename $i .fas) /usr/local/bioinfo/bin/trimal -in $dir/$prefix.fas -out $dir/$prefix.trim.aln -sident -sgt -automated1 done > ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.txt
Pour chaque alignement, extraction des % de sites avec des 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
awk 'BEGIN { min = max = "NaN" } { min = (NR==1 || $2<min ? $2 : min) max = (NR==1 || $2>max ? $2: max) } END { print min, max }' < ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.tab
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 meuilleurs alignements
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")
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
/home/formation/public_html/M2_Phylogenomique/scripts/trimal_selection.pl --alignment ~/work/ProchlorococcusSynechococcus/OG/alignment -gap 70 --identity 0.4
Phylogénie des souches
Super Matrice
Concaténation de 31 alignements
/home/formation/public_html/M2_Phylogenomique/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
/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 -redo -bb 1000 -alrt 1000 -nt AUTO
Input data: 16 sequences with 32355 nucleotide sites Number of constant sites: 10622 (= 32.8295% of all sites) Number of invariant (constant or ambiguous constant) sites: 10622 (= 32.8295% of all sites) Number of parsimony informative sites: 20377 Number of distinct site patterns: 16346 ModelFinder ----------- Best-fit model according to BIC: GTR+F+R4 Total CPU time used: 2936.33161 seconds (0h:48m:56s) Total wall-clock time used: 459.9323409 seconds (0h:7m:39s)
Enraciner l'arbre
~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31_nuc_GTR+F+R4_rooted.treefile
Phylogénie basée sur les séquences des ARNr
Nous utilisons 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 fichier de sortie ne sont pas vide!
ls -l /home/yquentin/work/*/prokka/Aaa*/Aaa*ssu*.rrna
Changer le motif pour obtenir les deux autres ARNr:
ls -l /home/yquentin/work/*/prokka/Aaa*/Aaa*lsu*.rrna ls -l /home/yquentin/work/*/prokka/Aaa*/Aaa*tsu*.rrna
Concaténer les fichiers:
mkdir ~/work/ProchlorococcusSynechococcus/rRNA cat /home/yquentin/work/*/prokka/Aaa*/Aaa*lsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/lsu.fas cat /home/yquentin/work/*/prokka/Aaa*/Aaa*ssu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/ssu.fas cat /home/yquentin/work/*/prokka/Aaa*/Aaa*tsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/tsu.fas
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 néanmoins la version mafft-qinsi.
ssu
srun --pty bash module load bioinfo/mafft-7.313 mafft-qinsi --maxiterate 1000 --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.fas > ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.aln trimal -in ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.aln -automated1 -phylip -out ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy
lsu
srun --pty bash module load bioinfo/mafft-7.313 mafft-qinsi --maxiterate 1000 --globalpair --thread -1 ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified.fas > ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified.aln trimal -in ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified.aln -automated1 -phylip -out ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified_trimed.phy
tsu
srun --pty bash module load bioinfo/mafft-7.313 mafft-qinsi --maxiterate 1000 --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified.fas > ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified.aln trimal -in ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified.aln -automated1 -phylip -out ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified_trimed.phy
Arbre SSU
Arbre SSU avec PhyML
qsub -V -b Y -N phymlRNAn -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/PhyML/current/PhyML-3.1_linux64 -i ~/work/ProchlorococcusSynechococcus/rRNA/ssu_trimed.phy -n 1 -b -1 -m GTR -v e -c 4 -a e -o n --quiet"
Vous pouvez observer que des génomes codent pour plusieurs copies du gènes codant pour l'ARN 16S. A l'aide de l'arbre et des scores de rnammer simplifier votre fichier. Il est également utilie de modifier le nom des séquences. Après ces modifications, vous relancez l'alignement, l'édition de l'alignement et le calcule de l'arbre avec des options plus poussées.
qsub -V -b Y -N phymlRNAtlr -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/PhyML/current/PhyML-3.1_linux64 -i ~/work/ProchlorococcusSynechococcus/rRNA/ssu_trimed.phy -n 1 -b -1 -m GTR -v e -c 4 -a e -o tlr --quiet"
Arbre SSU avec IQ-TREE
Pour seulement trouver le modèle le mieux adapté sans faire de reconstruction d'arbre, utilisez :qsub -V -b Y -N IQTreeSSUm -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy -m MF -redo -AIC"
Les résultats sont dans le fichier : ssu_renamed_simplified_trimed.phy.iqtree.
grep 'Best-fit model' ssu_renamed_simplified_trimed.phy.iqtree
lsu ssu GTR+F+R2 tsu K2P+G4
Évaluation des supports de branches avec approximation bootstrap ultra-rapide :
qsub -V -b Y -N ssuIQTreeSSUbb -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy -pre ssuGTRFR2bb1000bnni -m GTR+F+R2 -bb 1000 -redo -bnni -nt AUTO"
É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).
qsub -V -b Y -N ssuIQTreeSSUbbalrt -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.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 :
qsub -V -b Y -N ssuIQTreeSSUalrtb -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy -pre ssuGTRFR2alrtb -m GTR+F+R2 -alrt 1000 -b 100 -redo -nt AUTO"
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.
Copier l'arbre espèce de référence
Nous allons travailler dans un nouveau répertoire:
mkdir ~/work/ProchlorococcusSynechococcus/Ancestral_Characters cd ~/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 /home/formation/public_html/M2_Phylogenomique/scripts/gc_freq.pl --file $i done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt for i in ~/work/Synechococcus/prokka/Aaa*/Aaa*.fna do /home/formation/public_html/M2_Phylogenomique/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' 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)
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)
Ajouter les noms des colonnes et réordonner les lignes de strains_gc en fonction des feuilles de l'arbre tr.
colnames(data) <- c('Length', 'GC') data<- data[tr$tip.label,] plot(data$Length, data$GC, pch=20, ylim=c(0,1), main="Length/GC")
Existe-t-il un lien entre la taille du génome et son contenu en GC? Commentez.
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) 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, main ="Lengh")
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 dpour 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 obetnues 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") pairs(obj,pch=21,bg="grey",cex=1.5)
Que peut on conclure de cette analyse?
Tracé de l'arbre
plotTree(tr) nodelabels(pie=aceML$ace,cex=0.5)
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)
ln <- ln/10000 pic.gc <- pic(gc, tr) pic.ln <- pic(ln, tr) plot(tr) nodelabels(round(gc, 2), adj = c(0, -0.5), frame = "n") nodelabels(round(ln,2), adj = c(0, 1), frame = "n") tiplabels(round(gc, 2), adj = c(-2, -0.5), frame = "n") tiplabels(round(ln,2), adj = c(-1, 1), frame = "n")
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 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) eco_sym <- ace(ecotypes, tr, type="discrete", CI=T, model="SYM") eco_ard <- ace(ecotypes, tr, type="discrete", CI=T, model="ARD") The increase in likelihood with the additional parameter is not significant: 1 - pchisq(2*(eco_er$loglik - eco_sym$loglik), 1) In practice, there is an anova method for the class "ace" which simplifies the calculation of this test: anova(eco_er, eco_ard, eco_sym)
plot(tr) nodelabels(pie=eco_er$lik.anc,cex=0.7) </pre>
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
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
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
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) par(mar=c(5,10,1,1)) barplot(pat_sort, horiz=T, cex.names = 1, las=1)
Inférence des états ancestraux
Les différentes fonctions utilise 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ù toutes les taux de transitions sont égales,
- un modèle symétrique (SYM) dans lequel les taux de transitions dans les deux sens sont égales,
- 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
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 = ""))
Tracé MPR pour chaque motif
nb <- nb_start 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 }
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")
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)
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 a tous les arbres et tracer l'arbre est les probabilités postérieures aux nœuds.
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)
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é. 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")
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))
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 et 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érer 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
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)
Flux de gènes
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)
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