silico.biotoul.fr
 

M1 MABS BBS Math TD Proba

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Pseudo-comptage)
m (Pseudo-comptage)
(11 intermediate revisions not shown)
Line 1: Line 1:
-
= Recherche de site de fixation au ribosome =
+
= Recherche de sites de fixation au ribosome =
-
'''Rappel :'''
+
'''Rappel :''' Théorème de Bayes
-
<math>P(A/B) = \frac{P(A \cap B)}{P(B)} = \frac{P(B/A) \times P(A)}{P(B)}</math>
+
<math>P(A/B) = \frac{P(A \cap B)}{P(B)} = \frac{P(B/A) \times P(A)}{P(B)}</math> (pour plus de détails : [https://fr.wikipedia.org/wiki/Th%C3%A9or%C3%A8me_de_Bayes wikipédia])
Il s'agit dans la suite de détecter la présence de RBS (Ribosome Binding Site) sur une séquence génomique :
Il s'agit dans la suite de détecter la présence de RBS (Ribosome Binding Site) sur une séquence génomique :
Line 13: Line 13:
A l'aide du site [http://weblogo.threeplusone.com/ WebLogo], établir le Weblogo correspondant aux séquences de ''B. subtilis''. Où se situent les RBS ?
A l'aide du site [http://weblogo.threeplusone.com/ WebLogo], établir le Weblogo correspondant aux séquences de ''B. subtilis''. Où se situent les RBS ?
 +
== Calcul de <math>P(RBS)</math> ==
 +
 +
Pour la probabilité de présence d'un RBS, il est possible de l'approximer par le nombre de gènes (4,177) ''B. subtilis'' divisé par la taille de son génome (4,215,606 bp).
 +
 +
== Calcul de <math>P(sequence)</math> ==
 +
 +
Pour la probabilité d'observer une séquence, on pourra utiliser le produit des fréquences de chaque nucléotide dans le génome [[Media:BsubA.fas.gz]]. Après avoir récupéré et décompressé le fichier, pour lire le fichier FASTA, utilisez la librairie <tt>seqinr</tt> qui est normalement installée :
 +
<source lang='rsplus'>
 +
library(seqinr)
 +
# chargement
 +
genome=summary(read.fasta('BsubA.fas')$BsubA01)
 +
summary(genome)
 +
#            Length Class  Mode 
 +
#length      1      -none- numeric
 +
#composition 4      table  numeric
 +
#GC          1      -none- numeric
 +
genome$length
 +
#[1] 4215606
 +
genome$composition
 +
#
 +
#      a      c      g      t
 +
#1188073  919284  915112 1193137
 +
</source>
 +
 +
Vous pouvez donc calculer <tt>f_nt</tt>, la fréquence de chaque nucléotide dans le génome et écrire une fonction qui calcule la probabilité d'observer une séquence (on suppose l'indépendance à chaque position, ce qui est complètement faux) :
 +
<source lang='rsplus'>
 +
p_seq = function(s, freqs=f_nt ) {
 +
  # TO DO
 +
}
 +
# test : probabilité d'observer la séquence GGAGG
 +
p_seq(c('g','g','a','g','g'))
 +
# 0.0006258065
 +
</source>
== Calcul de <math>P(sequence/RBS)</math> ==
== Calcul de <math>P(sequence/RBS)</math> ==
Line 18: Line 51:
L'[[Media:Maths Proba Alignement RBS.fasta|alignement]] des séquences de RBS de ''B. subtilis'' vous est fournit au format FASTA. Il va vous servir à calculer P(sequence/RBS). Pour cela, vous devrez construire la matrice de fréquence de chaque résidu à chaque position (la fréquence observée sert à estimer la probabilité).
L'[[Media:Maths Proba Alignement RBS.fasta|alignement]] des séquences de RBS de ''B. subtilis'' vous est fournit au format FASTA. Il va vous servir à calculer P(sequence/RBS). Pour cela, vous devrez construire la matrice de fréquence de chaque résidu à chaque position (la fréquence observée sert à estimer la probabilité).
-
Pour lire le fichier FASTA, utilisez la librairie <tt>seqinr</tt> qui est normalement installée :
 
<source lang='rsplus'>
<source lang='rsplus'>
-
library(seqinr)
 
known_rbs = read.fasta('Maths_Proba_Alignement_RBS.fasta')
known_rbs = read.fasta('Maths_Proba_Alignement_RBS.fasta')
Line 59: Line 90:
Vous pouvez à présent définir une fonction qui calcule P(sequence/RBS) : la probabilité d'une séquence de 5 nucléotides sachant qu'elle correspond à un RBS (produit des probabilités des résidus observés à chaque position).
Vous pouvez à présent définir une fonction qui calcule P(sequence/RBS) : la probabilité d'une séquence de 5 nucléotides sachant qu'elle correspond à un RBS (produit des probabilités des résidus observés à chaque position).
-
Rappel: pour définir une fonction :
+
'''Rappel :''' pour définir une fonction :
<source lang='rsplus'>
<source lang='rsplus'>
Line 68: Line 99:
p_seq_rbs(c('g','g','a','g','g'))
p_seq_rbs(c('g','g','a','g','g'))
# 0.1117562
# 0.1117562
-
</source>
 
-
 
-
 
-
== Calcul de <math>P(RBS)</math> ==
 
-
 
-
Pour la probabilité de présence d'un RBS, il est possible de l'approximer par le nombre de gènes (4,177) ''B. subtilis'' divisé par la taille de son génome (4,215,606 bp).
 
-
 
-
== Calcul de <math>P(sequence)</math> ==
 
-
 
-
Pour la probabilité d'observer une séquence, on pourra utiliser le produit des fréquences de chaque nucléotide dans le génome [[Media:BsubA.fas.gz]]. Après avoir récupérer et décompressé le fichier :
 
-
<source lang='rsplus'>
 
-
# chargement
 
-
genome=summary(read.fasta('BsubA.fas')$BsubA01)
 
-
summary(genome)
 
-
#            Length Class  Mode 
 
-
#length      1      -none- numeric
 
-
#composition 4      table  numeric
 
-
#GC          1      -none- numeric
 
-
genome$length
 
-
#[1] 4215606
 
-
genome$composition
 
-
#
 
-
#      a      c      g      t
 
-
#1188073  919284  915112 1193137
 
-
</source>
 
-
 
-
Vous pouvez donc calculer <tt>f_nt</tt>, la fréquence de chaque nucléotide dans le génome et écrire une fonction qui calcule la probabilité d'observer une séquence (on suppose l'indépendance à chaque position, ce qui est complètement faux) :
 
