Atelier Phylogénomique
From silico.biotoul.fr
Matériel pédagogique
- Présentation 'Quest for Orthologs'.
- Références
- Présentation 'Reconstruction of ancestral traits'.
- Ancestral State Reconstruction with phytools
- Ancestral State Reconstruction with phytools2
Références
- Kettler et al., PLoS Genet. 2007 Dec;3(12):e231 Patterns and implications of gene gain and loss in the evolution of Prochlorococcus.
- Yan et al., Appl Environ Microbiol. 2018 Genome rearrangement shapes Prochlorococcus ecological adaptation.
- Biller et al., Nat. Rev. Microbiol. 2015 13(1) 13-27 Prochlorococcus: the structure and function of collective diversity.
- Chisholm Current Biology 27, R431–R510 Prochlorococcus.
Ressources informatiques
Nous allons utiliser les ressources de GenoToul.
Vous avec dans les FAQ, les réponses à toutes vos questions concernant l'utilisation de la ressource.
Ask for Install
trimal Gblocks PhyML IQ-tree
Sortcuts
Vous allez vous connecter avec un compte fleur:
ssh -Y <login>@genologin.toulouse.inra.fr
Vous pouvez ouvrir deux connections afin de pouvoir lancer qlogin de façon indépendante.
Échange de fichiers:
scp file login@genologin.toulouse.inra.fr:~/work
Logiciels disponibles
ou plus directement
ls /usr/local/bioinfo/src/
La documentation est disponible sur le site WEB et dans le répertoire correspondant au logiciel (fichiers How_to_use and/or Readme).
SLURM cluster softwares
soumission de jobs
1. First write a script (ex: myscript.sh) with the command line as following:
#!/bin/bash #SBATCH -J test #SBATCH -o output.out #SBATCH -e error.out #SBATCH -t 01:00:00 #SBATCH --mem=8G #SBATCH --mail-type=BEGIN,END,FAIL (the email address is automatically LDAP account's one) #Purge any previous modules module purge #Load the application module load bioinfo/ncbi-blast-2.2.29+ # My command lines I want to run on the cluster blastall ...
to submit the job, use the sbatch command line as following
- sbatch (qsub): submit a batch job to slurm (default workq partition()
- sbatch myscript.sh
- sarray (qarray): submit a batch job-array to slurm
- scancel (qdel): kill the specified job
INTERACTIVE
- srun [--pty bash] (qrsh): submit an interactive session with a compute node (default workq partition).
- runVisuSession.sh (qlogin): submit a TurboVNC / VirtualGL session with the graphical node (interq partition). Just for graphics jobs.
Pour controler les jobs
- sinfo (qhost): display nodes, partitions, reservations
- squeue (qstat): display jobs and state
- scontrol show : get informations on jobs, nodes, partitions
- sstat (qstat -j): show status of running jobs
- sview (qmon): graphical user interface
- sacct (qacct) : display accounting data
Utiliser R sur le cluster
Tutoriel expliquant comment utiliser R (et compiler des fichiers Rmd) sur le cluster (slurm) de la PF Bioinformatique de Toulouse (Gaëlle Lefort et Nathalie Vialaneix):
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
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 pair 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
genoplotR
Nous allons utiliser genoplotR pour visualiser les similarités entre les pairs de génomes.
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 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
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-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.8-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 porthomcl/0.input_faa/*.faa > panoct/peptides 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
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/matchtable2panoct.result.pl --results_dir results /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 indiquez sur la figure, le génome coeur, le génome accessoire, les singletons et la pan génome commentez les résultats obtenus. Quels sont les choix méthodologiques qui peuvent influencer la taille du génome coeur?
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
align_clusters.pl - create a multiple alignment file of the members in each cluster
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 des 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
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 porthomcl/0.input_faa/*.faa > panoct/peptides 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 /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
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 ~/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 /home/yquentin/work/ProchlorococcusSynechococcus/OG/alignment /home/formation/public_html/M2_Phylogenomique/scripts/aa_to_dna_aln.pl -dna /home/yquentin/work/ProchlorococcusSynechococcus/OG/DNA/OG_1471.fas --outdir /home/yquentin/work/ProchlorococcusSynechococcus/OG/alignment
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 /home/yquentin/work/ProchlorococcusSynechococcus/OG/DNA --outdir /home/yquentin/work/ProchlorococcusSynechococcus/OG/alignment --sample 10
Construction des arbres des groupes de gènes orthologues
Trimer les alignements avec trimal (pas sur genologin!)
/home/formation/public_html/M2_Phylogenomique/scripts/trimal_conservation.pl --alignment ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_1063.fas --erase 1
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"
Transformation en alignement multiples des séquences ADN des gènes
Créer un fichier avec les séquences nucléotidiques:
mkdir ~/work/ProchlorococcusSynechococcus/DNA cat ~/work/*/prokka/Aaa*/Aaa*.ffn > ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn grep -c '>' ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn
/home/formation/public_html/M2_Phylogenomique/scripts/muscle_loop.pl --directory /home/yquentin/work/ProchlorococcusSynechococcus/peptides_OG --sample 10 --outdir /home/yquentin/work/ProchlorococcusSynechococcus/OG/peptide_alignments
mkdir ~/work/ProchlorococcusSynechococcus/DNA_OG /home/formation/public_html/M2_Phylogenomique/scripts/index_fasta.pl --fasta ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn --matchtable ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt --outdir ~/work/ProchlorococcusSynechococcus/DNA_OG grep -c '>' /home/yquentin/work/ProchlorococcusSynechococcus/DNA_OG/*fas
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"
old stuff
Autres ressources
PATRIC
PATRIC est une ressource importante pour les bactéries patogènes: " The Pathosystems Resource Integration Center (PATRIC): the bacterial Bioinformatics Resource Center (Wattam et al., 2017)".
Joint Genome Institute
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
Ensembl Bacteria is a browser for bacterial and archaeal genomes.
Search for a genome: Prochlorococcus 20 entries Downloads, Filter Prochlorococcus 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": prokka home.
Citation: Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014 Jul 15;30(14):2068-9. 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.
/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
HMM profiles
Pour annoter les gènes MLST dans nos génomes, nous allons utiliser hmmer. La première étape est de construire des fichiers profiles pour les 5 gènes.
/home/formation/public_html/M2_Phylogenomique/scripts/fasta2hmm.pl --mlst_dir MLST --gene_list "gatB coxA hcpA ftsZ fbpA" --verbose 0 --erase 0
Vérifier que les fichiers MLST/*.hmm ne sont pas vides!
Identification des gènes dans les génomes
Nous allons extraire les gènes MLST et les concaténer pour chaque génome.
Annotation de chaque profile HMM
/usr/local/bioinfo/bin/nhmmer -E 1e-50 -o MLST/gatB_wAU.out -A MLST/gatB_wAu.stk --noali MLST/gatB.hmm prokka/wAu/wAu.fna;
Le fichier avec les alignements est en format STK (STOCKHOLM), nous devons le transformer en format FASTA
head MLST/gatB_wAu.stk
STK à FASTA
/usr/local/bioinfo/bin/esl-reformat --informat stockholm -o MLST/gatB_wAu.fa fasta MLST/gatB_wAu.stk;
Automatisation
Les deux étapes sont automatisées pour tous les gènes et tous les génomes.
/home/formation/public_html/M2_Phylogenomique/scripts/hmm2genes.pl --prokka_dir prokka --mlst_dir MLST --gene_list "gatB coxA hcpA ftsZ fbpA" --verbose 0 --erase 0
Les gènes concaténés sont dans le fichier :
MLST/concat_file.fa
Alignement et arbre
muscle:
/usr/local/bioinfo/bin/muscle -in MLST/concat_file.fa -out MLST/concat_file.mfa
readal (fasta à phylip):
/usr/local/bioinfo/bin/readal -in MLST/concat_file.mfa -out MLST/concat_file.phy -phylip;
PhyML:
/usr/local/bioinfo/src/PhyML/PhyML-3.1/PhyML-3.1_linux64 -i MLST/concat_file.phy -d nt -b -4 -m HKY85 -f m -t e -v e -c 4 -a e -s NNI -o tlr --quiet --no_memory_check
Visualisation de l'arbre
Lancer l'interface graphique sur votre poste de travail de la P0.
seaview MLST/concat_file.phy_phyml_tree.txt&
Le fichier ST_concatenated_reference.fas contient des séquence type de référence. Elles peuvent être ajoutées à votre alignement.
Comparer les arbres obtenus avec l'ARNr 16S et les gènes MLST.
Le fichier NCBI_Wolbachia.csv contient une colonne Clade_ID. Comparer cette annotation avec la topologie obtenue.
Groupes de gènes orthologues
OrthoMCL
Le guide utilisateurs sur le serveur: /usr/local/bioinfo/src/OrthoMCL/current/doc/OrthoMCLEngine/Main/UserGuide.txt
For alternate documentation online, please read Unit 6.12 of the Current Protocols of Bioinformatics available at:
OrthoMCL Protocols
For details on the orthomcl algorithm, please read the OrthoMCL Algorithm Document:
OrthoMCL Algorithm Document
En entrée, le programme prend un ensemble de protéomes.
La sortie de orthoMCL est un ensemble de fichiers :
pairs/ potentialOrthologs.txt potentialCoorthologs.txt potentialInparalogs.txt groups.txt
Les fichiers dans le répertoire pairs contiennent les paires de relations entre les protéines associées aux scores. Comme décrit dans l'article, on distingue :
- les orthologues potentiels
- les co-orthologues
- les inparalogues
Il y a trois grandes étapes:
- blastp tous_contre_tous
- l'identification des paires (répertoire pairs)
- le partitionnement du graphe des pairs avec le programme MCL program (groups.txt).
Les différentes étapes du l'analyse sont lancées indépendamment mais elles s'enchainent dans un ordre bien précis. Il est possible de modifier une étape pour tester des alternatives ou optimiser l'exécution. Toutefois, il est impératif de respecter les formats des fichiers.
à modifier!
orthoMCL pipe 1
Pour faciliter les premières étapes (Step 5 et 6), elles ont été intégrées dans le script prep_orthoMLC.pl.
/home/formation/public_html/M2_Phylogenomique/scripts/prep_orthoMCL.pl
Step 5: orthomclAdjustFasta
Input:
- fasta files as acquired from the genome resource.
Output:
- the my_orthomcl_dir/compliantFasta/ directory of orthomcl-compliant fasta files (see Step 6)
Use orthomclAdjustFasta to produce a compliant file from any input file that conforms to the following pattern (for other files, provide your own script to produce complaint fasta files):
- has one or more fields that are separated by white space or the '|' character (optionally surrounded by white space) - has the unique ID in the same field of every protein.
Create an empty my_orthomcl_dir/compliantFasta/ directory, and change to that directory.
Run orthomclAdjustFasta once for each input proteome fasta file. It will produce a compliant file in the new directory. Check each file to ensure that the proteins all have proper IDs.
Step 6: orthomclFilterFasta
Input:
- my_orthomcl_dir/compliantFasta/ - optionally a gene->protein mapping file
Output:
- my_orthomcl_dir/goodProteins.fasta - my_orthomcl_dir/poorProteins.fasta - report of suspicious proteomes (> 10% poor proteins)
This step produces a single goodProteins.fasta file to run BLAST on. It filters away poor-quality sequences (placing them in poorProteins.fasta). The filter is based on length and percent stop codons. You can adjust these values.
The input requirements are:
# a compliantFasta/ directory which contains all and only the proteome .fasta files, one file per proteome. # each .fasta file must have a name in the form 'xxxx.fasta' where xxxx is a three or four letter unique taxon code. For example: hsa.fasta or eco.fasta # each protein in those files must have a definition line in the following format: >xxxx|yyyyyyyy
where xxxx is the three or four letter taxon code and yyyyyyy is a sequence identifier unique within that taxon.
Change dir to my_orthomcl_dir/ and run orthomclFilterFasta.
orthoMCL pipe 2
Les blastp tous-contre-tous (Step 7) sont réalisé savec le script ava_BLAST_orthoMCL.pl.
Pour éviter de charger le cluster avec les blastp, nous allons les réaliser séparément pour chaque génome. Chaque compte fleur va lancer les blastp sur un ou deux genomes comme query et tous les autres génomes comme data base.
Association fleur/génome
anemone wAu wNo arome wCle aster wDim bleuet wHa camelia wMel capucine wNfe wMeP chardon wPiM clematite wViB cobee w Cam coquelicot wDci cyclamen wGmo wPpe dahlia wMun digitale wObr gerbera wRi geranium wBma wWb glaieul wCqu hortensia wLsi
Paramètres:
--query_list "wMun" --target_list "wAu wCle wDim wHa wMel wNfe wPiM wViB wCam wDci wGmo wMun wObr wRi wBma wCqu wLsi wMeP wNo wPpe wWb"
Copier les résultats dans un répertoire accessible à tous.
A la fin, concaténer tous les résultats blasp dans un seul fichier:
cat wolbachia/blastp/*fasta > wolbachia/blastp/wol_21_blt
Step 7: All-v-all BLAST
Input:
- goodProteins.fasta
Output:
- your_blast_results_in_tab_format
You must run your own BLAST.
We expect you to:
- use NCBI BLAST - run with the -m 8 option to provide tab delimited output required by Step 8 - for IMPORTANT DETAILS about other BLAST arguments, see: the OrthoMCL Algorithm Document
Blastp parameters:
-F 'm S': mask with Seg -v 100000: a "don't care" value -b 100000: a "don't care" value -z protein_database_size: the number of proteins in the set. S -e 1e-5: recommended e-value cutoff
Step 8: orthomclBlastParser
Input:
- your_blast_results_in_tab_format - my_orthomcl_dir/compliantFasta/
Output:
- my_orthomcl_dir/similarSequences.txt
This step parses NCBI BLAST -m 8 output into the format that can be loaded into the orthomcl database.
Use the orthomclBlastParser program for this. In addition to formatting, it computes the percent match of each hit.
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclBlastParser wolbachia//blastp/wol_21_blt wolbachia/my_orthomcl_dir/compliantFasta >> wolbachia/my_orthomcl_dir/similarSequences.txt
head wolbachia/my_orthomcl_dir/similarSequences.txt
orthoMCL mySQL DB
OrthoMCL nécessite une base de données (mySQL par exemple). Elle a été crée par Marie-Stephane Trotard (support genotoul)
Copier le fichier modèle
/usr/local/bioinfo/src/OrthoMCL/orthomclSoftware-v2.0.9/my_orthomcl_dir_template/orthomcl.config
sous votre work par exemple et personnaliser le ainsi
Pour créer le schema de la base:
/usr/local/bioinfo/src/OrthoMCL/orthomclSoftware-v2.0.9/bin/orthomclInstallSchema <path_vers_votre_orthomcl.config <path_verts_votre repertoire de travail>/install_schema.log
/usr/local/bioinfo/src/OrthoMCL/orthomclSoftware-v2.0.9/bin/orthomclInstallSchema wolbachia/my_orthomcl_dir/orthomcl.config wolbachia/my_orthomcl_dir/install_schema.log
DROP Tables
OrthoMCL expects a "clean" and empty database after orthoMCLInstallSchema, so in between runs it's easiest to run
mysql -h 147.99.108.14 -D orthomcl_yquentin -u ortho_yquentin -p
mysql>SHOW TABLES;
mysql>DROP TABLE IF EXISTS BestHit, BestInterTaxonScore, BestQueryTaxonScore, BetterHit, CoOrthNotOrtholog, CoOrtholog, CoOrthologAvgScore, CoOrthologCandidate, CoOrthologTaxon, CoOrthologTemp, InParalog, InParalog2Way, InParalogAvgScore, InParalogOrtholog, InParalogTaxonAvg, InParalogTemp, InplgOrthTaxonAvg, InplgOrthoInplg, InterTaxonMatch, Ortholog, Ortholog2Way, OrthologAvgScore, OrthologTaxon, OrthologTemp, OrthologUniqueId, SimilarSequences, UniqSimSeqsQueryId;
inside the mysql command line and then to do orthoMCLinstallSchema again.
Step 9: orthomclLoadBlast
Input:
- similarSequences.txt
Output:
- SimilarSequences table in the database
This step loads the BLAST results into the orthomcl database.
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclLoadBlast wolbachia/my_orthomcl_dir/orthomcl.config work/wolbachia/my_orthomcl_dir/similarSequences.txt
Step 10: orthomclPairs
Input:
- SimilarSequences table in the database
Output:
- PotentialOrthologs table - PotentialInParalogs table - PotentialCoOrthologs table
This is a computationally major step that finds protein pairs. It executes the algorithm described in the OrthoMCL Algorithm Document, using a relational database.
The program proceeds through a series of internal steps, each creating an intermediate database table or index. There are about 20 such tables created. Finally, it populates the output tables.
The cleanup= option allows you to control the cleaning up of the intermediary tables. T
- 'yes' option drops the intermediary tables once they are no longer needed.
- 'no' option keeps the intermediary tables in the database.
In total, they are expected to be about 50 percent of the SimilarSequences table. They are useful mostly for power users or developers who would like to query them. They can be removed afterwards with the 'only' or 'all' options. The latter also removes the final tables, and should only be done after Step 11 below has dumped them to files.
The startAfter= option allows you to pick up where the program left off, if it stops for any reason. Look in the log to find the last completed step, and use its tag as the value for startAfter=
Because this program will run for many hours, we recommend you run it using the UNIX 'screen' program, so that it does not abort in the middle. (If it does, use startAfter=).
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclPairs wolbachia/my_orthomcl_dir/orthomcl.config wolbachia/my_orthomcl_dir/orthomcl_pairs.log cleanup=no
Step 11: orthomclDumpPairsFiles
Input:
- database with populated pairs tables
Output
- pairs/ directory. - mclInput file
Run the orthomclDumpPairsFiles.
The pairs/ directory contains three files: ortholog.txt, coortholog.txt, inparalog.txt. Each of these has three columns:
- protein A - protein B - their normalized score (See the Orthomcl Algorithm Document).
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclDumpPairsFiles wolbachia/my_orthomcl_dir/orthomcl.config
head wolbachia/my_orthomcl_dir/pairs/coorthologs.txt head wolbachia_save/my_orthomcl_dir/pairs/inparalogs.txt head wolbachia_save/my_orthomcl_dir/pairs/orthologs.txt
Step 12: mcl
Input:
- mclInput file
Output:
- mclOutput file
/usr/local/bioinfo/bin/mcl wolbachia/my_orthomcl_dir/mclInput --abc -I 1.5 -o work/wolbachia/my_orthomcl_dir/mclOutput
Step 13: orthomclMclToGroups
Input:
- mclOutput file
Output:
- groups.txt
Change to my_orthomcl_dir and run:
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclMclToGroups wolb 1000 < mclOutput > groups.txt
- my_prefix is an arbitrary string to use as a prefix for your group IDs.
- 1000 is an arbitrary starting point for your group IDs.
Step 14: orthomclSingletons
Input:
- groups.txt - goodProteins.txt
Output:
- singletons.txt
Change to my_orthomcl_dir and run:
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclSingletons wolbachia/my_orthomcl_dir/good_proteins_file wolbachia/my_orthomcl_dir/groups.txt >> wolbachia/my_orthomcl_dir/singletons.txt
wc -l wolbachia/my_orthomcl_dir/singletons.txt 1639
Other utilities
Extract protein IDs from an orthomcl groups file.
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclExtractProteinIdsFromGroupsFile wolbachia/my_orthomcl_dir/groups.txt
Extract protein ID pairss from an orthomcl groups file.
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclExtractProteinPairsFromGroupsFile wolbachia/my_orthomcl_dir/groups.txt
</pre>
Perl scripts
PROKKA
- rename_ncbi_files.pl
- prokka_loop.pl
MLST
- fasta2hmm.pl
- hmm2genes.pl
orthMCL
- prep_orthoMCL.pl