Atelier Phylogénomique Analyses phylogénomiques
From silico.biotoul.fr
m (→Liens) |
m (→Comparaison des résultats de PanOCT avec OrthoFinder) |
||
(51 intermediate revisions not shown) | |||
Line 1: | Line 1: | ||
==Liens== | ==Liens== | ||
- | *[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Analyses_phylog.C3.A9nomiques | + | *retour à [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Analyses_phylog.C3.A9nomiques Atelier Phylogénomique] |
+ | ==Analyses phylogénomiques== | ||
+ | Comme dans [https://www.ncbi.nlm.nih.gov/pubmed/18159947 Kettler et al., 2007], nous allons utiliser quatre génomes de Synechococcus comme groupe externe dans nos analyses. | ||
+ | |||
+ | Mise à jour des scripts! | ||
+ | <pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | cp /home/formation/public_html/M2_Phylogenomique/scripts/* ~/work/scripts | ||
+ | </pre> | ||
+ | ===''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: [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Synechococcus/NCBI/accession_file.lst accession_file.lst] | ||
+ | |||
+ | Comme précédemment, les génomes sont renommés en utilisant un code à quatre lettres. | ||
+ | |||
+ | [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Synechococcus/NCBI/accession_labels_file.lst accession_labels_file.lst] | ||
+ | |||
+ | Les fichiers GenBank renommés sont disponibles dans le répertoire: [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Synechococcus/GenBank GenBank] | ||
+ | |||
+ | Les séquences ADN ont été extraites des fichiers GenBank : [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Synechococcus/DNA DNA] | ||
+ | |||
+ | Fichier avec les informations sur les souches: [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Synechococcus/NCBI/species_strain_names.txt species_strain_names.txt] | ||
+ | |||
+ | Nous allons utiliser la même procédure que pour les génomes de Prochlorococcus ([http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Prokka#Automatisation_des_annotations_prokka_sur_l.27ensemble_des_g.C3.A9nomes Automatisation des annotations prokka sur l'ensemble des génomes]) | ||
+ | |||
+ | <source lang='bash'> | ||
+ | mkdir -p ~/work/Synechococcus/prokka | ||
+ | cd ~/work/Synechococcus | ||
+ | ~/work/scripts/prokka_loop.pl --sample Synechococcus | ||
+ | |||
+ | squeue -l -u $USER | ||
+ | </source> | ||
+ | Une fois les jobs terminés, vérifiez que les fichiers de sorties de prokka existent et ne sont pas vides. | ||
+ | <source lang='bash'> | ||
+ | ls -l ~/work/Synechococcus/prokka/Aaa*/*.faa | ||
+ | </source> | ||
+ | |||
+ | ====''Synechococcus'' blastp All-All==== | ||
+ | Nous allons reproduire le même enchaînement de scipts utilisés avec ''Prochlorococcus'' ([http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Blast#Pr.C3.A9liminaires Atelier_Phylogénomique_Blast Préliminaires]) en prenant soin de remplacer Prochlorococcus par Synechococcus dans les noms des répertoires. | ||
+ | |||
+ | Copiez les fichier *.faa de prokka dans ~/work/Synechococcus/peptide/ | ||
+ | |||
+ | MSK | ||
+ | <source lang='bash'> | ||
+ | mkdir -p ~/work/Synechococcus/peptide | ||
+ | cd ~/work/Synechococcus/peptide | ||
+ | cp ~/work/Synechococcus/prokka/Aaa*/*.faa ~/work/Synechococcus/peptide/. | ||
+ | |||
+ | ls -l ~/work/Synechococcus/peptide | ||
+ | </source> | ||
+ | |||
+ | Créer la ''blast database''. | ||
+ | |||
+ | MSK | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | for i in ~/work/Synechococcus/peptide/*.faa; | ||
+ | do | ||
+ | echo "module load bioinfo/ncbi-blast-2.7.1+; makeblastdb -in $i -dbtype prot;" | ||
+ | done > makeblastdb.sh | ||
+ | |||
+ | cat makeblastdb.sh | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | <source lang='bash'> | ||
+ | sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 makeblastdb.sh | ||
+ | squeue -l -u $USER | ||
+ | </source> | ||
+ | |||
+ | Boucle sur les genomes de ''Synechococcus''. | ||
+ | |||
+ | MSK | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | mkdir ~/work/Synechococcus/BlastP | ||
+ | |||
+ | evalue=1e-5 | ||
+ | dbsize=100000000 | ||
+ | for i in ~/work/Synechococcus/peptide/*.faa; | ||
+ | do | ||
+ | ip=$(basename $i .faa) | ||
+ | for j in ~/work/Synechococcus/peptide/*.faa; | ||
+ | do | ||
+ | jp=$(basename $j .faa) | ||
+ | outfile="~/work/Synechococcus/BlastP/"$ip"_"$jp".tab" | ||
+ | echo "module load bioinfo/ncbi-blast-2.7.1+; blastp -query $i -db $j -seg yes -dbsize $dbsize -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile;" | ||
+ | done | ||
+ | done > blast_allall.sh | ||
+ | |||
+ | cat blast_allall.sh | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | <source lang='bash'> | ||
+ | sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh | ||
+ | squeue -l -u $USER | ||
+ | </source> | ||
+ | |||
+ | <!-- | ||
+ | <pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | module load bioinfo/ncbi-blast-2.7.1+ | ||
+ | for i in *.faa; | ||
+ | do | ||
+ | srun -n1 -l makeblastdb -in $i -dbtype prot; | ||
+ | done | ||
+ | </pre> | ||
+ | <pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | mkdir -p ~/work/Synechococcus/BlastP | ||
+ | cd ~/work/Synechococcus | ||
+ | ~/work/scripts/blastp_intra.pl | ||
+ | </pre> | ||
+ | <pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | ~/work/scripts/blastp_inter.pl | ||
+ | </pre> | ||
+ | --> | ||
+ | |||
+ | <source lang='bash'> | ||
+ | squeue -l -u $USER | ||
+ | ls ~/work/Synechococcus/BlastP | wc -l | ||
+ | </source> | ||
+ | |||
+ | ===''Prochlorococcus'' ''versus'' ''Synechococcus''=== | ||
+ | ====Liens symboliques sur les fichiers peptides==== | ||
+ | <source lang='bash'> | ||
+ | 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/. | ||
+ | </source> | ||
+ | |||
+ | ====Liens symboliques sur les fichiers blastp==== | ||
+ | <source lang='bash'> | ||
+ | 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/. | ||
+ | </source> | ||
+ | |||
+ | ====Compléter les paires de comparaisons==== | ||
+ | MSK | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | evalue=1e-5 | ||
+ | dbsize=100000000 | ||
+ | for i in ~/work/Synechococcus/peptide/*.faa; | ||
+ | do | ||
+ | ip=$(basename $i .faa) | ||
+ | for j in ~/work/Prochlorococcus/peptide/*.faa; | ||
+ | do | ||
+ | jp=$(basename $j .faa) | ||
+ | outfile="~/work/ProchlorococcusSynechococcus/BlastP/"$ip"_"$jp".tab" | ||
+ | echo "module load bioinfo/ncbi-blast-2.7.1+; blastp -query $i -db $j -seg yes -dbsize $dbsize -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile;" | ||
+ | |||
+ | outfile="~/work/ProchlorococcusSynechococcus/BlastP/"$jp"_"$ip".tab" | ||
+ | echo "module load bioinfo/ncbi-blast-2.7.1+; blastp -query $j -db $i -seg yes -dbsize $dbsize -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile;" | ||
+ | done | ||
+ | done > blast_allall.sh | ||
+ | |||
+ | cat blast_allall.sh | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | <source lang='bash'> | ||
+ | sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh | ||
+ | squeue -l -u $USER | ||
+ | </source> | ||
+ | Vérifiez que vous avez bien le nombre de fichiers attendus! | ||
+ | |||
+ | ===Groupes de gènes orthologues avec PanOCT=== | ||
+ | ====Préparation des fichiers ''combined''==== | ||
+ | <source lang='bash'> | ||
+ | srun --pty bash | ||
+ | mkdir -p ~/work/Synechococcus/panoct/results | ||
+ | mkdir -p ~/work/ProchlorococcusSynechococcus/panoct/results | ||
+ | </source> | ||
+ | '''combined.att''' | ||
+ | Créer un fichier avec les coordonnées, noms, fonction et souches des gènes. | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | cd ~/work/Synechococcus | ||
+ | for file in peptide/*.faa | ||
+ | do | ||
+ | prefix=$(basename $file .faa) | ||
+ | echo "~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output prokka/$prefix/$prefix.tab" | ||
+ | ~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output panoct/results/$prefix.tab | ||
+ | done | ||
+ | squeue -l -u $USER | ||
+ | ls panoct/results/*.tab | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | <source lang='bash'> | ||
+ | cat ~/work/Prochlorococcus/panoct/results/*.tab ~/work/Synechococcus/panoct/results/*.tab> ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att | ||
+ | head ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att | ||
+ | </source> | ||
+ | '''genomes.list''' | ||
+ | Liste des souches à analyser. | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | cd ~/work/ProchlorococcusSynechococcus/ | ||
+ | for i in peptide/*.faa | ||
+ | do | ||
+ | echo $(basename $i .faa) | ||
+ | done > panoct/results/genomes.list | ||
+ | more panoct/results/genomes.list | ||
+ | </syntaxhighlight> | ||
+ | '''combined.fasta''' | ||
+ | Concaténer les fichier peptides dans un seul fichier. | ||
+ | <source lang='bash'> | ||
+ | cat peptide/*.faa > panoct/results/combined.fasta | ||
+ | |||
+ | grep -c '>' panoct/results/combined.fasta | ||
+ | </source> | ||
+ | '''combined.blast''' | ||
+ | Concaténer les résultats des blastp dans un seul fichier. | ||
+ | <source lang='bash'> | ||
+ | cat BlastP/*.tab > panoct/results/combined.blast | ||
+ | |||
+ | head panoct/results/combined.blast | ||
+ | </source> | ||
+ | |||
+ | ====run panOCT Prochlorococcus vs Synechococcus==== | ||
+ | Exemple de script "[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/scripts/panoct_PS.csh panoct_PS.csh]" (les chemins sont à changer). | ||
+ | <source lang='bash'> | ||
+ | cd ~/work/ProchlorococcusSynechococcus/panoct | ||
+ | mkdir ~/work/ProchlorococcusSynechococcus/images | ||
+ | </source> | ||
+ | <source lang='bash'> | ||
+ | sbatch panoct_PS.csh | ||
+ | squeue -l -u $USER | ||
+ | </source> | ||
+ | <!-- | ||
+ | /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 | ||
+ | </pre> | ||
+ | --> | ||
+ | <pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | Question 3.1: | ||
+ | Résumez dans une figure les différentes étapes réalisées jusqu'à l’obtention des groupes de gènes orthologues. N'oubliez pas la légende! | ||
+ | </pre> | ||
+ | Heatmap des groupes de gènes orthologues avec le fichier ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable_0_1.txt | ||
+ | |||
+ | MSK | ||
+ | <syntaxhighlight lang="python"> | ||
+ | srun --pty bash # si nécessaire! | ||
+ | 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") | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | *[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Figures/matchtable_heatmap.pdf matchtable_heatmap.pdf] | ||
+ | |||
+ | <pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | Question 3.2: | ||
+ | Décrivez la figure obtenue. | ||
+ | Quelles informations vous apporte-t-elle? | ||
+ | </pre> | ||
+ | |||
+ | ====Gènes espèces spécifique==== | ||
+ | =====Gènes trouvés chez toutes les souches de ''Prochlorococcus'' mais dans aucune de ''Synechococcus'':===== | ||
+ | Vous allez utiliser un script ''awk'' pour réaliser ce filtrage. | ||
+ | |||
+ | *Editez le fichier [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Figures/matchtable.txt matchtable.txt] pour comprendre sa structure. | ||
+ | *Que représente chaque ligne? | ||
+ | *Comment sont organisées les colonnes? | ||
+ | *Identifiez le motif signalant l'absence d'un gène dans un OG. | ||
+ | *Pour chaque OG: | ||
+ | **Compter le nombre de gènes présents chez les ''Prochlorococcus''. | ||
+ | **Compter le nombre de gènes absents chez les ''Synechococcus''. | ||
+ | **Editer la ligne si elle remplit les conditions | ||
+ | |||
+ | MSK | ||
+ | <syntaxhighlight lang="awk"> | ||
+ | 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"); | ||
+ | } | ||
+ | }' < ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt > ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt | ||
+ | </syntaxhighlight> | ||
+ | |||
+ | Quelles-sont les fonctions des gènes retenus? | ||
+ | |||
+ | '''Prokka''' | ||
+ | <pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | rm ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt | ||
+ | |||
+ | awk '{ system("grep "$3" ~/work/Prochlorococcus/prokka/Aaaa/Aaaa.gff >> ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt")}' < ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt | ||
+ | |||
+ | cat ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt | ||
+ | </pre> | ||
+ | |||
+ | '''eggNOG''' | ||
+ | <pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | rm ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_eggNOG_annotations.txt | ||
+ | |||
+ | awk '{ system("grep "$3" /home/formation/work/Prochlorococcus/eggNOG/Aaaa.emapper.annotations >> ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_eggNOG_annotations.txt")}' < ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt | ||
+ | |||
+ | cat ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_eggNOG_annotations.txt | ||
+ | </pre> | ||
+ | |||
+ | =====Gènes trouvés chez toutes les souches de ''Synechococcus'' mais dans aucune de ''Prochlorococcus'':===== | ||
+ | |||
+ | MSK | ||
+ | <!-- | ||
+ | <syntaxhighlight lang="awk"> | ||
+ | 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==0 && sp ==4 ) { | ||
+ | h++; | ||
+ | printf("%d\t%d", h, $1); | ||
+ | for (i=2; i<=17; i++) { | ||
+ | if( $i !~ /-/ ) printf("\t%s", $i); | ||
+ | } | ||
+ | printf("\n"); | ||
+ | } | ||
+ | }' < ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt | ||
+ | </syntaxhighlight> | ||
+ | --> | ||
+ | <pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> | ||
+ | Question 3.3: | ||
+ | Analysez ces résultats et les confronter à ceux disponibles dans la littérature. | ||
+ | </pre> | ||
+ | Vous pouvez vous inspirer des articles de Kettler et al. (2007), Garczarek et al. 2007, et Partensky and Garczarek (2010). | ||
+ | |||
+ | ====Comparaison des résultats de PanOCT avec et sans les génomes de ''Synechococcus''==== | ||
+ | Ces pairs sont à extraire des fichiers: | ||
+ | *~/work/Prochlorococcus/panoct/results/matchtable.txt | ||
+ | *~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt | ||
+ | |||
+ | A laide d'un script ''awk'', pour les deux fichiers | ||
+ | *pour chaque OG, | ||
+ | **sélectionner les gènes de Prochlorococcus, | ||
+ | **les assembler par paires et les imprimer. | ||
+ | |||
+ | MSK | ||
+ | <syntaxhighlight lang="awk"> | ||
+ | awk ' | ||
+ | { | ||
+ | col[""]=0; | ||
+ | j=0; | ||
+ | for (i=2; i<=13; i++) { | ||
+ | if( $i !~ /-/ ) { | ||
+ | j++; | ||
+ | col[j]=$i; | ||
+ | } | ||
+ | } | ||
+ | ncol=j; | ||
+ | if ( ncol > 1 ) { | ||
+ | #printf("%d\n", $1); | ||
+ | for (k=1; k<=ncol-1; k++) { | ||
+ | for (m=k+1; m<=ncol; m++) { | ||
+ | printf("%s%s\n", col[k], col[m]); | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | delete col #printf("\n"); | ||
+ | }' ~/work/Prochlorococcus/panoct/results/matchtable.txt > ~/work/ProchlorococcusSynechococcus/Ppairs.txt | ||
+ | |||
+ | awk ' | ||
+ | { | ||
+ | col[""]=0; | ||
+ | j=0; | ||
+ | for (i=2; i<=13; i++) { | ||
+ | if( $i !~ /-/ ) { | ||
+ | j++; | ||
+ | col[j]=$i; | ||
+ | } | ||
+ | } | ||
+ | ncol=j; | ||
+ | if ( ncol > 1 ) { | ||
+ | #printf("%d\n", $1); | ||
+ | for (k=1; k<=ncol-1; k++) { | ||
+ | for (m=k+1; m<=ncol; m++) { | ||
+ | printf("%s%s\n", col[k], col[m]); | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | delete col | ||
+ | #printf("\n"); | ||
+ | }' ~/work/ProchlorococcusSynechococcus//panoct/results/matchtable.txt > ~/work/ProchlorococcusSynechococcus/PSpairs.txt | ||
+ | </syntaxhighlight> | ||
+ | La comparaison peut être réalisée graphiquement avec un digramme de Venn. | ||
+ | |||
+ | MSK | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | srun --pty bash | ||
+ | module load system/R-3.5.1 | ||
+ | R | ||
+ | library('gplots') | ||
+ | setwd('~/work/ProchlorococcusSynechococcus') | ||
+ | pdf_file <- 'ProchlorococcusSynechococcus_Venn.pdf' | ||
+ | Ppairs <- read.table('Ppairs.txt') | ||
+ | PSpairs <- read.table('PSpairs.txt') | ||
+ | input <- list(P=Ppairs, PS=PSpairs) | ||
+ | |||
+ | pdf(file=pdf_file, paper="a4r") | ||
+ | venn(input) | ||
+ | dev.off() | ||
+ | |||
+ | a<-venn(input, show.plot=FALSE) | ||
+ | intersections <- attr(a,"intersections") | ||
+ | intersections$P | ||
+ | intersections$PS | ||
+ | </syntaxhighlight> | ||
+ | *[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Figures/ProchlorococcusSynechococcus_Venn.pdf ProchlorococcusSynechococcus_Venn.pdf] | ||
+ | |||
+ | ====Comparaison des résultats de PanOCT avec OrthoFinder==== | ||
+ | <!-- | ||
+ | OrthoFinder File: ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/WorkingDirectory/OrthoFinder/Results_ProSyn/Phylogenetic_Hierarchical_Orthogroups/N0.tsv | ||
+ | --> | ||
+ | MSK | ||
+ | <!-- | ||
+ | cat ~/work/OrthoFinder/Prochlorococcus/OrthoFinder/Results_Pro/WorkingDirectory/OrthoFinder/Results_ProSyn/Phylogenetic_Hierarchical_Orthogroups/N0.tsv | awk -F $'\t' '{ | ||
+ | if ( $1~/N0.HOG/ ) { | ||
+ | j=0; | ||
+ | for (i=4; i<=NF; i++) { | ||
+ | if ( $i~/[A-Z]/ ) { | ||
+ | split($i, paral, ", "); | ||
+ | for ( k in paral ) { | ||
+ | if ( paral[k]~/[A-Z]/ ) { | ||
+ | j++; | ||
+ | col[j]=paral[k]; | ||
+ | } | ||
+ | } | ||
+ | delete paral; | ||
+ | } | ||
+ | } | ||
+ | |||
+ | ncol=j; | ||
+ | if ( ncol > 1 ) { | ||
+ | n=asort(col) | ||
+ | for (k=1; k<=ncol; k++) { | ||
+ | for (m=k+1; m<=ncol; m++) { | ||
+ | printf("%s-%s\n", col[k], col[m]); | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | delete col | ||
+ | } | ||
+ | }' > ~/work/OrthoFinder/Figures/ProSynHOGN0_pairs.txt | ||
+ | |||
+ | PorthoMCL | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | awk ' | ||
+ | { | ||
+ | j=0; | ||
+ | for (i=1; i<=NF; i++) { | ||
+ | j++; | ||
+ | col[j]=$i; | ||
+ | } | ||
+ | ncol=j; | ||
+ | if ( ncol > 1 ) { | ||
+ | for (k=1; k<=ncol; k++) { | ||
+ | for (m=k+1; m<=ncol; m++) { | ||
+ | printf("%s%s\n", col[k], col[m]); | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | delete col | ||
+ | }' ~/work/Prochlorococcus/PorthoMCL/8.all.ort.group > ~/work/Prochlorococcus/PorthoMCL/8.all.ort.pairs.txt | ||
+ | sed -i -r 's/\|//g' ~/work/Prochlorococcus/PorthoMCL/8.all.ort.pairs.txt | ||
+ | </syntaxhighlight> | ||
+ | --> | ||
+ | |||
+ | MSK | ||
+ | <!-- | ||
+ | <syntaxhighlight lang="bash"> | ||
+ | srun --pty bash | ||
+ | module load system/R-3.5.1 | ||
+ | |||
+ | R | ||
+ | library('gplots') | ||
+ | |||
+ | #PorthoMCL <- read.table('~/work/Prochlorococcus/PorthoMCL/8.all.ort.pairs.txt') | ||
+ | orthofinder <- read.table('~/work/OrthoFinder/Figures/ProSynHOGN0_pairs.txt') | ||
+ | panoct <- read.table('~/work/ProchlorococcusSynechococcus/PSpairs.txt') | ||
+ | input <- list(panOCT=panoct, OrthoFinder=orthofinder) | ||
+ | |||
+ | pdf_file <- '~/work/ProchlorococcusSynechococcus/panOCT_OrthoFinder_Venn_ProSyn.pdf' | ||
+ | pdf(file=pdf_file, paper="a4r") | ||
+ | venn(input) | ||
+ | dev.off() | ||
+ | |||
+ | orthofinder <- read.table('~/work/OrthoFinder/Figures/HOGN0_pairs.txt') | ||
+ | panoct <- read.table('~/work/ProchlorococcusSynechococcus/Ppairs.txt') | ||
+ | input <- list(panOCT=panoct, OrthoFinder=orthofinder) | ||
+ | |||
+ | pdf_file <- '~/work/ProchlorococcusSynechococcus/panOCT_OrthoFinder_Venn_Pro.pdf' | ||
+ | pdf(file=pdf_file, paper="a4r") | ||
+ | venn(input) | ||
+ | dev.off() | ||
+ | |||
+ | orthogroups <- read.table('~/work/OrthoFinder/Figures/Orthogroups_pairs.txt') | ||
+ | panoct <- read.table('~/work/ProchlorococcusSynechococcus/Ppairs.txt') | ||
+ | input <- list(panOCT=panoct, OrthoFinder=orthogroups) | ||
+ | |||
+ | pdf_file <- '~/work/ProchlorococcusSynechococcus/panOCT_OGOrthoFinder_Venn_Pro.pdf' | ||
+ | pdf(file=pdf_file, paper="a4r") | ||
+ | venn(input) | ||
+ | dev.off() | ||
+ | |||
+ | |||
+ | </syntaxhighlight> | ||
+ | --> | ||
+ | <!-- | ||
+ | *[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Figures/panOCT_OrthoFinder_Venn_ProSyn.pdf panOCT_OrthoFinder_Venn_ProSyn.pdf] | ||
+ | *[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Figures/panOCT_OrthoFinder_Venn_Pro.pdf panOCT_OrthoFinder_Venn_Pro.pdf] | ||
+ | *[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Figures/panOCT_OGOrthoFinder_Venn_Pro.pdf panOCT_OGOrthoFinder_Venn_Pro.pdf] | ||
+ | --> | ||
---- | ---- | ||
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Analyses_phylog.C3.A9nomiques Analyses phylogénomiques] | *[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Analyses_phylog.C3.A9nomiques Analyses phylogénomiques] |
Current revision as of 10:06, 7 December 2022
Contents
|
Liens
- retour à Atelier Phylogénomique
Analyses phylogénomiques
Comme dans Kettler et al., 2007, nous allons utiliser quatre génomes de Synechococcus comme groupe externe dans nos analyses.
Mise à jour des scripts!
cp /home/formation/public_html/M2_Phylogenomique/scripts/* ~/work/scripts
Synechococcus
Génomes et annotations
- accession "CP000435, CP000097, BX548020, CP000110"
Les fichiers GenBank peuvent être obtenus avec la commande wget en utilisant les chemins enregistrés dans le fichier: accession_file.lst
Comme précédemment, les génomes sont renommés en utilisant un code à quatre lettres.
Les fichiers GenBank renommés sont disponibles dans le répertoire: GenBank
Les séquences ADN ont été extraites des fichiers GenBank : DNA
Fichier avec les informations sur les souches: species_strain_names.txt
Nous allons utiliser la même procédure que pour les génomes de Prochlorococcus (Automatisation des annotations prokka sur l'ensemble des génomes)
mkdir -p ~/work/Synechococcus/prokka cd ~/work/Synechococcus ~/work/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 utilisés avec Prochlorococcus (Atelier_Phylogénomique_Blast Préliminaires) en prenant soin de remplacer Prochlorococcus par Synechococcus dans les noms des répertoires.
Copiez les fichier *.faa de prokka dans ~/work/Synechococcus/peptide/
MSK
mkdir -p ~/work/Synechococcus/peptide cd ~/work/Synechococcus/peptide cp ~/work/Synechococcus/prokka/Aaa*/*.faa ~/work/Synechococcus/peptide/. ls -l ~/work/Synechococcus/peptide
Créer la blast database.
MSK
for i in ~/work/Synechococcus/peptide/*.faa; do echo "module load bioinfo/ncbi-blast-2.7.1+; makeblastdb -in $i -dbtype prot;" done > makeblastdb.sh cat makeblastdb.sh
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 makeblastdb.sh squeue -l -u $USER
Boucle sur les genomes de Synechococcus.
MSK
mkdir ~/work/Synechococcus/BlastP evalue=1e-5 dbsize=100000000 for i in ~/work/Synechococcus/peptide/*.faa; do ip=$(basename $i .faa) for j in ~/work/Synechococcus/peptide/*.faa; do jp=$(basename $j .faa) outfile="~/work/Synechococcus/BlastP/"$ip"_"$jp".tab" echo "module load bioinfo/ncbi-blast-2.7.1+; blastp -query $i -db $j -seg yes -dbsize $dbsize -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile;" done done > blast_allall.sh cat blast_allall.sh
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh squeue -l -u $USER
squeue -l -u $USER ls ~/work/Synechococcus/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 paires de comparaisons
MSK
evalue=1e-5 dbsize=100000000 for i in ~/work/Synechococcus/peptide/*.faa; do ip=$(basename $i .faa) for j in ~/work/Prochlorococcus/peptide/*.faa; do jp=$(basename $j .faa) outfile="~/work/ProchlorococcusSynechococcus/BlastP/"$ip"_"$jp".tab" echo "module load bioinfo/ncbi-blast-2.7.1+; blastp -query $i -db $j -seg yes -dbsize $dbsize -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile;" outfile="~/work/ProchlorococcusSynechococcus/BlastP/"$jp"_"$ip".tab" echo "module load bioinfo/ncbi-blast-2.7.1+; blastp -query $j -db $i -seg yes -dbsize $dbsize -evalue $evalue -outfmt 6 -num_threads 1 -out $outfile;" done done > blast_allall.sh cat blast_allall.sh
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh squeue -l -u $USER
Vérifiez que vous avez bien le nombre de fichiers attendus!
Groupes de gènes orthologues avec PanOCT
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 "~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output prokka/$prefix/$prefix.tab" ~/work/scripts/prokkagff2panoct.pl --gffdir prokka/$prefix --output panoct/results/$prefix.tab done squeue -l -u $USER 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
Exemple de script "panoct_PS.csh" (les chemins sont à changer).
cd ~/work/ProchlorococcusSynechococcus/panoct mkdir ~/work/ProchlorococcusSynechococcus/images
sbatch panoct_PS.csh squeue -l -u $USER
Question 3.1: Résumez dans une figure les différentes étapes réalisées jusqu'à l’obtention des groupes de gènes orthologues. N'oubliez pas la légende!
Heatmap des groupes de gènes orthologues avec le fichier ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable_0_1.txt
MSK
srun --pty bash # si nécessaire! 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")
Question 3.2: Décrivez la figure obtenue. Quelles informations vous apporte-t-elle?
Gènes espèces spécifique
Gènes trouvés chez toutes les souches de Prochlorococcus mais dans aucune de Synechococcus:
Vous allez utiliser un script awk pour réaliser ce filtrage.
- Editez le fichier matchtable.txt pour comprendre sa structure.
- Que représente chaque ligne?
- Comment sont organisées les colonnes?
- Identifiez le motif signalant l'absence d'un gène dans un OG.
- Pour chaque OG:
- Compter le nombre de gènes présents chez les Prochlorococcus.
- Compter le nombre de gènes absents chez les Synechococcus.
- Editer la ligne si elle remplit les conditions
MSK
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"); } }' < ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt > ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt
Quelles-sont les fonctions des gènes retenus?
Prokka
rm ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt awk '{ system("grep "$3" ~/work/Prochlorococcus/prokka/Aaaa/Aaaa.gff >> ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt")}' < ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt cat ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt
eggNOG
rm ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_eggNOG_annotations.txt awk '{ system("grep "$3" /home/formation/work/Prochlorococcus/eggNOG/Aaaa.emapper.annotations >> ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_eggNOG_annotations.txt")}' < ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt cat ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_eggNOG_annotations.txt
Gènes trouvés chez toutes les souches de Synechococcus mais dans aucune de Prochlorococcus:
MSK
Question 3.3: Analysez ces résultats et les confronter à ceux disponibles dans la littérature.
Vous pouvez vous inspirer des articles de Kettler et al. (2007), Garczarek et al. 2007, et Partensky and Garczarek (2010).
Comparaison des résultats de PanOCT avec et sans les génomes de Synechococcus
Ces pairs sont à extraire des fichiers:
- ~/work/Prochlorococcus/panoct/results/matchtable.txt
- ~/work/ProchlorococcusSynechococcus/panoct/results/matchtable.txt
A laide d'un script awk, pour les deux fichiers
- pour chaque OG,
- sélectionner les gènes de Prochlorococcus,
- les assembler par paires et les imprimer.
MSK
awk ' { col[""]=0; j=0; for (i=2; i<=13; i++) { if( $i !~ /-/ ) { j++; col[j]=$i; } } ncol=j; if ( ncol > 1 ) { #printf("%d\n", $1); for (k=1; k<=ncol-1; k++) { for (m=k+1; m<=ncol; m++) { printf("%s%s\n", col[k], col[m]); } } } delete col #printf("\n"); }' ~/work/Prochlorococcus/panoct/results/matchtable.txt > ~/work/ProchlorococcusSynechococcus/Ppairs.txt awk ' { col[""]=0; j=0; for (i=2; i<=13; i++) { if( $i !~ /-/ ) { j++; col[j]=$i; } } ncol=j; if ( ncol > 1 ) { #printf("%d\n", $1); for (k=1; k<=ncol-1; k++) { for (m=k+1; m<=ncol; m++) { printf("%s%s\n", col[k], col[m]); } } } delete col #printf("\n"); }' ~/work/ProchlorococcusSynechococcus//panoct/results/matchtable.txt > ~/work/ProchlorococcusSynechococcus/PSpairs.txt
La comparaison peut être réalisée graphiquement avec un digramme de Venn.
MSK
srun --pty bash module load system/R-3.5.1 R library('gplots') setwd('~/work/ProchlorococcusSynechococcus') pdf_file <- 'ProchlorococcusSynechococcus_Venn.pdf' Ppairs <- read.table('Ppairs.txt') PSpairs <- read.table('PSpairs.txt') input <- list(P=Ppairs, PS=PSpairs) pdf(file=pdf_file, paper="a4r") venn(input) dev.off() a<-venn(input, show.plot=FALSE) intersections <- attr(a,"intersections") intersections$P intersections$PS
Comparaison des résultats de PanOCT avec OrthoFinder
MSK
MSK