simulation meiose

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

herve chapuis
Messages : 110
Enregistré le : 05 Déc 2008, 15:26

simulation meiose

Messagepar herve chapuis » 03 Mar 2017, 15:14

Bonjour à tous.

Je voudrais simuler la transmission d'un jeu de marqueurs dans une population diploïde : A un locus donné, un animal reçoit un des deux allèles de son père et un des deux de sa mère.
Je suis en générations séparées, c'est à dire que la connaissance du génotype de tous les individus de la génération précédente doit en principe me permettre de passer à la génération suivante.

J'ai d'abord créé un pedigree , le tableau ped0 qui contient les colonnes suivantes : animal, pere, mere, sexe, generation,ipar où ipar est un indicateur binaire qui vaut 1 si l'animal a été utilisé comme parent ou 0 sinon.

Le but de l'opération est de simuler des marqueurs pour les fondateurs puis de les transmettre successivement de générations en générations.

J'ai NANIM_GEN individus par générations, et je connais NPERE et NMERE le nombre de parents, qui est constant d'une génération à l'autre.
Les générations vont de 0 (les fondateurs) à NGEN0.
Les fondateurs sont au nombre de NPERE+NMERE.

Un jeu de marqueurs consiste en NLOCI marqueurs microsatellites, qui ont tous le même nombre d'allèles. NALL. Au départ l'initialisation consiste à tirer pour chaque animal 2*NLOCI numeros d'allèles compris entre 1 et NALL.

Code : Tout sélectionner

MICROSAT<-array(NA,dim=c(nrow(ped0),NLOCI*2))
for (ianim in 1:(NPERE+NMERE)){
  MICROSAT[ianim,]<-sample(x=1:NALL,size=NLOCI*2,replace=T)
}

Peut-on vectoriser sans boucler ? Comme j'ai peu de fondateurs cela ne me dérange pas de boucler, mais si on peut gagner du temps je suis preneur car il s'agit de simulations.

Là où les boucles deviennent problématiques c'est pour mimer la transmission :

Code : Tout sélectionner

gene=seq(1,NGEN0,1)
...
for (i in gene){
  tmp<-ped0[ped0$generation==i,]
  list_anim<-tmp$animal
  tmp_microsat<-MICROSAT[list_anim,]
 
#    cat("generation =",i,"\n")
# Transmission des marqueurs
  x1<-matrix(sample(c(1,2),size=NANIM_GEN*NLOCI,replace=T),ncol=NLOCI) #marqueurs reçus du père
  x2<-matrix(sample(c(1,2),size=NANIM_GEN*NLOCI,replace=T),ncol=NLOCI) #marqueurs reçus de la mère
 
for (j in 1:NANIM_GEN) {
    for (imarq in 1:NLOCI){
    tmp_microsat[j,2*(imarq-1)+1]<-MICROSAT[tmp[j,"pere"],2*(imarq-1) +x1[j,imarq]]
    tmp_microsat[j,2*(imarq-1)+2]<-MICROSAT[tmp[j,"mere"],2*(imarq-1) +x2[j,imarq]]
}
}
  ped0[ped0$generation==i,]<-tmp
  MICROSAT[list_anim,]<-tmp_microsat
  }


Pour un animal j au locus imarq, je veux remplir les colonnes 2*(imarq-1) +1 et 2*(imarq-1)+2 qui correspondent aux allèles respectivement reçus du père et de la mère. Et à chaque fois j'ai le choix entre les deux colonnes correspondantes chez le père et chez la mère.

J'aurais voulu vectoriser, mais je n'y suis pas parvenu. Même si x1 et x2 (indicateurs de l'héritage du père ou de la mère) sont obtenus hors d'une boucle.
Donc j'ai fait une boucle par animal puis une boucle par marqueur. Mais j'ai l'impression que c'est inefficient même si c'est efficace.
J'ai choisi de commencer par la boucle animal car j'en profite pour travailler sur deux tableaux de marqueurs : le tableau MICROSAT (15 locus de 8 allèles chacun) et un tableau DNASHIP (200 SNP avec 2 allèles) que je n'ai pas mentionné mais que je remplis de la même façon.

Comment puis-je obtenir le même résultat sans passer par toutes ces boucles ?
Je suis prêt à revoir complètement l'organisation des tableaux. Peut-être faut-il utiliser une "liste" dont chaque locus serait un élément ? Appliquer ensuite une fonction type lapply ? Ou un apply sur un tableau ? J'ai l'impression que c'est possible mais ma pratique est insuffisante. Toujours des réflexes de f90...

Merci à vous !

Hervé
Ingénieur de recherche INRAE Toulouse

herve chapuis
Messages : 110
Enregistré le : 05 Déc 2008, 15:26

Re: simulation meiose

Messagepar herve chapuis » 03 Mar 2017, 16:09

J'ai une autre question, reliée à la première.
A chaque génération j'ai un tableau tmp_microsat qui a NANIM_GEN lignes et 2 * LOCI colonnes.
Ces colonnes sont constituées de chiffres de 1 à NALL.

Je voudrais connaitre la façon la plus astucieuse pour obtenir un tableau des fréquences alléliques, sachant qu'il faut à chaque fois pooler les deux origines des marqueurs (celui hérité du père et celui hérité de la mère). Je dois donc calculer la fréquence à partir de 2 * NANIM_GEN observations
Si je fais

Code : Tout sélectionner

freq<-apply(tmp_microsat,2,table)

j'obtiens une liste de 30 éléments. Je voudrais une liste de 15 éléments et pour ensuite calculer les fréquences de chaque allèle pour un locus donné.

Quand je veux passer des effectifs aux fréquences (en supposant que j'aurai obtenu les pools de deux colonnes) j'ai d'ailleurs un message :

Code : Tout sélectionner

freq/lapply(freq,sum)
Error in freq/lapply(freq, sum) :
  argument non numérique pour un opérateur binaire

Pensant faire plus simple je divise chaque élément par 200 (qui est le nombre d'animaux)

Code : Tout sélectionner

> freq/200
Error in freq/200 : argument non numérique pour un opérateur binaire

Ca ne marche pas mieux...

Tout conseil est bienvenu.
Merci à vous.
Ingénieur de recherche INRAE Toulouse

Mickael Canouil
Messages : 1315
Enregistré le : 04 Avr 2011, 08:53
Contact :

Re: simulation meiose

Messagepar Mickael Canouil » 07 Mar 2017, 08:29

Bonjour,

sans exemple fonctionnel (même "lent"), il est très difficile de proposer des solutions (temps et code).

Votre boucle:

Code : Tout sélectionner

MICROSAT<-array(NA,dim=c(nrow(ped0),NLOCI*2))
for (ianim in 1:(NPERE+NMERE)){
  MICROSAT[ianim,]<-sample(x=1:NALL,size=NLOCI*2,replace=T)
}


Une version avec lapply, cela à l'avantage de pouvoir passer en parallèle facilement (parallel::mclapply)

Code : Tout sélectionner

MICROSAT <- do.call(rbind, lapply(seq_len(NPERE+NMERE), function (ianim) {
    sample(x=1:NALL,size=NLOCI*2,replace=TRUE)
}))


Cordialement,
Mickaël
mickael.canouil.fr | rlille.fr


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Google [Bot] et 1 invité