Atelier Phylogénomique

From silico.biotoul.fr

(Difference between revisions)

Jump to: navigation, search

Revision as of 07:51, 3 October 2018

1 Matériel pédagogique
2 Références
3 Ressources informatiques
4 Disponibilité des génomes
5 Annotation
- 5.1 Prokka
6 Séquence conservation entre souches de Prochlorococcus
- 6.1 Genome pairs
  - 6.1.1 genoplotR
  - 6.1.2 ACT
7 Groupes de gènes orthologues
8 old stuff
- 8.1 Autres ressources
9 Annotations des génomes
10 Annotation des souches
11 16 rRNA
- 11.1 Extraire la séquence des gènes d'ARNr 16S
- 11.2 Arbre des ARNr 16S
12 MLST
13 Groupes de gènes orthologues
14 Perl scripts

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.

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

software_faq

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

softwares

soumission de jobs

Seminar_cluster_SLURM

1. First write a script (ex: myscript.sh) with the command line as following:

#!/bin/bash
#SBATCH -J test
#SBATCH -o output.out
#SBATCH -e error.out
#SBATCH -t 01:00:00
#SBATCH --mem=8G
#SBATCH --mail-type=BEGIN,END,FAIL (the email address is automatically LDAP account's one)
#Purge any previous modules
module purge

#Load the application
module load bioinfo/ncbi-blast-2.2.29+

# My command lines I want to run on the cluster
blastall ...

to submit the job, use the sbatch command line as following

sbatch (qsub): submit a batch job to slurm (default workq partition()
sbatch myscript.sh
sarray (qarray): submit a batch job-array to slurm

scancel (qdel): kill the specified job

INTERACTIVE

srun [--pty bash] (qrsh): submit an interactive session with a compute node (default workq partition).
runVisuSession.sh (qlogin): submit a TurboVNC / VirtualGL session with the graphical node (interq partition). Just for graphics jobs.

Pour controler les jobs

sinfo (qhost): display nodes, partitions, reservations
squeue (qstat): display jobs and state
scontrol show : get informations on jobs, nodes, partitions
sstat (qstat -j): show status of running jobs
sview (qmon): graphical user interface
sacct (qacct) : display accounting data

Utiliser R sur le cluster

Tutoriel expliquant comment utiliser R (et compiler des fichiers Rmd) sur le cluster (slurm) de la PF Bioinformatique de Toulouse (Gaëlle Lefort et Nathalie Vialaneix):

tutoR_cluster.pdf

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.

prokka files.

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

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!

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

How to use SLURM MCL

Location: /usr/local/bioinfo/src/MCL

srun --pty bash 
#Load modules
module load bioinfo/mcl-14-137

cat 6.orthologs.-5.50/*.tsv >> 8.all.ort.-5.50.tsv

mcl 8.all.ort.-5.50.tsv --abc -I 1.5 -t 4 -o 8.all.ort.group.1.5-5.50

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

Pan-genome Ortholog Clustering Tool, is a program written in PERL for pan-genomic analysis of closely related prokaryotic species or strains. Unlike traditional graph-based ortholog detection programs, it uses micro synteny or conserved gene neighborhood (CGN) in addition to homology to accurately place proteins into orthologous clusters.

panoct project

Vous trouverez le package dans le repertoire suivant:

/home/formation/public_html/M2_Phylogenomique/PanOCT/panoct_v3.23/

Test:

cd /home/formation/public_html/M2_Phylogenomique/PanOCT/panoct_v3.23/example_dir
../bin/panoct.pl -b ../example_dir -t example_blast.txt -.pep -S Y -L 1 -M Y -H Y -V Y -N Y -F 1.33 -G y -c 0,25,50,75,100 -T

Préparation des fichiers d'entrée de PanOCT

srun --pty bash
mkdir -p ~/work/Prochlorococcus/panoct/results

gene.att

Créer un fichier avec les coordonnées, noms, fonction et souches des gènes.

cd ~/work/Prochlorococcus
for file in peptide/*.faa
do
 prefix=$(basename $file .faa)
 echo "/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.

Definitions

D'après Vernikos et al. 2015

Le pan-génome a été défini pour la première fois par dTettelin 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.

An open pan-genome indicates that a very large, undetermined number of additional genomes would be needed to identify all genes accessible to the species.

A closed pan-genome, indicates that additional genomes sequenced do not provide additional new genes to expand the pan-genome — the entire gene repertoire has been characterized, assuming that the sampling of strains sequenced is not biased. These observations are derived from extrapolations based on the current sample of bacterial genomes analyzed.

paralog_matchtable.pl

Takes a PanOCT matchtable_0_1 and a paralogs file (output from PanOCT) and outputs a parlogs matchtable where instead of 0 or 1 for presence or absence of the cluster a number is given indicating the number of 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

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

Afficher les courbes:

R
pangenome_size <- 'results/pangenome_size'
pdf_file <- '../images/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")

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

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

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

P_marinus_meta

@@ Line 618: / Line 618: @@
 =====Restrict view to OG without paralogs=====
 <pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
-mkdir ~/work/Prochlorococcus/images
 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

Atelier Phylogénomique

From silico.biotoul.fr

Revision as of 07:51, 3 October 2018

Contents

Matériel pédagogique

Références

Ressources informatiques

Sortcuts

SLURM cluster softwares

soumission de jobs

to submit the job, use the sbatch command line as following

INTERACTIVE

Pour controler les jobs

Utiliser R sur le cluster

Logiciels à installer sur vos postes de travail

Scripts

Disponibilité des génomes

Prochlorococcus

NCBI

GenBank

DNA

Annotation

Prokka

Exemple d'utilisation

Automatisation des annotations prokka sur l'ensemble des génomes

Visualisation des annotations

Séquence conservation entre souches de Prochlorococcus

Genome pairs

genoplotR

ACT

Groupes de gènes orthologues

Préliminaires

Blast All-All

make blast database

Intra genomes

Inter genomes

PorthoMCL

old stuff

Running MCL

Visualisation of OG

Global view

Restrict view to OG without paralogs

Restrict view to one OG

Tester l'effet de l'IF surla taille des OG

Tester l'effet de l'IF sur la paralogie

old stuff2

PanOCT

Préparation des fichiers d'entrée de PanOCT

gene.att

tags.txt

peptides

blast.txt

run panOCT

panoct.pl

Analyses pan-génomiques

Definitions

paralog_matchtable.pl

compute_pangenome.R

run_panoct.pl

compute_pangenome.R

old stuff

Autres ressources

PATRIC

Joint Genome Institute

EnsemblBacteria

Annotations des génomes

PROKKA

Exemples d'utilisations

Dirofilaria immitis wolbachia

Litomosoides sigmodontis wolbachia

Visualisation des annotations

Comparaison des génomes

Compatibilité avec MAUVE

Annotation des souches

16 rRNA

Extraire la séquence des gènes d'ARNr 16S

Arbre des ARNr 16S

MLST

MLST genes

HMM profiles