silico.biotoul.fr
 

Atelier Phylogénomique

From silico.biotoul.fr

Revision as of 13:42, 12 September 2018 by Quentin (Talk | contribs)
Jump to: navigation, search

Contents

Matériel pédagogique



Références


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>@genotoul.toulouse.inra.fr

Vous pouvez ouvrir deux connections afin de pouvoir lancer qlogin de façon indépendante.

Échange de fichiers: 

scp file login@genotoul.toulouse.inra.fr:~/work

Logiciels disponibles

ou plus directement 

ls /usr/local/bioinfo/src/

La documentation est disponible sur le site WEB et dans le répertoire correspondant au logiciel (fichiers How_to_use and/or Readme).

SLURM cluster softwares

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

2 - To submit the job, use the sbatch command line as following: 

sbatch myscript.sh

Pour controler les jobs 

q stat -u login : liste des jobs de l'utilisateur.
qstat -j job_id : informations sur le job spécifié.
qstat -s r : liste des job en cours d'exécution.

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

Prokka

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

Exemple d'utilisation: 

srun --pty bash
module load bioinfo/prokka-1.12
prokka  /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaaa.fas  --outdir /tmp/Aaaa --compliant --addgenes --prefix Aaaa  --locustag Aaaa.g --genus Prochlorococcus --species 'Prochlorococcus marinus subsp. marinus' --strain CCMP1375 --kingdom Bacteria --cpus 2

Le programme génère plusieurs fichiers pour chaque réplicon, dont:

Peptides

Pour chaque génome les peptides de chaque replicon sont assemblés dans un seul fichier: peptide.

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 -R /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/peptide .
make blast database
makeblastdb -in Aaaa.pep -dbtype prot

for i in *.pep; 
do      
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: 

#!/usr/bin/perl -w

use strict;
use Cwd;

my $hostname = $ENV{'HOSTNAME'};
my $node = '';

if ( defined $ENV{'SLURMD_NODENAME'} ) {
	$node = $ENV{'SLURMD_NODENAME'};
}
if ( $hostname !~ /genologin/ || $node ne '' ) {
	print STDERR "Error: log on $node, change for genologin server.\n";
	exit(1);
}

my $cwd = cwd();
if ( !-e '../BlastP' ){
	print STDERR "Error ../BlastP is expected!\n";
	exit;
}

# parameters
my $pattern = '.pep$';
my $evalue = 1e-5;

print "$hostname\t$node\n$cwd\n";

# look for files in the current directory
opendir(my $rdir, $cwd) or die "Cannot open $cwd: $!";
my @file_list = grep /$pattern/, sort(readdir $rdir);
closedir $rdir;
if ( scalar @file_list == 0 ) {
	print STDERR "Error: no file $pattern found in  $cwd!\n";
	exit(1);
}
print scalar @file_list, " $pattern entries in $cwd\n";

my($cmd, $prefix, $outfile, $workname, $sbatch_file, $sbatch_out, $sbatch_err);
foreach my $file (@file_list ) {
	print "$file\n";
	if ( $file =~ /(\w+)$pattern/ ) {
		$prefix = $1;
		my $outfile = '../BlastP/'. $prefix . '_' . $prefix . '.tab';
		if ( -e $outfile ) {
			print "skip: $outfile\n";
		} else {
			$cmd = "blastp -query $file -db $file -seg yes -dbsize 100000000  -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile";
			print "$cmd\n";
			$workname = 'blastp_' . $prefix;
			$sbatch_file = $prefix . '.sh';
			$sbatch_out  = $prefix . '.out';
			$sbatch_err  = $prefix . '.err';
			run_script($sbatch_file, $workname, $sbatch_out, $sbatch_err, $cmd);
			print "$sbatch_file\n";
			`sbatch $sbatch_file`;
		}
	} 
}

