Atelier Phylogénomique

From silico.biotoul.fr

Revision as of 09:53, 12 October 2018 by Quentin (Talk | contribs)

(diff) ← Older revision | Current revision (diff) | Newer revision → (diff)

Jump to: navigation, search

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 Analyses pan-génomiques
9 Analyses phylogénomiques
10 Reconstruction d'états ancestraux
11 old stuff
12 Annotations des génomes
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.

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

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

Installation du package genoPlotR

srun --pty bash
module load system/R-3.5.1
R
install.packages('genoPlotR')
Warning in install.packages("genoPlotR") :
  'lib = "/tools/R/R-3.5.1/lib64/R/library"' ne peut être ouvert en écriture
Would you like to use a personal library instead? (yes/No/cancel) yes
Would you like to create a personal library
‘~/R/x86_64-pc-linux-gnu-library/3.5’
to install packages into? (yes/No/cancel) yes
--- SVP sélectionner un miroir CRAN pour cette session ---
Secure CRAN mirrors
29: France (Lyon 2) [https]

* installing *source* package ‘ade4’ ...
* installing *source* package ‘genoPlotR’ ...

Les packages source téléchargés sont dans
	‘/tmp/RtmpLdNj99/downloaded_packages’
library(genoPlotR)

Logiciels à installer sur vos postes de travail

seaview : Multiplatform GUI for molecular phylogeny
mauve : Multiple genome alignment
Artemis : Genome browser and annotation tool
Artemis Comparison Tool : Java application for displaying pairwise comparisons between two or more DNA sequences

Scripts

Des scripts perl ont été écrits pour faciliter certaines étapes du TP. Ils sont disponibles dans le répertoire:

/home/formation/public_html/M2_Phylogenomique/scripts

Vous pouvez les copier dans votre espace de travail et les modifier à votre convenance.

Disponibilité des génomes

Caractéristiques des souches étudiées