-
<source lang='rsplus'>
 
-
p_seq =
 
-
p_seq = function(s, freqs=f_nt ) {
 
-
  # TO DO
 
-
}
 
-
# test : probabilité d'observer la séquence GGAGG
 
-
p_seq(c('g','g','a','g','g'))
 
-
# 0.0006258065
 
</source>
</source>
Line 109: Line 105:
Vous pouvez maintenant écrire une fonction qui calcule <math>P(RBS/ggagg) = \frac{P(ggagg/RBS) \times P(RBS)}{P(ggagg)}</math> :
Vous pouvez maintenant écrire une fonction qui calcule <math>P(RBS/ggagg) = \frac{P(ggagg/RBS) \times P(RBS)}{P(ggagg)}</math> :
<source lang='rsplus'>
<source lang='rsplus'>
-
p_rbs_seq = function(s) {
+
p_rbs_seq = function(s, frp=f_res_pos, fnt=f_nt) p_seq_rbs(s, frp) * p_RBS / p_seq(s, fnt)
-
  p_seq_rbs(s) * p_RBS / p_seq(s)
+
-
}
+
# test avec GGAGG
# test avec GGAGG
p_rbs_seq(c('g','g','a','g','g'))
p_rbs_seq(c('g','g','a','g','g'))
Line 153: Line 147:
Pour éviter cela, on utilise la méthode de pseudo-comptage qui consiste à attribuer une faible probabilité plutôt que 0. Ceci est obtenu en ajoutant une faible quantité à toutes les évènements observés. La (pseudo-)fréquence est donc <math>freq = (nb_{occurrences} + 1) / (nb_{elements} + 1)</math>
Pour éviter cela, on utilise la méthode de pseudo-comptage qui consiste à attribuer une faible probabilité plutôt que 0. Ceci est obtenu en ajoutant une faible quantité à toutes les évènements observés. La (pseudo-)fréquence est donc <math>freq = (nb_{occurrences} + 1) / (nb_{elements} + 1)</math>
-
En R :
+
 
 +
