Atelier Phylogénomique

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

for i in ~/work/Prochlorococcus/peptide/*.faa; 
do      
  echo "module load system/Python-2.7.2; module load bioinfo/eggnog-mapper-2.1.8; emapper.py -i $i --cpu 4 --output $i_NOG -m diamond --data_dir $EGGNOG_DATA_DIR"
done > eggnog-mapper.sh
 
cat eggnog-mapper.sh

sarray -J eggNOG -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=4 eggnog-mapper.sh
squeue -l -u $USER
 
=== OrthoFinder ===
'''OrthoFinder''' est une plateforme rapide, précise et complète pour la génomique comparative. Il trouve des '''orthogroupes''' et des orthologues, déduit des arbres de gènes enracinés pour tous les orthogroupes et identifie tous les événements de duplication de gènes dans ces arbres. 
 
Il déduit également un arbre des espèces enraciné pour l'espèce analysée et fait correspondre les événements de duplication de gènes des arbres génétiques aux branches de l'arbre des espèces. 
 
OrthoFinder fournit également des statistiques complètes pour les analyses génomiques comparatives. 
 
Suivre [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_OrthoFinder OrthoFinder]]!
 
===Préliminaires===
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.5:
Selon vous qu'est-ce qui guide le choix du type de séquences à utiliser dans les comparaisons (peptides ou nucléotidiques)?
</pre>
[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Blast Blast All-vs-All]
<!--
====Blast All-All====
Nous allons utiliser [http://bioinfo.genotoul.fr/index.php/how-to-use/?software=How_to_use_SLURM_NCBI_Blast%2B NCBI_Blast+].
 
Nous allons copier les fichiers peptides dans un répertoire de travail:
<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/Prochlorococcus/peptide
cp ~/work/Prochlorococcus/prokka/Aaa*/*.faa ~/work/Prochlorococcus/peptide/.
 
ls -l ~/work/Prochlorococcus/peptide
</pre>
=====make blast database=====
Exemple:
<pre style="color:green;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+
srun -n1 -l makeblastdb -in Aaaa.faa -dbtype prot
</pre>
Vous allez procéder comme précédemment, avec un script donné à ''sarray'', pour réaliser le makeblastdb sur tous les fichiers.
 