Table 1. General Characteristics of the Prochlorococcus and Synechococcus Isolates Used in ([https://www.ncbi.nlm.nih.gov/pubmed/18159947 Kettler et al., 2007).

Cyanobacterium	Isolate	 Light  Length(bp)      GC%     Number  Isol.   Region          Date           Accession
                         Adap.                          Genes   Depth

Prochlorococcus	MED4	 	HL(I)	1,657,990	30.8	1,929	5m	Med. Sea	Jan. 1989	BX548174
	        MIT9515	 	HL(I)	1,704,176	30.8	1,908	15m	Eq. Pacific	Jun. 1995	CP000552
	        MIT9301	 	HL(II)	1,642,773	31.4	1,907	90m	Sargasso Sea	Jul. 1993	CP000576
	        AS9601	 	HL(II)	1,669,886	31.3	1,926	50m	Arabian Sea	Nov. 1995	CP000551
	        MIT9215	 	HL(II)	1,738,790	31.1	1,989	5m	Eq. Pacific	Oct. 1992	CP000825
	        MIT9312	 	HL(II)	1,709,204	31.2	1,962	135m	Gulf Stream	Jul. 1993	CP000111
	        NATL1A	 	LL(I)	1,864,731	35.1	2,201	30m	N. Atlantic	Apr. 1990	CP000553
	        NATL2A	 	LL(I)	1,842,899	35	2,158	10m	N. Atlantic	Apr. 1990	CP000095
	        CCMP1375/SS120	LL(II)	1,751,080	36.4	1,925	120m	Sargasso Sea	May 1988	AE017126
	        MIT9211	 	LL(III)1,688,963	38	1,855	83m	Eq. Pacific	Apr. 1992	CP000878
	        MIT9303	 	LL(IV)	2,682,807	50.1	3,022	100m	Sargasso Sea	Jul. 1992	CP000554
	        MIT9313	 	LL(IV)	2,410,873	50.7	2,843	135m	Gulf Stream	Jul. 1992	BX548175
Synechococcus	CC9311	 	Syn.	2,606,748	52.5	3017	95m	Calif.		Current 1993    CP000435
	        CC9902	 	Syn.	2,234,828	54.2	2504	5m	Calif. 		Current 1999	CP000097
	        WH8102	 	Syn.	2,434,428	59.4	2787		Sargasso Sea	Mar. 1981	BX548020
	        CC9605	 	Syn.	2,510,659	59.2	2991	51m	Calif.		Current 1996    CP000110

Prochlorococcus

NCBI

Génomes de Prochlorococcus disponibles au NCBI browse

Accessions des génomes utilisés dans Kettler et al., PLoS Genet. 2007:

accession "BX548174,CP000552,CP000576,CP000551,CP000825,CP000111,CP000553,CP000095,AE017126,CP000878,CP000554,BX548175"

Les fichiers GenBank peuvent être obtenus avec la commande wget en utilisant les chemins enregistrés dans le fichier: accession_file.lst

Pour des raisons de compatibilité avec les programmes, les génomes sont renommés en utilisant un code à quatre lettres. accession_labels_file.lst

GenBank

Les fichiers GenBank renommés sont disponibles dans le répertoire: GenBank

DNA

Les séquences ADN ont été extraites des fichiers GenBank : DNA

Annotation

Prokka

Les réplicons des génomes sont annotés avec le logiciel prokka.

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 --sample Prochlorococcus

squeue -l -u <user>

Une fois les jobs terminés, vérifiez que les fichiers de sorties de prokka existent et ne sont pas vides.

ls -l ~/work/Prochlorococcus/prokka/Aaa*/*.faa

Les fichiers avec le suffix .err renferment la sortie standard de prokka. Si tout c'est bien passé, vous pouvez supprimer les fichiers .err et .sh.

Visualisation des annotations

Nous pouvons utiliser le logiciel art (Artemis) pour visualiser les annotations des génomes:

Il est fortement recommandé d'utiliser ce logiciel en local sur votre poste de travail.

runVisuSession.sh
art ~/work/Prochlorococcus/prokka/Aaaa/Aaaa.gff

Séquence conservation entre souches de Prochlorococcus

Genome pairs

Afin d'estimer les conservations entre les différents génomes, nous allons les comparer par paire de génomes dans l'ordre suivant, à l'aide de blastn:

('Aaab', 'Aaag', 'Aaaj', 'Aaaf', 'Aaak', 'Aaae', 'Aaai', 'Aaad', 'Aaaa', 'Aaah', 'Aaal', 'Aaac')

Nous allons utiliser l'option BLAST-2-Sequences de blastn. Exemple avec une pair

srun -n1 -l blastn -query /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaab.fas -subject /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaag.fas -evalue 1e-05 -outfmt 6 -num_threads 1 -out BlastN/Aaab_vs_Aaag.tab

mkdir ~/work/Prochlorococcus/BlastN

/home/formation/public_html/M2_Phylogenomique/scripts/blastn_genome_pair.pl

ls -l ~/work/Prochlorococcus/BlastN

genoplotR

Nous allons utiliser genoplotR pour visualiser les similarités entre les pairs de génomes.

Installation du package genoPlotR

genoplotR nécessite plusieurs objets:

dna_seg: un objet dna_seg est un ensemble de gènes ou d'éléments le long d'un génome, à représenter sur une carte. Nous allons utiliser les fichier en format gbk crées par prokka.
comparison: une comparaison est un ensemble de similitudes, représentant la comparaison entre deux segments d'ADN. Nous allons utiliser les résultats des blastn entre pairs de genomes.
annotation: Un objet d'annotation est utilisé pour annoter un segment d'ADN. Nous ne l'utilisons ici.
tree: un arbre au format Newick qui peut être analysée à l'aide du paquetage ade4. Nous ne l'utilisons plus tard!

mkdir ~/work/Prochlorococcus/images
srun --pty bash
module load system/R-3.5.1
Rscript /home/formation/public_html/M2_Phylogenomique/scripts/genoplot_blastn_links.R

evince ~/work/Prochlorococcus/images/genoplot_blastn_links.pdf

ACT

Il est également possible d'utiliser le logiciel act (documentation).

Copier les fichiers sur votre porte de travail en P0 et lancer :

act prokka/Aaab/Aaab.gbk BlastN/Aaab_vs_Aaag.tab prokka/Aaag/Aaag.gbk

Vous pouvez aussi utiliser les fichiers en format gff.

Groupes de gènes orthologues

Préliminaires

Blast All-All

Nous allons utiliser NCBI_Blast+.

Nous allons copier les fichiers peptides dans un répertoire de travail:

mkdir -p ~/work/Prochlorococcus/peptide
cd ~/work/Prochlorococcus/peptide
cp ~/work/Prochlorococcus/prokka/Aaa*/*.faa ~/work/Prochlorococcus/peptide/.

ls -l ~/work/Prochlorococcus/peptide

make blast database

cd ~/work/Prochlorococcus/peptide
module load bioinfo/ncbi-blast-2.6.0+
srun -n1 -l makeblastdb -in Aaaa.faa -dbtype prot

Boucle sur tous les fichiers

for i in *.faa; 
do      
srun -n1 -l makeblastdb -in $i  -dbtype prot; 
done

Intra genomes

blastp -query {input} -db {input} -seg {SEG} -dbsize 100000000  -evalue {EVALUE} -outfmt 6 -num_threads 1 -out {output}

Un script perl peut réaliser les blastp pour l'ensemble des génomes avec sbatch:

mkdir -p ~/work/Prochlorococcus/BlastP
cd ~/work/Prochlorococcus
/home/formation/public_html/M2_Phylogenomique/scripts/blastp_intra.pl

squeue -l -u <user>

ls BlastP | wc -l

Combien de fichiers attendez-vous?

Inter genomes

Nous allons utiliser un script similaire au précédant, mais avec une double boucle (sur -query et -db).

cd ~/work/Prochlorococcus
/home/formation/public_html/M2_Phylogenomique/scripts/blastp_inter.pl

squeue -l -u <user>

ls BlastP | wc -l

Combien de fichiers attendez-vous?

PorthoMCL

PorthoMCL: Parallel orthology prediction using MCL for the realm of massive genome availability Ehsan Tabari and Zhengchang Su Big Data Analytics 2017 2:4 DOI: 10.1186/s41044-016-0019-8

Il s'inspire d'orthoMCL

Disponible: /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master.

Le logiciel est basé sur l'enchainement de programmes python. Ils peuvent être lancés séparément ou automatiquement à l'aide d'un script shell.

Version automatique (que nous n'utiliserons pas car nous avons déjà réalisé les blastp):

 
mkdir -p ~/work/Prochlorococcus/PorthoMCL_auto12/0.input_faa
cd ~/work/Prochlorococcus/PorthoMCL_auto12
cp /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/prokka/A*/*.faa 0.input_faa/.
srun --nodes=12 --pty bash 
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomcl.sh -t 12 -e 7

Nous allons utiliser une version plus simple:

 
mkdir ~/work/Prochlorococcus/PorthoMCL
cd ~/work/Prochlorococcus/PorthoMCL
srun --pty bash /home/formation/public_html/M2_Phylogenomique/scripts/short_PorthoMCL.sh

old stuff

Running MCL

C'est exactement comme la dernière étape d'OrthoMCL. Il suffit d'exécuter MCL sur les sorties des étapes 6 et 7 :

Notez que -t définit le nombre de processeurs que vous devez utiliser pour MCL.

How to use SLURM MCL

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

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

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

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

Paralogs

cat 7.paralogs.-5.50/*.tsv >> 8.all.par.tsv

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

Visualisation of OG

srun --pty bash
module load system/R-3.5.1

Global view

Rscript /home/formation/public_html/M2_Phylogenomique/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_all.pdf 

evince ~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_all.pdf

Restrict view to OG without paralogs

Rscript /home/formation/public_html/M2_Phylogenomique/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_nopara.pdf --max_paralogs=0

evince ~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_nopara.pdf

Restrict view to one OG

Rscript /home/formation/public_html/M2_Phylogenomique/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_OG5.pdf --select_OG=5

evince ~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_OG5.pdf

Liste des membres du groupe 5

Aaaa.g_00266 Aaaa.g_00937 Aaaa.g_01196
Aaab.g_00239 Aaab.g_00820 Aaab.g_00825
Aaac.g_01231 Aaac.g_01233 Aaac.g_02086
Aaad.g_00743 Aaad.g_00747 Aaad.g_00909
Aaae.g_00242 Aaae.g_00786 Aaae.g_00791
Aaaf.g_00249 Aaaf.g_00782 Aaaf.g_00787
Aaag.g_00262 Aaag.g_00862 Aaag.g_00867
Aaah.g_00255 Aaah.g_00856 Aaah.g_01180
Aaai.g_00757 Aaai.g_00761 Aaai.g_00954
Aaaj.g_00252 Aaaj.g_00786 Aaaj.g_00791
Aaak.g_00251 Aaak.g_00834 Aaak.g_00839
Aaal.g_00941 Aaal.g_00943 Aaal.g_02214

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

Lancer MCL avec différents IF (respecter la syntaxe des noms!)

Le motif suivant est attendu: 8.all.ort.group.<IF>-5.50.

Rscript /home/formation/public_html/M2_Phylogenomique/scripts/MCL_partition_inflate.R --MCL_dir=~/work/Prochlorococcus/PorthoMCL --nb_strains=12 pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_IF.pdf

evince ~/work/Prochlorococcus/images/8.all.ort.group_IF.pdf

Tester l'effet de l'IF sur la paralogie

Rscript /home/formation/public_html/M2_Phylogenomique/scripts/MCL_OG_paralogy.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.1.5-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_paralogs.pdf
Rscript /home/formation/public_html/M2_Phylogenomique/scripts/MCL_OG_paralogy.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.8.0-5.50 --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group.8.0-5.50_paralogs.pdf

old stuff2

PanOCT with Prochlorococcus

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

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 peptide/*.faa > panoct/results/combined.fasta

head panoct/results/combined.fasta

blast.txt

Concaténer les résultats des blastp dans un seul fichier.

cat BlastP/*.tab > panoct/results/combined.blast

head panoct/results/combined.blast

run panOCT

panoct.pl

with default paramerters

 -t: name of btab (wublast-style or ncbi -m8 or -m9) input file [REQUIRED]
 -f: file containing unique genome identifier tags [REQUIRED]
 -g: gene attribute file (asmbl_id<tab>protein_identifier<tab>end5<tab>end3<tab>annotation<tab>genome_tag)
 -P: name of concatinated .pep file [REQUIRED to calc protein lengths]

cd panoct

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/panoct.pl -b results -t combined.blast -f genomes.list -g combined.att -P combined.fasta -S yes -L 1 -M Y -H Y -V Y -N Y -F 1.33 -G y -c 0,50,95,100 -T

Analyses pan-génomiques

Les analyses pan-génomiques fournissent un cadre pour déterminer la diversité génomique de l'ensemble des gènomes analysés, mais aussi pour prédire, par extrapolation, combien de séquences génomiques supplémentaires seraient nécessaires pour caractériser l'ensemble du pan-génome ou répertoire génétique.

Inside the Pan-genome - Methods and Software Overview

Définitions

D'après Vernikos et al. 2015

Le pan-génome a été défini pour la première fois par Tettelin et al., 2005.

le pan-génome englobe l'ensemble du répertoire de gènes accessibles au clade étudié ;
le génome coeur contient les gènes communs à toutes les souches du clade et comprend généralement des gènes responsables des aspects fondamentaux de la biologie du clade et de ses principaux traits phénotypiques ;
le génome accessoire est constitué des gènes communs à un sous-ensemble de souches et contribue à la diversité des espèces. Il peut coder des voies biochimiques supplémentaires et des fonctions qui ne sont pas essentielles à la croissance mais qui confèrent des avantages sélectifs, comme l'adaptation à différentes niches, la résistance antibiotique ou la colonisation d'un nouvel hôte (Medini et al.', 2005)
les gènes spécifiques d’une souche ou singletons désignent des gènes spécifiques à une souche n’ayant par d’orthologs dans les autres souches du clade.

Mise en oeuvre

Ce script permet de prendre en compte les paralogs (Chan et al., 2015).

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/paralog_matchtable.pl -M results/matchtable_0_1.txt -P results/paralogs.txt > results/matchtable_paralog.txt

Celui-ci permet de reformater la table matchtable.

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/matchtable2panoct.result.pl --results_dir results

Ce script calcule un ensemble de statistiques sur les résultats de panOCT.

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/pangenome_statistics.pl --data_file results/combined.att.dat --genomes_list results/genomes.list --method_result results/panoct.result -c results/centroids.fasta -f results/frameshifts.txt --fusion results/fragments_fusions.txt

Un résumé des résultats se trouve dans la fichier: ~/work/Prochlorococcus/panoct/overview_stats.txt

Commentez les résultats obtenus. Comparez ces résultats à ceux de [https://www.ncbi.nlm.nih.gov/pubmed/18159947 Kettler et ''al.'', 2007].

Tracé de l'histogramme

cat overview_stats.txt | perl -ne 'chomp; if ($p) {print "$_\n";} if (/^Cluster Size Breakdown/) {$p = 1;}' > core_cluster_histogram_data.txt

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/core_cluster_histogram.R -i core_cluster_histogram_data.txt

mv core_cluster_histogram.pdf ~/work/Prochlorococcus/images/.

Après avoir indiqué sur la figure, le génome cœur, le génome accessoire, les singletons et le pan génome commentez les résultats obtenus.

Quels sont les choix méthodologiques qui peuvent influencer la taille du génome cœur (argumentez)?

Taille du pan génome

La taille du pan génome à tendance à augmenter avec le nombre de génomes comparés. Estimer sa taille est un exemple d'une classe générale de mesures où, étant donné une collection d’entités et leurs attributs, le nombre d'attributs distincts observés est fonction du nombre d'entités considérées. C’est par exemple le cas de l’analyse du langage naturel, où les entités sont les textes et les attributs sont les mots. Dans ce contexte, l’augmentation du nombre n d'attributs distincts en fonction du nombre N d'entités considérées suit la loi de Heaps (https://en.wikipedia.org/wiki/Heaps%27_law Heaps 'law). Elle peut être représentée par la formule :

n=k*N-α,

où, dans un contexte génétique, n est le nombre attendu de gènes pour un nombre N de génomes et les paramètres k et α (α=1-γ) sont des paramètres libres qui sont déterminés empiriquement (Tettelin et al., 2008).

Appliqué à la découverte de nouveaux gènes et selon la loi de Heap, lorsque α > 1 (γ < 0), le pan-génome est considéré comme fermé, et l'ajout de nouveaux génomes n'augmentera pas significativement le nombre de nouveaux gènes. Par contre, lorsque α < 1 (0 < γ < 1), le pan-génome est ouvert, et pour chaque nouveau génome ajouté, le nombre de gènes augmente significativement (Tettelin et al., 2008).

Pour déterminer les paramètres k et α nous pouvons calculer toutes les combinaisons de 2 à N génomes et reporter la taille du pan génome pour chaque combinaison. Cependant, comme le nombre de combinaisons augmente très rapidement avec le nombre de génomes (C=N!/(n−1)!⋅(N−n), en pratique un échantillonnage des combinaisons possible est réalisé.

compute_pangenome.R

compute_pangenome.R usage:

 -i <input file name for tab delimited 0/1 pan-genome cluster table with column and row headers>
 -o <output file name for combinatorial counts>
 -p <percentage of genomes to be considered core - default 100>
 -q <percentage of genomes to be considered novel - default 0 but clearly at least one genome>
 -s <maximum number of combinatorial samples to generate>

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/compute_pangenome.R -i results/matchtable_paralog.txt -o results/pangenome_size -p 95 -q 0 -s 250

Ajout de la taille des génomes en première colonne:

for i in results/*.tab;  do       awk 'END{ printf "1\t"NR"\t0"NR"\t"NR"\t"NR"\t"NR"\t\n"}' $i; done > results/unique.txt
cat results/pangenome_size results/unique.txt > results/pangenome_size_comp

Afficher les courbes

Pan génome et génome coeur

Rscript /home/formation/public_html/M2_Phylogenomique/scripts/pan_core_plot.R
evince ~/work/Prochlorococcus/images/pangenome.pdf

Log x Log plot des nouveaux gènes

Rscript /home/formation/public_html/M2_Phylogenomique/scripts/new_genes_plot.R
evince ~/work/Prochlorococcus/images/newgenome.pdf

run_panoct.pl

--genome_list_file,-g : File containing list of genome names
--att_dir, -A : Directory containing genome .att files
--gene_att_file. -a : Use this gene_att_file instead of creating one from    att_dir
--blast_file,-b : Path to a file of all-vs-all BLAST results. Must be -m 8     or -m 9 tab delimited output.

--working_dir, -w : Path to be consedered the working directory for the     pipeline. Defaults to cwd ('./').
 --no_grid : Run all subtasks locally instead of in parallel on an     SGE/Univa grid.

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/run_panoct.pl --genome_list_file tags.txt --gene_att_file gene_att  --blast_file combined.blast --no_grid

results/R.plots/core_cluster_histogram.pdf

compute_pangenome.R

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/compute_pangenome.R -i /share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/matchtable_paralog.txt -o /share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/pangenome_size -p 95 -q 0 -s 250

Afficher les courbes:

R
pangenome_size <- '/share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/pangenome_size'
pdf_file <- '/share/gsi2/data/www/soap-data/M2BBS/Phylogenomic/Prochlorococcus/PanGenomePipeline/results/pangenome.pdf'

data <- read.table(file=pangenome_size, header=T)
data

ymax <- max(data$Genome, data$Pan_genome)
ymin <- 0

pdf(file=pdf_file, paper="a4r")
plot(data$Genome, data$Pan_genome, col='blue',xlim=c(2, 12), xlab='Number of genomes', ylab='Number of genes', ylim=c(ymin, ymax), pch=20)
points(data$Genome, data$Core, col='red', pch=20)
dev.off()
cat(pdf_file, "\n")

Analyses phylogénomiques

Comme dans Kettler et al., 2007, nous allons utiliser quatre génomes de Synechococcus comme groupe externe dans nos analyses.

Synechococcus

Génomes et annotations

accession "CP000435, CP000097, BX548020, CP000110"

Les fichiers GenBank peuvent être obtenus avec la commande wget en utilisant les chemins enregistrés dans le fichier: accession_file.lst

Comme précédemment, les génomes sont renommés en utilisant un code à quatre lettres.

accession_labels_file.lst

Les fichiers GenBank renommés sont disponibles dans le répertoire: GenBank

Les séquences ADN ont été extraites des fichiers GenBank : DNA

Fichier avec les informations sur les souches: species_strain_names.txt

mkdir -p ~/work/Synechococcus/prokka
cd ~/work/Synechococcus
/home/formation/public_html/M2_Phylogenomique/scripts/prokka_loop.pl --sample Synechococcus

squeue -l -u <user>

Une fois les jobs terminés, vérifiez que les fichiers de sorties de prokka existent et ne sont pas vides.

ls -l ~/work/Synechococcus/prokka/Aaa*/*.faa

Synechococcus blastp All-All

Nous allons reproduire le même enchaînement de scipts:

mkdir -p ~/work/Synechococcus/peptide
cd ~/work/Synechococcus/peptide
cp ~/work/Synechococcus/prokka/Aaa*/*.faa ~/work/Synechococcus/peptide/.

ls -l ~/work/Synechococcus/peptide

module load bioinfo/ncbi-blast-2.6.0+
for i in *.faa; 
do      
srun -n1 -l makeblastdb -in $i  -dbtype prot; 
done

mkdir -p ~/work/Synechococcus/BlastP
cd ~/work/Synechococcus
/home/formation/public_html/M2_Phylogenomique/scripts/blastp_intra.pl

/home/formation/public_html/M2_Phylogenomique/scripts/blastp_inter.pl

squeue -l -u <user>
ls BlastP | wc -l

Prochlorococcus versus Synechococcus

Liens symboliques sur les fichiers peptides

mkdir -p ~/work/ProchlorococcusSynechococcus/peptide
ln -s ~/work/Prochlorococcus/peptide/*.faa* ~/work/ProchlorococcusSynechococcus/peptide/.
ln -s ~/work/Synechococcus/peptide/*.faa* ~/work/ProchlorococcusSynechococcus/peptide/.

ls -l ~/work/ProchlorococcusSynechococcus/peptide/.

Liens symboliques sur les fichiers blastp

mkdir -p ~/work/ProchlorococcusSynechococcus/BlastP
ln -s ~/work/Prochlorococcus/BlastP/*.tab ~/work/ProchlorococcusSynechococcus/BlastP/.
ln -s ~/work/Synechococcus/BlastP/*.tab ~/work/ProchlorococcusSynechococcus/BlastP/.

ls -l ~/work/ProchlorococcusSynechococcus/BlastP/.

Compléter les pairs de comparisons

cd ~/work/ProchlorococcusSynechococcus
/home/formation/public_html/M2_Phylogenomique/scripts/blastp_inter.pl

squeue -l -u <user>

ls BlastP | wc -l

Groupes de gènes orthologues

Préparation des fichiers combined

srun --pty bash
mkdir -p ~/work/Synechococcus/panoct/results
mkdir -p ~/work/ProchlorococcusSynechococcus/panoct/results

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

cd ~/work/Synechococcus
for file in peptide/*.faa
do
 prefix=$(basename $file .faa)
 echo "/home/formation/public_html/M2_Phylogenomique/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output prokka/$prefix/$prefix.tab"
 /home/formation/public_html/M2_Phylogenomique/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output panoct/results/$prefix.tab
done
ls panoct/results/*.tab

cat ~/work/Prochlorococcus/panoct/results/*.tab ~/work/Synechococcus/panoct/results/*.tab> ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att
head ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att

genomes.list Liste des souches à analyser.

cd ~/work/ProchlorococcusSynechococcus/
for i in peptide/*.faa
do      
echo $(basename $i .faa)
done > panoct/results/genomes.list
more panoct/results/genomes.list

combined.fasta Concaténer les fichier peptides dans un seul fichier.

cat peptide/*.faa > panoct/results/combined.fasta

grep -c '>' panoct/results/combined.fasta

combined.blast Concaténer les résultats des blastp dans un seul fichier.

cat BlastP/*.tab > panoct/results/combined.blast

head panoct/results/combined.blast

run panOCT Prochlorococcus vs Synechococcus

cd ~/work/ProchlorococcusSynechococcus/panoct
mkdir ~/work/ProchlorococcusSynechococcus/images

/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/panoct.pl -b results -t combined.blast -f genomes.list -g combined.att -P combined.fasta -S yes -L 1 -M Y -H Y -V Y -N Y -F 1.33 -G y -c 0,50,95,100 -T

Heatmap des groupes de gènes orthologues

srun --pty bash
module load system/R-3.5.1;
R
pdf_file <- '~/work/ProchlorococcusSynechococcus/images/matchtable_heatmap.pdf'
strain_file <- '~/work/ProchlorococcusSynechococcus/panoct/results/genomes.list'
matchtable_0_1 <- '~/work/ProchlorococcusSynechococcus/panoct/results/matchtable_0_1.txt'

strains <- read.delim(file=strain_file, header=FALSE)
data <- read.delim(file=matchtable_0_1, , sep="\t", header=FALSE, row.names=1)
colnames(data) <- t(strains)

pdf(file=pdf_file, paper="a4r")
heatmap(t(as.matrix(data)), scale='none', col=c('white', 'darkblue'), xlab="Strains", labCol=NA)
dev.off()
cat(pdf_file, "\n")

Gènes espèces spécifique

Gènes trouvés chez toutes les souches de Prochlorococcus mais dans aucune de Synechococcus:

awk 'BEGIN {
	h=0;
}
{
	mp=0
	for (i=2; i<=13; i++) {
	if( $i !~ /-/ ) mp++
	}
	sp=0
	for (i=14; i<=17; i++) {
	if( $i !~ /-/ ) sp++
	}
	if ( mp==12 && sp ==0 ) {
	h++;
	printf("%d\t%d", h, $1);
	for (i=2; i<=17; i++) {
	if( $i !~ /-/ ) printf("\t%s", $i);
	}
	printf("\n");

	}
	
}' < /home/formation/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt > /home/formation/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt

Quelles sont les fonctions de ces gènes?

awk '{ system("grep "$3"  ../Prochlorococcus/prokka/Aaaa/Aaaa.gff >> Prochlorococcus_specific_annotations.txt")}' < /home/formation/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt

Gènes trouvés chez toutes les souches de Synechococcus mais dans aucune de Prochlorococcus:

A vous de jouer!

Extraction des séquences nucléotidiques des gènes orthologues

Créer un fichier avec toutes les séquences nucléotidiques:

mkdir ~/work/ProchlorococcusSynechococcus/DNA
cat ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.ffn ~/work/Synechococcus/prokka/Aaa*/Aaa*.ffn > ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn

grep -c '>' ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn

Extraire les groupes de gènes orthologues

mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA
/home/formation/public_html/M2_Phylogenomique/scripts/extract_sequences_from_matchtable.pl --fasta ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn --matchtable ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt --outdir ~/work/ProchlorococcusSynechococcus/OG/DNA --quorum 16

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

grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas

Transformation en alignement multiples des séquences ADN des gènes

Le script aa_to_dna_aln.pl prend en entrée un fichier fasta contenant des séquences d'ADN. Il traduit des séquences, réalise un alignement multiple avec muscle et retraduit en ADN cet alignement multiple (BioPerl).

mkdir ~/work/ProchlorococcusSynechococcus/OG/alignment

/home/formation/public_html/M2_Phylogenomique/scripts/aa_to_dna_aln.pl -dna ~/work/ProchlorococcusSynechococcus/OG/DNA/OG_1471.fas --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment

La traduction en peptides ajoute un '*' à la fin des séquences. Vous pouvez ignorer le message:

*** WARNING *** Invalid character '*' in FASTA sequence data, ignored

En sortie:

pep_OG_1471.fas       peptides
ali_pep_OG_1471.fas   alignement des peptides
ali_dna_OG_1471.fas   alignement des nucléotides

Le programme suivant réalise le traitement pour un sous ensemble de sample fichiers tirés au hasard.

/home/formation/public_html/M2_Phylogenomique/scripts/aa_to_dna_aln_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100

Construction des arbres des groupes de gènes orthologues

Edition des alignements avec trimal

Trimer les alignements avec trimal (pas sur genologin!). Nous pouvons utiliser les résultats de trimal pour identifier des alignements de mauvaise qualité.

ssh -Y <user>@genotoul.toulouse.inra.fr

for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_*.fas
do      
dir=$(dirname $i)
prefix=$(basename $i .fas)
/usr/local/bioinfo/bin/trimal -in $dir/$prefix.fas -out $dir/$prefix.trim.aln -sident -sgt -automated1
done > ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.txt

Pour chaque alignement, extraction des % de sites avec des indels et de la conservation moyenne des alignements

awk 'BEGIN {
	i=1;
	col[""]=0;
}
{ 
	if ( $5==0) 
	{
	col[i "," 1]=$2
	} 
	else 
	if ( $2=="AverageIdentity" && $3!="Average") 
	{
	col[i "," 2]=$3
	i = i+1;
	} 
} 
END {
	max = i
	i=1;
	while (i < max) {
	print i "\t" col[i "," 1] "\t" col[i "," 2];
	i = i+1;
	}
}' < ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.txt > ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.tab

awk 'BEGIN { min = max = "NaN" }
{
    min = (NR==1 || $2<min ? $2 : min)
    max = (NR==1 || $2>max ? $2: max)
}
END { print min, max }' < ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.tab

Tracé de la distribution des deux paramètres

Tracé de la distribution des deux paramètres pour définir des seuils permettant de retenir les meuilleurs alignements

R
pdf_file <- '~/work/ProchlorococcusSynechococcus/images/alignment_statistics.pdf'
data <- read.table('~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.tab', row.names=1)
pdf(file=pdf_file, paper="a4r")
plot(data, xlim=c(1,100), ylim=c(0,1), pch=20, xlab="% without gap", ylab="AverageIdentity")
text(data, labels=row.names(data), pos=4)
dev.off()
cat(pdf_file, "\n")

Sélection des meilleurs alignements

Sélection des meilleurs alignements à partir de la distribution du nombre de sites sans indels et de la conservation moyenne des alignements

/home/formation/public_html/M2_Phylogenomique/scripts/trimal_selection.pl --alignment ~/work/ProchlorococcusSynechococcus/OG/alignment -gap 70 --identity 0.4

Phylogénie des souches

Super Matrice

Concaténation de 31 alignements

/home/formation/public_html/M2_Phylogenomique/scripts/concat_aligments.pl  --alignments ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70.lst  --outfile ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 -nb_ali 31

IQ-TREE

/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 -redo  -bb 1000 -alrt 1000 -nt AUTO

Input data: 16 sequences with 32355 nucleotide sites
Number of constant sites: 10622 (= 32.8295% of all sites)
Number of invariant (constant or ambiguous constant) sites: 10622 (= 32.8295% of all sites)
Number of parsimony informative sites: 20377
Number of distinct site patterns: 16346

ModelFinder
-----------

Best-fit model according to BIC: GTR+F+R4

Total CPU time used: 2936.33161 seconds (0h:48m:56s)
Total wall-clock time used: 459.9323409 seconds (0h:7m:39s)

Enraciner l'arbre

~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31_nuc_GTR+F+R4_rooted.treefile

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

Nous utilisons rnammer pour annoter les ARNr (lsu, ssu, tsu) dans les génomes.

/home/formation/public_html/M2_Phylogenomique/scripts/rnammer_loop.pl --prokka_dir ~/work/Synechococcus/prokka --model ssu
/home/formation/public_html/M2_Phylogenomique/scripts/rnammer_loop.pl --prokka_dir ~/work/Prochlorococcus/prokka --model ssu

Vérifiez que les fichier de sortie ne sont pas vide!

ls -l /home/yquentin/work/*/prokka/Aaa*/Aaa*ssu*.rrna

Changer le motif pour obtenir les deux autres ARNr:

ls -l /home/yquentin/work/*/prokka/Aaa*/Aaa*lsu*.rrna
ls -l /home/yquentin/work/*/prokka/Aaa*/Aaa*tsu*.rrna

Concaténer les fichiers:

mkdir ~/work/ProchlorococcusSynechococcus/rRNA
cat /home/yquentin/work/*/prokka/Aaa*/Aaa*lsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/lsu.fas
cat /home/yquentin/work/*/prokka/Aaa*/Aaa*ssu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/ssu.fas
cat /home/yquentin/work/*/prokka/Aaa*/Aaa*tsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/tsu.fas

Alignements des ARNr

Mafft comporte deux options, Q-INS-i et X-INS-i, dans lesquelles les informations de structure secondaire de l'ARN sont prises en compte. Ces méthodes sont adaptées à un alignement global de séquences d'ARNc très divergentes. Pour les ARN relativement conservés, tels que les ARNr SSU et LSU, l'avantage de ces méthodes est faible (Katoh et al., 2103). Nous utilisons néanmoins la version mafft-qinsi.

ssu

srun --pty bash
module load bioinfo/mafft-7.313
mafft-qinsi  --maxiterate 1000 --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.fas > ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.aln
trimal -in ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.aln -automated1 -phylip -out ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy

lsu

srun --pty bash
module load bioinfo/mafft-7.313
mafft-qinsi  --maxiterate 1000 --globalpair --thread -1 ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified.fas > ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified.aln
trimal -in ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified.aln -automated1 -phylip -out ~/work/ProchlorococcusSynechococcus/rRNA/lsu_renamed_simplified_trimed.phy

tsu

srun --pty bash
module load bioinfo/mafft-7.313
mafft-qinsi  --maxiterate 1000 --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified.fas > ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified.aln
trimal -in ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified.aln -automated1 -phylip -out ~/work/ProchlorococcusSynechococcus/rRNA/tsu_renamed_simplified_trimed.phy

Arbre SSU

Arbre SSU avec PhyML

qsub -V -b Y -N phymlRNAn -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/PhyML/current/PhyML-3.1_linux64 -i ~/work/ProchlorococcusSynechococcus/rRNA/ssu_trimed.phy -n 1 -b -1 -m GTR -v e -c 4 -a e -o n --quiet"

Vous pouvez observer que des génomes codent pour plusieurs copies du gènes codant pour l'ARN 16S. A l'aide de l'arbre et des scores de rnammer simplifier votre fichier. Il est également utilie de modifier le nom des séquences. Après ces modifications, vous relancez l'alignement, l'édition de l'alignement et le calcule de l'arbre avec des options plus poussées.

qsub -V -b Y -N phymlRNAtlr -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/PhyML/current/PhyML-3.1_linux64 -i ~/work/ProchlorococcusSynechococcus/rRNA/ssu_trimed.phy -n 1 -b -1 -m GTR -v e -c 4 -a e -o tlr --quiet"

Arbre SSU avec IQ-TREE

Pour seulement trouver le modèle le mieux adapté sans faire de reconstruction d'arbre, utilisez :

qsub -V -b Y -N IQTreeSSUm -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy  -m MF -redo -AIC"

Les résultats sont dans le fichier : ssu_renamed_simplified_trimed.phy.iqtree.

grep 'Best-fit model' ssu_renamed_simplified_trimed.phy.iqtree

lsu ssu GTR+F+R2 tsu K2P+G4

dna-models

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

qsub -V -b Y -N ssuIQTreeSSUbb -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy  -pre ssuGTRFR2bb1000bnni -m GTR+F+R2 -bb 1000 -redo -bnni -nt AUTO"

Évaluer les supports de branche avec des tests de branche simple :

IQ-TREE propose le test du rapport approximatif de vraisemblance de type SH (Guindon et al., 2010).

qsub -V -b Y -N ssuIQTreeSSUbbalrt -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy -pre ssuGTRFR2bbalrt -m GTR+F+R2 -bb 1000 -alrt 1000 -redo -nt AUTO"

Évaluation des supports de branche avec un bootstrap non paramétrique standard :

qsub -V -b Y -N ssuIQTreeSSUalrtb -l h_vmem=20G -l mem=18G "/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_trimed.phy -pre ssuGTRFR2alrtb -m GTR+F+R2 -alrt 1000 -b 100 -redo -nt AUTO"

Reconstruction d'états ancestraux

mkdir ~/work/ProchlorococcusSynechococcus/Ancestral_Characters 
cd ~/work/ProchlorococcusSynechococcus/Ancestral_Characters

Copier l'arbre espèce de référence

cp ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31_nuc_GTR+F+R4_rooted.treefile ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph

Composition en GC des génomes

for i in ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.fna
do      
/home/formation/public_html/M2_Phylogenomique/scripts/gc_freq.pl --file $i
done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt

for i in ~/work/Synechococcus/prokka/Aaa*/Aaa*.fna
do      
/home/formation/public_html/M2_Phylogenomique/scripts/gc_freq.pl --file $i
done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Synechococcus_gc.txt

cat ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Synechococcus_gc.txt > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt

Arbre espèces

R

library(ape)
treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 

strains_gc_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt"

tr <- read.tree(treefile)

Vérifier que l'arbre tr est enraciné:

is.rooted(tr)

plot(tr)
plot(tr <- root(tr, interactive = TRUE))

Informations sur les souches

strains_info_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt"
strains_info <- read.table(file=strains_info_file, header=T,sep="\t", row.names=1)
strains_info <- strains_info[tr$tip.label,]

Changement des labels de l'arbre

tr$tip.label <- as.character(strains_info$Strain)
tipcol <- c()
tipcol[1:4] <- 'red'
tipcol[5:16] <- 'blue'

plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol)
nodelabels(tr$node, bg='white', col='purple', frame='n', adj=c(1,-1), cex=0.5)
tiplabels(strains_info$Light, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
add.scale.bar(length=0.1)

Ecotypes

Codage des ecotypes HL et LL

ecotypes <- c()
ecotypes[grep('Syn', strains_info$Light)] <- 'Syn'
ecotypes[grep('LL', strains_info$Light)] <- 'LL'
ecotypes[grep('HL', strains_info$Light)] <- 'HL'
names(ecotypes)<-tr$tip.label
ecotypes

eco_er <- ace(ecotypes, tr, type="discrete", CI=T, model="ER")
plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol)
tiplabels(strains_info$Light, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=eco_er$lik.anc,cex=0.5)
add.scale.bar(length=0.1)

eco_sym <- ace(ecotypes, tr, type="discrete", CI=T, model="SYM")

eco_ard <- ace(ecotypes, tr, type="discrete", CI=T, model="ARD")

The increase in likelihood with the additional parameter
is not significant:

1 - pchisq(2*(eco_er$loglik - eco_sym$loglik), 1)

In practice, there is an anova method for the class "ace" which simplifies the
calculation of this test:

anova(eco_er, eco_ard, eco_sym)

plot(tr) nodelabels(pie=eco_er$lik.anc,cex=0.7) </pre>

Flux de gènes

matchtable

matchtable_file <- "~/work/ProchlorococcusSynechococcus/panoct/results//matchtable_0_1.txt"
genomes.list_file <- "~/work/ProchlorococcusSynechococcus/panoct/results//genomes.list"

matchtable <- read.table(matchtable_file, sep="\t", row.names=1)
genomes.list <- read.table(genomes.list_file)
colnames(matchtable) <- t(genomes.list)
head(matchtable)

Simplifier la table

R

library(ape)
treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 

strains_gc_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt"

tr <- read.tree(treefile)

Transposer la matrice et ordonner les lignes comme les feuilles de l'arbre

t_matchtable <- t(matchtable)
t_matchtable <- t_matchtable[tr$tip.label,]

nb_strains <- nrow(t_matchtable)

Compiler les différents motifs observés dans une liste

motifs <- list()
for ( cluster in 1:ncol(t_matchtable) ) {
   pattern <- paste(t_matchtable[,cluster], sep='', collapse='')
   nb <- length(motifs[[pattern]]$occurences)
   if ( nb == 0 ) {
      motifs[[pattern]]$sring <- t_matchtable[,cluster]
   }
   motifs[[pattern]]$occurences[[nb+1]] <- cluster
   cat(cluster, pattern, nb, sep="\t", "\n")
}

Traiter chaque motifs

Inférence des états ancestraux

nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2)
nb <- 400
nbpattern <- 1
#while (nb <= length(motifs) ) {
while (nb < 401 ) {
   pattern <- names(motifs)[nb]
   if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) {
      er <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ER")
      motifs[[pattern]]$er <- er
   }
   nb <- nb+1
}

