silico.biotoul.fr
 

Atelier Phylogénomique Arbre espèces

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Super-arbres)
m (Evaluation des alignements avec t_coffee)
 
(19 intermediate revisions not shown)
Line 1: Line 1:
==Liens==
==Liens==
-
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces:_pr.C3.A9paration_des_fichiers Arbre espèces]
+
*retour à [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces:_pr.C3.A9paration_des_fichiers Atelier Phylogénomique]
==Introduction==
==Introduction==
Line 10: Line 10:
Si le répertoire n'existe pas déjà, le créer.
Si le répertoire n'existe pas déjà, le créer.
-
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
+
<source lang='bash'>
mkdir -p ~/work/ProchlorococcusSynechococcus/DNA
mkdir -p ~/work/ProchlorococcusSynechococcus/DNA
cat ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.ffn ~/work/Synechococcus/prokka/Aaa*/Aaa*.ffn > ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn
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
grep -c '>' ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn
-
</pre>
+
</source>
===Extraire les groupes de gènes orthologues===
===Extraire les groupes de gènes orthologues===
Exemple de la création d'un script et lancement du job avec sbatch
Exemple de la création d'un script et lancement du job avec sbatch
-
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
+
<source lang='bash'>
mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA
mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA
-
</pre>
+
</source>
Créer le fichier suivant ''commande.sh''.
Créer le fichier suivant ''commande.sh''.
-
<pre>
+
<source lang='bash'>
#!/usr/bin/bash
#!/usr/bin/bash
~/work/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
~/work/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
-
</pre>
+
</source>
-
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">sbatch commande.sh
+
<source lang='bash'>
-
</pre>
+
sbatch commande.sh
 +
</source>
Vérifier que les fichiers contiennent au moins 16 séquences:
Vérifier que les fichiers contiennent au moins 16 séquences:
-
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
+
<source lang='bash'>
grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
-
</pre>
+
</source>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.1:
Question 5.1:
Line 39: Line 40:
Et quand les espèces étudiées sont très éloignées phylogénétiquement parlant ?
Et quand les espèces étudiées sont très éloignées phylogénétiquement parlant ?
</pre>
</pre>
 +
