IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
TP : Analyse de séquence avec R
Objectifs :
L’objectif de ce TP est :
présentera les différents types d'analyse de séquence qui peuvent être effectués avec R et les
bibliothèques et packages R correspondants.
I. Récupérer une séquence
Pour commencer la tâche d'analyse de séquence, la première chose dont nous avons besoin est une
séquence d'ADN ou de protéine.
Ces données de séquence peuvent être récupérées soit en visitant la page d'hébergement de la
base de données via un navigateur, soit en accédant au formulaire du magasin de données dans R
via les commandes/fonctions correspondantes. Dans ce TP, nous introduisons quelques méthodes
pour récupérer une séquence.
La récupération de séquences à partir de bases de données nécessite les éléments suivants :
Accès à la base de données via Internet
Un ID de séquence ou un mot-clé pour la recherche
Les bibliothèques biomaRt et seqinr installées et chargées dans la session R
Pour récupérer une séquence de données via R, procédez comme suit :
1. Obtenez une séquence des données. Le moyen le plus simple d'obtenir une séquence
d'intérêt est de visiter la page Web de la base de données et de la télécharger manuellement.
2. ou de l'obtenir via le package seqinr. Pour utiliser le package seqinr, chargez d'abord la
bibliothèque avec la commande suivante :
> [Link]("seqinr")
> library(seqinr)
NB. Utilisez la fonction lseqinr()[1:9] pour lister les 9 premières fonctions de ce package.
3. Choisissez la banque de données dont vous avez besoin pour récupérer la séquence. Tout
d'abord, vous pouvez vérifier les banques de données disponibles, puis en choisir une. Par
exemple, dans notre cas, nous choisissons genbank comme suit :
> choosebank()
> choosebank("genbank")
4. La banque de données sélectionnée est interrogée pour l'ensemble des termes de recherche
d'entrée (K pour mot clé et SP pour espèce) et toutes les séquences disponibles pour les
attributs de recherche sont affectées à l'objet défini qui a été nommé et renvoyé (dans notre
cas, BRCA1). Une fois la banque de données ouverte, interrogez la banque avec la recherche
qui vous intéresse à l'aide de la commande de requête suivante (notez que cette étape prend
du temps):
> query("BRCA1", "SP=Homo sapiens AND K=BRCA1")
En dehors de ces listes, nous pouvons filtrer les résultats non pertinents en examinant les
attributs de la liste, tels que le nom de la séquence. Dans notre code présenté, nous avons
utilisé la position de la séquence pertinente dans la liste. La commande getsequence récupère
ensuite la séquence à partir de toutes les informations disponibles lors de la récupération.
5. Examinez la valeur renvoyée par la commande de requête en examinant les différents
[Link] Page 1/8 IABM S3
IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
composants de l'objet comme suit :
> attributes(BRCA1)
6. Pour vérifier l'ensemble de toutes les séquences récupérées, saisissez la commande suivante
et obtenez le nom, la longueur et d'autres attributs pour chaque séquence disponible dans la
requête :
> BRCA1$req
7. Pour récupérer une séquence spécifique par son numéro d'accession, utilisez l'attribut AC
dans la commande de requête comme suit :
> query ("BRCA1", "SP=Homo sapiens ET AC=U61268")
8. Pour extraire une séquence spécifique de l'objet de requête, utilisez la commande
getsequence suivante (dans notre cas, nous l'utilisons pour la première séquence de l'objet
BRCA) :
> myseq <- getSequence(BRCA1$req[[1]])
Extraire l'attribut name des séquences comme suit :
> mynames <- getName(BRCA1)
9. Examinez l'objet qui a été créé avec la commande getsequence suivante :
> myseq
10. Vous pouvez récupérer l'annotation (les autres attributs d'une séquence) avec la commande
suivante :
> annots <- getAnnot(BRCA1$req[[1]])
11. Il est toujours conseillé de fermer la banque de données après la requête afin d'éviter
plusieurs banques ouvertes. Faites cela avec la fonction closebank suivante :
> closebank()
12. Lecture et écriture du fichier FASTA
Essayez d'écrire les données de séquence récupérées contenant des informations dans
différents champs sous forme de séquences multiples (ceci n'est applicable que si plus d'une
séquence a été récupérée, sinon vous recevrez une seule séquence) dans un fichier FASTA
nommé [Link] comme suit :
> [Link](myseqs, mynames, [Link] = "[Link]")
Lisez le fichier FASTA à partir de son emplacement (fournissez le chemin complet s'il ne se trouve
pas dans le répertoire de travail actuel) comme suit :
> masequence <- [Link](file = "[Link]")
II. Obtenir le détail d'une composition de séquence
Une fois que nous avons récupéré une séquence, nous devons en savoir plus à son sujet ; par
exemple, nous devons connaître la fréquence des nucléotides ou des acides aminés, et les bases
nucléotidiques de la guanine et de la cytosine (teneur en GC). Pour illustrer, si nous avons des bases
répétitives, disons AAAAAA, la séquence n'est pas très intéressante en raison du faible contenu en
informations. Le contenu d'une séquence aide également à déterminer certaines propriétés de la
molécule entière, par exemple, l'hydrophobicité de la basicité acide dans les protéines en fonction
des acides aminés présents dans la séquence.
[Link] Page 2/8 IABM S3
IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
Un aspect intéressant est la teneur en GC dans un nucléotide. Il fait référence à la fraction de
Guanine (G) et de Cytosine (C) dans la séquence, et certains génomes, notamment chez les
bactéries, montrent une différence significative à cette échelle et des variations en termes de
régions génomiques. Pour illustrer, certaines Actinobactéries peuvent avoir plus de 70 % de GC,
tandis que certaines Protéobactéries peuvent avoir moins de 20 % de GC. De plus, la teneur en GC
est également utilisée pour prédire la température d'hybridation de la séquence lors des
expériences de PCR. Ces aspects rendent l'analyse du contenu d'une séquence importante. Dans
cette recette, nous expliquons comment extraire ces informations des données de séquence.
Pour calculer certains détails de la composition de la séquence, procédez comme suit :
1. Dans notre cas, nous allons chercher deux séquences d'ARN polymérase (sous-unités bêta) de
GenBank (notez que les séquences appartiennent à deux espèces différentes, Mycobacterium
tuberculosis (Actinobacterium) et Escherichia coli (Proteobacterium) on utilisant
choosebank("genbank") et extraire les séquences comme suit :
# Gets the [Link] sequences
> query(listname = "actino", query="SP=Mycobacterium tuberculosis AND K=rpoB")
# Gets the [Link] sequences
> query(listname = "proteo", query="SP=Escherichia coli AND K=rpoB")
2. Observer le résumé des séquences, telles que les noms et la longueur, à l'aide des
commandes suivantes :
> actino$req # 694 Séquences au 17 avril 2014
> proteo$req # 171 Séquences au 17 avril 2014
3. Parmi ces séquences, choisissez une séquence de chaque espèce (M. tuberculosis et E. coli)
qui sera comparée. Dans notre cas, nous allons avec JX303316 dans le premier et
[Link] dans la dernière espèce (les numéros d'index de séquence 644 et 1,
respectivement, dans les objets qwa correspondants) comme suit :
> myActino <- getSequence(actino$req[[644]])
> myProteo <- getSequence(proteo$req[[1]])
NB. Maintenant fermez la banque en utilisant closebanc()
4. Pour calculer le nombre de chaque base dans la séquence, utilisez la fonction de table
générique comme suit :
> table(myActino)
> table(myProteo)
5. Pour obtenir la fraction pour chaque base, faites-le sur une ligne en divisant le résultat de la
fonction de table par la longueur de la séquence (multipliez-le par 100 pour obtenir le résultat
en pourcentage) comme indiqué dans la commande suivante :
> table(myProteo)/length(myProteo)
6. Pour calculer les fréquences individuelles des nucléotides, utilisez la fonction de table
suivante de la même manière que précédemment :
> table(myseqf)
7. Passant aux aspects intéressants du contenu GC des séquences, utilisez simplement la
commande GC du package seqinr pour les deux séquences qui ont été récupérées comme
suit :
> GC(myActino)
[Link] Page 3/8 IABM S3
IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
> GC(myProteo)
8. Créez un graphique à barres simple tout en comparant différentes séquences dans le contenu
GC avec une commande sur une ligne pour visualiser les résultats, comme indiqué dans la
commande suivante :
> barplot(c(Actinobacteria= GC(myActino), Proteobacteria= GC(myProteo)), main=
"Contenu en GC dans différentes bactéries")= GC(myProteo), main="Contenu en GC
dans différentes bactéries")
9. Pour connaître la fréquence de chaque paire possible de nucléotides dans la séquence, utilisez
la fonction count, comme indiqué dans l'exemple suivant, pour chaque paire de caractères
(vous pouvez également le faire pour les triplets et ainsi de suite en choisissant la bonne
valeur pour l’argument de taille de mots):
> seqinr ::count(myActino,wordsize=2)
10. Si vous décidez de copier manuellement la séquence depuis la page Web de la base de
données, elle sera considérée comme un objet personnage ; il doit être converti en un
vecteur de caractères comme suit :
> myseq2 <- "AAAATGCAGTAACCCATGCCAAAATGCAGTAA"
> myseq2 <- strsplit(myseq, "")
> myseq2 <- unlist(myseq)
11. Nous pouvons le faire manuellement en utilisant les valeurs de C (index 2) et G (index 3)
comme suit :
> GC(myseq2)
> myGC <- sum(table(myseq)[2], table(myseq)[3])/sum(table(myseq))
> myGC
Pour rechercher un motif dans une séquence, procédez comme suit :
1. La recherche des codons de départ et d'arrêt dans R peut être considérée comme un simple
problème de correspondance de motifs et peut être résolue avec la bibliothèque Biostrings.
Chargez cette bibliothèque avec la fonction suivante :
> library(Biostrings)
2. Ensuite, créez la séquence à analyser comme suit (dans notre cas, nous utilisons une
séquence aléatoire, mais nous pouvons également utiliser les séquences extraites de
GenBank) :
> monnucléotide <- DNAString("aacataatgcagtagaacccatgagccc")
3. Recherchez un motif qui vous intéresse, tel qu'un codon d'initiation ATG, comme illustré dans
l'exemple suivant :
> matchPattern(DNAString("ATG"), monnucléotide)
4. De même, recherchez le modèle pour les codons d'arrêt, tels que TAA ou d'autres codons
d'arrêt, comme indiqué dans l'exemple suivant :
> matchPattern("TAA", monnucléotide)
III. Alignement de séquences par paires
Comme mentionné dans la section d'introduction de ce chapitre, lors de l'analyse de séquences de
[Link] Page 4/8 IABM S3
IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
gènes ou de protéines, nous avons souvent besoin de les comparer pour connaître leurs similitudes
et leurs différences. Cela sert des objectifs à partir de diverses perspectives, comme pour les études
évolutives et pour comprendre la structure et la fonction d'une nouvelle séquence en le comparant
à la séquence connue. Pour savoir si les deux séquences gène/protéine que nous étudions sont
similaires ou différentes au niveau quantitatif, nous mesurons leurs similarités.
L'alignement de séquence nécessite les éléments suivants :
La bibliothèque Biostrings installée et chargée dans la session R
Séquences qui nous intéressent
La bibliothèque Biostrings peut être installée et chargée comme les autres bibliothèques
Bioconductor. Nous aurons besoin que les séquences soient alignées.
L'alignement calcule également un score d'alignement, pour lequel nous devons définir le système
de notation qui dépend du type de séquence (nucléotide ou protéine). Cela sera discuté dans les
prochaines parties de TP.
Pour effectuer un alignement de séquences par paires, procédez comme suit :
1. Commencez par installer et charger la bibliothèque Biostrings en tapant les commandes
suivantes :
> if (!requireNamespace("BiocManager", quietly = TRUE))
[Link]("BiocManager")
BiocManager::install("Biostrings")
> library(Biostrings)
2. Lisez les deux séquences. Dans l'exemple suivant, nous avons montré des séquences saisies
manuellement pour plus de commodité, mais la même méthode peut également être
utilisée pour d'autres types de séquences :
> sequence1 <- "GAATTCGGCTA"
> sequence2 <- "GATTACCTA"
3. Attribuez des pénalités pour la non-concordance et l'écart afin d'obtenir une matrice de
notation pour les nucléotides comme suit :
> myScoringMat <- nucleotideSubstitutionMatrix(match = 1, mismatch = -1, baseOnly =
TRUE)
> myScoringMat
4. Attribuez des pénalités d'écart pour les alignements comme suit :
> gapOpen <- 2
> gapExtend <- 1
5. Exécutez la fonction pairwiseAlignment suivante pour effectuer un alignement global des
séquences :
> myAlignment <- pairwiseAlignment(sequence1, sequence2, substitutionMatrix =
myScoringMat, gapOpening = gapOpen, gapExtension = gapExtend, type="global",
scoreOnly = FALSE)
> show(myAlignment)
6. Dans le cas des séquences protéiques, nous avons des matrices de notation appelées
matrices de substitution. Pour les voir, utilisez la fonction suivante :
> data(package="Biostrings")
[Link] Page 5/8 IABM S3
IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
7. Attribuez/sélectionnez l'une de ces matrices de substitution pour aligner les séquences
protéiques comme suit :
> data (BLOSUM62)
> subMat <- "BLOSUM62"
8. Maintenant, effectuez l'alignement des séquences protéiques avec la matrice de
substitution sélectionnée comme suit :
> sequence3 <- "PAWHEAE"
> sequence4 <- "HEAGAWGHE"
> myAlignProt <- pairwiseAlignment(sequence3, sequence4, substitutionMatrix =
subMat, gapOpening = gapOpen, gapExtension =gapExtend, type="global", scoreOnly =
FALSE)
> show(myAlignProt)
9. Vous pouvez faire un alignement local en changeant le type en mettant "local"
IV. Alignement de séquences multiples
Il y a des occasions dans l'analyse de séquence où nous devons comparer plusieurs séquences les
unes par rapport aux autres. Un exemple général d'un tel cas est l'analyse phylogénétique, qui sera
discutée dans la recette suivante. Une telle analyse révèle à quel point les séquences sont similaires
ou dissemblables les unes par rapport aux autres. Dans cette recette, nous expliquons comment
faire un alignement de séquences multiples (MSA) en utilisant R. Le package msa fournit une
interface R unifiée aux algorithmes d'alignement de séquences multiples ClustalW, ClustalOmega et
MUSCLE. Le package ne nécessite aucun package logiciel supplémentaire et fonctionne sur toutes
les principales plates-formes pour MSA.
L’alignement de séquences multiples nécessite les éléments suivants :
Une nouvelle bibliothèque nommé msa doit être installé et chargé dans la session R. Les
algorithmes ClustalW, ClustalOmega et MUSCLE sont très bien établis et fréquemment
utilisés par la bioinformatique.
Pour les séquences qui nous intéressent, nous utiliserons les séquences protéiques du
cytochrome oxydase (COX-2) de cinq espèces différentes.
Suivez les étapes suivantes pour effectuer l'alignement de séquences multiples :
1. Commencez par installer et charger les packages requis en saisissant les commandes
suivantes :
> if (!require("BiocManager", quietly = TRUE))
[Link]("BiocManager")
BiocManager::install("msa")
> library (msa)
2. Ensuite, obtenez toutes les séquences à aligner. Ils peuvent être lus directement à partir
d'un fichier FASTA comme suit en utilisant la fonction readAAStringSet de Biostrings library:
> myseq3 <- readAAStringSet("[Link]")
Examinez les séquences du fichier qui a été lu, comme suit :
[Link] Page 6/8 IABM S3
IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
> myseq3
3. Maintenant, vous devez récupérer la séquence dans les référentiels et les écrire dans un
fichier FASTA avant de les utiliser comme expliqué dans les recettes précédentes.
4. Maintenant, utilisez la fonction msa pour effectuer l'alignement avec la commande
suivante :
> myAlignment<- msa(myseq3)
> myAlignment
5. Pour imprimer l'alignement à l'écran, utilisez la fonction suivante :
> print(myAlignment, show = "complete")
6. Enregistrez l'alignement se forma pdf, utilisez la fonction suivante :
> msaPrettyPrint(myAlignment, file=tmpFile, output="pdf",
showNames="none", showNumbering="none", showLogo="none",
verbose=FALSE, askForOverwrite=FALSE)
V. Analyse phylogénétique et arborescence
L'analyse phylogénétique consiste à trouver la relation évolutive entre les espèces (organismes),
dans notre cas, sur la base des données de séquence. Une fois que nous avons un ensemble de
séquences provenant de différentes sources, il devient très intéressant de comprendre à quel point
elles sont proches ou éloignées en termes d'évolution moléculaire. À la suite des mutations au
cours de l'évolution, des différences émergent au niveau de la séquence. Ces différences peuvent
être représentées en termes de mesures de distance (consultez la section Voir aussi à la fin de cette
recette). Ces mesures peuvent ensuite être utilisées pour estimer les relations évolutives entre les
espèces, souvent représentées sous forme d'arbres phylogénétiques. Cette relation est très
souvent décrite en termes d'arbre phylogénétique. Jusqu'à présent, nous avons examiné les
différents aspects de la récupération, de l'alignement et de l'analyse des séquences. Cette section a
pour but d'entraîner la recette pour effectuer une analyse phylogénétique sur les données de
séquence.
Les conditions préalables sont les suivantes :
Une session R avec le package ape installé et chargé
Séquences à analyser (au moins leurs identifiants)
Pour effectuer l'analyse phylogénétique sur les séquences de votre choix, suivez les étapes
suivantes :
1. Vous pouvez créer un arbre phylogénétique après un alignement multiple, par exemple, en
utilisant l'alignement multiple effectué précédemment (myAlignment), nous le convertirons
au format seqinr en utilisant msaconvert comme suit :
> alignment <- msaConvert(myAlignment, "seqinr::alignment")
2. Calculez la matrice de distance pour les séquences avec la fonction [Link] comme
suivante :
> distMatrix <- [Link](alignment, "similarity")
3. Regrouper la matrice de distance à l'aide de hclust
> clustering <- hclust(distMatrix)
4. Tracer les distances
[Link] Page 7/8 IABM S3
IABM Bioinformatique / 2021-2022 FMPM MARRAKECH
> plot(clustering)
5. Et pour faire des arbres plus pratiques avec des colores et ainsi de suite nous utiliserons la
bibliothèque ape, suivez les étapes suivantes:
6. Installez et chargez le package ape en saisissant la commande suivante :
> if (!requireNamespace("BiocManager", quietly = TRUE))
[Link]("BiocManager")
BiocManager::install("ape")
> library(ape)
7. Définissez les séquences qui vous intéressent en tant qu'ensemble en termes d'identifiants
comme suit :
> myset <- c("U15717", "U15718", "U15719", "U15720", "U15721",
"U15722", "U15723", "U15724")
8. Récupérez les séquences qui vous intéressent, comme indiqué dans l'exemple suivant :
> myseqs <- [Link](myset)
9. Calculez la matrice de distance pour les séquences avec la fonction [Link] suivante :
> mydist <- [Link](myseqs)
> mydist
10. Pour obtenir l'objet phylo pour les arbres phylogénétiques, utilisez la fonction triangMtd
suivante qui utilise la méthode des triangles pour reconstruire l'arbre :
> myphylo <- triangMtd(mydist)
11. A ce stade, nous pouvons faire un bootstrap sur l'objet phylo, si besoin, en utilisant le boot.
fonction phylo. Pour des informations détaillées, tapez ?[Link] dans la session R.
12. Créez les différents types d'arbres phylogénétiques pour votre analyse comme suit :
> plot(myphylo, type="phylogram", col="red", cex=1, lwd=1,main="(A) Phylogramme")
> plot(myphylo, type="cladogram", col="red", cex=1, lwd=1, main="(B) Cladogramme")
[Link] Page 8/8 IABM S3