Atelier Phylogénomique
From silico.biotoul.fr
m (→Tracé pour chaque motif) |
m (→Tracé pour chaque motif) |
||
Line 1,621: | Line 1,621: | ||
bgcol <- unlist(lapply(vec, FUN= function(x) mycolors[x+1])) | bgcol <- unlist(lapply(vec, FUN= function(x) mycolors[x+1])) | ||
tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5) | tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5) | ||
- | nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, | + | nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors) |
- | edgelabels(text=st$direction, edge=st$transition, frame="r", | + | 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) | add.scale.bar(length=0.1) | ||
Line 1,630: | Line 1,630: | ||
} | } | ||
nb <- nb+1 | nb <- nb+1 | ||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
} | } | ||
Revision as of 13:42, 12 October 2018
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.
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
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
mkdir ~/work/ProchlorococcusSynechococcus/Ancestral_Characters cd ~/work/ProchlorococcusSynechococcus/Ancestral_Characters
Copier l'arbre espèce de référence
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
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
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érifier que l'arbre tr est enraciné: is.rooted(tr) plot(tr) plot(tr <- root(tr, interactive = TRUE))
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.5, 0.5), cex=1) add.scale.bar(length=0.1)
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>
Flux de gènes
matchtable
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)
Simplifier la table
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)
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 <- 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") }
Traiter chaque motifs
Inférence des états ancestraux
nb_start <- 1 nb_end <- 10 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 } nb <- nb+1 }
Fonction pour annoter les états des feuilles de l'arbre
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
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 ) { cat(i, parent , parent_states, '->', child_states, child, "\n") transition <- c(transition, i) if ( parent_states == 0 ) { direction <- c(direction, '+') transcolor <- c(transcolor, 'blue') } else { direction <- c(direction, '-') transcolor <- c(transcolor, 'red') } } } translist <- list(transition=transition, direction=direction, transcolor=transcolor) return(translist) }
Exemple d'appel de la fonction
states_transitions(tr, edge_states)
cat ("Press [enter] to continue") line <- readline() plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol) tiplabels(motifspattern$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1) edgelabels(text=direction, edge=transition, frame="r", bg=transcolor, adj=c(0,-0.5),cex=0.2) </pre>
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, tip.color=tipcol) 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 } ==old stuff== <!-- Lien sur les génomes de ''Prochlorococcus'' du NCBI: *Répertoire avec les [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Prochlorococcus/NCBI souches] du NCBI. Tables avec des information sur les projets de séquençages: *[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Prochlorococcus/NCBI/Prochlorococcus.csv Prochlorococcus.csv] Nous avons édité le fichier pour associer un identifiant de moins de 5 lettres pour chaque génomes (compatibilité avec certains logiciels utilisés). ====Autres souches==== Les génomes des souches suivantes sont absentes du NCBI. *''Wolbachia'' endosymbiont of ''Dirofilaria immitis'' wDi.2.2 [http://dirofilaria.nematod.es/ http://dirofilaria.nematod.es/] *''Wolbachia'' endosymbiont of ''Litomosoides sigmodontis'' wLs.2.0 [http://litomosoides.nematod.es/ http://litomosoides.nematod.es/] [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/Additional Autres souches] téléchargées. ====Groupes externes==== Les souches suivantes peuvent être utilisées pour enraciner un arbre réalisé sur les souches de Wolbachia: *Ace Anaplasma centrale str. Israel — 1 1 206 806 [https://www.ncbi.nlm.nih.gov/nuccore/NC_013532.1?report=fasta NC_013532] *Aph Anaplasma phagocytophilum HZ — 1 1 471 282 [https://www.ncbi.nlm.nih.gov/nuccore/NC_007797.1?report=fasta NC_007797] *Ech Ehrlichia chaffeensis str. Arkansas — 1 1 176 248 [https://www.ncbi.nlm.nih.gov/nuccore/NC_007799.1?report=fasta NC_007799] *Eru Ehrlichia ruminantium str. Gardel — 1 1 499 920 [https://www.ncbi.nlm.nih.gov/nuccore/NC_006831.1?report=fasta NC_006831] *Nri Neorickettsia risticii str. Illinois — 1 879 977 [https://www.ncbi.nlm.nih.gov/nuccore/NC_013009.1?report=fasta NC_013009] *Nse Neorickettsia sennetsu str. Miyayama — 1 859 006 [https://www.ncbi.nlm.nih.gov/nuccore/NC_007798.1?report=fasta NC_007798] *Ccr Caulobacter crescentus CB15 — 1 4 016 947 [https://www.ncbi.nlm.nih.gov/nuccore/NC_002696.1?report=fasta NC_002696] [http://saba.ibcg.biotoul.fr/soap-data/Atelier_Phylogenomique/Wolbachia/Outgroup_strains Outgroup_strains] --> <!-- ===Autres ressources=== ====PATRIC==== PATRIC est une ressource importante pour les bactéries patogènes: " The Pathosystems Resource Integration Center ([https://www.patricbrc.org/ PATRIC]): the bacterial Bioinformatics Resource Center ([https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkw1017 Wattam et al., 2017])". ====Joint Genome Institute==== [https://jgi.doe.gov/ JGI] is a DOE Office of Science User Facility managed by Lawrence Berkeley National Laboratory© 1997-2017 The Regents of the University of California. ====EnsemblBacteria==== [http://bacteria.ensembl.org/index.html Ensembl Bacteria] is a browser for bacterial and archaeal genomes. Search for a genome: Prochlorococcus 20 entries Downloads, Filter Prochlorococcus 16 entries --> <!-- ==Genome resources== The Pathosystems Resource Integration Center ([https://www.patricbrc.org/ PATRIC])is the bacterial Bioinformatics Resource Center ([https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkw1017 Wattam et al., 2017]). All Data Types: wolbachia Genomes: 46 [https://www.patricbrc.org/view/GenomeList/?keyword(wolbachia)#view_tab=genomes Genome List View] ====Download genomes==== Select all genomes (in Genome Name) Use DWNLD with More Options to download files in FASTA format. Genomes to add: *''Wolbachia'' endosymbiont of ''Dirofilaria immitis'' wDi.2.2 [http://dirofilaria.nematod.es/ http://dirofilaria.nematod.es/] *''Wolbachia'' endosymbiont of ''Litomosoides sigmodontis'' wLs.2.0 [http://litomosoides.nematod.es/ http://litomosoides.nematod.es/] Outgroup strains *Ace Anaplasma centrale str. Israel — 1 1 206 806 NC_013532 *Aph Anaplasma phagocytophilum HZ — 1 1 471 282 NC_007797 *Ech Ehrlichia chaffeensis str. Arkansas — 1 1 176 248 NC_007799 *Eru Ehrlichia ruminantium str. Gardel — 1 1 499 920 NC_006831 *Nri Neorickettsia risticii str. Illinois — 1 879 977 NC_013009 *Nse Neorickettsia sennetsu str. Miyayama — 1 859 006 NC_007798 *Ccr Caulobacter crescentus CB15 — 1 4 016 947 NC_002696 ===Joint Genome Institute=== [https://jgi.doe.gov/ JGI] is a DOE Office of Science User Facility managed by Lawrence Berkeley National Laboratory© 1997-2017 The Regents of the University of California. ===EnsemblBacteria=== [http://bacteria.ensembl.org/index.html Ensembl Bacteria] is a browser for bacterial and archaeal genomes. Search for a genome: wolbachia 20 entries Downloads, Filter wolbachia 16 entries --> ==Annotations des génomes== <!-- Nous allons annoter tous les génomes avec la même stratégie pour avoir une annotation homogène. ===PROKKA=== "Prokka is a software tool for the rapid annotation of prokaryotic genomes": [http://www.vicbioinformatics.com/software.prokka.shtml prokka home]. Citation: Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014 Jul 15;30(14):2068-9. [http://www.ncbi.nlm.nih.gov/pubmed/24642063 PMID:24642063] Nous allons utiliser la version prokka 1.10 (current) qui est disponible sur le serveur. ====Exemples d'utilisations==== ====''Dirofilaria immitis wolbachia''==== --outdir répertoire avec les fichiers résultats. --compliant forcer la compatibilité avec GenBank/EMBL Créer un répertoire : mkdir /home/<user_name>/work/wolbachia cd /home/<user_name>/work/wolbachia en suivant les recommendations du fichier How_to_use (head /usr/local/bioinfo/src/PROKKA/How_to_use) module load bioinfo/prokka-1.10 car la version current correspond à prokka-1.10. Exécuter le programme à l'aide de qsub (fichier commande) ou connectez-vous sur un nœud du cluster avec qlogin. <pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> /usr/local/bioinfo/src/PROKKA/current/bin/prokka --force --outdir prokka/wDim --addgenes --prefix wDim --locustag wDim --species 'Wolbachia endosymbiont' --strain wDim --compliant /home/formation/public_html/M2_Phylogenomique/data/Wolbachia/Additional/Dirofilaria_immitis_wolbachia_2.2.fna
Output Files
.gff This is the master annotation in GFF3 format, containing both sequences and annotations. It can be viewed directly in Artemis or IGV. .gbf This is a standard Genbank file derived from the master .gff. If the input to prokka was a multi-FASTA, then this will be a multi-Genbank, with one record for each sequence. .fna Nucleotide FASTA file of the input contig sequences. .faa Protein FASTA file of the translated CDS sequences. .ffn Nucleotide FASTA file of all the prediction transcripts (CDS, rRNA, tRNA, tmRNA, misc_RNA) .sqn An ASN1 format "Sequin" file for submission to Genbank. It needs to be edited to set the correct taxonomy, authors, related publication etc. .fsa Nucleotide FASTA file of the input contig sequences, used by "tbl2asn" to create the .sqn file. It is mostly the same as the .fna file, but with extra Sequin tags in the sequence description lines. .tbl Feature Table file, used by "tbl2asn" to create the .sqn file. .err Unacceptable annotations - the NCBI discrepancy report. .log Contains all the output that Prokka produced during its run. This is a record of what settings you used, even if the --quiet option was enabled. .txt Statistics relating to the annotated features found. .tsv Tab-separated file of all features: locus_tag,ftype,gene,EC_number,product
more prokka/wDim/wDim.txt
organism: Genus wolbachia endosymbiont wDim contigs: 2 bases: 921012 rRNA: 3 gene: 793 CDS: 755 tRNA: 34 tmRNA: 1
Litomosoides sigmodontis wolbachia
/usr/local/bioinfo/src/PROKKA/current/bin/prokka --force --outdir prokka/wLsi --addgenes --prefix wLsi --locustag wLsi --species 'Wolbachia endosymbiont' --strain wLsi --compliant /home/formation/public_html/M2_Phylogenomique/data/Wolbachia/Additional/Litomosoides_sigmodontis_wolbachia_2.0.fna
more prokka/wLsi/wLsi.txt
organism: Genus wolbachia endosymbiont wLsi contigs: 10 bases: 1048936 rRNA: 3 gene: 918 CDS: 880 tRNA: 34 tmRNA: 1
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.
Transférer les fichiers.
art prokka/wDim/wDim.gff
Comparaison des génomes
Nous allons utiliser le logiciel act (documentation).
Il est nécessaire de réaliser au préalable une comparison des deux génome (BLAST ou MEGABLAST).
Exemple d'utilisation de blastall
formatdb -i prokka/wDim/wDim.fna -p F
Créer un nouveau répertoire: act
blastall -p blastn -i prokka/wLsi/wLsi.fna -d prokka/wDim/wDim.fna -m 8 -o act/wLsi-wDim.blt
Copier les fichiers sur votre porte de travail en P0 et lancer :
act prokka/wLsi/wLsi.gbf act/wLsi-wDim.blt prokka/wDim/wDim.gbf&
Vous pouvez aussi utiliser les fichiers en format gff.
Commentez les résultats!
Compatibilité avec MAUVE
Il peut y avoir un problème lors du chargement des fichiers gbk dans MAUVE. Une solution à été proposée:
"Mauve uses BioJava to parse GenBank files, and it is very picky about Genbank files. It does not like long contig names, like those from Velvet or Spades. One solution is to use --centre XXX in Prokka and it will rename all your contigs to be NCBI (and Mauve) compliant. It does not like the ACCESSION and VERSION strings that Prokka produces via the "tbl2asn" tool. The following Unix command will fix them:"
egrep -v '^(ACCESSION|VERSION)' prokka.gbk > mauve.gbk
Annotation des souches
Nous allons retenir les souches dont le génome est 'Complete Genome' ou 'Scaffold'.
Un script perl va automatiser la procédure:
/home/formation/public_html/M2_Phylogenomique/scripts/rename_ncbi_files.pl --prokka_dir prokka --verbose 0
Les fichiers fasta sont dans /home/<user_name>/work/wolbachia/prokka/
Le script suivant permet de lancer prokka sur un ensemble de génomes en utilisant qsub.
/home/formation/public_html/M2_Phylogenomique/scripts/prokka_loop.pl --directory /home/formation/public_html/M2_Phylogenomique/data/Wolbachia/NCBI --prokka_dir prokka --genome_list "wMel wMun" --verbose 0 --erase 0
16 rRNA
Vérifier la présence d'une annotation pour le gène d'ARNr 16S dans les différents génomes.
Extraire la séquence des gènes d'ARNr 16S
Exemple pour le génome wAu, le nom du gène est wAu_01218.
perl -ne 'if(/^>(\S+)/){$c=grep{/^$1$/}qw(wAu_01218)}print if $c' prokka/wAu/wAu.ffn
Arbre des ARNr 16S
L'ARNr 16S est très souvent utilisé pour inférer la phylogénies des bactéries. Nous allons utiliser les séquences extraites pour réaliser un arbre souches sur nos données.
La (les) méthode(s) à mettre en oeuvre est (sont) laissée(s) à votre choix!
MLST
Nous allons utiliser la base de données MLST : Wolbachia PubMLST database (wolbachia) pour extraire les gènes MLST utilisés pour classer les souches de Wolbachia.
MLST genes
Les fichiers sont dans le répertoire : /home/formation/public_html/M2_Phylogenomique/data/Wolbachia//MLST