########################################################################
sub run_script {
	my($script_file, $workname, $out, $err, $cmd) = @_;
	
	open (my $SCRTFILE, ">$script_file") || die "Blast Error Unable to create temporary file: $script_file: $!";

	print $SCRTFILE '#!/bin/bash' . "\n";
	print $SCRTFILE '#SBATCH -J t' . "$workname\n";
	print $SCRTFILE '#SBATCH -o ' .  "$out\n";
	print $SCRTFILE '#SBATCH -o ' .  "$err\n";
	print $SCRTFILE '#SBATCH --time=00:5:00' .  "\n";
	print $SCRTFILE '#SBATCH --cpus-per-task=1' .  "\n";
	print $SCRTFILE "module purge\n";	 
	print $SCRTFILE "module load bioinfo/ncbi-blast-2.6.0+\n";	 
	print $SCRTFILE "$cmd\n";	 
	close($SCRTFILE);
}

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

squeue -l -u <user>
Inter genomes

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

#!/usr/bin/perl -w

use strict;
use Cwd;

my $hostname = $ENV{'HOSTNAME'};
my $node = '';

if ( defined $ENV{'SLURMD_NODENAME'} ) {
	$node = $ENV{'SLURMD_NODENAME'};
}
if ( $hostname !~ /genologin/ || $node ne '' ) {
	print STDERR "Error: log on $node, change for genologin server.\n";
	exit(1);
}

my $cwd = cwd();
if ( !-e '../BlastP' ){
	print STDERR "Error ../BlastP is expected!\n";
	exit;
}

# parameters
my $pattern = '.pep$';
my $evalue = 1e-5;


print "$hostname\t$node\n$cwd\n";

# look for files in the current directory
opendir(my $rdir, $cwd) or die "Cannot open $cwd: $!";
my @file_list = grep /$pattern/, sort(readdir $rdir);
closedir $rdir;
if ( scalar @file_list == 0 ) {
	print STDERR "Error: no file $pattern found in  $cwd!\n";
	exit(1);
}
print scalar @file_list, " $pattern entries in $cwd\n";

my($cmd, $qprefix, $dbprefix, $prefix, $outfile, $workname, $sbatch_file, $sbatch_out, $sbatch_err);
foreach my $query (@file_list ) {
	if ( $query =~ /(\w+)$pattern/ ) {
			$qprefix = $1;
	} else {
		next;
	}
	foreach my $db (@file_list ) {
		next if ($query eq  $db );
		if ( $db =~ /(\w+)$pattern/ ) {
			$dbprefix = $1;
			$prefix = $qprefix . '_' . $dbprefix;
			my $outfile = '../BlastP/'. $qprefix . '_' . $dbprefix . '.tab';
			if ( -e $outfile ) {
				print "skip: $outfile\n";
			} else {
				$cmd = "blastp -query $query -db $db -seg yes -dbsize 100000000  -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile";
				print "$cmd\n";
				$workname = 'blastp_' . $prefix;
				$sbatch_file = $prefix . '.sh';
				$sbatch_out  = $prefix . '.out';
				$sbatch_err  = $prefix . '.err';
				run_script($sbatch_file, $workname, $sbatch_out, $sbatch_err, $cmd);
				print "$sbatch_file\n";
				`sbatch $sbatch_file`;
			}
		} 
	}
}
########################################################################
sub run_script {
	my($script_file, $workname, $out, $err, $cmd) = @_;
	
	open (my $SCRTFILE, ">$script_file") || die "Blast Error Unable to create temporary file: $script_file: $!";

	print $SCRTFILE '#!/bin/bash' . "\n";
	print $SCRTFILE '#SBATCH -J t' . "$workname\n";
	print $SCRTFILE '#SBATCH -o ' .  "$out\n";
	print $SCRTFILE '#SBATCH -o ' .  "$err\n";
	print $SCRTFILE '#SBATCH --time=00:5:00' .  "\n";
	print $SCRTFILE '#SBATCH --cpus-per-task=1' .  "\n";
	print $SCRTFILE "module purge\n";	 
	print $SCRTFILE "module load bioinfo/ncbi-blast-2.6.0+\n";	 
	print $SCRTFILE "$cmd\n";	 
	close($SCRTFILE);
}



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

squeue -l -u <user>

PorthoMCL

PorthoMCL

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

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