Refaire donc le calcul pour obtenir le nombre d’occurrences à chaque position +1 :
<source lang='rsplus'>
<source lang='rsplus'>
-
# OCCURENCE A CHAQUE POSITION
+
   [,1] [,2] [,3] [,4] [,5]
-
alphabet = c('a','c','g','t')
+
a  36  12  42  45    7
-
motif_length = length(getSequence(known_rbs[1])[[1]])
+
c    1    1    1    1    1
-
pf_res_pos = matrix(nrow = length(alphabet), ncol = motif_length)
+
g  69  94  63  58  93
-
rownames(pf_res_pos) = alphabet
+
t    2    1    2    4    7
-
for (res in alphabet)
+
</source>
-
for (pos in 1:motif_length) {
+
-
pf_res_pos[res, pos] = length( which( alignment[,pos]==res) ) + 1
+
-
}
+
-
pf_res_pos
+
-
#   [,1] [,2] [,3] [,4] [,5]
+
-
# a  36  12  42  45    7
+
-
# c    1    1    1    1    1
+
-
# g  69  94  63  58  93
+
-
# t    2    1    2    4    7
+
-
 
+
-
# on obtient ensuite les frequences:
+
-
pf_res_pos = pf_res_pos / (nrow(alignment)+1)
+
 +
Et donc les fréquences suivantes :
 +
<source lang='rsplus'>
round( pf_res_pos ,2)
round( pf_res_pos ,2)
-
#   [,1] [,2] [,3] [,4] [,5]
+
   [,1] [,2] [,3] [,4] [,5]
-
# a 0.34 0.11 0.40 0.43 0.07
+
a 0.33 0.11 0.39 0.42 0.06
-
# c 0.01 0.01 0.01 0.01 0.01
+
c 0.01 0.01 0.01 0.01 0.01
-
# g 0.66 0.90 0.60 0.55 0.89
+
g 0.64 0.87 0.58 0.54 0.86
-
# t 0.02 0.01 0.02 0.04 0.07
+
t 0.02 0.01 0.02 0.04 0.06
</source>
</source>
Line 204: Line 189:
[[Image:Maths Proba Predits Positifs.png]]
[[Image:Maths Proba Predits Positifs.png]]
 +
Intervalle de confiance sur la moyenne de la probabilité d'un RBS dans le jeu d'apprentissage :
 +
<source lang='rsplus'>
 +
# INTERVALLE DE CONFIANCE A 95%
 +
m = mean(true_pos_probas) # 0.140906
 +
abline(h=m,  lwd=3)
 +
s = sd(true_pos_probas) # 0.1069476
 +
error = qnorm(0.975)*s/sqrt(nrow(alignment))
 +
error
 +
abline(h=m-error,  lwd=2)
 +
</source>
 +
 +
[[Image:Maths Proba Predits Positifs et intervalle de confiance.png]]
 +
 +
Attention, un tel seuil engendrera de nombre faux négatifs :
 +
<source lang='rsplus'>
 +
> which(true_pos_probas<m-error)
 +
[1]  3  4  5  6  8  12  13  14  17  18  21  24  26  27  28  29  31  32  34  35  36  40  41  44  45  46  47  49  50  53  58  59  60  61  63
 +
[36]  72  75  76  77  78  81  82  84  86  93  94  97 103
 +
</source>
 +
 +
<!--
= Loi hypergéométrique =
= Loi hypergéométrique =
<math>p=\sum^{min(q,t)}_{k=c}\frac{C^k_t\times C^{q-k}_{g-t}}{C^q_g}</math>
<math>p=\sum^{min(q,t)}_{k=c}\frac{C^k_t\times C^{q-k}_{g-t}}{C^q_g}</math>
Line 226: Line 232:
* Quelle est le plus grand nombre <math>x</math> pour lequel vous pouvez calculer <math>x!</math> ?
* Quelle est le plus grand nombre <math>x</math> pour lequel vous pouvez calculer <math>x!</math> ?
 +
-->

Revision as of 09:04, 23 October 2017

Contents

Recherche de sites de fixation au ribosome

Rappel : Théorème de Bayes

P(A/B) = \frac{P(A \cap B)}{P(B)} = \frac{P(B/A) \times P(A)}{P(B)} (pour plus de détails : wikipédia)

Il s'agit dans la suite de détecter la présence de RBS (Ribosome Binding Site) sur une séquence génomique : P(RBS/sequence) = \frac{P(sequence/RBS) \times P(RBS)}{P(sequence)}