===Transformation en alignement multiples des séquences ADN des gènes===
===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 ([https://bioperl.org/ BioPerl]).
+
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 (merci [https://bioperl.org/ BioPerl]).
-
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
+
<source lang='bash'>
mkdir ~/work/ProchlorococcusSynechococcus/OG/alignment
mkdir ~/work/ProchlorococcusSynechococcus/OG/alignment
-
</pre>
+
</source>
Attention à faire un script et le lancer en sbatch ou en srun --pty bash
Attention à faire un script et le lancer en sbatch ou en srun --pty bash
-
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> ~/work/scripts/aa_to_dna_aln.pl -dna ~/work/ProchlorococcusSynechococcus/OG/DNA/OG_1471.fas --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment
+
<source lang='bash'>
-
</pre>
+
~/work/scripts/aa_to_dna_aln.pl -dna ~/work/ProchlorococcusSynechococcus/OG/DNA/OG_1471.fas --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment
 +
</source>
La traduction en peptides ajoute un '*' à la fin des séquences. Vous pouvez ignorer le message:
La traduction en peptides ajoute un '*' à la fin des séquences. Vous pouvez ignorer le message:
Line 62: Line 65:
Pour monitorer vos jobs : squeue -u $USER
Pour monitorer vos jobs : squeue -u $USER
-
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
+
<source lang='bash'>
~/work/scripts/aa_to_dna_aln_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100
~/work/scripts/aa_to_dna_aln_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100
-
</pre>
+
</source>
Comptez les nombre d'alignements obtenus.  
Comptez les nombre d'alignements obtenus.  
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Line 74: Line 77:
</pre>
</pre>
 +
MSK
 +
<!--
==Construction des arbres des groupes de gènes orthologues==
==Construction des arbres des groupes de gènes orthologues==
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Line 81: Line 86:
Connaissez-vous d'autres logiciels réalisant un traitement comparable des alignements multiples?
Connaissez-vous d'autres logiciels réalisant un traitement comparable des alignements multiples?
</pre>
</pre>
-
===Evaluation des alignements avec t_coffee===
+
-->
 +
===Evaluation des alignements avec ''t_coffee''===
 +
*[https://tcoffee.readthedocs.io/en/latest/tcoffee_main_documentation.html?highlight=score#transitive-consistency-score-tcs Transitive-Consistency-Score]
 +
TCS is an alignment evaluation score that makes it possible to identify the most correct positions in an MSA. It has been shown that these positions are the most likely to be structuraly correct and also the most informative when estimating phylogenetic trees. The TCS '''evaluation''' and '''filtering procedure''' can be used to evaluate and filter any MSA.
 +
sample_seq1.score_ascii displays the score of the MSA, the sequences and the residues.
 +
sample_seq1.score_html displays a colored version score of the MSA, the sequences and the residues.
 +
TCS allows you to filter out from your alignment regions that appears unreliable according to the consistency score; the filtering can be made at the residue level or the column level:
 +
sample_seq1.tcs_residue_filter3 All residues with a TCS score lower than 3 are filtered out
 +
sample_seq1.tcs_column_filter3 All columns with a TCS score lower than 3 are filtered out
 +
 
<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
-
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_*.fas  
+
original_alignments=~/work/ProchlorococcusSynechococcus/OG/alignment/
 +
for i in $original_alignments/ali_dna_OG_*.fas  
do  
do  
   ip=$(basename $i .aln)     
   ip=$(basename $i .aln)     
Line 98: Line 113:
Sélection alignements avec un SCORE=1000.
Sélection alignements avec un SCORE=1000.
<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
 +
goodDNA_alignments=~/work/ProchlorococcusSynechococcus/OG/goodDNA_alignments/
 +
goodPEP_alignments=~/work/ProchlorococcusSynechococcus/OG/goodPEP_alignments/
 +
original_alignments=~/work/ProchlorococcusSynechococcus/OG/alignment/
 +
mkdir $goodDNA_alignments $goodPEP_alignments
 +
 +
goodDNA_alignment_list=~/work/ProchlorococcusSynechococcus/OG/goodDNA_alignments/DNASCORE1000.lst
 +
goodPEP_alignment_list=~/work/ProchlorococcusSynechococcus/OG/goodPEP_alignments/PEPSCORE1000.lst
 +
 +
if [[ -f $goodDNA_alignment_list ]]; then
 +
rm $goodDNA_alignment_list
 +
fi
 +
 +
if [[ -f $goodPEP_alignment_list ]]; then
 +
rm $goodPEP_alignment_list
 +
fi
 +
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/*.ali;  
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/*.ali;  
do  
do  
-
  ip=$(basename $i .aln)   
+
  if [[ $(grep -H 'SCORE=1000' $i) ]]; then
-
grep -H 'SCORE=1000' $i
+
  ip=$(basename $i .ali)
 +
  oriDNA="$original_alignments$ip"
 +
  goodDNA="$goodDNA_alignments$ip"
 +
  cp $oriDNA $goodDNA
 +
  echo $oriDNA >> $goodDNA_alignment_list
 +
 
 +
  PEP=${ip:11:15}
 +
  oriPEP="${original_alignments}ali_pep_OG_${PEP}"
 +
  goodPEP="${goodPEP_alignments}ali_pep_OG_${PEP}"
 +
  cp $oriPEP $goodPEP
 +
  echo $oriPEP >> $goodPEP_alignment_list
 +
fi
done
done
 +
cat $goodDNA_alignment_list
</syntaxhighlight>
</syntaxhighlight>
 +
Edition des alignements avec trimal: MSK
 +
<!--
===Edition des alignements avec trimal===
===Edition des alignements avec trimal===
-
MSK
 
-
<!--
 
Trimer les alignements avec trimal. Nous pouvons utiliser les résultats de trimal pour identifier des alignements de mauvaise qualité.
Trimer les alignements avec trimal. Nous pouvons utiliser les résultats de trimal pour identifier des alignements de mauvaise qualité.
Line 206: Line 249:
-->
-->
-
==Super-arbres==
 
<!--
<!--
 +
==Super-arbres==
 +
Nous allons utiliser les arbres individuels protéiques ainsi que celui reconstruit à partir de la petite sous-unité de l’rRNA.
Nous allons utiliser les arbres individuels protéiques ainsi que celui reconstruit à partir de la petite sous-unité de l’rRNA.
Sur genologin le script ssuTree.sh était celui-là :
Sur genologin le script ssuTree.sh était celui-là :
Line 317: Line 361:
-->
-->
 +
<!--
==Comparaison des arbres==
==Comparaison des arbres==
Line 335: Line 380:
Commentez les résultats. Quel est l’arbre le plus différent des autres ? Qu’est-ce qui pourrait l’expliquer ?
Commentez les résultats. Quel est l’arbre le plus différent des autres ? Qu’est-ce qui pourrait l’expliquer ?
</pre>
</pre>
-
 
+
-->
----
----
-
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces:_pr.C3.A9paration_des_fichiers Arbre espèces]
+
*retour à [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Arbre_esp.C3.A8ces:_pr.C3.A9paration_des_fichiers Atelier Phylogénomique]

Current revision as of 14:49, 8 December 2022

Contents

Liens

Introduction

Support de cours : supports

Nous allons utiliser un sous ensemble de gènes concervés chez Prochlorococcus et Synechococcus pour expérimenter les différentes méthodes de reconstruction phylogénomiques. Nous nous initierons à la comparaison d’arbres.

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

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

Si le répertoire n'existe pas déjà, le créer.

mkdir -p ~/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

Exemple de la création d'un script et lancement du job avec sbatch

mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA

Créer le fichier suivant commande.sh.

#!/usr/bin/bash
~/work/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
sbatch commande.sh

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

grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
Question 5.1:
A quoi correspond ce quorum. Pourquoi utiliser un seuil de 16? 
Combien de fichiers avez-vous obtenus?
Est-ce pertinent dans les situations où vous avez un grand nombre de souches?
Et quand les espèces étudiées sont très éloignées phylogénétiquement parlant ?

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 (merci BioPerl).

mkdir ~/work/ProchlorococcusSynechococcus/OG/alignment

Attention à faire un script et le lancer en sbatch ou en srun --pty bash

 ~/work/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.

Il soumet sur le cluster donc, merci de vous remettre sur genologin pour lancer la commande ci-dessous.

Pour monitorer vos jobs : squeue -u $USER

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

Comptez les nombre d'alignements obtenus.

Question 5.2:
Quelle est la raison de passer par un alignement des peptides pour obtenir l'alignement en nucléotides?
Pour quelle raison allons nous travailler sur un sous ensemble d'alignements? 
Est-il pertinent de réaliser un échantillonnage des alignements par tirage aléatoire? Quelles autres approches pouvez-vous envisager?
Serait-il pertinent de réaliser plusieurs tirages? Quel usage pourriez-vous en faire?

MSK

Evaluation des alignements avec t_coffee

TCS is an alignment evaluation score that makes it possible to identify the most correct positions in an MSA. It has been shown that these positions are the most likely to be structuraly correct and also the most informative when estimating phylogenetic trees. The TCS evaluation and filtering procedure can be used to evaluate and filter any MSA.

sample_seq1.score_ascii displays the score of the MSA, the sequences and the residues.
sample_seq1.score_html displays a colored version score of the MSA, the sequences and the residues.

TCS allows you to filter out from your alignment regions that appears unreliable according to the consistency score; the filtering can be made at the residue level or the column level:

sample_seq1.tcs_residue_filter3 All residues with a TCS score lower than 3 are filtered out
sample_seq1.tcs_column_filter3 All columns with a TCS score lower than 3 are filtered out
original_alignments=~/work/ProchlorococcusSynechococcus/OG/alignment/
for i in $original_alignments/ali_dna_OG_*.fas 
do 
  ip=$(basename $i .aln)     
  outfile="~/work/ProchlorococcusSynechococcus/OG/alignment/"$ip".ali"
  echo "module load  bioinfo/T-COFFEE_11.00.8cbe486; t_coffee -infile $i -output score_ascii, aln, score_html -outfile $outfile;" 
done > OG_t_coffee_TCS.sh
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 OG_t_coffee_TCS.sh
squeue -l -u $USER

Sélection alignements avec un SCORE=1000.

goodDNA_alignments=~/work/ProchlorococcusSynechococcus/OG/goodDNA_alignments/
goodPEP_alignments=~/work/ProchlorococcusSynechococcus/OG/goodPEP_alignments/
original_alignments=~/work/ProchlorococcusSynechococcus/OG/alignment/
mkdir $goodDNA_alignments $goodPEP_alignments
 
goodDNA_alignment_list=~/work/ProchlorococcusSynechococcus/OG/goodDNA_alignments/DNASCORE1000.lst
goodPEP_alignment_list=~/work/ProchlorococcusSynechococcus/OG/goodPEP_alignments/PEPSCORE1000.lst
 
if [[ -f $goodDNA_alignment_list ]]; then
 rm $goodDNA_alignment_list
fi
 
if [[ -f $goodPEP_alignment_list ]]; then
 rm $goodPEP_alignment_list
fi
 
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/*.ali; 
do 
 if [[ $(grep -H 'SCORE=1000' $i) ]]; then
   ip=$(basename $i .ali)
   oriDNA="$original_alignments$ip"
   goodDNA="$goodDNA_alignments$ip"
   cp $oriDNA $goodDNA
   echo $oriDNA >> $goodDNA_alignment_list
 
   PEP=${ip:11:15}
   oriPEP="${original_alignments}ali_pep_OG_${PEP}"
   goodPEP="${goodPEP_alignments}ali_pep_OG_${PEP}"
   cp $oriPEP $goodPEP
   echo $oriPEP >> $goodPEP_alignment_list
 fi
done
cat $goodDNA_alignment_list

Edition des alignements avec trimal: MSK