MSK
<syntaxhighlight lang="bash">
for i in ~/work/Prochlorococcus/peptide/*.faa; 
do      
  echo "module load bioinfo/ncbi-blast-2.7.1+; makeblastdb -in $i  -dbtype prot;" 
done > makeblastdb.sh
 
cat makeblastdb.sh
</syntaxhighlight>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 makeblastdb.sh
squeue -l -u <user>
</pre>
=====Intra genomes=====
Nous vous proposons un script perl pour réaliser les ''blastp'' de l'ensemble des génomes avec '''sbatch''':
<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/Prochlorococcus/BlastP
 
~/work/scripts/blastp_intra.pl
 
squeue -l -u <user>
 
ls BlastP | wc -l
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.6:
Explicitez les paramètres de blast utilisés dans le script blastp_intra.pl.
Combien de fichiers attendez-vous?
</pre>
 
=====Paire de genomes=====
<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/Prochlorococcus/BlastP
</pre>
Exemple : 
<pre style="color:green;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+
srun -n1 -l blastp -query ~/work/Prochlorococcus/peptide/Aaaa.faa -db ~/work/Prochlorococcus/peptide/Aaaa.faa -seg yes -dbsize 100000000 -evalue 1e-5 -outfmt 6 -num_threads 1 -out ~/work/Prochlorococcus/BlastP/Aaaa_Aaaa.tab
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.6:
Explicitez et justifiez les paramètres de blast utilisés dans votre script.
</pre>
 
=====Boucle sur les genomes=====
Utilisez ''sarray'' pour réaliser les ''blast'' en toutes les paires de génomes.
 
MSK
<syntaxhighlight lang="bash">
evalue=1e-5
dbsize=100000000
for i in ~/work/Prochlorococcus/peptide/*.faa; 
do 
  ip=$(basename $i .faa)     
  for j in ~/work/Prochlorococcus/peptide/*.faa;   
  do
    jp=$(basename $j .faa)     
    outfile="~/work/Prochlorococcus/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>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh
squeue -l -u <user>
</pre>
 
Pour les blastp inter génomes, nous allons utiliser un script similaire au précédent, mais avec une double boucle (sur -query et -db).
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cd ~/work/Prochlorococcus
~/work/scripts/blastp_inter.pl
 
squeue -l -u <user>
 
ls BlastP | wc -l
</pre>
vérifiez que les fichiers de sorties de blast sont présents et non vides.
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.7:
Combien de fichiers attendez-vous?
</pre>
-->
 
===PorthoMCL===
MSK
<!--
[https://github.com/etabari/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 [http://bdataanalytics.biomedcentral.com/articles/10.1186/s41044-016-0019-8 DOI: 10.1186/s41044-016-0019-8]
 
Il s'inspire d'[https://www.ncbi.nlm.nih.gov/pmc/articles/PMC403725/ orthoMCL]
 
Suivre [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_PorthoMCL PorthoMCL]
-->
<!--
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -o-pre-wrap"> 
Question 1.8:
Rappelez les différentes étapes réalisées par le logiciel. 
Précisez les paramètres utilisés pour identifier les paires de gènes orthologues.
</pre>
 
Le logiciel est disponible ici: /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master. Il est basé sur l'enchainement de programmes ''python''. Ils peuvent être lancés séparément ou automatiquement à l'aide d'un script shell. 
 
MSK
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> 
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
</pre>
 
Nous allons utiliser une version plus simple:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> 
mkdir ~/work/Prochlorococcus/PorthoMCL
cd ~/work/Prochlorococcus/PorthoMCL
srun --pty bash ~/work/scripts/short_PorthoMCL.sh
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> 
Question 1.9:
Editez le script et donnez une description des différentes tâches réalisées.
Que proposeriez-vous pour l'améliorer?
</pre>
====Step 4: Parse BLAST results====
Comme les différents blastp ont déjà été réalisés, nous passons directement à l'étape 4.
 
Mais nous devons:
=====Concaténer les fichiers de résultats blastp par souche:=====
<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/Prochlorococcus/PorthoMCL/3.blastres
</pre>
 
MSK
<syntaxhighlight lang="bash">
for i in ~/work/Prochlorococcus/peptide/*.faa; 
do 
  ip=$(basename $i .faa)     
  outfile=~/work/Prochlorococcus/PorthoMCL/3.blastres/$ip.tab;   
  if [ -f $outfile ]; then
    echo "rm $outfile"
    rm $outfile
  fi
  for j in ~/work/Prochlorococcus/peptide/*.faa; 
  do
    jp=$(basename $j .faa)     
    infile=~/work/Prochlorococcus/BlastP/$ip"_"$jp".tab"
    if [ -f $infile ]; then
      echo "$infile"
      cat $infile >> $outfile
    else
      echo "Error $infile not found"
    fi
  done
done 
</syntaxhighlight>
 
=====Modifier le nom des protéines dans les fichiers blast et fasta=====
Pour que les programmes suivants identifient le nom de souche et le nom de gène, il est nécessaire de rajouter un séparateur ('|') dans les noms des protéines: Aaaa.g_00001 devient Aaaa|.g_00001. Pour cela, nous pouvons utiliser la commande sed.
 
Dans les fichiers ~/work/Prochlorococcus/PorthoMCL/3.blastres/*.tab
 
MSK
<syntaxhighlight lang="bash">
for file in ~/work/Prochlorococcus/PorthoMCL/3.blastres/*.tab
do
  prefix=$(basename "$file")
  sed -i -r 's/([A-Z][a-z]{3})(\w*\.)/\1|\2/g' $file
done
</syntaxhighlight>
<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/Prochlorococcus/PorthoMCL/compliantFasta
</pre>
 
et dans les fichiers  ~/work/Prochlorococcus/peptide/*.faa
 
MSK
<syntaxhighlight lang="bash">
for file in ~/work/Prochlorococcus/peptide/*.faa
do
    prefix=$(basename $file .pep)
    sed -r 's/([A-Z][a-z]{3})/\1|/g' $file > ~/work/Prochlorococcus/PorthoMCL/compliantFasta/$prefix.fasta
done
</syntaxhighlight>
 
=====Changer le format des fichier blastp=====
<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/Prochlorococcus/PorthoMCL/4.splitSimSeq
</pre>
Exemple pour une souche:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclBlastParser /home/yquentin/work/Prochlorococcus/PorthoMCL/3.blastres/Aaaa.tab ~/work/Prochlorococcus/PorthoMCL/compliantFasta > ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq/Aaaa.ss.tsv
</pre>
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"> 
head ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq/Aaaa.ss.tsv
</pre>
'''Format de sortie''':
 query hit evalueMant evalueExp percentMatch
*evalueMant et evalueExp : mantisse (partie décimale) et exposant de la e-value
*percentMatch = longueur de l'alignement / longueur de la séquence la plus courte (query ou hit) *100
Appliquez à toutes les souches.
 
MSK
<syntaxhighlight lang="bash">
for file in ~/work/Prochlorococcus/peptide/*.faa
do
  prefix=$(basename $file .faa)
  /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclBlastParser ~/work/Prochlorococcus/PorthoMCL/3.blastres/$prefix.tab ~/work/Prochlorococcus/PorthoMCL/compliantFasta > ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq/$prefix.ss.tsv
done
</syntaxhighlight>
 
Vérifiez les fichiers obtenus.
 
====Step 5: Finding Best Hits====
Les paralogues sont identifiés et un score non normalisé leurs aient assignés. Les scores seront normalisés à l'étape 5.3 pour qu'ils soient comparables entre les génomes.
<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/Prochlorococcus/PorthoMCL/5.paralogTemp
mkdir -p ~/work/Prochlorococcus/PorthoMCL/5.besthit
</pre>
Créer un fichier avec la liste des taxon
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
ls -1 ~/work/Prochlorococcus/peptide/*faa | sed -r 's/.+([A-Z][a-z]{3}).faa/\1/' > ~/work/Prochlorococcus/PorthoMCL/taxon_list
</pre>
Test avec une souche:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsBestHit.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -s ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq -b ~/work/Prochlorococcus/PorthoMCL/5.besthit -q ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp -x 1
</pre>
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
head ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp/Aaaa.pt.tsv
head ~/work/Prochlorococcus/PorthoMCL/5.besthit/Aaaa.bh.tsv
</pre>
Boucle sur les souches:
 
MSK
<syntaxhighlight lang="bash">
for num in {1..12}
do
  /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsBestHit.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -s ~/work/Prochlorococcus/PorthoMCL/4.splitSimSeq -b ~/work/Prochlorococcus/PorthoMCL/5.besthit -q ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp -x $num
done
</syntaxhighlight>
 
====Step 6: Finding Orthologs====
A la sortie de cette étape nous obtenons l'ensemble des paires de gènes orthologues. Les fichiers de sortie sont concaténés. Le fichier concaténé est donnée en entrée à MCL.
<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/Prochlorococcus/PorthoMCL/6.orthologs
</pre>
Test avec une souche:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsOrthologs.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -b ~/work/Prochlorococcus/PorthoMCL/5.besthit -o ~/work/Prochlorococcus/PorthoMCL/6.orthologs -x 1
</pre>
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
head ~/work/Prochlorococcus/PorthoMCL/6.orthologs/Aaaa.ort.tsv
</pre>
Boucle sur les souches:
 
MSK
<syntaxhighlight lang="bash">
for num in {1..12}
do
  /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsOrthologs.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -b ~/work/Prochlorococcus/PorthoMCL/5.besthit -o ~/work/Prochlorococcus/PorthoMCL/6.orthologs -x $num
done
</syntaxhighlight>
 
====Step 7: Finding Paralogs====
Les scores calculés à l'étape 5 sont normalisés.  Pour la normalisation, PorthoMCL utilise la moyenne des scores des paralogues qui ont des orthologues.
 
Pour trouver les paralogues, nous devons extraire tous les gènes qui ont une relation orthologique dans chaque génome. Ceci peut être réalisé en utilisant un ensemble de commandes bash comme suit.
 
<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/Prochlorococcus/PorthoMCL/7.ogenes
cd ~/work/Prochlorococcus/PorthoMCL
</pre>
genes in the second column
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
awk -F'[|\t]' '{print $4 >> ("7.ogenes/"$3".og.tsv")}' 6.orthologs/*.ort.tsv
</pre>
genes in the first column
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
awk -F'[|\t]' '{print $2 >> ("7.ogenes/"$1".og.tsv")}' 6.orthologs/*.ort.tsv
</pre>
Avec ces fichiers (fichiers og) en main, nous pouvons normaliser les relations inparalogiques :
 
Le résultat de cette étape est l'ensemble des gènes paralogues avec des scores normalisés.
 
<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/Prochlorococcus/PorthoMCL/7.paralogs
</pre>
Test avec la première souche:
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsInParalogs.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -q ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp -o ~/work/Prochlorococcus/PorthoMCL/7.ogenes -p ~/work/Prochlorococcus/PorthoMCL/7.paralogs -x 1
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
head ~/work/Prochlorococcus/PorthoMCL/7.paralogs/Aaaa.par.tsv
</pre>
Boucle sur les souches:
 
MSK
<syntaxhighlight lang="bash">
for num in {1..12}
do
  /home/formation/public_html/M2_Phylogenomique/PorthoMCL-master/porthomclPairsInParalogs.py -t ~/work/Prochlorococcus/PorthoMCL/taxon_list -q ~/work/Prochlorococcus/PorthoMCL/5.paralogTemp -o ~/work/Prochlorococcus/PorthoMCL/7.ogenes -p ~/work/Prochlorococcus/PorthoMCL/7.paralogs -x $num
done
</syntaxhighlight>
 
====Step 8: MCL====
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.10:
Rappelez le principe de MCL et les paramètres utilisés.
Quel est l'effet de ces paramètres?
</pre>
C'est exactement comme la dernière étape d'OrthoMCL. Il suffit d'exécuter MCL sur les sorties du script :
 
Notez que -t définit le nombre de processeurs que vous devez utiliser pour MCL.
 
[http://bioinfo.genotoul.fr/index.php/how-to-use/?software=How_to_use_SLURM_MCL How to use SLURM MCL]
 
Location: /usr/local/bioinfo/src/MCL
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash 
#Load modules
module load bioinfo/mcl-14-137
 
cat ~/work/Prochlorococcus/PorthoMCL/6.orthologs/*.tsv >> ~/work/Prochlorococcus/PorthoMCL/8.all.ort.tsv
 
mcl ~/work/Prochlorococcus/PorthoMCL/8.all.ort.tsv --abc -I 1.5 -t 4 -o ~/work/Prochlorococcus/PorthoMCL/8.all.ort.group
</pre>
Taille des OG:
 
Vous pouvez utiliser [http://www.shellunix.com/awk.html awk] avec NF qui est le nombre de champs de l'enregistrement courant.
 
MSK
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
awk 'BEGIN {
 n=1;
 printf("class\tsize\n");
}
{ 
 printf(n"\t"NF"\n");
 n = n+1;
}' ~/work/Prochlorococcus/PorthoMCL/8.all.ort.group > ~/work/Prochlorococcus/PorthoMCL/8.all.ort.group_class_size.tab
</pre>
 
Taux de paralogie (avec ''awk'' mais en vous inspirant du [http://silico.biotoul.fr/p/M1_BBS_Graphes_TP_Recherche_de_communaut%C3%A9s_dans_les_graphes TP5] de M1):
 
MSK
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
awk '{
  j=0;
  for (i=1; i<=NF; i++) {
    str=substr($i,1,4);
    strlist[str]++;
    j++;
  }
  if ( j > 12 ) {
    nb=0;
    for (str in strlist ) {nb++}
    par=(NF-nb)/NF
    printf("%i\t%i\t%5.3f", NF, nb, par)
    for (str in strlist ) {
      printf("\t%s:%s", str, strlist[str]);
    }
    printf("\n")
  }
  delete strlist;
}' ~/work/Prochlorococcus/PorthoMCL/8.all.ort.group
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Que pensez-vous de la distribution des taille des OG? Quelle-est la taille maximum attendue?
</pre>
Paralogs
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">cat 7.paralogs/*.tsv >> 8.all.par.tsv
 
mcl 8.all.par.tsv  --abc -I 1.5 -t 4 -o 8.all.par.group
</pre>
 
====Visualisation des OG====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
module load system/R-3.5.1
</pre>
=====Vision globale avec genoplotR=====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Rscript ~/work/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_all.pdf 
 
evince work/Prochlorococcus/images/8.all.ort.group_all.pdf 
</pre>
[http://silico.biotoul.fr/enseignement/m2BBS/Phylogenomic/figures/8.all.ort.group_all.pdf 8.all.ort.group_all.pdf]
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.11:
Analysez la figure obtenue.
Comparez-là avec celle réalisée avec les blastn.
</pre>
 
=====Restreindre aux OG sans paralogues=====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Rscript ~/work/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_nopara.pdf --max_paralogs=0
 
evince work/Prochlorococcus/images/8.all.ort.group.1.5-5.50_nopara.pdf 
</pre>
[http://silico.biotoul.fr/enseignement/m2BBS/Phylogenomic/figures/8.all.ort.group_nopara.pdf 8.all.ort.group_nopara.pdf]
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.12:
Comparez les deux figures, avec et sans les paralogues.
Comment se répartissent les régions peu denses en gènes orthologues?
</pre>
 
=====Restreindre à un seul OG=====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Rscript ~/work/scripts/genoplot_OG_links.R --MCL_file=~/work/Prochlorococcus/PorthoMCL/8.all.ort.group --pdf_file=~/work/Prochlorococcus/images/8.all.ort.group_OG5.pdf --select_OG=5
 
evince work/Prochlorococcus/images/8.all.ort.group_OG5.pdf 
</pre>
[http://silico.biotoul.fr/enseignement/m2BBS/Phylogenomic/figures/8.all.ort.group_all.pdf 8.all.ort.group_OG5.pdf]
 
Liste des membres du groupe 5
<pre>
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
</pre>
Extraction du sous-graphe correspondant à ces gènes:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
for i in 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 
do
 grep $i ~/work/Prochlorococcus/PorthoMCL/8.all.ort.tsv
done | sort | uniq > ~/work/Prochlorococcus/PorthoMCL/G5.gr
</pre>
Tracez le graphe de ces gènes avec la librairie ''igraph'' (c.f. [http://silico.biotoul.fr/p/M1_BBS_Graphes_TP_Recherche_de_communaut%C3%A9s_dans_les_graphes TP5] de M1).
 
MSK
<syntaxhighlight lang="python">
library(igraph)
gr <- read.graph("work/Prochlorococcus/PorthoMCL/G5.gr",format="ncol")
plot(gr, vertex.label.cex=0.8, vertex.size=5, vertex.shape="sphere", edge.width=E(gr)$weight, edge.color="red")
</syntaxhighlight>
[http://silico.biotoul.fr/enseignement/m2BBS/Phylogenomic/figures/graph_OG5.png graph_OG5.png]
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.13:
Commentez ces résultats.
Pouvez vous formuler des hypothèses sur l'évolution de ces gènes?
Comment vérifier ces hypothèses?
</pre>
=====Tester l'effet de l'IF sur la 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. 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Rscript ~/work/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
</pre>
Code couleurs:
<pre>
bleu :   tous les OG
rouge :  OG sans paralogues
vert :   OG constitués uniquement d'orthologues 1:1
orange : OG avec au moins un paralogue(s)
</pre>
=====Tester l'effet de l'IF sur la paralogie=====
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Rscript ~/work/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 ~/work/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
</pre>
====old stuff2====
 
read partition as list in R
[http://bioinfo.genotoul.fr/index.php/how-to-use/?software=How_to_use_SLURM_R ow_to_use_SLURM_R]
<div style="height: 150px; overflow:auto; border: thin solid blue; background: grey">
<pre>
partition <- '/home/yquentin/work/Prochlorococcus/PorthoMCL/8.all.ort.group.5.50'
i <- 1
con = file(partition, "r")
while ( TRUE ) {
  line = readLines(con, n = 1)
  if ( length(line) == 0 ) {
    break
  }    
  list[[i]] <- (strsplit(line, "\t", fixed = FALSE, perl = TRUE))
  #print(list[[i]])
  i <- i + 1
}
close(con)
 
</pre>
</div>
 
<div style="height: 150px; overflow:auto; border: thin solid blue; background: grey">
<pre>
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()
</pre>
</div>
-->
 
=== Panaroo ===
[https://gtonkinhill.github.io/panaroo panaroo]
 
===PanOCT===
Pan-genome Ortholog Clustering Tool, est un programme écrit en PERL pour l'analyse pan-génomique d'espèces ou de souches procaryotes étroitement apparentées. Contrairement aux programmes traditionnels de détection d'orthologues basés sur des graphes, il utilise la micro-synténie ou le voisinage de gènes conservés (CGN) en plus de l'homologie pour placer avec précision les protéines dans des groupes orthologues.
 
[https://sourceforge.net/projects/panoct/ panoct project]
 
Suivre [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_PanOCT PanOCT]
<!--
Vous trouverez le package dans le repertoire suivant:
 /home/formation/public_html/M2_Phylogenomique/PanOCT/panoct_v3.23/
Test:
 <pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">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
</pre>
====Préparation des fichiers d'entrée de PanOCT====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
mkdir -p ~/work/Prochlorococcus/panoct/results
</pre>
=====gene.att=====
Créer un fichier avec les coordonnées, noms, fonction et souches des gènes. 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cd ~/work/Prochlorococcus
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
 
ls panoct/results/*.tab
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cat panoct/results/*.tab > panoct/results/combined.att
head panoct/results/combined.att
</pre>
 
=====tags.txt=====
Liste des souches à analyser.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
for i in peptide/*.faa
do      
echo $(basename $i .faa)
done > panoct/results/genomes.list
more panoct/results/genomes.list
</pre>
 
=====peptides=====
Concaténer les fichier peptides dans un seul fichier.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cat peptide/*.faa > panoct/results/combined.fasta
 
head panoct/results/combined.fasta
</pre>
 
=====blast.txt=====
Concaténer les résultats des blastp dans un seul fichier.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cat BlastP/*.tab > panoct/results/combined.blast
 
head panoct/results/combined.blast
</pre>
 
====run panOCT====
panoct.pl: avec les paramètres par défaut.
<pre>
 -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]
</pre>
La commande:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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>
Lancer avec sbatch et un script du type "[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/scripts/panoct_P.csh panoct_P.csh]" (les chemins sont à changer).
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sbatch panoct_P.csh
squeue -l -u <user>
</pre>
-->
 
==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. 
 
[https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4765519/ Inside the Pan-genome - Methods and Software Overview]
 
Suivre: [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Analyses_pan-g%C3%A9nomiques Analyses pan-génomiques]
<!--
===Définitions===
D'après [https://www.ncbi.nlm.nih.gov/pubmed/25483351 Vernikos ''et al.'' 2015]
 
Le pan-génome a été défini pour la première fois par [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1216834/ 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 ([https://www.sciencedirect.com/science/article/pii/S0959437X05001759 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’orthologues dans les autres souches du clade.
 
===Mise en oeuvre===
Ce script permet de prendre en compte les paralogs ([https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4507327/ Chan et al., 2015]). 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
Celui-ci permet de reformater la table ''matchtable''.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/matchtable2panoct.result.pl --results_dir results
</pre>
Ce script calcule un ensemble de statistiques sur les résultats de panOCT.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
Un résumé des résultats se trouve dans la fichier: ~/work/Prochlorococcus/panoct/overview_stats.txt
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.14:
Commentez les résultats obtenus. Comparez ces résultats à ceux de Kettler et al., 2007.
</pre>
 
Tracé de l'histogramme
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">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/.
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.15:
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)?
</pre>
 
===Taille du pan génome===
La taille du pan génome à tendance à augmenter avec le nombre de génomes comparés (pour une revue: [https://www.ncbi.nlm.nih.gov/pubmed/19086349/ Tettelin et al., 2008]). 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 ([https://www.ncbi.nlm.nih.gov/pubmed/19086349/ 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 ([https://www.ncbi.nlm.nih.gov/pubmed/19086349/ 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:
<pre>
 -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>
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
Ajout de la taille des génomes en première colonne:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">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
</pre>
====Pan génome et génome coeur====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Rscript /home/formation/public_html/M2_Phylogenomique/scripts/pan_core_plot.R
 
evince ~/work/Prochlorococcus/images/pangenome.pdf 
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.16:
Décrire la figure obtenue. 
Comment sont obtenus les différents points gris? 
A quoi correspondent les triangles et les carrés de couleurs et les deux courbes associées?
Comparez cette figure à la Figure 1 de Kettler et al., 2007. 
</pre>
 
=====Log x Log plot des nouveaux gènes =====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Rscript /home/formation/public_html/M2_Phylogenomique/scripts/new_genes_plot.R
evince ~/work/Prochlorococcus/images/newgenome.pdf
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.17:
Décrire cette figure. 
Que représente le 'a' (alpha)?
Comparez la valeur obtenue avec celle reportée dans Tettelin et al., 2008.
Que pourrait être l'effet de l'ajout de nouveaux génomes sur la taille du pan-génome (voir Biller et al., 2015)?
Que pouvez-vous en conclure?
</pre>
====run_panoct.pl====
<pre>
--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.
</pre>
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
 
 results/R.plots/core_cluster_histogram.pdf
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1.18:
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)?
</pre>
 
===compute_pangenome.R===
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/PanGenomePipeline/PanGenomePipeline-master/pangenome/bin/compute_pangenome.R -i ~/work/Prochlorococcus/panoct/results/matchtable_paralog.txt -o ~/work/Prochlorococcus/panoct/results/pangenome_size -p 95 -q 0 -s 250
</pre>
 
Afficher les courbes:
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
R
pangenome_size <- '~/work/Prochlorococcus/panoct/results/pangenome_size'
pdf_file <- '~/work/Prochlorococcus/images/pangenome_size.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")
</pre>
 
===PanOCT classification des proteins de l'OG 5 de PorthoMCL===
Nous avons vu que le OG 5 de PorthoMCL renfermait 36 séquences. Comment se distribuent ces protéines avec panOCT?
 
MSK
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
if [ -f ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt ]; 
then
  echo "rm ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt"
fi
 
for i in 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 
do
 awk -v var=$i '{
  pattern=var;
  if (match($0,pattern)) {
   printf("%s\t%d\n", var,$1);
  }
 }' ~/work/Prochlorococcus/panoct/all_clusters_members.txt >> ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt
done
 
sort -nk 2 ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">Question 1.18:
Commentez le résultat obtenu. 
Que pouvez-vous en conclure?
</pre>
 
Nous pouvons annoter le graphe de la classe 5 trouvée avec PorthoMCL avec les classes trouvées avec panoct. Comme les noms des séquences ne suivent pas la même règle, nous devons modifier le fichier classes pour ajouter un '|' entre le nom de souche et le nom de gène. 
 
MSK
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sed -i -r 's/([A-Z][a-z]{3})(\w*\.)/\1|\2/g'  ~/work/Prochlorococcus/panoct/OG5_clusters_members.txt
</pre>
 
MSK
<syntaxhighlight lang="python">
library(igraph)
gr <- read.graph("work/Prochlorococcus/PorthoMCL/G5.gr", format="ncol")
cl <- read.table("work/Prochlorococcus/panoct/OG5_clusters_members.txt", row.names=1)
colors <- rainbow(5)[as.factor(cl[V(gr)$name,])]
plot(gr, vertex.label.cex=0.8, vertex.size=10, vertex.shape="sphere", vertex.color=colors, edge.width=E(gr)$weight, edge.color="red")
</syntaxhighlight>
-->
 
==Alignement et comparaison de génomes complets==
Suivre : [http://silico.biotoul.fr/site/index.php?title=Atelier_Phylog%C3%A9nomique_Alignement_Genomes Alignement Genomes] 
<!--
:'''Jeu de données'''
 
Vous pouvez retrouver les informations sur les 12 génomes de Prochlorococcus [http://www.m2p-bioinfo.ups-tlse.fr/p/Atelier_Phylog%C3%A9nomique#Caract.C3.A9ristiques_des_souches_.C3.A9tudi.C3.A9es ici] et les données dans le répertoire: /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/ 
 
Copiez les 12 génomes de Prochlorococcus dans un répertoire de votre ~/work, par exemple genomes_prochlo:
<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/Alignement_genomes/genomes_prochlo
cp  /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/*.fas ~/work/Alignement_genomes/genomes_prochlo/
ls ~/work/Alignement_genomes/genomes_prochlo/
</pre>
 
===Exploration de la diversité génomique à partir de l’ANI et des distances Mash===
 
=====Diversité génomique basée sur l’ANI=====
:'''Calculer l’ANI entre toutes les paires de génomes en utilisant la version basée sur [https://mummer4.github.io/ Mummer].''' http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/2018_supports/CoursAligntGenomes2018.pdf
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
module load system/Python-3.6.3
module load bioinfo/mummer-4.0.0beta2
average_nucleotide_identity.py -h
 
average_nucleotide_identity.py -v -i ~/work/Alignement_genomes/genomes_prochlo/ -o  ~/work/Alignement_genomes/genomes_ANIm_output/  --gformat png,pdf,eps,svg --write_excel -m ANIm
</pre>
 
Exemple de script "[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/scripts/RunSLURM_ANI.csh RunSLURM_ANI.csh]" (les chemins sont à changer):
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sbatch RunSLURM_ANI.csh
squeue -l -u <user>
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.1:
Regardez les différents fichiers résultats.
Regardez la couverture et le pourcentage d’identité des alignements et commentez les valeurs obtenues.
Qu’en concluez-vous sur la diversité des génomes de Prochlorococcus ?
</pre>
 
:'''Construire un arbre de Neighbor-Joining basé sur le ANI (ANIm_percentage_identity et ANIm_alignment_coverage) avec le logiciel de votre choix'''
Vous pourrez par exemple utiliser la fonction ''nj'' du package R ''ape''. Notez que la commande ''nj'' prend en entrée une matrice de distance. La fonction ''heatmap'' ([https://www.r-graph-gallery.com/index.html r-graph-gallery]) peut être utile pour visualiser les relations entre les souches. 
 
<syntaxhighlight lang="python">
id_file <- "work/Alignement_genomes/genomes_ANIm_output/ANIm_percentage_identity.tab"
id_data <- read.table(file=id_file, header=T, row.names=1)
heatmap(as.matrix(id_data), scale="none", symm=T, main="ANIm_percentage_identity")
 
co_file <- "work/Alignement_genomes/genomes_ANIm_output/ANIm_alignment_coverage.tab"
co_data <- read.table(file=co_file, header=T, row.names=1)
heatmap(as.matrix(co_data), scale="none", symm=T, main="ANIm_alignment_coverage")
 
id_nj <- nj(as.matrix(1-id_data))
plot(id_nj, main="ANIm_percentage_identity")
 
co_nj <- nj(as.matrix(1-co_data))
plot(co_nj, main="ANIm_alignment_coverage")
</syntaxhighlight>
 
On pourra par exemple utiliser le site [http://www.phylogeny.fr/one_task.cgi?task_type=bionj phylogeny.fr] dans lequel on importera le fichier genomes_ANIm_output/ANIm_percentage_identity.tab (modifié en enlevant la 1ère ligne et en la remplaçant par le nombre de séquences 12).
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.2:
Interprétez les résultats.
</pre>
 
:'''Sélectionnez les génomes en ne gardant que le sous-groupe de 6 génomes qui ont au moins 28% de couverture avec tous les autres génomes (pour cela regardez le fichier ANIm_alignement_coverage.tab)'''
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.3:
Citez les.
</pre>
 
=====Distance Mash entre les génomes=====
Passez à l'étape suivante.
:'''Calculer la distance Mash entre toutes les paires de génomes'''
 
Documentation : [https://mash.readthedocs.io/en/latest/tutorials.html Mash]
 
En mode intéractif:
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
module load system/R-3.4.3 compiler/gcc-7.2.0
module load bioinfo/Mash-2.1
~/work/scripts/Mash_sketch.sh  ~/work/Alignement_genomes/genomes_prochlo/
</pre>
 
Les résultats se trouvent dans le répertoire data_MashSketches/.
 
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
~/work/scripts/Mash_dist_allpairs.sh data_MashSketches/
</pre>
 
Le résultat se trouve dans le fichier mash_dist.out
 
NB:
* 1ère colonne: Jaccard index (~fraction de kmers partagés)
* 2ème colonne : p-valeur
* 3ème colonne : distance Mash (estimation du taux de mutation selon un modèle d'évolution simple)
 
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.4:
Interprétez les résultats. 
Comparez les résultats de distance MASH à ceux de l’ANI.
</pre>
===Alignements Mauve et ProgressiveMauve===
<pre style="color:brown;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
NB : Commencez par lancer l’alignement ProgressiveMauve (environ 50-60 minutes de temps d’execution) avant de faire la question sur l'alignement Mauve !!!
</pre>
 
=====Alignements Mauve d'un sous-ensemble de 6 génomes=====
<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/Alignement_genomes/cat_genomes_prochlo
</pre>
======Concaténer les 6 génomes sélectionnés à la question précédente dans un fichier multifasta======
<syntaxhighlight lang="bash">
rm ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.fna
rm ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.gbk
rm ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.gff
 
for i in Aaab Aaag Aaaj Aaaf Aaak Aaae 
do
 echo $i
 cat ~/work/Alignement_genomes/genomes_prochlo/$i.fas >> ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.fna
 cat ~/work/Prochlorococcus/prokka/$i/$i.gbk >> ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.gbk
 cat ~/work/Prochlorococcus/prokka/$i/$i.gff >> ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.gff
done
</syntaxhighlight>
MSK
<syntaxhighlight lang="bash">
for i in Aaax  Aaay Aaaz
do
 echo $i
 cat ~/work/Alignement_genomes/genomes_prochlo/$i.fas >> ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.fna
 cat ~/work/Prochlorococcus/prokka/$i/$i.gbk >> ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.gbk
done
 
grep -c '>' ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.fna
grep -c 'LOCUS' ~/work/Alignement_genomes/cat_genomes_prochlo/6GC_Prochlorococcus.gbk
</syntaxhighlight>
 
======Lancement de l’alignement des 6 génomes sur le cluster SLURM======
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mkdir ~/work/Alignement_genomes/Mauve
cd ~/work/Alignement_genomes
</pre>
Exemple de script "[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/scripts/RunSLURM_Mauve_6GProchlo.csh RunSLURM_Mauve_6GProchlo.csh]" (les chemins sont à changer):
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sbatch RunSLURM_Mauve_6GProchlo.csh
squeue -l -u <user>
</pre>
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mkdir ~/work/Alignement_genomes/Mauve
cd ~/work/Alignement_genomes
srun --pty bash
Exemple de la commande à lancer avec le fichier en format gbk
<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/mauve_2.4.0
mauveAligner --output=Mauve/6GC_Prochlorococcus_gbk.mauve_def --permutation-matrix-output=Mauve/6GC_Prochlorococcus_gbk.permutation_matrix --output-guide-tree=Mauve/6GC_Prochlorococcus_gbk.tree --output-alignment=Mauve/6GC_Prochlorococcus_gbk_mauve.xmfa cat_genomes_prochlo/6GC_Prochlorococcus.gbk
</pre>
Soumission du job (vérifiez que vous êtes bien sur '''genologin2'''!):
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sbatch --chdir=~/work/Alignement_genomes RunSLURM_Mauve_6GProchlo.csh
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
#!/bin/bash
#SBATCH --time=02:00:00 #job time limit
#SBATCH -J testjob
#SBATCH -o RunSLURM_Mauve_6GProclo.out
#SBATCH -e RunSLURM_Mauve_6GProclo.err
#SBATCH --mem=8G
#SBATCH --cpus-per-task=4 #ncpu on the same node
#SBATCH --mail-type=BEGIN,END,FAIL (email address is LDAP accounts) 
#Purge any previous modules
module purge
#Load the application
module load bioinfo/mauve_2.4.0
# My command lines I want to run on the cluster
mauveAligner --output=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.mauve_def --permutation-matrix-output=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.permutation_matrix --output-guide-tree=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.tree --output-alignment=/home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus_mauve.xmfa  /home/hchiapello/work/TP_M2_BIOINFO/TP_2018/6GC_Prochlorococcus.fna
</pre>
 
======Analyser et interpréter les résultats en les visualisant via l’interface Mauve (commande Mauve)======
Remarques: 
*dans le fichier ''Mauve/6GC_Prochlorococcus_gbk_mauve.xmfa'', le chemin du fichier gbk et relatif, penser à lancer Mauve dans le bon répertoire pour avoir les annotations des gènes.
*lien entre le code et le nom de souche: [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Prochlorococcus/NCBI/species_strain_names.txt species_strain_names.txt]
<pre>
#FormatVersion Mauve1
#Sequence1File	cat_genomes_prochlo/6GC_Prochlorococcus.gbk
#Sequence1Entry	1
#Sequence1Format	GenBank
#Annotation1File	cat_genomes_prochlo/6GC_Prochlorococcus.gbk
...
</pre>
 
======Exploration du contexte génomique======
L'outil ''Sequence Navigator'' (les jumelles) permet de rechercher un ou plusieurs gènes sur différents critères. Nous allons utiliser cette fonctionnalité pour analyser le contexte génomique des gènes suivants. Les noms des gènes sont accessibles  par ''locus tag''. En vous plaçant sur le gène, vous avez ses annotations avec ''View Genbank...''. En quittant les jumelles, vous pouvez analyser la conservation du contexte à différentes échelles. 
<pre>
Aaab.g_00239 Aaab.g_00820 Aaab.g_00825 
Aaag.g_00262 Aaag.g_00862 Aaag.g_00867 
Aaaj.g_00252 Aaaj.g_00786 Aaaj.g_00791 
Aaaf.g_00249 Aaaf.g_00782 Aaaf.g_00787 
Aaak.g_00251 Aaak.g_00834 Aaak.g_00839 
Aaae.g_00242 Aaae.g_00786 Aaae.g_00791 
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.5:
Combien y’a-t-il de LCB dans l’alignement ? Quel est leur poids minimal ? Y’a–t-il des réarrangements globaux dans l’alignement et si oui lesquels ? Décrire la structure de l’alignement. Que se passe-t-il si on fait varier le poids des LCB ?
Qu'avez-vous appris de l'analyse du contexte génomique des gènes.
</pre>
 
=====Alignement Progressive Mauve de l’ensemble complet des 6 génomes=====
Afin de comparer les logiciels Mauve et ProgressiveMauve nous allons analyser l'ensemble de 6 génomes avec  ProgressiveMauve.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
ls
mkdir ~/work/Alignement_genomes/ProgressiveMauve
cd ~/work/Alignement_genomes
</pre>
Créer un ficher .csh en prenant pour exemple le fichier "[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/scripts/RunSLURM_PMauve_6GProchlo.csh RunSLURM_PMauve_6GProchlo.csh]"
avec comme ligne de commande:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
progressiveMauve --output=ProgressiveMauve/6GC_Prochlorococcus_PMauve.xmfa cat_genomes_prochlo/6GC_Prochlorococcus.gbk
</pre>
[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Mauve Mauve]
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap"
>mkdir ~/work/Alignement_genomes/ProgressiveMauve
srun --pty bash
module load bioinfo/mauve_2.4.0
progressiveMauve --output=ProgressiveMauve/6GC_Prochlorococcus_PMauve.xmfa cat_genomes_prochlo/6GC_Prochlorococcus.gbk
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.6:
Comparez et interprétez les résultats obtenus.
</pre>
 
=====Alignement Progressive Mauve de l’ensemble complet des 12 génomes=====
Concaténer les 12 génomes dans un fichier multifasta
 
MSK
<syntaxhighlight lang="bash">
rm ~/work/Alignement_genomes/cat_genomes_prochlo/12GC_Prochlorococcus.fna
rm ~/work/Alignement_genomes/cat_genomes_prochlo/12GC_Prochlorococcus.gbk
for i in Aaaa Aaab Aaac Aaad Aaae Aaaf Aaag Aaah Aaai Aaaj Aaak Aaal
do
 echo $i
 cat ~/work/Alignement_genomes/genomes_prochlo/$i.fas >> ~/work/Alignement_genomes/cat_genomes_prochlo/12GC_Prochlorococcus.fna
 cat ~/work/Prochlorococcus/prokka/$i/$i.gbk >> ~/work/Alignement_genomes/cat_genomes_prochlo/12GC_Prochlorococcus.gbk
done
 
grep -c '>' ~/work/Alignement_genomes/cat_genomes_prochlo/12GC_Prochlorococcus.fna
grep -c 'LOCUS' ~/work/Alignement_genomes/cat_genomes_prochlo/12GC_Prochlorococcus.gbk
</syntaxhighlight>
======Lancer l’alignement ProgressiveMauve des 12 génomes sur le cluster SLURM======
Exemple de script "[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/scripts/RunSLURM_PMauve_12GProchlo.csh RunSLURM_PMauve_12GProchlo.csh]" (les chemins sont à changer).
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sbatch RunSLURM_PMauve_12GProchlo.csh
squeue -l -u <user>
</pre>
[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/Mauve Mauve]
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mkdir ~/work/Alignement_genomes/ProgressiveMauve
cd ~/work/Alignement_genomes
srun --pty bash
module load bioinfo/mauve_2.4.0
progressiveMauve --output=ProgressiveMauve/12GC_Prochlorococcus_PMauve.xmfa genomes_prochlo/12GC_Prochlorococcus.fna
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sbatch RunSLURM_PMauve_12GProchlo.csh
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.7:
Analyser et interpréter les résultats en les visualisant via l’interface Mauve
 
Si vous avez l'annotation des gènes dans le résultat de ''Progressive Mauve'', vous pouvez visualiser la conservation du contexte de ces gènes de l'OG 5 de PorthoMCL et proposer une interprétation des liens complexes existant entre ces gènes homologues.
</pre>
L'article de Yan et al., 2018 Genome rearrangement shapes ''Prochlorococcus'' ecological adaptation. Appl Environ Microbiol 84:e01178-18. https://doi.org/10.1128/AEM.01178-18 peut vous aider dans l'interprétation des résultats.
 
===Reconstruction de l’histoire évolutive des réarrangements et des états ancestraux===
masqué
:'''Exporter le fichier de permutation des LCBs de l’alignement généré à la question précédente avec ProgressiveMauve en cliquant dans l’interface Mauve sur Tools => Export => Export Permutation. Ainsi vous pouvez créer les fichiers permutations.'''
*6GC_Prochlorococcus_pmauve.permutation (pour se familiariser sur les résultats!)
*12GC_Prochlorococcus_pmauve.permutation
 
NB1: Il n’est pas possible de lancer ProgressiveMauve directement avec l’option –permutation-matrix (ne fonctionne qu’avec Mauve).
 
NB2: Il est nécessaire d'éditer le fichier 12GC_Prochlorococcus_pmauve.permutation pour le mettre au format MLGO : ajouter une ligne « >id genome » dans le même ordre que le fichier d’entrée 12GC_Prochlorococcus.fna et remplacer toutes les virgules par des espaces. Ajouter un espace avant le $ de chaque fin de ligne
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
 grep ">" 12GC_Prochlorococcus.fna
</pre>
 
Le fichier doit avoir le format suivant 
<pre>
>Aaaa
1 56 55 54 26 5 52 53 27 $
>Aaab
1 56 55 54 26 5 52 53 27 $
Etc...
</pre>
 
:'''Utiliser le logiciel MLGO pour inférer l’histoire évolutive des réarrangements à partir de la matrice de permutation modifiée (n'oubliez pas de calculer les bootstraps ).'''
 
Lien vers [http://www.geneorder.org/server.php MLGO]
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2.6:
Visualiser les arbres produits par MLGO en utilisant l’outil de votre choix (par exemple FigTree ou iTOL) et interpréter les résultats en fonction des écotypes.
</pre>
 
NB1: Le fichier gene_order.tree.PMAG contient l’arbre avec les nœuds ancestraux (A1,A2,…etc) et le fichier gene_order.out contient l’ordre des gènes (ici LCBs) dans ces génomes ancestraux.
 
NB2: Le fichier de sortie MLGO gene_order_bs.tree qui contient les valeurs de bootstrap entre [] est non standard et n’est pas lisible par FigTree. Pour pouvoir le lire effectuer la transformation suivante:
 
<pre>
more gene_order_bs.tree
 
(Aaab:0.00038261697365918421,(((Aaal:0.00110779153484151733,Aaac:0.00045702512028280982):0.01128794608013189253[100],((Aaah:0.00195679558538791641,Aaaa:0.00249248732070895488):0.00170192911132310224[98],(Aaad:0.00000001574704787426,Aaai:0.00000001574704787426):0.00384727219137918474[100]):0.00019414433073579879[65]):0.00328248641555709160[100],(Aaaj:0.00000001574704787426,(Aaaf:0.00000001574704787426,(Aaae:0.00000001574704787426,Aaak:0.00000001574704787426):0.00000001574704787426[34]):0.00000001574704787426[33]):0.00056868823569016999[98]):0.00018460659686999882[50],Aaag:0.00040691005233194506);
 
Est à transformer en:
 
(Aaab:0.00038261697365918421,(((Aaal:0.00110779153484151733,Aaac:0.00045702512028280982)100:0.01128794608013189253,((Aaah:0.00195679558538791641,Aaaa:0.00249248732070895488)98:0.00170192911132310224,(Aaad:0.00000001574704787426,Aaai:0.00000001574704787426)100:0.00384727219137918474)65:0.00019414433073579879)100:0.00328248641555709160,(Aaaj:0.00000001574704787426,(Aaaf:0.00000001574704787426,(Aaae:0.00000001574704787426,Aaak:0.00000001574704787426)34:0.00000001574704787426)33:0.00000001574704787426)98:0.00056868823569016999)50:0.00018460659686999882,Aaag:0.00040691005233194506);
</pre>
Ou plus simplement, vous pouvez utiliser [https://itol.embl.de/ iTOL] pous visualiser l'arbre directement avec les valeurs de bootstraps et éventuellement l'annoter.
-->
 
==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. 
 
Suivre [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Analyses_phylog%C3%A9nomiques Analyses phylogénomiques]
<!--
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]
 
<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/prokka
cd ~/work/Synechococcus
~/work/scripts/prokka_loop.pl --sample Synechococcus
 
squeue -l -u <user>
</pre>
Une fois les jobs terminés, vérifiez que les fichiers de sorties de prokka existent et ne sont pas vides.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
ls -l ~/work/Synechococcus/prokka/Aaa*/*.faa
</pre>
 