Le logiciel est basé sur l'enchainement de programmes python. Nous avons écrit un script permettant d'automatiser l'exécution de ces programmes dans le cadre de notre analyse. 

 
rm -R ~/work/Prochlorococcus/PorthoMCL

/home/formation/public_html/M2_Phylogenomique/scripts/short_PorthoMCL.sh

Le package renferme un script (porthomcl.sh) pour automatiser la procédure pour une application plus général.

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 

#Load modules
module load bioinfo/mcl-14-137

rm ~/work/Prochlorococcus/PorthoMCL/8.all.ort.5.50.tsv
rm ~/work/Prochlorococcus/PorthoMCL/8.all.par.5.50.tsv

cat ~/work/Prochlorococcus/PorthoMCL/6.orthologs.5.50/*.tsv >> ~/work/Prochlorococcus/PorthoMCL/8.all.ort.5.50.tsv

mcl ~/work/Prochlorococcus/PorthoMCL/8.all.ort.5.50.tsv --abc -I 1.5 -t 4 -o ~/work/Prochlorococcus/PorthoMCL/8.all.ort.group.5.50

cat ~/work/Prochlorococcus/PorthoMCL/7.paralogs/*.tsv >> ~/work/Prochlorococcus/PorthoMCL/8.all.par.tsv

mcl ~/work/Prochlorococcus/PorthoMCL/8.all.par.tsv  --abc -I 1.5 -t 4 -o ~/work/Prochlorococcus/PorthoMCL/8.all.par.group.I1.5

directory <- '/home/yquentin/work/Prochlorococcus/PorthoMCL/4.splitSimSeq'
filenames <- list.files(directory, pattern = "*..ss.tsv*", full.names = TRUE)
nbfile <- length(filenames)

allintra <- NULL
for (i in 1:nbfile){ 
	cat(filenames[i], "\n")
	data <- read.table(filenames[i], header=F, stringsAsFactors = FALSE)
	intra <- subset(data, data[,1]!=data[,2])
	allintra <- rbind(allintra, intra)
}
nrow(allintra)

png(pdffile, width = 480, height = 480, units = "px", pointsize = 12, bg = "white")
hist(allintra$V5, nclass=200, xlim=c(0,100), xlab='percentMatch', col='lightblue', freq=F)
hist(allintra$V4, nclass=200, xlab='evalueExponent', col='lightblue', freq=F)
dev.off()

old stuff

Autres ressources

PATRIC

PATRIC est une ressource importante pour les bactéries patogènes: " The Pathosystems Resource Integration Center (PATRIC): the bacterial Bioinformatics Resource Center (Wattam et al., 2017)".

Joint Genome Institute

JGI is a DOE Office of Science User Facility managed by Lawrence Berkeley National Laboratory© 1997-2017 The Regents of the University of California.

EnsemblBacteria

Ensembl Bacteria is a browser for bacterial and archaeal genomes. 

Search for a genome: Prochlorococcus
20 entries
Downloads, Filter Prochlorococcus
16 entries

Annotations des génomes

Nous allons annoter tous les génomes avec la même stratégie pour avoir une annotation homogène.

PROKKA

"Prokka is a software tool for the rapid annotation of prokaryotic genomes": prokka home.

Citation: Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014 Jul 15;30(14):2068-9. PMID:24642063

Nous allons utiliser la version prokka 1.10 (current) qui est disponible sur le serveur.

Exemples d'utilisations

Dirofilaria immitis wolbachia

--outdir    répertoire avec les fichiers résultats.
--compliant forcer la compatibilité avec GenBank/EMBL

Créer un répertoire : 

mkdir /home/<user_name>/work/wolbachia
cd /home/<user_name>/work/wolbachia

en suivant les recommendations du fichier How_to_use (head /usr/local/bioinfo/src/PROKKA/How_to_use) 

module load bioinfo/prokka-1.10

car la version current correspond à prokka-1.10.

Exécuter le programme à l'aide de qsub (fichier commande) ou connectez-vous sur un nœud du cluster avec qlogin. 

/usr/local/bioinfo/src/PROKKA/current/bin/prokka --force --outdir prokka/wDim --addgenes --prefix wDim --locustag wDim --species 'Wolbachia endosymbiont' --strain wDim --compliant /home/formation/public_html/M2_Phylogenomique/data/Wolbachia/Additional/Dirofilaria_immitis_wolbachia_2.2.fna

Output Files 

.gff 	This is the master annotation in GFF3 format, containing both sequences and annotations. It can be viewed directly in Artemis or IGV.
.gbf 	This is a standard Genbank file derived from the master .gff. If the input to prokka was a multi-FASTA, then this will be a multi-Genbank, with one record for each sequence.
.fna 	Nucleotide FASTA file of the input contig sequences.
.faa 	Protein FASTA file of the translated CDS sequences.
.ffn 	Nucleotide FASTA file of all the prediction transcripts (CDS, rRNA, tRNA, tmRNA, misc_RNA)
.sqn 	An ASN1 format "Sequin" file for submission to Genbank. It needs to be edited to set the correct taxonomy, authors, related publication etc.
.fsa 	Nucleotide FASTA file of the input contig sequences, used by "tbl2asn" to create the .sqn file. It is mostly the same as the .fna file, but with extra Sequin tags in the sequence description lines.
.tbl 	Feature Table file, used by "tbl2asn" to create the .sqn file.
.err 	Unacceptable annotations - the NCBI discrepancy report.
.log 	Contains all the output that Prokka produced during its run. This is a record of what settings you used, even if the --quiet option was enabled.
.txt 	Statistics relating to the annotated features found.
.tsv 	Tab-separated file of all features: locus_tag,ftype,gene,EC_number,product

more prokka/wDim/wDim.txt

organism: Genus wolbachia endosymbiont wDim 
contigs: 2
bases: 921012
rRNA: 3
gene: 793
CDS: 755
tRNA: 34
tmRNA: 1

Litomosoides sigmodontis wolbachia

/usr/local/bioinfo/src/PROKKA/current/bin/prokka --force --outdir prokka/wLsi --addgenes --prefix wLsi --locustag wLsi --species 'Wolbachia endosymbiont' --strain wLsi --compliant /home/formation/public_html/M2_Phylogenomique/data/Wolbachia/Additional/Litomosoides_sigmodontis_wolbachia_2.0.fna

more prokka/wLsi/wLsi.txt

organism: Genus wolbachia endosymbiont wLsi contigs: 10
bases: 1048936
rRNA: 3
gene: 918
CDS: 880
tRNA: 34
tmRNA: 1

Visualisation des annotations

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

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

Transférer les fichiers. 

art prokka/wDim/wDim.gff

Comparaison des génomes

Nous allons utiliser le logiciel act (documentation).

Il est nécessaire de réaliser au préalable une comparison des deux génome (BLAST ou MEGABLAST).

Exemple d'utilisation de blastall 

formatdb -i prokka/wDim/wDim.fna -p F

Créer un nouveau répertoire: act 

blastall -p blastn -i prokka/wLsi/wLsi.fna -d prokka/wDim/wDim.fna -m 8 -o act/wLsi-wDim.blt

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

act prokka/wLsi/wLsi.gbf act/wLsi-wDim.blt prokka/wDim/wDim.gbf&

Vous pouvez aussi utiliser les fichiers en format gff.

Commentez les résultats!

Compatibilité avec MAUVE

Il peut y avoir un problème lors du chargement des fichiers gbk dans MAUVE. Une solution à été proposée:

"Mauve uses BioJava to parse GenBank files, and it is very picky about Genbank files. It does not like long contig names, like those from Velvet or Spades. One solution is to use --centre XXX in Prokka and it will rename all your contigs to be NCBI (and Mauve) compliant. It does not like the ACCESSION and VERSION strings that Prokka produces via the "tbl2asn" tool. The following Unix command will fix them:" 

egrep -v '^(ACCESSION|VERSION)' prokka.gbk > mauve.gbk

Annotation des souches

Nous allons retenir les souches dont le génome est 'Complete Genome' ou 'Scaffold'.

Un script perl va automatiser la procédure: 

/home/formation/public_html/M2_Phylogenomique/scripts/rename_ncbi_files.pl --prokka_dir prokka --verbose 0

Les fichiers fasta sont dans /home/<user_name>/work/wolbachia/prokka/

Le script suivant permet de lancer prokka sur un ensemble de génomes en utilisant qsub. 

/home/formation/public_html/M2_Phylogenomique/scripts/prokka_loop.pl --directory /home/formation/public_html/M2_Phylogenomique/data/Wolbachia/NCBI --prokka_dir prokka --genome_list "wMel wMun" --verbose 0 --erase 0

16 rRNA

Vérifier la présence d'une annotation pour le gène d'ARNr 16S dans les différents génomes.

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

Exemple pour le génome wAu, le nom du gène est wAu_01218. 

perl -ne 'if(/^>(\S+)/){$c=grep{/^$1$/}qw(wAu_01218)}print if $c' prokka/wAu/wAu.ffn

Arbre des ARNr 16S

L'ARNr 16S est très souvent utilisé pour inférer la phylogénies des bactéries. Nous allons utiliser les séquences extraites pour réaliser un arbre souches sur nos données.

La (les) méthode(s) à mettre en oeuvre est (sont) laissée(s) à votre choix!

MLST

Nous allons utiliser la base de données MLST : Wolbachia PubMLST database (wolbachia) pour extraire les gènes MLST utilisés pour classer les souches de Wolbachia.

MLST genes

Les fichiers sont dans le répertoire : /home/formation/public_html/M2_Phylogenomique/data/Wolbachia//MLST

HMM profiles

Pour annoter les gènes MLST dans nos génomes, nous allons utiliser hmmer. La première étape est de construire des fichiers profiles pour les 5 gènes. 

/home/formation/public_html/M2_Phylogenomique/scripts/fasta2hmm.pl --mlst_dir MLST --gene_list "gatB coxA hcpA ftsZ fbpA" --verbose 0 --erase 0

Vérifier que les fichiers MLST/*.hmm ne sont pas vides!

Identification des gènes dans les génomes

Nous allons extraire les gènes MLST et les concaténer pour chaque génome.

Annotation de chaque profile HMM

/usr/local/bioinfo/bin/nhmmer -E 1e-50 -o MLST/gatB_wAU.out -A MLST/gatB_wAu.stk --noali MLST/gatB.hmm prokka/wAu/wAu.fna;

Le fichier avec les alignements est en format STK (STOCKHOLM), nous devons le transformer en format FASTA 

head MLST/gatB_wAu.stk

STK à FASTA

/usr/local/bioinfo/bin/esl-reformat --informat stockholm -o MLST/gatB_wAu.fa fasta MLST/gatB_wAu.stk;

Automatisation

Les deux étapes sont automatisées pour tous les gènes et tous les génomes. 

/home/formation/public_html/M2_Phylogenomique/scripts/hmm2genes.pl --prokka_dir prokka --mlst_dir MLST --gene_list "gatB coxA hcpA ftsZ fbpA" --verbose 0 --erase 0

Les gènes concaténés sont dans le fichier : 

MLST/concat_file.fa

Alignement et arbre

muscle: 

/usr/local/bioinfo/bin/muscle -in MLST/concat_file.fa -out MLST/concat_file.mfa

readal (fasta à phylip): 

/usr/local/bioinfo/bin/readal -in MLST/concat_file.mfa -out MLST/concat_file.phy -phylip;

PhyML: 

/usr/local/bioinfo/src/PhyML/PhyML-3.1/PhyML-3.1_linux64 -i MLST/concat_file.phy -d nt -b -4 -m HKY85 -f m -t e -v e -c 4 -a e -s NNI -o tlr --quiet --no_memory_check

Visualisation de l'arbre

Lancer l'interface graphique sur votre poste de travail de la P0. 

seaview MLST/concat_file.phy_phyml_tree.txt&

Le fichier ST_concatenated_reference.fas contient des séquence type de référence. Elles peuvent être ajoutées à votre alignement.

Comparer les arbres obtenus avec l'ARNr 16S et les gènes MLST.

Le fichier NCBI_Wolbachia.csv contient une colonne Clade_ID. Comparer cette annotation avec la topologie obtenue.

Groupes de gènes orthologues

OrthoMCL

Le guide utilisateurs sur le serveur: /usr/local/bioinfo/src/OrthoMCL/current/doc/OrthoMCLEngine/Main/UserGuide.txt

For alternate documentation online, please read Unit 6.12 of the Current Protocols of Bioinformatics available at: 

 OrthoMCL Protocols

For details on the orthomcl algorithm, please read the OrthoMCL Algorithm Document: 

 OrthoMCL Algorithm Document

En entrée, le programme prend un ensemble de protéomes.

La sortie de orthoMCL est un ensemble de fichiers : 

  pairs/
     potentialOrthologs.txt
     potentialCoorthologs.txt
     potentialInparalogs.txt
  groups.txt

Les fichiers dans le répertoire pairs contiennent les paires de relations entre les protéines associées aux scores. Comme décrit dans l'article, on distingue :

  • les orthologues potentiels
  • les co-orthologues
  • les inparalogues

Il y a trois grandes étapes:

  • blastp tous_contre_tous
  • l'identification des paires (répertoire pairs)
  • le partitionnement du graphe des pairs avec le programme MCL program (groups.txt).

Les différentes étapes du l'analyse sont lancées indépendamment mais elles s'enchainent dans un ordre bien précis. Il est possible de modifier une étape pour tester des alternatives ou optimiser l'exécution. Toutefois, il est impératif de respecter les formats des fichiers.

à modifier!

orthoMCL pipe 1

Pour faciliter les premières étapes (Step 5 et 6), elles ont été intégrées dans le script prep_orthoMLC.pl.

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

Step 5: orthomclAdjustFasta

Input: 

 - fasta files as acquired from the genome resource.

Output: 

 - the my_orthomcl_dir/compliantFasta/ directory of orthomcl-compliant fasta files (see Step 6)

Use orthomclAdjustFasta to produce a compliant file from any input file that conforms to the following pattern (for other files, provide your own script to produce complaint fasta files): 

 - has one or more fields that are separated by white space or the '|' character (optionally surrounded by white space)
 - has the unique ID in the same field of every protein.


Create an empty my_orthomcl_dir/compliantFasta/ directory, and change to that directory.

Run orthomclAdjustFasta once for each input proteome fasta file. It will produce a compliant file in the new directory. Check each file to ensure that the proteins all have proper IDs.

Step 6: orthomclFilterFasta

Input: 

 - my_orthomcl_dir/compliantFasta/ 
 - optionally a gene->protein mapping file

Output: 

 - my_orthomcl_dir/goodProteins.fasta
 - my_orthomcl_dir/poorProteins.fasta
 - report of suspicious proteomes (> 10% poor proteins)

This step produces a single goodProteins.fasta file to run BLAST on. It filters away poor-quality sequences (placing them in poorProteins.fasta). The filter is based on length and percent stop codons. You can adjust these values.

The input requirements are: 

 # a compliantFasta/ directory which contains all and only the proteome .fasta files, one file per proteome.
 # each .fasta file must have a name in the form 'xxxx.fasta' where xxxx is a three or four letter unique taxon code.  For example: hsa.fasta or eco.fasta
 # each protein in those files must have a definition line in the following format:
    >xxxx|yyyyyyyy

where xxxx is the three or four letter taxon code and yyyyyyy is a sequence identifier unique within that taxon.

Change dir to my_orthomcl_dir/ and run orthomclFilterFasta.

orthoMCL pipe 2

Les blastp tous-contre-tous (Step 7) sont réalisé savec le script ava_BLAST_orthoMCL.pl.

Pour éviter de charger le cluster avec les blastp, nous allons les réaliser séparément pour chaque génome. Chaque compte fleur va lancer les blastp sur un ou deux genomes comme query et tous les autres génomes comme data base.

Association fleur/génome 

anemone 	wAu 	wNo
arome 	wCle
aster 	wDim
bleuet 	wHa
camelia 	wMel
capucine 	wNfe 	wMeP
chardon 	wPiM
clematite 	wViB
cobee w	Cam
coquelicot 	wDci
cyclamen 	wGmo 	wPpe
dahlia 	wMun
digitale 	wObr
gerbera 	wRi
geranium 	wBma 	wWb 
glaieul 	wCqu
hortensia 	wLsi

Paramètres: 

--query_list "wMun" 
--target_list "wAu wCle wDim wHa wMel wNfe wPiM wViB wCam wDci wGmo wMun wObr wRi wBma wCqu wLsi wMeP wNo  wPpe wWb"

Copier les résultats dans un répertoire accessible à tous.

A la fin, concaténer tous les résultats blasp dans un seul fichier: 

cat wolbachia/blastp/*fasta > wolbachia/blastp/wol_21_blt

Step 7: All-v-all BLAST

Input: 

 - goodProteins.fasta

Output: 

 - your_blast_results_in_tab_format

You must run your own BLAST.

We expect you to: 

 - use NCBI BLAST
 - run with the -m 8 option to provide tab delimited output required by Step 8
 - for IMPORTANT DETAILS about other BLAST arguments, see:
   the OrthoMCL Algorithm Document

Blastp parameters: 

-F  'm S': mask with Seg
-v 100000:  a "don't care" value
-b 100000:  a "don't care" value
-z protein_database_size:   the number of proteins in the set.  S
-e 1e-5:  recommended e-value cutoff

Step 8: orthomclBlastParser

Input: 

 - your_blast_results_in_tab_format
 - my_orthomcl_dir/compliantFasta/

Output: 

 - my_orthomcl_dir/similarSequences.txt

This step parses NCBI BLAST -m 8 output into the format that can be loaded into the orthomcl database.

Use the orthomclBlastParser program for this. In addition to formatting, it computes the percent match of each hit. 

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclBlastParser wolbachia//blastp/wol_21_blt wolbachia/my_orthomcl_dir/compliantFasta >> wolbachia/my_orthomcl_dir/similarSequences.txt

head wolbachia/my_orthomcl_dir/similarSequences.txt

orthoMCL mySQL DB

OrthoMCL nécessite une base de données (mySQL par exemple). Elle a été crée par Marie-Stephane Trotard (support genotoul)

Copier le fichier modèle 

/usr/local/bioinfo/src/OrthoMCL/orthomclSoftware-v2.0.9/my_orthomcl_dir_template/orthomcl.config

sous votre work par exemple et personnaliser le ainsi

Pour créer le schema de la base: 

/usr/local/bioinfo/src/OrthoMCL/orthomclSoftware-v2.0.9/bin/orthomclInstallSchema
<path_vers_votre_orthomcl.config <path_verts_votre repertoire de
travail>/install_schema.log

/usr/local/bioinfo/src/OrthoMCL/orthomclSoftware-v2.0.9/bin/orthomclInstallSchema wolbachia/my_orthomcl_dir/orthomcl.config wolbachia/my_orthomcl_dir/install_schema.log

DROP Tables

OrthoMCL expects a "clean" and empty database after orthoMCLInstallSchema, so in between runs it's easiest to run 

mysql -h 147.99.108.14 -D orthomcl_yquentin -u ortho_yquentin -p

mysql>SHOW TABLES;

mysql>DROP TABLE IF EXISTS BestHit,  BestInterTaxonScore, BestQueryTaxonScore, BetterHit, CoOrthNotOrtholog, CoOrtholog, CoOrthologAvgScore, CoOrthologCandidate, CoOrthologTaxon, CoOrthologTemp, InParalog, InParalog2Way, InParalogAvgScore, InParalogOrtholog, InParalogTaxonAvg, InParalogTemp, InplgOrthTaxonAvg, InplgOrthoInplg, InterTaxonMatch, Ortholog, Ortholog2Way, OrthologAvgScore, OrthologTaxon, OrthologTemp, OrthologUniqueId, SimilarSequences, UniqSimSeqsQueryId;

inside the mysql command line and then to do orthoMCLinstallSchema again.

Step 9: orthomclLoadBlast

Input: 

 - similarSequences.txt

Output: 

 - SimilarSequences table in the database

This step loads the BLAST results into the orthomcl database. 

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclLoadBlast wolbachia/my_orthomcl_dir/orthomcl.config work/wolbachia/my_orthomcl_dir/similarSequences.txt

Step 10: orthomclPairs

Input: 

 - SimilarSequences table in the database

Output: 

 - PotentialOrthologs table
 - PotentialInParalogs table
 - PotentialCoOrthologs table

This is a computationally major step that finds protein pairs. It executes the algorithm described in the OrthoMCL Algorithm Document, using a relational database.

The program proceeds through a series of internal steps, each creating an intermediate database table or index. There are about 20 such tables created. Finally, it populates the output tables.

The cleanup= option allows you to control the cleaning up of the intermediary tables. T

  • 'yes' option drops the intermediary tables once they are no longer needed.
  • 'no' option keeps the intermediary tables in the database.

In total, they are expected to be about 50 percent of the SimilarSequences table. They are useful mostly for power users or developers who would like to query them. They can be removed afterwards with the 'only' or 'all' options. The latter also removes the final tables, and should only be done after Step 11 below has dumped them to files.

The startAfter= option allows you to pick up where the program left off, if it stops for any reason. Look in the log to find the last completed step, and use its tag as the value for startAfter=

Because this program will run for many hours, we recommend you run it using the UNIX 'screen' program, so that it does not abort in the middle. (If it does, use startAfter=). 

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclPairs wolbachia/my_orthomcl_dir/orthomcl.config wolbachia/my_orthomcl_dir/orthomcl_pairs.log cleanup=no

Step 11: orthomclDumpPairsFiles

Input: 

 - database with populated pairs tables

Output 

 - pairs/ directory.  
 - mclInput file

Run the orthomclDumpPairsFiles.

The pairs/ directory contains three files: ortholog.txt, coortholog.txt, inparalog.txt. Each of these has three columns: 

  - protein A
  - protein B
  - their normalized score (See the Orthomcl Algorithm Document).

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclDumpPairsFiles wolbachia/my_orthomcl_dir/orthomcl.config

head wolbachia/my_orthomcl_dir/pairs/coorthologs.txt
head wolbachia_save/my_orthomcl_dir/pairs/inparalogs.txt
head wolbachia_save/my_orthomcl_dir/pairs/orthologs.txt

Step 12: mcl

Input: 

 - mclInput file

Output: 

 - mclOutput file

/usr/local/bioinfo/bin/mcl wolbachia/my_orthomcl_dir/mclInput --abc -I 1.5 -o work/wolbachia/my_orthomcl_dir/mclOutput

Step 13: orthomclMclToGroups

Input: 

 - mclOutput file

Output: 

 - groups.txt

Change to my_orthomcl_dir and run: 

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclMclToGroups wolb 1000 < mclOutput > groups.txt
  • my_prefix is an arbitrary string to use as a prefix for your group IDs.
  • 1000 is an arbitrary starting point for your group IDs.

Step 14: orthomclSingletons

Input: 

 - groups.txt
 - goodProteins.txt

Output: 

 - singletons.txt

Change to my_orthomcl_dir and run: 

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclSingletons wolbachia/my_orthomcl_dir/good_proteins_file wolbachia/my_orthomcl_dir/groups.txt >> wolbachia/my_orthomcl_dir/singletons.txt

wc -l wolbachia/my_orthomcl_dir/singletons.txt
1639

Other utilities

Extract protein IDs from an orthomcl groups file. 

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclExtractProteinIdsFromGroupsFile wolbachia/my_orthomcl_dir/groups.txt

Extract protein ID pairss from an orthomcl groups file. 

/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclExtractProteinPairsFromGroupsFile wolbachia/my_orthomcl_dir/groups.txt

</pre>

Perl scripts

PROKKA

  • rename_ncbi_files.pl
  • prokka_loop.pl

MLST

  • fasta2hmm.pl
  • hmm2genes.pl

orthMCL

  • prep_orthoMCL.pl
Retrieved from "http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique"