Compte rendu stage

semaine 1 Installation + Bibliographie

Markdown is an easy to read and write text format:

semaine 2 Bibliographie + Implémentation du Peeling Recursif.

Algorithme récursif

Premiere implémentation d’un calcul scientifique en python: ##### Fonction Anterior

\[ a_i(u_i) =\sum_{u_m}\left\{a_m(u_m)g(y_m|u_m)\prod_{\stackrel{j\in S_m}{j\ne f}}p_{mj}(u_m) \times \sum_{u_f}\left\{a_f(u_f)g(y_f|u_f)\prod_{\stackrel{j\in S_f}{j\ne m}}p_{fj}(u_f)\times tr(u_i|um,uf)) \times \prod_{\stackrel{j\in C_mf}{j\ne i}}\Bigg[\sum_{u_j}tr(u_j|u_m ,u_f)g(y_i|u_i) \prod_{k \in S_i}p_{kj}(u_k)\Bigg]\right\}\right\} \]

def AnteriorLoop(self,indi:PedNode,ui) -> float:
        #Stop condition
        if indi in self.d.ped.founders(): #indi is a founder -> no anterior can be calculate
            return self.alleleFreq

        (indf, indm) = (
            self.d.ped.get_father(indi.indiv),
            self.d.ped.get_mother(indi.indiv),
        )
        ## Méthode avec boucle >>> Horrible...  objectif faire la meme chose avec des tenseurs
        res=0
        for um in range(3):
            res += Anterior(indm,ui)*Penetrance(indm,ui)*PosteriorPi(indm,indf)
            sum_uf=0
            for uf in range(3):
                sum_uf += Anterior(indm,ui)*Penetrance(indm,ui)*PosteriorPi(indm,indf)*transmission(ui,um,uf)
                res_sib = 1
                for j in fullsib(indi):
                    sum_uj=0
                    for uj in range(3):
                        sum_uj+= transmission(uj,um,uf)*Penetrance(indm,uj)*PosteriorPi(j)
                    res_sib*sum_uj
                sum_uf*res_sib
            res*sum_uf
        return res

Meme fonction avec numpy et utilisation de tenseur:

def Anterior(self,indiv:PedNode) -> np.ndarray:
        ## Stop Condition: if indiv is a founder return AlleleFreq
        if indiv in self.d.ped.founders:
            return self.d.alleleFreq
        mother=indiv.mother
        father=indiv.father
        sibilist=self.getfullsib(indiv)

        ## if indiv got no sibling then Pi_probajoin = 1 
        if len(sibilist) ==0:
            sumF = np.sum(self.Peeling(father,mother.indiv)*self.d.transmission,axis=2)
        else:
            sumF = np.sum(self.Peeling(father,mother.indiv)*self.d.transmission*self.Pi_Probajoin(indiv,sibilist),axis=2) 
        
        sumM = np.sum(self.Peeling(mother,father.indiv)*sumF,axis=1)
        return sumM #shape =(3,) pour chaqu'un des génotype de indiv. a vérifier

Penetrance + Matrice d’erreur

La matrice d’erreur est la matrice qui traduit les erreurs de génotypage pour chaque individu. Elle est construite de manière a ce qu’elle associe une probabilitée de chaque génotype par individu par locus.

Ainsi si un individu est génotypé AA (homozygote référent) alors on renvoi \((1-ε)^2\) comme Probabilitée que se soit effectivement un individu AA mais donne aussi les probabilitées pour aa: \(ε^2\) et pour Aa: \((1-ε)\times ε\). Cela correspond donc a la probabilitée conditionnelle d’observer le vrais génotype sachant le génotypage observé. ie. \(P(y_i|u_i)\). Ce qui correspond exactement a notre fonction de Penetrance

Erreur

Stockage des donnée

self.Pi_post ={}
        # J'aurais pu stocker les donnée pour chaque couple p_ij mais compliquer a générer des clé pour le dictionnaire (ou alors utiliser une matrice symétrique)

En réalité c’est une fausse bonne idée de faire ca: par fois je vais calculer le posterior en ne prenant pas compte un des individus. donc cela peut me faussé les résultats si je l’enregistre le produit de Posterieur de tous les couples de i mais quen suite je veux le récupérer en excluant un couple. on pourrait diviser le résultat par le Posterior(i,excluded)

Solution: faire un dictionaire pour les Posterieurs et mettre en clé un frozenset {i,j} (les sets ne sont pas ordonner ie. {i,j} = {j,i}) ou un tuple trié

Test de l’algorythme récursif sur tous les locus du jeux de deonnée de Fernando (25 individus, 10 loci)

En gardant l’information totale des génotype (sans aucun NA). On obtiens un résultat inattendu a savoir 249 bonnes inférences et 1 erreur. Il sagit du troisieme locus du deuxieme individu. En effet il est inférer en tant que Hétérozygote alors qu’il est Homozygote référant. Cette erreur est peut-être dû aux hyperparametres a savoir:

  • le taux d’erreur ici fixé a 1%
  • Une fréquence allélique des fondateurs fixé a 33% pour chaque Haplotypes

Apres modification du script une amélioration a été ajouter aux fréquences des haplotypes, celle ci se fait désormé par locus en prend en compte la fréquence d’alléle chez les fondateur. Puis on calcule la probabilitée de chaque aplotype en avec la loi de Hary-Weinberg qui fait l’hypothése de panmixie:

\[ \stackrel{AA = p² }{\stackrel{Aa = 2pq}{aa = q²}} \]

Et il y a aussi eu un ajustement du taux d’erreur de séquencage a 0.001%.