Pour se convaincre de la présence d'un motif conservé, les séquences en amont du site d’initiation de la traduction d'une centaine de séquences de Bacillus subtilis sont fournies (Séquences au format fasta). Elles sont alignées sur le codon start.

A l'aide du site WebLogo, établir le Weblogo correspondant aux séquences de B. subtilis. Où se situent les RBS ?

Calcul de P(RBS)

Pour la probabilité de présence d'un RBS, il est possible de l'approximer par le nombre de gènes (4,177) B. subtilis divisé par la taille de son génome (4,215,606 bp).

Calcul de P(sequence)

Pour la probabilité d'observer une séquence, on pourra utiliser le produit des fréquences de chaque nucléotide dans le génome Media:BsubA.fas.gz. Après avoir récupéré et décompressé le fichier, pour lire le fichier FASTA, utilisez la librairie seqinr qui est normalement installée :

library(seqinr)
# chargement
genome=summary(read.fasta('BsubA.fas')$BsubA01)
summary(genome)
#            Length Class  Mode   
#length      1      -none- numeric
#composition 4      table  numeric
#GC          1      -none- numeric
genome$length
#[1] 4215606
genome$composition
#
#      a       c       g       t 
#1188073  919284  915112 1193137

Vous pouvez donc calculer f_nt, la fréquence de chaque nucléotide dans le génome et écrire une fonction qui calcule la probabilité d'observer une séquence (on suppose l'indépendance à chaque position, ce qui est complètement faux) :

p_seq = function(s, freqs=f_nt ) {
   # TO DO
}
# test : probabilité d'observer la séquence GGAGG
p_seq(c('g','g','a','g','g'))
# 0.0006258065

Calcul de P(sequence / RBS)

L'alignement des séquences de RBS de B. subtilis vous est fournit au format FASTA. Il va vous servir à calculer P(sequence/RBS). Pour cela, vous devrez construire la matrice de fréquence de chaque résidu à chaque position (la fréquence observée sert à estimer la probabilité).

known_rbs = read.fasta('Maths_Proba_Alignement_RBS.fasta')
 
# affichage
known_rbs
 
# premiere sequence (objet sequence)
known_rbs[1]
 
# premiere sequence (sequence en nucleotides)
attr(known_rbs[1], 'name') # identifiant
getSequence(known_rbs[1]) # objet liste
getSequence(known_rbs[1])[[1]] # 1ere sequence sous forme vecteur de caracteres

Manipulez l'objet known_rbs afin d'avoir une matrice alignment comme suit :

> head(alignment)
     [,1] [,2] [,3] [,4] [,5]
[1,] "g"  "g"  "a"  "g"  "g" 
[2,] "g"  "g"  "a"  "g"  "g" 
[3,] "a"  "g"  "a"  "g"  "g" 
[4,] "g"  "g"  "g"  "t"  "g" 
[5,] "g"  "a"  "g"  "g"  "t" 
[6,] "a"  "g"  "g"  "t"  "g"

A partir de cette matrice, calculer la matrice f_res_pos qui donne la fréquence d'un résidu à une position. Vous devriez obtenir :

> round(f_res_pos,2)
  [,1] [,2] [,3] [,4] [,5]
a 0.34 0.11 0.39 0.42 0.06
c 0.00 0.00 0.00 0.00 0.00
g 0.65 0.89 0.60 0.55 0.88
t 0.01 0.00 0.01 0.03 0.06

Vous pouvez à présent définir une fonction qui calcule P(sequence/RBS) : la probabilité d'une séquence de 5 nucléotides sachant qu'elle correspond à un RBS (produit des probabilités des résidus observés à chaque position).

Rappel : pour définir une fonction :

p_seq_rbs = function(s, f = f_res_pos) { # s: sequence à utiliser, f: frequence des résidus à chaque position
   # TO DO
}
# test de la fonction: probabilité d'observer la sequence GGAGG sachant que c'est un RBS
p_seq_rbs(c('g','g','a','g','g'))
# 0.1117562

Calcul de P(RBS / sequence)

Vous pouvez maintenant écrire une fonction qui calcule P(RBS/ggagg) = \frac{P(ggagg/RBS) \times P(RBS)}{P(ggagg)} :

p_rbs_seq = function(s, frp=f_res_pos, fnt=f_nt) p_seq_rbs(s, frp) * p_RBS / p_seq(s, fnt)
# test avec GGAGG
p_rbs_seq(c('g','g','a','g','g'))
# 0.1769441