Fonction pour annoter les états des feuilles de l'arbre

tip_states  <- function(pattern) {
   leaves <- matrix(c(0), nrow=nb_strains, ncol=2)
   for ( i in 1:nb_strains ) {
      leaves[i,1] <- 0
      leaves[i,2] <- 0
      if (pattern[i] == 0 ) {
         leaves[i,1]  <- 1
      } else {
         leaves[i,2] <- 1
      }
      #cat(pattern[i], leaves[i,1], leaves[i,2], "\n")
   }
   return(leaves)
}

Exemple d'appel de la fonction

leaves <- tip_states(motifs[[pattern]]$sring)

Concaténation des deux matrices (feuilles et noeuds)

edge_states <- rbind(leaves, er$lik.anc)

Fonction pour annoter les transitions d'états sur les branches de l'arbre

states_transitions <- function(tr, edge_states) {
   
   nb_edges <- nrow(tr$edge)
   direction <- c()
   transition <- c()
   transcolor <- c()

   for ( i in 1:nb_edges ) {
      child <- tr$edge[i,2]
      parent <- tr$edge[i,1]
      if ( edge_states[child, 1] > 0.5 ) {
         child_states <- 0 
      } else {
         child_states <- 1
      } 
      if ( edge_states[parent, 1] > 0.5 ) {
         parent_states <- 0 
      } else {
         parent_states <- 1
      } 
      if ( parent_states != child_states ) {
         cat(i, parent , parent_states, '->', child_states, child, "\n")
         transition <- c(transition, i)
         if ( parent_states == 0 ) {
            direction <- c(direction, '+')
            transcolor <- c(transcolor, 'blue')
         } else {
            direction <- c(direction, '-')
            transcolor <- c(transcolor, 'red')
        }
      }
   }
   translist <- list(transition=transition, direction=direction, transcolor=transcolor)
   return(translist)
}