====''Synechococcus'' blastp All-All====
Nous allons reproduire le même enchaînement de scipts utilisés avec ''Prochlorococcus'' ([[Atelier_Phylogénomique#Blast_All-All]]) en prenant soin de remplacer Prochlorococcus par Synechococcus dans les noms des répertoires.
 
Copiez les fichier *.faa de prokka dans ~/work/Synechococcus/peptide/
 
<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/peptide
cd ~/work/Synechococcus/peptide
cp ~/work/Synechococcus/prokka/Aaa*/*.faa ~/work/Synechococcus/peptide/.
 
ls -l ~/work/Synechococcus/peptide
</pre>
 
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>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 makeblastdb.sh
squeue -l -u <user>
</pre>
 
Boucle sur les genomes de ''Synechococcus''.
 
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/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>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh
squeue -l -u <user>
</pre>
<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>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
squeue -l -u <user>
ls BlastP | wc -l
</pre>
 
===''Prochlorococcus'' ''versus'' ''Synechococcus''===
====Liens symboliques sur les fichiers peptides====
<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/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/.
</pre>
 
====Liens symboliques sur les fichiers blastp====
<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/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/.
</pre>
====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>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 blast_allall.sh
squeue -l -u <user>
</pre>
Vérifiez que vous avez bien le nombre de fichiers attendus!
 
===Groupes de gènes orthologues avec PanOCT===
====Préparation des fichiers combined====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
mkdir -p ~/work/Synechococcus/panoct/results
mkdir -p ~/work/ProchlorococcusSynechococcus/panoct/results
</pre>
'''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>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cat ~/work/Prochlorococcus/panoct/results/*.tab ~/work/Synechococcus/panoct/results/*.tab> ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att
head ~/work/ProchlorococcusSynechococcus/panoct/results/combined.att
</pre>
'''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.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cat peptide/*.faa > panoct/results/combined.fasta
 
grep -c '>' panoct/results/combined.fasta
</pre>
'''combined.blast'''
Concaténer les résultats des blastp dans un seul fichier.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cat BlastP/*.tab > panoct/results/combined.blast
 
head panoct/results/combined.blast
</pre>
 
====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).
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cd ~/work/ProchlorococcusSynechococcus/panoct
mkdir ~/work/ProchlorococcusSynechococcus/images
</pre>
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sbatch panoct_PS.csh
squeue -l -u <user>
</pre>
 
/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>
 
<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 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==11 && 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?
 
MSK
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
awk '{ system("grep "$3"  ~/work/Prochlorococcus/prokka/Aaaa/Aaaa.gff >>  ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific_annotations.txt")}' < ~/work/ProchlorococcusSynechococcus/Prochlorococcus_specific.txt
</pre>
 
=====Gènes trouvés chez toutes les souches de ''Synechococcus'' mais dans aucune de ''Prochlorococcus'':=====
 
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==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 lang="awk">
</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; 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; 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">
library('gplots')
 
Ppairs <- read.table('work/ProchlorococcusSynechococcus/Ppairs.txt')
PSpairs <- read.table('work/ProchlorococcusSynechococcus/PSpairs.txt')
input <- list(P=Ppairs, PS=PSpairs)
venn(input)
 
a<-venn(input, show.plot=FALSE)
intersections <- attr(a,"intersections")
intersections$P
intersections$PS
</syntaxhighlight>
 
====Comparaison des résultats de PanOCT avec PorthoMCL====
MSK
<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">
library('gplots')
 
PorthoMCL <- read.table('work/Prochlorococcus/PorthoMCL/8.all.ort.pairs.txt')
panoct <- read.table('work/ProchlorococcusSynechococcus/Ppairs.txt')
input <- list(PorthoMCL=PorthoMCL, panoct=panoct)
venn(input)
</syntaxhighlight>
-->
 
==Phylogénie basée sur les séquences des ARNr==
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 4.1:
Quel-est l’intérêt de réaliser des arbres avec les séquences de l'ARNr? Quels-sont les ARNr présents dans les génomes de procaryotes? A quelle(s) sous-unité(s) ribosomique sont-ils associés?
</pre>
<!--
    La grande sous-unité ribosomique 50S des procaryotes contient les ARNr suivants :
        ARNr 23S (2904 nucléotides chez E. coli1) ;
        ARNr 5S ; il n'est pas lié à l'ARNr 23S.
 
    La petite sous-unité ribosomique 30S contient l'ARNr suivant :
        ARNr 16S (1541 nucléotides chez E. coli2).
-->
 
Suivre : [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Phylog%C3%A9nie_ARNr Phylogénie ARNr]
 
<!--
===Annotation des ARNr===
Nous utilisons le logiciel ''rnammer'' pour annoter les ARNr (lsu, ssu, tsu) dans les génomes.
 
<pre style="color:green;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
search_module rnammer
srun --pty bash
module load bioinfo/rnammer-1.2
rnammer -S bac -m ssu -f ~/work/Prochlorococcus/prokka/Aaaa/Aaaa_ssu.rrna < /home/formation/public_html/M2_Phylogenomique/data/Prochlorococcus/DNA/Aaaa.fas
</pre>
Vous allez procéder comme précédemment, avec un script donné à ''sarray'', pour réaliser le rnammer sur tous les fichiers et les trois types d'ARNr.
 
MSK
<syntaxhighlight lang="bash">
for s in Prochlorococcus Synechococcus
do
  for t in ssu lsu tsu
  do   
    for i in /home/formation/public_html/M2_Phylogenomique/data/$s/DNA/*.fas 
    do   
      genome=$(basename "$i" .fas)
      output="~/work/$s/prokka/"$genome"/"$genome"_"$t".rrna"
      echo "module load bioinfo/rnammer-1.2; rnammer -S bac -m $t -f $output < $i;" 
    done
  done
done > rnammer.sh
cat rnammer.sh
</syntaxhighlight>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
sarray -J mkdb -o %j.out -e %j.err -t 01:00:00 --cpus-per-task=1 rnammer.sh
squeue -l -u <user>
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
 
Vérifiez que les fichiers de sortie ne sont pas vide!
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
ls -l ~/work/*/prokka/Aaa*/Aaa*su*.rrna
</pre>
Concaténer les fichiers: 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mkdir ~/work/ProchlorococcusSynechococcus/rRNA
cat ~/work/Prochlorococcus/prokka/Aaa*/Aaa*lsu.rrna ~/work/Synechococcus/prokka/Aaa*/Aaa*lsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/lsu.fas
cat ~/work/Prochlorococcus/prokka/Aaa*/Aaa*ssu.rrna ~/work/Synechococcus/prokka/Aaa*/Aaa*ssu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/ssu.fas
cat ~/work/Prochlorococcus/prokka/Aaa*/Aaa*tsu.rrna ~/work/Synechococcus/prokka/Aaa*/Aaa*tsu.rrna > ~/work/ProchlorococcusSynechococcus/rRNA/tsu.fas
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 4.2:
Combien de gènes codant pour les gènes d'ARNr sont prédits dans les différentes souches?
Commentez.
</pre>
 
===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 ([https://academic.oup.com/mbe/article/30/4/772/1073398 Katoh ''et al.,'' 2103]). Nous utilisons la version ''mafft'' pour des raisons de rapidités.
 
ssu
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
module load bioinfo/mafft-7.313
mafft  --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/ssu.fas > ~/work/ProchlorococcusSynechococcus/rRNA/ssu.aln
</pre>
lsu
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
module load bioinfo/mafft-7.313
mafft  --globalpair --thread 1 ~/work/ProchlorococcusSynechococcus/rRNA/lsu.fas > ~/work/ProchlorococcusSynechococcus/rRNA/lsu.aln
</pre>
tsu
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
module load bioinfo/mafft-7.313
mafft --globalpair ~/work/ProchlorococcusSynechococcus/rRNA/tsu.fas > ~/work/ProchlorococcusSynechococcus/rRNA/tsu.aln
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 4.3:
Pensez-vous que les alignements auraient été de meilleure qualité avec mafft-qinsi et l'option --maxiterate 1000?
</pre>
 
===Arbre avec seaview===
Utilisez le logiciel seaview pour calculer les arbres avec les trois types ARNr.
 
Expérimentez plusieurs méthodes avec différents paramètres.
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 4.4:
Comparez les résultats obtenus.
</pre>
Éditez les fichiers pour ne retenir qu'une seule copie de chaque gènes par souche.
Renommer les séquences par le code à quatre lettres.
 
Concaténez les trois types d'ARNr et calculer l'arbre avec la méthode de votre choix.
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Discutez ces résultats.
</pre>
Code R pour obtenir une illustration des réarrangements présents entre deux arbres (source: [http://blog.phytools.org/2016/08/finding-association-between-two-trees.html phytools blog]).
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
library('phytools')
ta <-read.tree(file='all_mod-PhyML_tree.ph')
tl <-read.tree(file='lsu_mod-PhyML_tree.ph')
ts <-read.tree(file='ssu_mod-PhyML_tree.ph')
 
plot.cophylo(cophylo(ta,tl,rotate=TRUE),fsize=0.7, link.type="curved", link.col="blue")
plot.cophylo(cophylo(ta,ts,rotate=TRUE),fsize=0.7, link.type="curved", link.col="blue")
</pre>
 
===Arbre SSU avec IQ-TREE===
[http://www.iqtree.org/doc/ IQ-tree] doc.
 
IQ-TREE utilise ModelFinder ([https://www.nature.com/articles/nmeth.4285 Kalyaanamoorthy et al., 2017]) pour sélectionner le meilleur modèle adaptés aux données. 
 
Pour seulement trouver le modèle le mieux adapté sans faire de reconstruction d'arbre, utilisez :
<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/iqtree-1.6.7
iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy  -m MF -redo -AIC
</pre>
Les résultats sont dans le fichier : ssu_renamed_simplified.phy.iqtree.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
grep 'Best-fit model' ssu_renamed_simplified.phy.iqtree
</pre>
lsu 
ssu GTR+F+R2
tsu K2P+G4
 
[http://www.iqtree.org/doc/Substitution-Models#dna-models dna-models]
 
====Évaluation des supports de branches avec approximation bootstrap ultra-rapide (UFBoot):====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy  -pre ssuGTRFR2bb1000bnni -m GTR+F+R2 -bb 1000 -redo -bnni -nt AUTO"
</pre>
NOTE: les valeurs de support de l'UFBoot ont des interprétations différentes de celles du bootstrap non paramétrique. Suivez le lien [http://www.iqtree.org/doc/Frequently-Asked-Questions#how-do-i-interpret-ultrafast-bootstrap-ufboot-support-values UFBoot support values interpretation] pour plus d'information.
 
====Évaluer les supports de branche avec des tests de branche simple :====
IQ-TREE propose le test du rapport approximatif de vraisemblance de type SH ([https://academic.oup.com/sysbio/article/59/3/307/1702850 Guindon et al., 2010]). 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy -pre ssuGTRFR2bbalrt -m GTR+F+R2 -bb 1000 -alrt 1000 -redo -nt AUTO"
</pre>
 
====Évaluation des supports de branche avec un bootstrap non paramétrique standard :====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
iqtree -s ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy -pre ssuGTRFR2alrtb -m GTR+F+R2 -alrt 1000 -b 100 -redo -nt AUTO"
</pre>
 
===Arbre SSU avec FastTree===
[http://www.microbesonline.org/fasttree/ FastTree] doc.
<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/FastTree-2.1.10
fasttree -nt -gtr < ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified.phy > ~/work/ProchlorococcusSynechococcus/rRNA/ssu_renamed_simplified_fasttree.ph 
</pre>
Comparez et commentez les résultats obtenus avec IQ-TREE et FastTree.
-->
 
==Arbre espèces: préparation des fichiers==
'''Support de cours :''' [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/2020_supports/ 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.
 
Suivre : [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Arbre_esp%C3%A8ces Arbre espèces]
<!--
===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.
<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/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
</pre>
====Extraire les groupes de gènes orthologues====
 
Exemple de la création s'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">
mkdir -p ~/work/ProchlorococcusSynechococcus/OG/DNA
</pre>
Créer le fihier suivant commande.sh
<pre>
#!/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
</pre>
<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
</pre>
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">
grep -c '>' ~/work/ProchlorococcusSynechococcus/OG/DNA/*fas
</pre>
<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:
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 ?
</pre>
====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]).
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mkdir ~/work/ProchlorococcusSynechococcus/OG/alignment
</pre>
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
</pre>
 
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:
<pre>
pep_OG_1471.fas       peptides
ali_pep_OG_1471.fas   alignement des peptides
ali_dna_OG_1471.fas   alignement des nucléotides
</pre>
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 "Mon_Login"
<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_loop.pl --directory ~/work/ProchlorococcusSynechococcus/OG/DNA --outdir ~/work/ProchlorococcusSynechococcus/OG/alignment --sample 100
</pre>
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">
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?
</pre>
 
===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">
Question 5.3:
Quels sont les traitements réalisés par le logiciel trimal? 
Quelles sont les options possibles pour réaliser ces traitements ? Ne pas lister toutes les options, les regrouper par type de traitement!
Connaissez-vous d'autres logiciels réalisant un traitement comparable des alignements multiples?
</pre>
====Evaluation des alignements avec t_coffee====
<syntaxhighlight lang="bash">
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
</syntaxhighlight>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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
</pre>
 
Sélection alignements avec un SCORE=1000.
<syntaxhighlight lang="bash">
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/*.ali; 
do 
 ip=$(basename $i .aln)     
 grep -H 'SCORE=1000' $i
done
</syntaxhighlight>
 
====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é.
 
genologin softwares : [http://bioinfo.genotoul.fr/index.php/resources-2/softwares/?searchll=trimal trimal] 
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
module load bioinfo/trimal-1.4.1
n=1
for i in ~/work/ProchlorococcusSynechococcus/OG/alignment/ali_dna_OG_*.fas
do
echo -e "Nb:" $n
echo -e "Alignement:" $i "\n"
dir=$(dirname $i)
prefix=$(basename $i .fas)
trimal -in $dir/$prefix.fas -out $dir/$prefix.trim.aln -sident -sgt
((n++))
done > ~/work/ProchlorococcusSynechococcus/OG/alignment/statistics.txt
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.4:
Editez le fichier de résultats et analysez son contenu.
</pre>
=====Pour chaque alignement, extraction des % de sites sans indels et de la conservation moyenne des alignements=====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.5:
Vérifiez que vous avez extrait les informations attendues.
</pre>
 
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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
</pre>
=====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 meilleurs alignements
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
module load system/R-3.6.1
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")
quit()
</pre>
 
Vous pouvez utiliser filezilla (sous Windows) ou un scp (sous linux) pour télécharger le fichier pdf en local.
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.6:
Éditez les alignements de plus mauvaises qualités.
Commentez.
</pre>
 
=====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.
 
Choisir --gap dans {0,100} et --identity dans {0, 1}:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
~/work/scripts/trimal_selection.pl --alignment ~/work/ProchlorococcusSynechococcus/OG/alignment --gap {0,100} --identity {0, 1}
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.7:
Combien d'alignements avez-vous retenu?
</pre>
-->
 
==Super-alignement (Jeudi)==
 
Comme dans l’article publié en 2018 de Yan et al., nous utiliserons 31 gènes du core genome tirés au hasard.
Nous partirons des arbres effectués sur le super-alignement protéique de ces gènes et sur celui retranscrits en nucléotides.
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.8: 
Proposez une méthode pour obtenir le super-alignement protéique de ces 31 gènes concaténés.
</pre>
<!--
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.9 : 
De la même façon comment obtenons nous la « rétro-traduction » en codons ? Dans quels cas cela peut-il être intéressant d’effectuer l’alignement en nucléotides de cette manière ?
</pre>
-->
====Concaténation de 31 alignements====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/work/scripts/concat_aligments.pl  --alignments /home/formation/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70.lst  --outfile ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 -nb_ali 31
</pre>
Liste des alignements retenus:
 /home/formation/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31.lst
Aller regarder le fichier de sortie. Est-il conforme à l'attendu ?
 
====IQ-TREE====
genologin softwares : [http://bioinfo.genotoul.fr/index.php/resources-2/softwares/?searchll=iq-tree iq-tree] 
 
FAQ : http://www.iqtree.org/doc/Frequently-Asked-Questions
 
Documentation : http://www.iqtree.org/doc/
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mkdir ~/work/ProchlorococcusSynechococcus/phyloG/
cd ~/work/ProchlorococcusSynechococcus/phyloG/
cp ~/work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31 .
</pre>
Nous allons inférer un arbre à partir du super-alignement en codons:
Créer le fichier condTree.sh contenant les lignes suivantes. Lancez le, puis répondez aux questions suivantes car cela va être assez long, vous y reviendrez ensuite. Surtout faites un sbatch --cpus-per-task=4 à partir de genologin pour le lancer car nous avons demandé 4 slots de calcul.
N'oubliez pas de spécifier le modèle dans la ligne de commande vous gagnerez 8 heures de calcul ;-)
 
Bonne pratique pour quand vous travaillerez : prenez l'habitude de vérifier la version la plus récente installée sur genologin. Pensez à vérifier si elle doit être mise à jour, si oui, demander la mise à jour. Si non, utilisez la dernière version.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
#!/bin/bash
module load bioinfo/iqtree-2.0.6
iqtree -s ~/work/ProchlorococcusSynechococcus/phyloG/alignment_0.4_70_31 -redo -bb 1000 -alrt 1000 -st CODON -nt 4 -m KOSI07+F+R6
</pre>
<!--
[http://www.iqtree.org/doc/Substitution-Models Substitution-Models]
<pre>
Akaike Information Criterion:           KOSI07+F+R7
Corrected Akaike Information Criterion: KOSI07+F+R7
Bayesian Information Criterion:         KOSI07+F+R6
Best-fit model: KOSI07+F+R6 chosen according to BIC
</pre>
Empirical codon model ([https://academic.oup.com/mbe/article/24/7/1464/986344 Kosiol et al., 2007]).
-->
Pour vous, nous avons inféré un arbre à partir du super-alignement protéique, la ligne de commande lancée était : 
 iqtree -s all_pep.fas  -bb 1000 -alrt 1000 -nt 4
 
Vous trouverez les fichiers générés ici : /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/all_pep.fas.*. Vous pouvez les copier dans votre répertoire ~/work/ProchlorococcusSynechococcus/phyloG/.
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.10: 
Pouvez-vous m’expliquer ce qu’on a demandé à IQTREE dans les deux cas ? Quels ont été les modèles choisis ? Pouvez-vous les expliquer ?
Ceci peut vous aider : 
http://www.iqtree.org/doc/Substitution-Models
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.11:
Comparez les deux arbres (*.treefile) visuellement avec figtree par exemple (ou avec des outils en ligne tels que : phylo.io, PhyD3, iTOL...).
Pour afficher les valeurs de bootstraps, il est nécessaire que figtree les charge comme « labels », ainsi dans la rubrique « node labels » il vous faudra sélectionner « labels » dans le menu déroulant « display ».
Mettez l'arbre généré dans le rapport de TP.
Quelles sont les principales différences ? Que pensez-vous de ces arbres ? Leurs supports ? Leurs congruences ? Les grands clades connus sont-ils retrouvés ? Comparer avec les arbres des articles ou revues suivant(e)s :
La revue de Biller et al : Prochlorococcus  : the structure and function of collective diversity
L’article de Kettler et al. : Patterns and implications of gene gain and loss in the evolution of Prochlorococcus et celui de Yan et al. : Genome rearrangement shapes Prochlorococcus ecological adaptation.
N’oubliez pas de commenter les valeurs de support des nœuds. Pensez à enraciner les arbres en utilisant l’outgroup Synechococcus.
</pre>
 
NB: Les noms des espèces et les informations des clades correspondants sont disponibles dans le fichier /home/formation/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt ou dans le tableau [http://www.m2p-bioinfo.ups-tlse.fr/p/Atelier_Phylog%C3%A9nomique#Disponibilit.C3.A9_des_g.C3.A9nomes ici].
 
===Super-arbres===
 
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à :
 #!/bin/bash
 iqtree -s ./ssu_renamed_simplified.aln -nt 1 -AIC -bb 1000 -alrt 1000 -redo
 
Et on a tapé :
 module load bioinfo/iqtree-2.0.6
 sbatch --mem=20G ssuTree.sh
 
Voici l’arbre obtenu : /home/formation/work/ProchlorococcusSynechococcus/phyloG/ssu_renamed_simplified.aln.treefile
 
'''A vous :''' 
 
Loguez-vous sur genologin.
 
Lançons les arbres protéiques. Pour cela faire un seul script que vous appellerez ind_pep_trees.sh et qui écrira toutes les commandes en bouclant sur chaque fichier d’entrée. Le but étant d’obtenir une ligne par commande iqtree sur un fichier d’alignement protéique.  
Les fichiers d’input sont /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/*_renamed.fas.
 
Il vous faut les copier dans votre work car sinon iqtree écrira les fichiers d'output sur le répertoire du compte de formation et pas dans le votre...
 
Pour vous aider inspirez vous de la réponse à la question "How to generate an sarray command file with bash for single (fastq) file ?" sur la page http://bioinfo.genotoul.fr/index.php/faq/bioinfo_tips_faq/.
 
'''Attention à bien spécifier -nt 1. Si vous souhaitez utiliser plus d’un CPU il faut le réserver avec l’option --cpus-per-task dans la commande sbatch / sarray. Pas besoin d’augmenter la RAM pour les protéines. Pensez aussi à mettre -AIC comme critère de sélection de modèle.'''
 
Regarder le contenu de ind_pep_trees.sh. Est-il correct ?
 
Rajouter la première ligne obligatoire sur genologin avec vi par exemple : 
 #!/bin/bash 
 
Pour le lancer en parallèle sur plusieurs nœuds du cluster :
<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/iqtree-2.0.6
sarray ind_pep_trees.sh
</pre>
 
Pour monitorer votre job :
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
squeue -l -u <login>
</pre>
 
Pour le killer si besoin : 
 scancel jobid
 
Quand tous vos jobs sont terminés, vérifiez que les fichiers de sorties ne soient pas vides et que ça s’est bien passé en faisant par exemple :
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
tail *_renamed.fas.log
</pre>
 
Vous pouvez aussi regarder les modèles sélectionnés :
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
grep 'Best-fit model:' *_renamed.fas.log
</pre>
 
Concaténer tous les arbres (les 31 arbres protéiques et l’arbre ARNr obtenu hier) avec la commande cat.
Nommez le fichier alltrees.tree.
 
'''Test de plusieurs méthodes de super-arbres :'''
 
Commençons par la méthode la plus répandue : le MRP.
 
Pour aller sur un nœud :
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
</pre>
puis taper :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
module load system/R-3.5.1
R
library(phytools)
trees=read.tree("./alltrees.tree")
supertrees<-mrp.supertree(trees,rearrangements="SPR", start="NJ")
</pre>
Vous avez obtenu les super-arbres les plus parcimonieux. Sauvez-les en utilisant la fonction write.tree de R.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
write.tree(supertrees, file = "./superTrees.tree")
quit()
</pre>
Nous sommes sortis de R.
 
Dans notre cas, nous avons un seul arbre le plus parcimonieux, mais nous aurions pu en obtenir plusieurs.
 
Utiliser iqtree pour obtenir le consensus des 32 arbres avec la règle majoritaire étendue.
 
Pour cela restez sur le nœud et tapez :
<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/iqtree-2.0.6
iqtree -con -t alltrees.tree -nt 1
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.12:
D’après vous que signifie les valeurs aux nœuds sur cet arbre consensus ?
</pre>
 
Lancez aussi l’arbre consensus en réseau : 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
iqtree -net -t alltrees.tree -nt 1
</pre>
 
Et visualisez-le avec splitstree en local. 
 
https://software-ab.informatik.uni-tuebingen.de/download/splitstree5/welcome.html
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.13:
Commentez ce réseau par rapport aux autres arbres obtenus. Qu’en pensez-vous ?
</pre>
 
===Comparaison des arbres===
 
Concaténez maintenant les deux arbres de super-matrice (celui sur les codons que vous avez lancé toute à l'heure et celui sur les protéines) ainsi que les deux super-arbres (le consensus et le MRP).
 
Si vous avez eu des difficultés à obtenir l'arbre à partir de l'alignement en codon il est disponible ici : /home/formation/work/ProchlorococcusSynechococcus/phyloG/alignment_0.4_70_31.fas*.
 
Attention marquez quelques part l'ordre avec lequel vous les avez concaténés pour créer le fichier d'arbres à comparer afin de vous en souvenir ensuite.
Lancer ensuite le calcul de la distance de Robinson and Foulds sur ce fichier avec iqtree entre les 4 arbres 2 à 2.
 
Attention à faire le module load bioinfo/iqtree-2.0.6
 
Et à rajouter -nt 1 dans les options.
Restez sur le nœud ou faites un sbatch, mais ne lancez rien sur le nœud maître.
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5.14: 
Commentez les résultats. Quel est l’arbre le plus différent des autres ? Qu’est-ce qui pourrait l’expliquer ?
</pre>
 
==Analyse des liens phylogénétiques des gènes du OG5 de PortMCL==
<!--
Éditez un fichier avec la liste des gènes à étudier.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
~/work/scripts/extract_sequences_from_list.pl -fasta ~/work/ProchlorococcusSynechococcus/DNA/combined_dna.ffn --outname ~/work/ProchlorococcusSynechococcus/OG5/OG5_genes.faa --list ~/work/ProchlorococcusSynechococcus/OG5_gene_list.txt 
 
~/work/scripts/aa_to_dna_aln.pl -dna ~/work/ProchlorococcusSynechococcus/OG5/OG5_genes.faa --outdir ~/work/ProchlorococcusSynechococcus/OG5
</pre>
Réalisez l'analyse phylogénétique avec le logiciel de votre choix. Reportez sur l'arbre obtenu les différents évènements de gain/pert/HT de gènes.
-->
 
==Reconstruction d'états ancestraux==
===Introduction===
Si certaines combinaisons de caractères sont systématiquement associées à plusieurs espèces, cela pourrait suggérer que des forces évolutives, comme la sélection, ont façonné ces associations. Toutefois, les associations non aléatoires de certains traits chez certaines espèces peuvent être dues à l'héritage commun de leurs ancêtres et, par conséquent, les changements concomitants dans le temps ne peuvent être déduits. 
 
Si les caractères ont évolué de façon aléatoire sans association, les espèces plus étroitement apparentées sont plus susceptibles d'être semblables que les autres, créant ainsi des relations apparentes entre les caractères.
 
Il est donc nécessaire de considérer les relations phylogénétiques entre les espèces lors de l'analyse de leurs caractères. 
 
Deux questions peuvent être abordées lors de l'intégration de la phylogénie dans les données comparatives :
* prise en compte de la non-indépendance inter-espèces dans l'étude des caractères et de leurs relations,
* estimation des paramètres d'évolution des caractères.
Ces deux questions sont étroitement liées. En effet l'impact de la phylogénie sur la distribution des caractères dépend non seulement de la phylogénie mais aussi de la façon dont ces caractères évoluent.
 
Suivre : [http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique_Etats_ancestraux Etats ancestraux]
<!--
===Répertoire===
Nous allons travailler dans un nouveau répertoire: 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mkdir ~/work/ProchlorococcusSynechococcus/Ancestral_Characters 
cd ~/work/ProchlorococcusSynechococcus/Ancestral_Characters
cp /home/formation/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt ~/work/ProchlorococcusSynechococcus/Ancestral_Characters
</pre>
 
===Copier l'arbre espèce de référence===
Choisissez l'arbre que vous allez utiliser comme référence et copiez le dans le nouveau répertoire.
 
Comme alternative, vous pouvez aussi copier celui-ci:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cp /home/formation//work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31_nuc_GTR+F+R4_rooted.treefile ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph
</pre>
 
===Composition en GC des génomes===
Comme caractère continue, nous allons étudier l'évolution de la composition en G+C des génomes. Le taux de GC est calculé avec ''gc_freq.pl''. 
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
for i in ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.fna
do      
~/work/scripts/gc_freq.pl --file $i
done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt
 
for i in ~/work/Synechococcus/prokka/Aaa*/Aaa*.fna
do      
~/work/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
</pre>
 
===Arbre espèces===
Nous allons utiliser le logiciel [https://www.r-project.org/ ''R''] pour réaliser nos analyses et les librairies ''ape'' et ''phytools''.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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)
</pre>
Vérifiez que l'arbre tr est enraciné:
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
is.rooted(tr)
 
plot(tr)
plot(tr <- root(tr, interactive = TRUE)) # pour changer la racine de façon interactive
</pre>
Informations sur les souches
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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,]
</pre>
Changement des labels de l'arbre
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
tr$tip.label <- as.character(strains_info$Strain)
tipcol <- c()
tipcol[1:4] <- 'red'
tipcol[5:16] <- 'blue'
 
pdf(file="~/work/ProchlorococcusSynechococcus/images/species_tree.pdf", paper="a4r")
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.2, 0.5), cex=1)
add.scale.bar(length=0.1)
dev.off()
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
evince ~/work/ProchlorococcusSynechococcus/images/species_tree.pdf
</pre>
 
===GC et taille des génomes===
Lire les fichiers arbre et données:
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
library(phytools)
 
treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 
tr <- read.tree(treefile)
 
strains_gc_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt"
data <- read.table(file=strains_gc_file, header=F,sep="\t", row.names=1)
</pre>
Ajout des noms des colonnes et réordonnement des lignes de ''strains_gc'' en fonction des feuilles de l'arbre ''tr''.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/length_GC.pdf", paper="a4r")
colnames(data) <- c('Length', 'GC')
data<- data[tr$tip.label,]
plot(data$Length, data$GC, pch=20, ylim=c(0,1), main="Length/GC")
dev.off()
</pre>
Number of genes
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
for file in ~/work/Prochlorococcus/prokka/Aaa*/*.txt ~/work/Synechococcus/prokka/Aaa*/*.txt
do
 prefix=$(basename $file .txt)
 awk -v pf=$prefix '{
  if (match($0,"gene:")) {
   printf("%s\t%d\n", pf, $2);
  }
 }' $file >> ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_genes.txt
done
</pre>
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
strains_genes_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_genes.txt"
pdf(file="~/work/ProchlorococcusSynechococcus/images/length_genes.pdf", paper="a4r")
genes <- read.table(file=strains_genes_file, header=F,sep="\t", row.names=1)
genes<- genes[tr$tip.label,]
plot(data$Length, genes, pch=20, ylim=c(0,3000), main="Length/#genes", xlab="Length", ylab="#genes")
dev.off()
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.1:
Existe-t-il un lien entre la taille du génome et son contenu en GC? Commentez.
Existe-t-il un lien entre la taille du génome et son contenu en gènes?
</pre>
 
====Cartographie de l'évolution d'un caractère continue sur l'arbre====
La fonction trace l'arbre sur lequel est reporté l'évolution d'un  caractère continu. La cartographie est réalisée en estimant les états aux nœuds internes à l'aide du ML avec la fonction ''fastAnc'', et l'interpolation des états le long de chaque branche (Felsenstein,  1985). 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
gc <- data$GC
names(gc) <- rownames(data)
ln <- data$Length
names(ln) <- rownames(data)
 
pdf(file="~/work/ProchlorococcusSynechococcus/images/length_GC_length.pdf", paper="a4r")
XX <- contMap(tr, gc, sig=2, outline=F, lwd=5, lims=c(0.3,0.6))
XX <- contMap(tr, ln, sig=2, outline=F, lwd=5)
dev.off()
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.2:
Commentez les figures obtenues.
</pre>
 
====Reconstruction par Maximum de vraisemblance====
Il existe plusieurs fonctions R utilisant le ML pour estimer les états ancestraux aux nœuds internes pour les caractères continus.: ''ace(method="ML")'', ''fastAnc'', et ''anc.ML''. Les trois méthodes effectuent l'estimation du ML en supposant un modèle brownien pour modéliser l'évolution du caractère. La seule différence est la méthode d'optimisation. fastAnc, par exemple, utilise une méthode de ré-enracinement dans laquelle l'arbre est ré-enraciné à chaque nœud interne et l'algorithme de contraste est utilisé pour obtenir l'état ML pour ce noeud. ''anc.ML'' utilise une optimisation numérique, mais initie la recherche avec les valeurs obtenues avec ''fastAnc''.
 
Comparaison des trois méthodes sur notre exemple (pour chaque méthode '<method>$ace' est un vecteur contenant les états des nœuds internes)
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
aceML <- ace(gc,tr,type="continuous",method="ML")
fAnc <- fastAnc(tr,gc,vars=TRUE,CI=TRUE)
ancML <- anc.ML(tr,gc)
obj<-cbind(aceML$ace,fAnc$ace,ancML$ace)
colnames(obj)<-c("ace(method=\"ML\")","fastAnc","anc.ML")
pdf(file="~/work/ProchlorococcusSynechococcus/images/ML_fastAnc_anc.pdf", paper="a4r")
pairs(obj,pch=21,bg="grey",cex=1.5)
dev.off()
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.3:
Que peut on conclure de cette analyse?
</pre>
Tracé de l'arbre
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_aceML.pdf", paper="a4r")
plotTree(tr)
nodelabels(pie=aceML$ace,cex=0.5)
dev.off()
</pre>
 
===Lien entre caractères===
Felsenstein a proposé une méthode qui prend en compte la phylogénie dans l'analyse des données comparatives. L'idée qui sous-tend sa méthode des '''contrastes''' est que, si nous supposons qu'un trait continu évolue de façon aléatoire dans n'importe quelle direction (modèle de mouvement brownien), le contraste entre deux espèces devrait avoir une distribution normale avec une moyenne nulle et une variance proportionnelle au temps écoulé depuis la divergence. Si les contrastes sont mis à l'échelle avec ces derniers, alors ils ont une variance égale à un.
 
La fonction pic (phylogenetically independent contrasts)
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
lns <- ln/10000000
pic.gc <- pic(gc, tr)
pic.lns <- pic(lns, tr)
 
pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_contrasts.pdf", paper="a4r")
plot(tr)
nodelabels(round(gc, 2), adj = c(0, -0.5), frame = "n")
nodelabels(round(lns,2), adj = c(0, 1), frame = "n")
tiplabels(round(gc, 2), adj = c(-1, -0.5), frame = "n")
tiplabels(round(lns,2), adj = c(-1, 1), frame = "n")
dev.off()
</pre>
Lien entre les pairs de contrastes
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">plot(pic.lns, pic.gc)
pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_pairs_contrasts.pdf", paper="a4r")
plot(pic.lns, pic.gc, xlim=c(-0.15,0.15), ylim=c(-0.15,0.15))
abline(a = 0, b = 1, lty = 2) # x = y line
 
cor(pic.gc, pic.lns)
lm(pic.lns ~ pic.gc)
dev.off()
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.4:
Existe-t-il un lien entre la taille du génome et son contenu en GC (coefficients de la régression)? Commentez.
</pre>
Remarque, faire la régression entre les PICs en passant par l'origine est justifiée si les caractères évoluent sous un modèle de mouvement brownien et s'il existe une relation linéaire entre eux.
 
===Autocorrélation phylogénétique===
masquée
 
Comme les espèces ne sont pas indépendantes par leurs relations phylogénétiques,
nous pouvons utiliser cette dépendance pour quantifier l'association entre les variables observées au niveau des espèce. Gittleman et Kot[109] propose une méthode basée sur
une approche d'autocorrélation. Elle utilise l'indice d'autocorrélation I[208] de Moran qui fait intervenir ''wij'' un poids qui quantifie la proximité phylogénétique entre espèces ''i'' et ''j''. Avec un arbre phylogénétique, nous avons ''wij'' = 1/''dij'', où ''dij'' sont les distances mesurées sur un arbre, et ''wii'' = 0.
 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">w <- 1/cophenetic(tr)
diag(w) <- 0
Moran.I(gc, w)
</pre>
Le résultat est une liste avec quatre éléments : 
* la valeur observée de I (observée),
* sa valeur attendue sous l'hypothèse nulle d'aucune corrélation,
*l'écart-type de l'I (sd) observé, 
*la P-valeur de l'hypothèse nulle.
 
La fonction ''gearymoran'' d'''ade4'' calcule le coefficient de Moran et test sa signification à l'aide d'une procédure de randomisation. L'option ''nrepet'' spécifie le nombre de réplications de la randomisation (999 par défaut). 
 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
library(ade4)
gearymoran(w, data.frame(gc, ln))
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.5:
Existe-t-il un lien entre l'évolution des caractères (taille du génome et contenu en GC) et l'évolution des souches? Commentez.
</pre>
 
 
===Habitats===
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.6:
Rappelez d'ou proviennent les différents isolats analysés et expliquez comment ont été défini les écotypes.
</pre>
 
Nous allons maintenant nous intéresser à la profondeur à laquelle les organismes vivent. Cette information est absente pour la souche aaap. Nous allons donc supprimer cette feuille à l'arbre.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
tr_d <- drop.tip(tr, "Aaap")
strains_info_d <- strains_info[rownames(strains_info) != 'Aaap', ]
dp <- strains_info_d$Depth
names(dp) <- rownames(strains_info_d)
 
pdf(file="~/work/ProchlorococcusSynechococcus/images/Habitats_profondeur.pdf", paper="a4r")
XX <- contMap(tr_d, dp)
dev.off()
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.7:
Existe-t-il un lien entre profondeur et l'évolution des souches (voir Biller et al., 2015)? Commentez.
 
</pre>
 
===Ecotypes===
Codage des ecotypes HL et LL
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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
</pre>
Reconstruction des états du caractère aux nœuds de l'arbre:
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/Ecotypes.pdf", paper="a4r")
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)
dev.off()
</pre>
Vous pouvez utiliser toutes l'information des ecotype et refaire l'inférence des états du caractère aux nœuds de l'arbre.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
ecotypes <- strains_info$Light
</pre>
 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/Ecotypes_full.pdf", paper="a4r")
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)
dev.off()
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.8:
Sur quelles branches de l'arbre placez vous les transitions d'écotypes? Commentez.
</pre>
 
===Profils des gènes===
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
R
library(ape)
</pre>
Lecture de l'arbre espèces
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 
tr <- read.tree(treefile)
</pre>
 
====Matrice de présence/absence====
Nous allons utiliser les groupes de gènes orthologues calculés par panOCT. 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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)
</pre>
=====Transposer la matrice et ordonner les lignes comme les feuilles de l'arbre=====
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
t_matchtable <- t(matchtable)
t_matchtable <- t_matchtable[tr$tip.label,]
 
nb_strains <- nrow(t_matchtable)
</pre>
=====Compiler les différents motifs observés dans une liste=====
Plusieurs groupes de gènes orthologues peuvent avoir le même profil phylogénétique. Il n'est donc pas nécessaire de faire les calculs sur la totalité.
<pre>
motifs[[pattern]]$sring : motif associé aux souches
motifs[[pattern]]$occurences : occurence du motif 
</pre>
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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")
}
</pre>
 
====Définir un ensemble de motifs====
Pour la mise en main des scripts, il peut être utile de travailler avec un sous ensemble de motifs!
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
nb_start <- 1
nb_end <- 10
if ( nb_end > length(motifs)) nb_end <- length(motifs) 
</pre>
ou
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">nb_start <- 1
nb_end <- length(motifs) 
</pre>
=====Taille des patterns trouvés=====
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
nb <- nb_start
pat_size <- c()
pat_name <- c()
i <- 1
 
while (nb <= nb_end ) {
   pattern <- names(motifs)[nb]
   cat(pattern, "\n")
   size <- length(motifs[[pattern]]$occurences)
   motifs[[pattern]]$size <- size
   if ( size > 10 ) {
      pat_size[i] <- size
      pat_name[i] <- pattern
      #cat(pattern, size, "\n")
      i <- i+1
   }
   nb <- nb+1
}
names(pat_size) <- pat_name
pat_sort <- sort(pat_size, decreasing=T)
 
pdf(file="~/work/ProchlorococcusSynechococcus/images/taille_des_motifs.pdf", paper="a4r")
par(mar=c(5,10,1,1))
barplot(pat_sort, horiz=T, cex.names = 0.6, las=1, cex=1)
dev.off()
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.9:
Quel sera l'usage de ce calcul?
</pre>
 
===Inférence des états ancestraux===
Les différentes fonctions utilisent des paramètres similaires.
 
Le nombre d'états est déduit des données.
 
Le modèle est spécifié avec une matrice d'entiers représentant les indices des paramètres : 1 représente le premier paramètre, 2 le second, et ainsi de suite. 
 
Vous devrez choisir la matrice qui spécifie les probabilités de transition entre les états de votre caractère discret. Ace propose des raccourcis pour 
* un modèle à un paramètre (ER) où tous les taux de transitions sont égaux,
* un modèle symétrique (SYM) dans lequel les taux de transitions dans les deux sens sont égaux, 
* un modèle où tous les taux sont différents (ARD), toutes les transitions possibles entre les états reçoivent des paramètres distincts. 
Vous pouvez également spécifier votre propre matrice. Pour indiquer qu'une transition est impossible, un zéro doit être donné dans la cellule appropriée de la matrice.
 
====La reconstruction la plus parcimonieuse====
Cette fonction (Most Parsimonious Reconstruction our MPR) réalise la reconstruction des caractères ancestraux par parcimonie, comme décrit dans Hanazawa et al. (1995) et modifié par Narushima et Hanazawa (1997). Ils ont utilisé l'approche proposée par Farris (1970) et Swofford and Maddison’s (1987) pour reconstruire les états ancestraux. Le caractère est supposé prendre des valeurs entières. L'algorithme recherche les ensembles de valeurs pour chaque nœud sous forme d'intervalles avec des valeurs inférieures et supérieures.
 
L'arbre doit être non enraciné et entièrement dichotomique.
 
Exemple
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/plus_parcimonieuse.pdf", paper="a4r")
nb <- 3
pattern <- names(motifs)[nb]
mrp <- MPR(motifs[[pattern]]$sring, unroot(tr), "Aaao")
 
plot(unroot(tr), label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(paste("[", mrp[, 1], ",", mrp[, 2], "]", sep = ""))
dev.off()
</pre>
 
=====Tracé MPR pour chaque motif=====
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/plus_parcimonieuse_all.pdf", paper="a4r")
nb <- nb_start
nb_end < 10
while (nb < nb_end ) {
 
   pattern <- names(motifs)[nb]
   if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) {
      mrp <- MPR(motifs[[pattern]]$sring, unroot(tr), "Aaao")
 
      plot(unroot(tr), label.offset=0.2, show.node.label=F, cex=1, main=paste("MRP", nb))
      tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame="n", adj=c(-0.5, 0.5), cex=1)
      nodelabels(paste("[", mrp[, 1], ",", mrp[, 2], "]", sep = ""), bg="white", col="purple")
 
      #cat ("Press [enter] to continue")
      #line <- readline()
 
  }
   nb <- nb+1
}
dev.off()
nb_end <- length(motifs) 
</pre>
 
====Marquage stochastique des caractères ====
La fonction ''make.simmap'' de ''phytools'' permet de réaliser un marquage stochastique des caractères aux nœuds de l'arbre.
* nsim: number of simulations.
* pi: gives the prior distribution on the root node of the tree, options are equal, estimated, or a vector with the frequencies. 
 
Exemple
<pre>
motif 3 ("0000001100101111")
</pre>
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/marquage_stochastique.pdf", paper="a4r")
nb <- 3
pattern <- names(motifs)[nb]
ms_er_trees <- make.simmap(tr, motifs[[pattern]]$sring, model="ER", nsim=100, pi="estimated")
 
cols<-setNames(c("blue","red"),c(0,1))
plotSimmap(ms_er_trees[[1]], cols)
dev.off()
</pre>
Calculons maintenant les fréquences d'état à partir des marquages stochastiques pour chaque noeud interne. Ce sont nos probabilités postérieures pour les noeuds.
 
Fonction pour calculer les états:
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
map_states <- function(x){
   y <- sapply(x$maps,function(x) names(x)[1])
   names(y) <- x$edge[,1]
   y <- y[as.character(length(x$tip)+1:x$Nnode)]
   return(y)
}
</pre>
Apliquer la fonction à tous les arbres et tracer l'arbre et les probabilités postérieures aux nœuds.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/probabilites_posterieures.pdf", paper="a4r")
AA <- sapply(ms_er_trees, map_states)
piesA <- t(apply(AA,1,function(x,levels,Nsim) summary(factor(x,levels))/Nsim,levels=c(0,1), Nsim=100))
 
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=piesA, cex=0.5, piecol=cols)
dev.off()
</pre>
 
====Maximum de vraisemblance====
La commande ace de la librairie ''ape'' peut reconstruire des états ancestraux pour des caractères discrets en utilisant le maximum de vraisemblance (Pagel 1994 Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters. Proceedings of the Royal Society of London. Series B. Biological Sciences, 255, 37–45.). 
 
Par défaut, la probabilité que les différents états ancestraux des caractères discrets sont calculés à l'aide d'une procédure d'estimation conjointe (joint) en utilisant une procédure semblable à celle décrite dans Pupko et al. (2000).  Si "marginal = VRAI", une procédure d'estimation marginale est utilisée. Avec cette méthode, les valeurs de vraisemblance à un nœud donné sont calculées en utilisant uniquement les informations provenant des feuilles (et des branches) descendantes de ce nœud. 
Avec l'estimation conjointe, toutes les informations sont utilisées pour chaque nœud. 
La mise en œuvre actuelle de l'estimation conjointe utilise un algorithme "en deux passes" qui est beaucoup plus rapide que l'approche stochastique alors que les estimations des deux méthodes sont très proches.
 
Si l'option CI = TRUE est utilisée, alors la probabilité de chaque état ancestral est retournée pour chaque noeud dans une matrice appelée lik.anc. 
For discrete characters only maximum likelihood estimation is available
 
obj<-anc.Bayes(tree,x,ngen=10000) # sample ancestral states
print(obj,printlen=20) ## estimates
 
Remarque:  la fonction ''fitDiscrete'' de ''geiger'' utilise la même syntaxe que ace et donne des résultats similaires, mais les valeurs de log-vraisemblance sont mises à l'échelle différemment cependant le test du rapport de vraisemblance reste le même.
 
Exemple 
motif 3 ("0000001100101111")
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/ML3.pdf", paper="a4r")
nb <- 3
pattern <- names(motifs)[nb]
er <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ER")
sym <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="SYM")
ard <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ARD")
 
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=er$lik.anc, cex=0.3, adj=c(0.5, 0.5))
nodelabels(pie=sym$lik.anc, cex=0.3, adj=c(0.53, 0.5))
nodelabels(pie=ard$lik.anc, cex=0.3, adj=c(0.56, 0.5))
dev.off()
</pre>
La sortie du modèle ER correspond à celle du modèle SYM. Pour un caractère binaire, ces modèles sont identiques, bien qu'ils soient distincts pour les caractères à trois états ou plus. Pour un caractère à trois états, ER est un modèle à un paramètre, SYM un modèle à trois paramètres et ARD un modèle à six paramètres. 
 
Le modèle ARD donne généralement la probabilité la plus élevée. Cependant, il inclut aussi plus de paramètres que le modèle ER. Il est bien connu que l'ajout de paramètres à un modèle augmente généralement sa probabilité. Pour déterminer si l'utilisation du modèle ARD est appropriée, vous devez exécuter un test de vraisemblance.
 
======Test de vraisemblance======
La différence de logarithme de probabilité entre deux modèles est distribuée sous forme de chi2, avec des degrés de liberté égaux au nombre de paramètres qui diffèrent entre les modèles. Ainsi, un test de vraisemblance de ces modèles demande si la différence de vraisemblance est suffisamment grande pour se situer dans la queue la plus à droite de la distribution du chi2.
 
La fonction ''pchisq(value, df)'' donne le pourcentage de la fonction de distribution cumulative pour chi2 qui se trouve à gauche de la valeur donnée pour un degré de liberté souhaité(df). Ainsi,
 
  1-pchisq(2*abs(er$loglik - ard$loglik), 1)
 
donne le pourcentage de la distribution du chi2 qui se trouve à droite de la différence de vraisemblance observée pour les modèles ER et ARD (qui diffèrent d'un paramètre).
======Boucle sur tous les motifs======
Par curiosité nous allons réaliser le test pour tous les motifs et imprimer les valeurs les plus significatives.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2, )
nb <- nb_start
nbpattern <- 1
#while (nb <= length(motifs) ) {
while (nb <= nb_end ) {
   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
      ard <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ARD")
      motifs[[pattern]]$er <- ard
      ml_test <- 1-pchisq(2*abs(er$loglik - ard$loglik), 1)
      if ( ml_test < 0.01 ) {
         cat(nb, pattern, er$loglik, ard$loglik, ml_test, "\n")
      }
   }
   nb <- nb+1
}
</pre>
 
=====Fonction pour annoter les états des feuilles de l'arbre=====
Les inférences des états du caractère au nœuds sont calculées sur les noeuds interne de l'arbre, nous ajoutons ici les noeuds feuilles.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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)
}
</pre>
Exemple d'appel de la fonction
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
leaves <- tip_states(motifs[[pattern]]$sring)
</pre>
Concaténation des deux matrices (feuilles et noeuds)
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
edge_states <- rbind(leaves, er$lik.anc)
</pre>
 
=====Fonction pour annoter les transitions d'états sur les branches de l'arbre=====
Une transition d'état est identifiée sur une branche (''edge'') si l'état du caractère présent au nœud parent est différent de l'état du caractère présent au nœud fils.
*En entrée l'arbre ''tr'' et les états du caractère inférés aux noeuds (''edge_states''). 
*En sortie une liste avec les listes de transitions (type de tansition, direction et une couleur).
 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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 ) {
         transition <- c(transition, i)
         if ( parent_states == 0 ) {
            direction <- c(direction, '+')
            transcolor <- c(transcolor, 'blue')
         } else {
            direction <- c(direction, '-')
            transcolor <- c(transcolor, 'red')
        }
        cat(i, 'parent', parent , parent_states, '->', child_states, 'child', child, direction, "\n")     }
   }
   translist <- list(transition=transition, direction=direction, transcolor=transcolor)
   return(translist)
}
</pre>
Exemple d'appel de la fonction
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
states_transitions(tr, edge_states)
</pre>
=====Tracé pour chaque motif=====
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2)
mycolors <- c('red', 'blue')
 
nb <- nb_start
nbpattern <- 1
#while (nb <= length(motifs) ) {
while (nb < nb_end ) {
 
   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]]$er$lik.anc)
      st <- states_transitions(tr, edge_states)
 
      plot(tr, label.offset=0.2, show.node.label=F, cex=1)
      vec <- as.vector(motifs[[pattern]]$sring)
      bgcol <- unlist(lapply(vec, FUN= function(x) mycolors[x+1]))
      tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5)
      nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors)
      edgelabels(text=st$direction, edge=st$transition, frame="r", bg=st$transcolor, adj=c(0,-0.5),cex=0.2)
      add.scale.bar(length=0.1)
 
      cat ("Press [enter] to continue")
      line <- readline()
 
  }
   nb <- nb+1
}
</pre>
=====Flux de gènes=====
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2)
mycolors <- c('red', 'blue')
 
nb <- nb_start
nbpattern <- 1
gain <- vector (mode='numeric',length=nrow(tr$edge))
loss <- vector (mode='numeric',length=nrow(tr$edge))
 
while (nb < nb_end ) {
 
   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]]$er$lik.anc)
      st <- states_transitions(tr, edge_states)
      nb_trans <- length(st$transition)
      if ( nb_trans > 0 ) {
         for ( j in 1:nb_trans) {
           edge <- st$transition[j]
           cat(nb_trans, j, edge, st$direction[j], "\n")
           if ( st$direction[j] == "+" ) {
            cat('gain', nb_trans, j, st$direction[j], "\n")
              gain[edge] <- gain[edge] + motifs[[pattern]]$size
            } else if ( st$direction[j] == "-" ) {
               loss[edge] <- loss[edge] - motifs[[pattern]]$size
            }
          }
      } else {
          cat("no transition, unexpected!\n")
      }
  }
   nb <- nb+1
}
</pre>
======Gains et pertes de gènes======
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/gains_pertes_genes.pdf", paper="a4r")
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5)
nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors)
edgelabels(text=gain, frame="n", col='blue', bg='white', adj=c(0,-0.5),cex=0.8)
edgelabels(text=loss, frame="n", col='red', bg='white', adj=c(0,1.5),cex=0.8)
add.scale.bar(length=0.1)
dev.off()
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
evince ~/work/ProchlorococcusSynechococcus/images/gains_pertes_genes.pdf
</pre>
 
======Flux de gènes======
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
pdf(file="~/work/ProchlorococcusSynechococcus/images/flux_genes.pdf", paper="a4r")
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
flux <- gain + loss
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5)
nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors)
edgelabels(text=flux, frame="n", col='black', bg='white', adj=c(0,-0.5),cex=0.8)
add.scale.bar(length=0.1)
dev.off()
</pre>
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6.10:
Commentez les figures obtenues.
Comparez ces résultats avec ceux de Sun et Blanchard, 2014)
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
evince ~/work/ProchlorococcusSynechococcus/images/flux_genes.pdf
</pre>
 
===Inférence Bayésienne ===
====BayesTraits====
[http://www.evolution.rdg.ac.uk/BayesTraitsV3.0.1/BayesTraitsV3.0.1.html BayesTraits] is a computer package for performing analyses of trait evolution among groups of species for which a phylogeny or sample of phylogenies is available. This new package incoporates our earlier and separate programes Multistate, Discrete and Continuous. BayesTraits can be applied to the analysis of traits that adopt a finite number of discrete states, or to the analysis of continuously varying traits. Hypotheses can be tested about models of evolution, about ancestral states and about correlations among pairs of traits. 
 
[http://www.evolution.rdg.ac.uk/BayesTraitsV3.0.1/Files/BayesTraitsV3.Manual.pdf BayesTraitsV3.Manual.pdf]
 
La mise en œuvre de ce logiciel est complexe et ne sera pas traitée dans le cadre de cd TP.
-->
 
==Liens utiles==
===Mesquite===
[http://www.mesquiteproject.org/ Mesquite]: A modular system for evolutionary analysis
 
<!--
==TP Phylogénomique du 25 octobre==
 
'''Support de cours :''' [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/2018_supports/ supports]
 
Aujourd’hui nous utiliserons certains gènes de Prochlorococcus et de Synechococcus pour expérimenter les différentes méthodes de reconstruction phylogénomiques. Nous nous initierons à la comparaison d’arbres.
 
 
'''Attention : sur le cluster genotoul :''' toujours utiliser qsub ou qrsh et configurer les jobs pour qu'ils écrivent dans le /work de votre compte fleur.
 
'''Sur le cluster genologin :''' toujours utiliser sbatch ou srun et écrire sur le /work de votre compte fleur.
 
Rappel des comptes fleur : mdp : f1o2r3 !
 
anemone arome aster bleuet camelia capucine chardon clematite cobee coquelicot cosmos cyclamen dahlia digitale geranium gerbera glaieul hortensia
 
 
'''Lien utile pour l’utilisation du cluster :''' [http://bioinfo.genotoul.fr/index.php/faq/ FAQ]
 
'''Lien vers la doc de IQTREE :''' [http://www.iqtree.org/doc/ IQTREE]
 
===I)	Super-alignement===
 
Comme dans l’article publié en 2018 de Yan et al., nous utiliserons 31 gènes du core genome tirés au hasard.
Nous partirons des arbres effectués sur le super-alignement protéique de ces gènes et sur celui retranscrits en nucléotides.
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 1 : Proposez une méthode pour obtenir l’alignement protéique de ces 31 gènes concaténés 
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2 : De la même façon comment obtenons nous la « rétro-traduction » en codons ? Dans quels cas cela peut-il être intéressant d’effectuer l’alignement en nucléotides de cette manière ?
</pre>
 
Vous trouverez l’arbre inféré à partir du super-alignement protéique ici :
/home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/all_pep.fas.treefile.
 
La ligne de commande lancée était : 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s all_pep.fas  -bb 1000 -alrt 1000 -nt 4
</pre>
Et celui inféré à partir du super-alignement en codons ici :
/home/formation/work/ProchlorococcusSynechococcus/phyloG/alignment_0.4_70_31.fas.treefile.
 
La ligne de commande lancée était : 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/IQ-TREE/iqtree-1.6.3-Linux/bin/iqtree -s alignment_0.4_70_31  -bb 1000 -alrt 1000 -st codon -nt 4
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 2 : Pouvez-vous m’expliquer ce qu’on a demandé à IQTREE dans les deux cas ? Quels ont été les modèles choisis ? Pouvez-vous les expliquer ?
</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 : comparez les deux arbres visuellement avec figtree par exemple (ou avec des outils en ligne tels que : phylo.io, PhyD3, iTOL...).
Pour utiliser figtree sur le cluster pensez à se loguer en -XY. Sinon il s’installe assez facilement en local du moment qu’on a le plugin java. Pour afficher les valeurs de bootstraps, il est nécessaire que figtree les charge comme « labels », ainsi dans la rubrique « node labels » il vous faudra sélectionner « labels » dans le menu déroulant « display ».
Les noms des espèces et les informations des clades correspondantes sont disponibles dans le fichiers /home/formation/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt.
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 4 : Quelles sont les principales différences ? Que pensez-vous de ces arbres ? Leurs supports ? Leurs congruences ? Les grands clades connus sont-ils retrouvés ? Comparer avec les arbres des articles ou revues suivant(e)s :
La revue de Biller et al : Prochlorococcus  : the structure and function of collective diversity
L’article de Kettler et al. : Patterns and implications of gene gain and loss in the evolution of Prochlorococcus et celui de Yan et al. : Genome rearrangement shapes Prochlorococcus ecological adaptation.
N’oubliez pas de commenter les valeurs de support des nœuds. Et pensez à raciner les arbres en utilisant l’outgroup Synechococcus.
</pre>
 
===II)	Super-arbres===
 
Nous allons utiliser les arbres individuels protéiques ainsi que celui reconstruits à partir de la petite sous-unité de l’rRNA.
Sur genologin le script était celui-là :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
#!/bin/bash
iqtree -s ./ssu_renamed_simplified.aln -nt 1 -AIC -bb 1000 -alrt 1000 -redo
</pre>
Et on a tapé :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
module load bioinfo/iqtree-1.6.7
sbatch --mem=20G ssuTree.sh
</pre>
 
Voici l’arbre obtenu : /home/formation/work/ProchlorococcusSynechococcus/phyloG /ssu_renamed_simplified.aln.treefile
 
'''A vous :'''  
Loguez-vous sur genologin (et pas genotoul).
 
Lançons les arbres protéiques. Pour cela faire un seul script que vous appellerez ind_pep_trees.sh et qui écrira toutes les commandes en bouclant sur chaque fichier d’entrée. Le but étant d’obtenir une ligne par commande iqtree sur un fichier d’alignement protéique.  
Les fichiers d’input sont /home/formation/work/ProchlorococcusSynechococcus/phyloG/proteine/*_renamed.fas
 
'''Attention à bien spécifier -nt 1. Si vous souhaitez utiliser plus d’un CPU il faut le réserver avec l’option --cpus-per-task dans la commande sbatch / sarray. Pas besoin d’augmenter la RAM pour les protéines. Pensez aussi à mettre -AIC comme critère de sélection de modèle.'''
 
Regarder le contenu de ind_pep_trees.sh. Est-il correct ?
 
Rajouter la première ligne obligatoire sur genologin avec vi par exemple : 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
#!/bin/bash 
</pre>
 
Pour le lancer en parallèle sur plusieurs nœuds du cluster :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
module load bioinfo/iqtree-1.6.7
sarray ind_pep_trees.sh
</pre>
 
Pour monitorer votre job :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
squeue –u login
</pre>
 
Pour le killer si besoin : 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
scancel jobid
</pre>
 
Quand tous vos jobs sont terminés, vérifiez que les fichiers de sorties ne soient pas vides et que ça s’est bien passé en faisant par exemple :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
tail *_renamed.fas.log
</pre>
 
Vous pouvez aussi regarder les modèles sélectionnés :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
grep 'Best-fit model:' *_renamed.fas.log
</pre>
 
Concaténer tous les arbres (les 31 arbres protéiques et l’arbre ARNr) avec la commande cat.
Nommez le fichier alltrees.tree.
 
'''Test de plusieurs méthodes de super-arbres :'''
 
Commençons par la méthode la plus répandue : le MRP.
 
Pour aller sur un nœud :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
srun --pty bash
</pre>
puis taper :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
module load system/R-3.5.1
R
library(phytools)
trees=read.tree("MonPath/alltrees.tree ")
supertrees<-mrp.supertree(trees,rearrangements="SPR", start="NJ")
</pre>
Vous avez obtenu les super-arbres les plus parcimonieux. Sauvez-les en utilisant la fonction write.tree de R.
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
write.tree(supertrees, file = MonPath/superTrees.tree")
quit()
</pre>
Nous sommes sortis de R.
 
Dans notre cas, nous avons un seul arbre le plus parcimonieux, mais nous aurions pu en obtenir plusieurs.
 
Utiliser iqtree pour obtenir le consensus des 32 arbres avec la règle majoritaire étendue.
 
Pour cela restez sur le nœud et tapez :
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
module load bioinfo/iqtree-1.6.7
iqtree -con -t alltrees.tree -nt 1
</pre>
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 5 : D’après vous que signifie les valeurs aux nœuds sur cet arbre consensus ?
</pre>
 
Lancez aussi l’arbre consensus en réseau : 
<pre style="color:purple;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
iqtree -net -t alltrees.tree -nt 1
</pre>
 
Et visualisez-le avec splitstree en local. 
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 6 : Commentez ce réseau par rapport aux autres arbres obtenus. Qu’en pensez-vous ?
</pre>
 
===III)	Comparaison des arbres===
 
Concaténez maintenant les deux arbres de super-matrice (celui sur les codons et celui sur les protéines) ainsi que les deux super-arbres (le consensus et le MRP).
 
Lancer ensuite le calcul de la distance de Robinson and Foulds sur ce fichier avec iqtree.
 
Attention à faire le module load bioinfo/iqtree-1.6.7
 
Et à rajouter -nt 1 dans les options.
 
<pre style="color:red;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
Question 7 : Commentez les résultats. Quel est l’arbre le plus différent des autres ? Qu’est-ce qui pourrait l’expliquer ?
</pre>
 
Fin du TP du 25 octobre 2018
-->
 
 
 
 
 
 
 
 
 
 
 
 
 
<!--
==old stuff==
 
Lien sur les génomes de ''Prochlorococcus'' du NCBI:
*Répertoire avec les [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Prochlorococcus/NCBI souches] du NCBI.
 
Tables avec des information sur les projets de séquençages:
 
*[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Prochlorococcus/NCBI/Prochlorococcus.csv Prochlorococcus.csv]
 
 
Nous avons édité le fichier pour associer un identifiant de moins de 5 lettres pour chaque génomes (compatibilité avec certains logiciels utilisés).
 
====Autres souches====
Les génomes des souches suivantes sont absentes du NCBI.
*''Wolbachia'' endosymbiont of ''Dirofilaria immitis'' wDi.2.2  [http://dirofilaria.nematod.es/ http://dirofilaria.nematod.es/]
*''Wolbachia'' endosymbiont of ''Litomosoides sigmodontis'' wLs.2.0  [http://litomosoides.nematod.es/ http://litomosoides.nematod.es/]
 
[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/Additional Autres souches] téléchargées.
 
====Groupes externes====
Les souches suivantes peuvent être utilisées pour enraciner un arbre réalisé sur les souches de Wolbachia:
*Ace	Anaplasma centrale str. Israel	—	1	1 206 806	[https://www.ncbi.nlm.nih.gov/nuccore/NC_013532.1?report=fasta NC_013532]
*Aph	Anaplasma phagocytophilum HZ	—	1	1 471 282	[https://www.ncbi.nlm.nih.gov/nuccore/NC_007797.1?report=fasta NC_007797]
*Ech	Ehrlichia chaffeensis str. Arkansas	—	1	1 176 248	[https://www.ncbi.nlm.nih.gov/nuccore/NC_007799.1?report=fasta NC_007799]
*Eru	Ehrlichia ruminantium str. Gardel	—	1	1 499 920	[https://www.ncbi.nlm.nih.gov/nuccore/NC_006831.1?report=fasta NC_006831]
*Nri	Neorickettsia risticii str. Illinois	—	1	879 977	[https://www.ncbi.nlm.nih.gov/nuccore/NC_013009.1?report=fasta NC_013009]
*Nse	Neorickettsia sennetsu str. Miyayama	—	1	859 006	[https://www.ncbi.nlm.nih.gov/nuccore/NC_007798.1?report=fasta NC_007798]
*Ccr	Caulobacter crescentus CB15	—	1	4 016 947	[https://www.ncbi.nlm.nih.gov/nuccore/NC_002696.1?report=fasta NC_002696]
 
[http://saba.ibcg.biotoul.fr/soap-data/Atelier_Phylogenomique/Wolbachia/Outgroup_strains Outgroup_strains]
-->
<!--
===Autres ressources===
====PATRIC====
PATRIC est une ressource importante pour les bactéries patogènes: "
The Pathosystems Resource Integration Center ([https://www.patricbrc.org/ PATRIC]): the bacterial Bioinformatics Resource Center ([https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkw1017 Wattam et al., 2017])".
 
====Joint Genome Institute====
[https://jgi.doe.gov/ 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====
[http://bacteria.ensembl.org/index.html Ensembl Bacteria] is a browser for bacterial and archaeal genomes. 
 
 Search for a genome: Prochlorococcus
 20 entries
 Downloads, Filter Prochlorococcus
 16 entries
-->
<!--
==Genome resources==
The Pathosystems Resource Integration Center ([https://www.patricbrc.org/ PATRIC])is the bacterial Bioinformatics Resource Center ([https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkw1017 Wattam et al., 2017]).
 
 All Data Types: wolbachia
 Genomes: 46
 
[https://www.patricbrc.org/view/GenomeList/?keyword(wolbachia)#view_tab=genomes Genome List View]
 
====Download genomes====
 Select all genomes (in Genome Name)
 Use DWNLD with More Options to download files in FASTA format.
 
Genomes to add:
*''Wolbachia'' endosymbiont of ''Dirofilaria immitis'' wDi.2.2  [http://dirofilaria.nematod.es/ http://dirofilaria.nematod.es/]
*''Wolbachia'' endosymbiont of ''Litomosoides sigmodontis'' wLs.2.0  [http://litomosoides.nematod.es/ http://litomosoides.nematod.es/]
 
Outgroup strains
*Ace	Anaplasma centrale str. Israel	—	1	1 206 806	NC_013532
*Aph	Anaplasma phagocytophilum HZ	—	1	1 471 282	NC_007797
*Ech	Ehrlichia chaffeensis str. Arkansas	—	1	1 176 248	NC_007799
*Eru	Ehrlichia ruminantium str. Gardel	—	1	1 499 920	NC_006831
*Nri	Neorickettsia risticii str. Illinois	—	1	879 977	NC_013009
*Nse	Neorickettsia sennetsu str. Miyayama	—	1	859 006	NC_007798
*Ccr	Caulobacter crescentus CB15	—	1	4 016 947	NC_002696
 
===Joint Genome Institute===
[https://jgi.doe.gov/ 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===
[http://bacteria.ensembl.org/index.html Ensembl Bacteria] is a browser for bacterial and archaeal genomes. 
 
 Search for a genome: wolbachia
 20 entries
 Downloads, Filter wolbachia
 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": [http://www.vicbioinformatics.com/software.prokka.shtml prokka home].
 
Citation: 
Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014 Jul 15;30(14):2068-9. [http://www.ncbi.nlm.nih.gov/pubmed/24642063 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. 
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
Output Files
<pre>
.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
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
more prokka/wDim/wDim.txt
</pre>
<pre>
organism: Genus wolbachia endosymbiont wDim 
contigs: 2
bases: 921012
rRNA: 3
gene: 793
CDS: 755
tRNA: 34
tmRNA: 1
</pre>
 
====''Litomosoides sigmodontis wolbachia''====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
more prokka/wLsi/wLsi.txt
</pre>
<pre>
organism: Genus wolbachia endosymbiont wLsi contigs: 10
bases: 1048936
rRNA: 3
gene: 918
CDS: 880
tRNA: 34
tmRNA: 1
</pre>
 
===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.
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
art prokka/wDim/wDim.gff
</pre>
 
===Comparaison des génomes===
Nous allons utiliser le logiciel act ([ftp://ftp.sanger.ac.uk/pub/resources/software/act/act.pdf 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
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
formatdb -i prokka/wDim/wDim.fna -p F
</pre>
Créer un nouveau répertoire: act
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
blastall -p blastn -i prokka/wLsi/wLsi.fna -d prokka/wDim/wDim.fna -m 8 -o act/wLsi-wDim.blt
</pre>
 
Copier les fichiers sur votre porte de travail en P0 et lancer :
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
act prokka/wLsi/wLsi.gbf act/wLsi-wDim.blt prokka/wDim/wDim.gbf&
</pre>
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:
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/scripts/rename_ncbi_files.pl --prokka_dir prokka --verbose 0
</pre>
 
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. 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
 
==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 ([http://pubmlst.org/wolbachia 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
 
*[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/MLST/coxA.fas coxA.fas] 
*[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/MLST/fbpA.fas fbpA.fas] 
*[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/MLST/ftsZ.fas ftsZ.fas]
*[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/MLST/gatB.fas gatB.fas]
*[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/MLST/hcpA.fas hcpA.fas]
<!--
*[http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/MLST/concat_file.mfa  concat_file.mfa]
-->
<!--
===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.
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/home/formation/public_html/M2_Phylogenomique/scripts/fasta2hmm.pl --mlst_dir MLST --gene_list "gatB coxA hcpA ftsZ fbpA" --verbose 0 --erase 0
</pre>
 
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====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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;
</pre>
Le fichier avec les alignements est en format STK (STOCKHOLM), nous devons le transformer en format FASTA
 head MLST/gatB_wAu.stk
 
====STK à FASTA====
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/bin/esl-reformat --informat stockholm -o MLST/gatB_wAu.fa fasta MLST/gatB_wAu.stk;
</pre>
 
====Automatisation====
Les deux étapes sont automatisées pour tous les gènes et tous les génomes.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
 
Les gènes concaténés sont dans le fichier : 
 MLST/concat_file.fa
 
====Alignement et arbre====
muscle:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/bin/muscle -in MLST/concat_file.fa -out MLST/concat_file.mfa
</pre>
readal (fasta à phylip):
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/bin/readal -in MLST/concat_file.mfa -out MLST/concat_file.phy -phylip;
</pre>
PhyML:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
 
====Visualisation de l'arbre====
Lancer l'interface graphique sur votre poste de travail de la P0.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
seaview MLST/concat_file.phy_phyml_tree.txt&
</pre>
 
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 [http://genoweb.toulouse.inra.fr/~formation/M2_Phylogenomique/data/Wolbachia/NCBI/NCBI_Wolbachia.csv 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:
  [http://onlinelibrary.wiley.com/doi/10.1002/0471250953.bi0612s35/full OrthoMCL Protocols]
 
For details on the orthomcl algorithm, please read the OrthoMCL Algorithm Document: 
  [https://docs.google.com/document/d/1RB-SqCjBmcpNq-YbOYdFxotHGuU7RK_wqxqDAMjyP_w/pub OrthoMCL Algorithm Document]
 
En entrée, le programme prend un ensemble de protéomes.
 
La sortie de orthoMCL est un ensemble de fichiers :
<pre>
  pairs/
     potentialOrthologs.txt
     potentialCoorthologs.txt
     potentialInparalogs.txt
  groups.txt
</pre>
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
<pre>
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
</pre>
 
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:
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
cat wolbachia/blastp/*fasta > wolbachia/blastp/wol_21_blt
</pre>
 
=== 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 [https://docs.google.com/document/d/1RB-SqCjBmcpNq-YbOYdFxotHGuU7RK_wqxqDAMjyP_w/pub 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.
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclBlastParser wolbachia//blastp/wol_21_blt wolbachia/my_orthomcl_dir/compliantFasta >> wolbachia/my_orthomcl_dir/similarSequences.txt
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
head wolbachia/my_orthomcl_dir/similarSequences.txt
</pre>
 
===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:
<pre>
/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
</pre>
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/orthomclSoftware-v2.0.9/bin/orthomclInstallSchema wolbachia/my_orthomcl_dir/orthomcl.config wolbachia/my_orthomcl_dir/install_schema.log
</pre>
====DROP Tables====
OrthoMCL expects a "clean" and empty database after orthoMCLInstallSchema, so in between runs it's easiest to run 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
mysql -h 147.99.108.14 -D orthomcl_yquentin -u ortho_yquentin -p
</pre>
<pre>
mysql>SHOW TABLES;
</pre>
<pre style="color:grey;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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;
</pre>
 
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. 
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclLoadBlast wolbachia/my_orthomcl_dir/orthomcl.config work/wolbachia/my_orthomcl_dir/similarSequences.txt
</pre>
 
=== 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=).
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclPairs wolbachia/my_orthomcl_dir/orthomcl.config wolbachia/my_orthomcl_dir/orthomcl_pairs.log cleanup=no
</pre>
 
=== 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).
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclDumpPairsFiles wolbachia/my_orthomcl_dir/orthomcl.config
</pre>
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
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
</pre>
 
=== Step 12: mcl ===
Input:
  - mclInput file
Output:
  - mclOutput file
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/bin/mcl wolbachia/my_orthomcl_dir/mclInput --abc -I 1.5 -o work/wolbachia/my_orthomcl_dir/mclOutput
</pre>
 
=== Step 13: orthomclMclToGroups ===
Input:
  - mclOutput file
Output:
  - groups.txt
 
Change to my_orthomcl_dir and run:
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclMclToGroups wolb 1000 < mclOutput > groups.txt
</pre>
*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:
 
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/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
</pre>
 wc -l wolbachia/my_orthomcl_dir/singletons.txt
 1639
 
===Other utilities===
Extract protein IDs from an orthomcl groups file.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclExtractProteinIdsFromGroupsFile wolbachia/my_orthomcl_dir/groups.txt 
</pre>
Extract protein ID pairss from an orthomcl groups file.
<pre style="color:blue;white-space: pre-wrap;white-space: -moz-pre-wrap;white-space: -pre-wrap;white-space: -o-pre-wrap">
/usr/local/bioinfo/src/OrthoMCL/current/bin/orthomclExtractProteinPairsFromGroupsFile wolbachia/my_orthomcl_dir/groups.txt</pre>
</pre>
-->
 
<!--
==Perl scripts==
 
PROKKA
* rename_ncbi_files.pl
* prokka_loop.pl
MLST
* fasta2hmm.pl
* hmm2genes.pl
orthMCL
*prep_orthoMCL.pl
 
 
[http://pangenome.de/P_marinus_meta P_marinus_meta]
-->