Manipulation des alignements

Multiple Sequence Alignment :

import : from Bio import AlignIO
lecture d'alignement :
- pour un seul alignement : align = AlignIO.read('align.txt', 'clustal')
- pour plusieurs alignements : for align in AlignIO.read('align.txt', 'clustal')
- retourne un objet de type Bio.Align.MultipleSeqAlignment
ecriture d'alignements : AlignIO.write(align, 'myFile.txt', 'clustal')
propriétés et fonctions :
- len(align) : le nombre de séquences alignées.
- align.get_alignment_length() : la taille de l'alignement.
- align.format('clustal') : alignement formatté.
- align[0].id : l'identifiant de la première séquence alignée.
- align[0:2] : les 2 premières séquences alignées.
- align[:,0:5] : les 5 premières colonnes de l'alignement, sous forme d'objet MultipleSeqAlignment.
- align[:, 0] : la première colonne sous forme d'un chaîne de caractères avec la lettre de la première séquence en première position, la lettre de la deuxième séquence en deuxième position, etc ... (chaîne de longueur le nombre de séquences).
- align[2].seq[6] : la 7ème base de la 3ème séquence alignée.
- on peut boucler sur les alignements :
```
for alignedSeq in align:
    print(alignedSeq.id)
      
```
- lancement d'un programme d'alignement multiple :
  - pour avoir la liste, faire : from Bio.Align.Applications import ClustalwCommandline; help(ClustalwCommandline)
  - exemple de lancement :
```
com = ClustalwCommandline('clustalw', infile='myfile.fas')
myStdout, myStderr = com()
	
```
    si statut de retour non 0, une exception est levée. Le stdout se trouve dans la variable myStdout
  - autres exemples classiques: MuscleCommandline, TCoffeeCommandline
- utilisation de water (Smith-Waterman, local) et needle (Needleman-Wunsch, global) de EMBOSS :
  - Alignement avec needle : from Bio.Emboss.Applications import NeedleCommandline; com = NeedleCommandline(asequence='seq1.fas', bsequence='seq2.fas', gapopen=10, gapextend=0.5, outfile='needle.txt'); myStdout, myStderr = com(); align = AlignIO.read('needle.txt', 'emboss')
  - Alignement avec water : from Bio.Emboss.Applications import WaterCommandline; com = WaterCommandline(asequence='seq1.fas', bsequence='seq2.fas', gapopen=10, gapextend=0.5, outfile='needle.txt'); myStdout, myStderr = com(); align = AlignIO.read('needle.txt', 'emboss')
il semble que les positions ne soient pas parsées.

Blast (en local) :

biopython offre un lanceur et parseur pour le blast du NCBI (pas celui de Washington University). Attention, ce n'est pas pour la vieille version blastall (écrite en C), mais la nouvelle version (écrite en C++ en 2009)
from Bio.Blast.Applications import NcbiblastnCommandline pour blastn (pour blastx par exemple, c'est NcbiblastxCommandline)
com = NcbiblastnCommandline(query = 'seqFile.fas', db = 'blastDatabase', evalue = 0.01, outfmt = 5, out = 'blast2.xml') : définit la recherche sur la base blast blastDatabase (fichiers blastDatabase.*), avec un format de sortie en xml, mais ne la lance pas !
si on veut utiliser une seule séquence plutôt qu'une base blast, utiliser subject = 'seqFile2.fas' au lieu de db = 'blastDatabase'.
pour la lancer, il faut faire ensuite : myStdout, myStderr = com()
on peut mettre outfmt = 0 pour avoir un format txt traditionnel, ou outfmt = 6 pour le format tableau.
si l'outil blastn (par exemple) n'est pas dans le path, utiliser en plus l'argument cmd = '/myPath/blastn'
2 types de parseurs de blast : parseur du format xml, largement recommandé car beaucoup plus stable, et parseur plain text (mais format plain text évolue souvent selon les versions de blast)

Parsing de blast :

import : from Bio.Blast import NCBIXML (il faut avoir lancé le blast avec l'option outfmt=5 (pour avoir le xml).
si une seule query (et donc un seul résultat blast) : fh = open('blastOut.xml'); blastRecord = NCBIXML.read(fh) : résultat de type Bio.Blast.Record.Blast
si plusieurs queries : blastRecords = NCBIXML.parse(fh) : renvoie un itérateur.

lecture du contenu de blastRecord :

for alignment in blastRecord.alignments:
    for hsp in alignment.hsps:
        print(hsp.expect)
        print(hsp.query)

blastRecord.query : la query avec sa description, comme dans le fichier blast
blastRecord.query_length : la longueur de la query
blastRecord a un champ description contenant la liste de toutes les descriptions (avec les champs title, score, e, num_alignments) : blastRecord.description[0].title
blastRecord a un champ alignments contenant la liste des tous les alignements avec les champs title, length, hsps (liste des hsps) : blastRecord.alignments[0].title
blastRecord.alignments[0].hit_def : la description de la première séquence matchée (identifiant est au début)
blastRecord.alignments[0].length : la longueur de la première séquence matchée.
blastRecord.alignments[0].hsps[0] est de type Bio.Blast.Record.HSP et contient les champs :
- align_length
- bits, expect, score
- frame : pour du blastn, (1, 1) si brin +, (1, -1) si brin - (car c'est la subject qui est reverse-complémentée si nécessaire). Pour du blastx, le frame est (n, 0) avec n dans {-1, -2, -3, 1, 2, 3} et les frames négatifs sont à partir de la fin (on reverse complément la séquence et on regarde le frame : s'il faut n, alors le frame final sera -n). Avec blastp, le frame renvoyé est (0, 0)
- gaps : nombre de gaps
- identities : nombre de lettres identiques
- match : la ligne des | d'homologie
- query, sbjct : les séquences gappées pour la query et le subject
- query_start, query_end, sbjct_start, sbjct_end : les positions, mais ATTENTION : si l'alignement est sur le brin moins, hsp.sbjct_end < hsp.sbjct_start, c'est à dire que hsp.sbjct_start est toujours en face de hsp.query_start, quelque soit le brin (différent de bioperl !). Ce n'est pas vrai pour du blastx pour lequel on a toujours start < end pour la query et le subject.
sur un hsp, quand on fait str(hsp), on a le HSP au format qu'on veut, mais tronqué au milieu si longueur de plus de 50 !

Blast via le web :

from Bio.Blast import NCBIWWW; result = CBIWWW.qblast('blastn', 'nt', seqRecord.format('fasta')) : par défault, le résultat est sous format xml.
le type retourné est cStringIO.StringI.
on peut faire une seule fois result.read() ou alors aussi result.getvalue()

programmer en python, tutoriel python, graphes en python, Aymeric Duclert