Exemple d'appel de la fonction

states_transitions(tr, edge_states)

cat ("Press [enter] to continue") line <- readline() plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol) tiplabels(motifspattern$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1) edgelabels(text=direction, edge=transition, frame="r", bg=transcolor, adj=c(0,-0.5),cex=0.2) </pre>

Tracé pour chaque motif

nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2)
nb <- 400
nbpattern <- 1
#while (nb <= length(motifs) ) {
while (nb < 401 ) {

   pattern <- names(motifs)[nb]
   if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) {
      leaves <- tip_states(motifs[[pattern]]$sring)
      edge_states <- rbind(leaves, motifs[[pattern]]$sring$er$lik.anc)
      states_transitions(tr, edge_states)

      plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol)
      tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
      nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5)
      add.scale.bar(length=0.1)

      cat ("Press [enter] to continue")
      line <- readline()

  }
   nb <- nb+1
}

plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol)
tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=nodesflux/nbpattern,cex=0.5)
nodelabels()
add.scale.bar(length=0.1)

   tr$edge

Atelier Phylogénomique

From silico.biotoul.fr

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

Installation du package genoPlotR

Logiciels à installer sur vos postes de travail

Scripts

Disponibilité des génomes

Caractéristiques des souches étudiées

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 with Prochlorococcus

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

gene.att

tags.txt

peptides

blast.txt

run panOCT

panoct.pl

Analyses pan-génomiques

Définitions

Mise en oeuvre

Taille du pan génome

compute_pangenome.R

Afficher les courbes

Pan génome et génome coeur

Log x Log plot des nouveaux gènes

run_panoct.pl

compute_pangenome.R

Analyses phylogénomiques

Synechococcus

Génomes et annotations

Synechococcus blastp All-All

Prochlorococcus versus Synechococcus

Liens symboliques sur les fichiers peptides

Liens symboliques sur les fichiers blastp

Compléter les pairs de comparisons

Groupes de gènes orthologues

Préparation des fichiers combined

run panOCT Prochlorococcus vs Synechococcus

Gènes espèces spécifique

Gènes trouvés chez toutes les souches de Prochlorococcus mais dans aucune de Synechococcus:

Gènes trouvés chez toutes les souches de Synechococcus mais dans aucune de Prochlorococcus:

Extraction des séquences nucléotidiques des gènes orthologues