silico.biotoul.fr
 

Atelier Phylogénomique

From silico.biotoul.fr

Revision as of 16:40, 25 October 2018 by Desousa (Talk | contribs)
Jump to: navigation, search

Contents

Matériel pédagogique

Support de cours : Quest for Orthologs

Support de cours : Alignements de génomes et Phylogénomique

Références

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		BIRBES		CLEMENT
arome		BOUCHIBA	YOUNES
aster		BURNARD		CALLUM
bleuet		DECAESTECKER	JULIEN
camelia		JEREMIE		MARIE
capucine	KURYLO		CYRIL
chardon		LATOUR		JUSTINE
clematite	LEGUET		MANON
cobee		LUDWIG		CORALIE
coquelicot	MARICOURT	PIERRE-LOUIS
cosmos		PÜZIN		WILLIAM
cyclamen	QUELLERY	ELISABETH
dahlia		SIMON		CHRIS
digitale	TIHIC		EDI
geranium	VAN BLERK	SEBASTIAN
gerbera		YLMAZ		FULYA
glaieul 
hortensia
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

Seminar_cluster_SLURM

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):

tutoR_cluster.pdf

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)

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

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.

prokka files.

Question 1.1:
Pourquoi pensez-vous qu'il soit nécessaire d'annoter les génomes?
Quelles sont les annotations réalisées?
Quels sont les logiciels utilisés?

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:

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.

Question 1.2:
Comparez globalement les résultats obtenus avec ceux reportés dans la publication.
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.

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 paire:

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

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

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.

How to use SLURM 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 ~/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.
Restrict view to OG without paralogs
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?
Restrict view to one 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
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.

panoct project

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?

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 peut être utile pour visualiser les relations entre les souches.

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
Calculer la distance Mash entre toutes les paires de génomes

Documentation : Mash

En mode intéractif:

srun --pty bash
module load system/R-3.4.3 compiler/gcc-7.2.0
module load bioinfo/Mash-2.1
~/work/scripts/Mash_sketch.sh  ~/work/Alignement_genomes/genomes_prochlo/

Les résultats se trouvent dans le répertoire data_MashSketches/.

~/work/scripts/Mash_dist_allpairs.sh data_MashSketches/

Le résultat se trouve dans le fichier mash_dist.out

