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é