Recherche de RBS sur une séquence

Séquence test

Ecrire une fonction qui calcule la probabilité de présence d'un RBS à chaque position d'une séquence donnée. Appliquée à la séquence test, on obtient le résultat suivant :

# chargement de la séquence test
newseq=as.vector(t(read.table('sequence_test.txt',colClasses='character')[1,]))
newseq
 
# recherche de RBS
probas = search_rbs(newseq)
probas
 
# affichage
plot(1:length(probas), probas, type="b", xaxt='n', xlab='sequence', ylab='P(RBS/sequence)')
axis(1, at=1:length(probas), labels=newseq[1:length(probas)])

Image:Maths Proba Sequence test res.png

Pseudo-comptage

Lors du calcul des fréquences de chaque résidu à chaque position dans l'alignement de départ, on observe que certaine fréquences sont nulles. Cela a pour effet d'attribuer la probabilité 0 à une position dès lors qu'un seul résidu a la malchance de tomber sur une de ces fréquences nulles.

Contenu de f_res_pos :

f_res_pos
         [,1]      [,2]        [,3]       [,4]       [,5]
a 0.336538462 0.1057692 0.394230769 0.42307692 0.05769231
c 0.000000000 0.0000000 0.000000000 0.00000000 0.00000000
g 0.653846154 0.8942308 0.596153846 0.54807692 0.88461538
t 0.009615385 0.0000000 0.009615385 0.02884615 0.05769231

On obtient donc une probabilité nulle, si la séquence contient un C (quelle que soit la position) ou un T en deuxième position.

Pour éviter cela, on utilise la méthode de pseudo-comptage qui consiste à attribuer une faible probabilité plutôt que 0. Ceci est obtenu en ajoutant une faible quantité à toutes les évènements observés. La (pseudo-)fréquence est donc freq = (nboccurrences + 1) / (nbelements + 1)

Refaire donc le calcul pour obtenir le nombre d’occurrences à chaque position +1 :

  [,1] [,2] [,3] [,4] [,5]
a   36   12   42   45    7
c    1    1    1    1    1
g   69   94   63   58   93
t    2    1    2    4    7

Et donc les fréquences suivantes :

round( pf_res_pos ,2)
  [,1] [,2] [,3] [,4] [,5]
a 0.33 0.11 0.39 0.42 0.06
c 0.01 0.01 0.01 0.01 0.01
g 0.64 0.87 0.58 0.54 0.86
t 0.02 0.01 0.02 0.04 0.06

Seuil pour la prédiction

A l'aide de la fonction apply, calculer les probabilités obtenues sur le jeu de séquences d'apprentissage (l'alignement). Vous devriez obtenir ceci :

> summary(true_pos_probas)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001209 0.043780 0.162500 0.140900 0.189800 0.355000

Soit un seuil de 0.001209 pour détecter tous les RBS avérés. Malheureusement, ce seuil relativement bas induit des faux positifs sur la séquence tests :

> threshold = min(true_pos_probas)
> which(probas2 >= threshold)
[1] 16 18 19 20 21 22
 
plot(1:length(probas2), probas2, type="l", xaxt='n', xlab='sequence', ylab='P(RBS/sequence)')
axis(1, at=1:length(probas2), labels=newseq[1:length(probas2)])
positives = which(probas2 >= threshold)
points(positives, probas2[positives] , pch=16, col='green')

Image:Maths Proba Predits Positifs.png

Intervalle de confiance sur la moyenne de la probabilité d'un RBS dans le jeu d'apprentissage :

# INTERVALLE DE CONFIANCE A 95%
m = mean(true_pos_probas) # 0.140906
abline(h=m,  lwd=3)
s = sd(true_pos_probas) # 0.1069476
error = qnorm(0.975)*s/sqrt(nrow(alignment))
error
abline(h=m-error,  lwd=2)

Image:Maths Proba Predits Positifs et intervalle de confiance.png

Attention, un tel seuil engendrera de nombre faux négatifs :

> which(true_pos_probas<m-error)
 [1]   3   4   5   6   8  12  13  14  17  18  21  24  26  27  28  29  31  32  34  35  36  40  41  44  45  46  47  49  50  53  58  59  60  61  63
[36]  72  75  76  77  78  81  82  84  86  93  94  97 103