Atelier Phylogénomique Arbre espèces
From silico.biotoul.fr
Contents |
Liens
- retour à [1]
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
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/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.
good_alignments=~/work/ProchlorococcusSynechococcus/OG/good_alignments/ mkdir $good_alignments for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/*.ali; do ip=$(basename $i) good="$good_alignments$ip" if [[ $(grep -H 'SCORE=1000' $i) ]]; then echo $good cp $i $good fi done
Edition des alignements avec trimal
MSK