NB:

  • 1ère colonne: Jaccard index (~fraction de kmers partagés)
  • 2ème colonne : p-valeur
  • 3ème colonne : distance Mash (estimation du taux de mutation selon un modèle d'évolution simple)
Question 2.4:
Interprétez les résultats. 
Comparez les résultats de distance MASH à ceux de l’ANI.

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
cat ~/work/Alignement_genomes/genomes_prochlo/<génome>.fas >> ~/work/Alignement_genomes/genomes_prochlo/6GC_Prochlorococcus.fna
Lancement de l’alignement des 6 génomes sur le cluster SLURM

Exemple de script "RunSLURM_Mauve_6GProchlo.csh" (les chemins sont à changer) :

#!/bin/bash
#SBATCH --time=02:00:00 #job time limit
#SBATCH -J testjob
#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 --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=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.mauve_def --permutation-matrix-output=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.permutation_matrix --output-guide-tree=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.tree --output-alignment=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus_mauve.xmfa  /home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.fna
Analyser et interpréter les résultats en les visualisant via l’interface Mauve (commande Mauve)
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 12 génomes

Concaténer les 12 génomes dans un fichier multifasta

cat /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/*.fas >> 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 Pmauve
#SBATCH -o RunSLURM_PMauve_12GProclo.out
#SBATCH -e RunSLURM_PMauve_12GProclo.err
#SBATCH --mem=8G
#SBATCH --cpus-per-task=4 #ncpu on the same node
#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=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/12GC_Prochlorococcus_PMauve.xmfa /home/hchiapello/work/TP_M2_BIOINFO/TP_2018/12GC_Prochlorococcus.fna
sbatch RunSLURM_PMauve_12GProchlo.csh
Analyser et interpréter les résultats en les visualisant via l’interface Mauve
Mauve 12GC_Prochlorococcus_pmauve.xmfa

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 le fichier 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 l’outil 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) 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);

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.

accession_labels_file.lst

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

Quelles sont les fonctions de ces gènes?

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!

Phylogénie basée sur les séquences des ARNr

Question 3.4:
Quel-est l’intérêt de réaliser des arbres avec les séquences de l'ARNr?

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 3.5:
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 3.6:
Pensez-vous que les alignements auraient été de meilleure qualité avec mafft-qinsi et l'option --maxiterate 1000?

Arbre ARNr

Arbre avec seaview

Question 3.7:
Utilisez le logiciel seaview pour calculer les arbres avec les trois types ARNr.
Expérimentez plusieurs méthodes avec différents paramètres.
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.

Arbre SSU avec IQ-TREE

IQ-tree doc.

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

dna-models

Évaluation des supports de branches avec approximation bootstrap ultra-rapide :

iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.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).

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 espèces

Support de cours : supports

Aujourd’hui nous utiliserons certains gènes de Prochlorococcus et de 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:

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

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

grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
Question 4.1:
A quoi correspond ce quorum. Pourquoi utiliser un seuil de 16? 
Combien d'alignements avez-vous obtenu?
Est-ce pertinent dans les situations où vous avez un grand nombre de souches?

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

~/work/scripts/aa_to_dna_aln_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100
Question 4.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 4.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

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

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")
Question 4.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 4.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 4.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

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

Nous allons inférer un arbre à partir du super-alignement en codons:

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


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 l'arbre ici : /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/all_pep.fas.treefile. Vous pouvez le copier dans votre répertoire ~/work/ProchlorococcusSynechococcus/phyloG/.

Question 4.9: 
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 ?
Question 4.10:
Comparez les deux arbres 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 ».
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 (et pas genotoul).

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

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.

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

Question 4.12:
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).

Lancer ensuite le calcul de la distance de Robinson and Foulds sur ce fichier avec iqtree.

Attention à faire le module load bioinfo/iqtree-1.6.7

Et à rajouter -nt 1 dans les options.

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

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.

Attention le nom de votre arbre de référence est différent pour chacun d'entre vous et si vous ne l'avez pas déjà fais enracinez le puis enregistrez !!

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.

Sur génotoul:

qrsh

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)

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")
Question 5.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)

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")
Question 5.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")
pairs(obj,pch=21,bg="grey",cex=1.5)
Question 5.3:
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)

lns <- ln/10000000
pic.gc <- pic(gc, tr)
pic.lns <- pic(lns, tr)

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")

Lien entre les pairs de contrastes

plot(pic.lns, pic.gc)
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)
Question 5.4:
Existe-t-il un lien entre la taille du génome et son contenu en GC? 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

Comme les espèces ne sont pas indépendantes par leurs relations phylogénétiques, nous pouvons utiliser cette dépendance pour quantifier l'association entre les variables observées au niveau des espèce. Gittleman et Kot[109] propose une méthode basée sur une approche d'autocorrélation. Elle utilise l'indice d'autocorrélation I[208] de Moran qui fait intervenir wij un poids qui quantifie la proximité phylogénétique entre espèces i et j. Avec un arbre phylogénétique, nous avons wij = 1/dij, où dij sont les distances mesurées sur un arbre, et wii = 0.

w <- 1/cophenetic(tr)
diag(w) <- 0
Moran.I(gc, w)

Le résultat est une liste avec quatre éléments :

  • la valeur observée de I (observée),
  • sa valeur attendue sous l'hypothèse nulle d'aucune corrélation,
  • l'écart-type de l'I (sd) observé,
  • la P-valeur de l'hypothèse nulle.

La fonction gearymoran d'ade4 calcule le coefficient de Moran et teste sa signification à l'aide d'une procédure de randomisation. L'option nrepet spécifie le nombre de réplications de la randomisation (999 par défaut).

library(ade4)
gearymoran(w, data.frame(gc, ln))
Question 5.5:
Existe-t-il un lien entre l'évolution des caractères (taille du génome et contenu en GC) et l'évolution des souches? Commentez.

Habitats

Question 5.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)

XX <- contMap(tr_d, dp)
Question 5.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:

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)

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

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 à tous les arbres et tracer l'arbre et 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é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")

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")
      # amotifs[[pattern]]$er <- ard
      # aml_test <- 1-pchisq(2*abs(er$loglik - ard$loglik), 1)
      # aif ( ml_test < 0.01 ) {
      # a   cat(nb, pattern, er$loglik, ard$loglik, ml_test, "\n")
      # a}
   }
   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()
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 5.10:
Commentez les figures obtenues.
Comparez ces résultats avec ceux de Sun et Blanchard, 2014)

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.

BayesTraitsV3.Manual.pdf

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