m_qtlmap_haplotype_V2
NAME
m_qtlmap_haplotype_V2
DESCRIPTION
NOTES
SEE ALSO
MAX_CAS_ALLOC
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Constants ]
NAME
MAX_CAS_ALLOC
DESCRIPTION
The maximum number of marker to use the pded_V5 analysis
NOTES
Constante pour eviter les segmentation fault du au overflow 2 < MAX_CAS_ALLOC < 50
SOURCE
51 integer , private, parameter :: MAX_CAS_ALLOC = 22
MAX_MARKER
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Constants ]
NAME
MAX_MARKER
DESCRIPTION
The maximum number of marker to use the pded_V5 analysis
NOTES
there are a lot of allocation 2**nmk ==> overflow if nmk is too big
SOURCE
40 integer, private, parameter :: MAX_MARKER = 39
ordre
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Variables ]
NAME
ordre
DESCRIPTION
NOTES
dim nr, lx
ordref
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Variables ]
NAME
ordref
DESCRIPTION
NOTES
dim nm, lx
ordrem
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Variables ]
NAME
ordrem
DESCRIPTION
NOTES
dim nm, lx
ordrep
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Variables ]
NAME
ordrep
DESCRIPTION
NOTES
Ordre 12 ou 21... : dimension (np,lx)
ancetre
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
ancetre
DESCRIPTION
NOTES
Mise a zero des probabilites des phases impossibles pour les peres et pour les meres d apres les phenotypes des grands parents
SOURCE
303 subroutine ancetre () 304 ! 305 ! Divers 306 logical typgp,typgm 307 integer(kind=KIND_PHENO) :: m1,m2,mc1,mc2,md1,md2 308 ! 309 integer ll,igp,jgm,ip,jm,kr,c 310 integer ngm1,ngm2,nr1,nr2 311 ! 312 ! Recherche des haplotypes impossibles pour les parents 313 ! Recherche des phenotypes des grands parents 314 do c=1,nchr 315 do ll=1,nmk(c) 316 do 10 igp=1,ngp 317 typgp=.false. 318 if(corregp(igp).ne.INT_NOT_DEFINED) then 319 if(pheno(c,ll,corregp(igp),1).ne.nmanque) then 320 typgp=.true. 321 m1=pheno(c,ll,corregp(igp),1) 322 m2=pheno(c,ll,corregp(igp),2) 323 end if 324 end if 325 ngm1=ngmgp(igp)+1 326 ngm2=ngmgp(igp+1) 327 do 10 jgm=ngm1,ngm2 328 typgm=.false. 329 if(corregm(jgm).ne.INT_NOT_DEFINED)then 330 if (pheno(c,ll,corregm(jgm),1).ne.nmanque) then 331 typgm=.true. 332 mc1=pheno(c,ll,corregm(jgm),1) 333 mc2=pheno(c,ll,corregm(jgm),2) 334 end if 335 end if 336 ! 337 ! Passage en revue des parents 338 nr1=nrgm(jgm)+1 339 nr2=nrgm(jgm+1) 340 do 10 kr=nr1,nr2 341 ordre(c,kr,ll)=0 342 if(correr(kr).eq.INT_NOT_DEFINED) go to 10 343 if(pheno(c,ll,correr(kr),1).eq.nmanque) go to 10 344 md1=pheno(c,ll,correr(kr),1) 345 md2=pheno(c,ll,correr(kr),2) 346 ! 347 ! Parent heterozygote 348 ! Un allele different chez un des grands parents 349 if((typgp).and.(md1.ne.m1.and.md1.ne.m2)) ordre(c,kr,ll)=21 350 if((typgp).and.(md2.ne.m1.and.md2.ne.m2)) ordre(c,kr,ll)=12 351 if((typgm).and.(md1.ne.mc1.and.md1.ne.mc2)) ordre(c,kr,ll)=12 352 if((typgm).and.(md2.ne.mc1.and.md2.ne.mc2)) ordre(c,kr,ll)=21 353 ! 354 10 continue 355 end do 356 ! 357 ! Stockage des probabilites pour les peres et les meres 358 ! La gestion des valeurs de phasp et phasm (d�crite ici en commentaires) est 359 ! d�port�e dans get_phasp et get_phasm 360 do ll=1,nmk(c) 361 do ip=1,np 362 ! phasp(c,ip)=.false. 363 ordrep(c,ip,ll)=0 364 if(reppere(ip).ne.INT_NOT_DEFINED) ordrep(c,ip,ll)=ordre(c,reppere(ip),ll) 365 ! if(ordrep(ip,ll).ne.0)phasp(c,ip)=.true. 366 end do 367 do jm=1,nm 368 !phasm(c,jm)=.false. 369 ordrem(c,jm,ll)=0 370 ordref(c,repfem(jm),ll)=0 371 if(repmere(jm).ne.INT_NOT_DEFINED) then 372 ordrem(c,jm,ll)=ordre(c,repmere(jm),ll) 373 ordref(c,repfem(jm),ll)=ordre(c,repmere(jm),ll) 374 end if 375 ! if(ordrem(jm,ll).ne.0)phasm(c,jm)=.true. 376 end do 377 end do 378 end do 379 end subroutine ancetre
calcul_phases_symmax2sat_dam
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
probin
DESCRIPTION
NOTES
SOURCE
4792 subroutine calcul_phases_symmax2sat_dam 4793 4794 integer :: c,nm1,nm2,m1,m2,ip,imark,kkd,nd1,nd2,ik,ik2,ll,jm,kd,m,kdIfem,ifem,ll1,ll2 4795 real(kind=dp) :: p1,p2,wmax 4796 real(kind=dp) , dimension(nd,maxval(nmk)+1) :: trans 4797 integer , dimension(maxval(nmk)) :: H 4798 integer , dimension(maxval(nmk),maxval(nmk)) :: nbp,nbm 4799 real(kind=dp) , dimension(maxval(nmk),maxval(nmk)) :: W 4800 integer , dimension(nfem) :: ndf,corref 4801 integer , dimension(nfem,ndm(nm+1)) :: repdes 4802 integer , dimension(nm,maxval(nmk)) :: phm 4803 real(kind=dp) , dimension(maxval(nmk),maxval(nmk)) :: r 4804 logical :: ok 4805 4806 trans = 0 4807 4808 repdes=0 4809 ndf=0 4810 W=0.d0 4811 do jm=1,nm 4812 ifem=repfem(jm) 4813 do kd=ndm(jm)+1,ndm(jm+1) 4814 ndf(ifem)=ndf(ifem)+1 4815 repdes(ifem,ndf(ifem))=kd 4816 end do 4817 corref(ifem)=correm(jm) 4818 end do 4819 4820 !****************************************************************************** 4821 ! 4822 ! on traite les femelles (et non les meres) une a une 4823 ! 4824 ! nombfem est le nombre de femelles dont on peut chercher les phases 4825 ! 4826 do c=1,nchr 4827 do ik=1,nmk(c)-1 4828 do ik2=ik+1,nmk(c) 4829 r(ik,ik2) = log( (1.d0-rm(c,ik,ik2))/rm(c,ik,ik2) ) 4830 end do 4831 end do 4832 4833 do ifem=1,nfem 4834 if(estfem(ifem))then 4835 kkd=0 4836 do kdIfem=1,ndf(ifem) 4837 ! 4838 !****************************************************************************** 4839 ! 4840 ! Si le nombre de pleins freres est trop faible, le genotype de la mere ne sera pas considere 4841 ! Dans ce cas le premier genotype rencontre possible est considere comme le bon 4842 4843 !if(ndf(ifem).lt.ndmin) then 4844 ! force(jm)=.true. 4845 4846 ! 4847 ! creation du tableau des transmissions observees 4848 ! 4849 kd = repdes(ifem,kdIfem) 4850 kkd=kkd+1 4851 do ll=1,nmk(c) 4852 trans(kkd,ll)=0 4853 p1=prot(c,ll,kd,1)+prot(c,ll,kd,3) 4854 if(p1.eq.1.d0)trans(kkd,ll)=1 4855 p2=prot(c,ll,kd,2)+prot(c,ll,kd,4) 4856 if(p2.eq.1.d0)trans(kkd,ll)=2 4857 ! if ( ll == 1 .and. ifem == 4 ) then 4858 ! 4859 ! print *,'kkd:',kkd,' kd:',kd,'p1:',p1,' p2:',p2,'=>',trans(kkd,ll),' prot:',prot(c,ll,kd,:),animal(kd),& 4860 ! get_pheno(c,pheno(c,ll,corred(kd),1)),get_pheno(c,pheno(c,ll,corred(kd),2)) 4861 ! end if 4862 end do ! fin ll 4863 end do ! fin kdIfem 4864 !******************************* 4865 4866 nbp=0 4867 nbm=0 4868 kkd=0 4869 4870 do kdIfem=1,ndf(ifem) 4871 kd = repdes(ifem,kdIfem) 4872 kkd=kkd+1 4873 ik=1 4874 do while (ik<nmk(c)) 4875 do while (ik<=(nmk(c)-1) .and. trans(kkd,ik) == 0 ) 4876 ik=ik+1 4877 end do 4878 if ( ik >= nmk(c)) exit 4879 ik2=ik+1 4880 do while (ik2<=nmk(c) .and. trans(kkd,ik2) == 0 ) 4881 ik2=ik2+1 4882 end do 4883 if ( ik2 > nmk(c)) exit 4884 if ( trans(kkd,ik) == trans(kkd,ik2) ) then 4885 nbp(ik,ik2) = nbp(ik,ik2) + 1 4886 else 4887 nbm(ik,ik2) = nbm(ik,ik2) + 1 4888 end if 4889 4890 ik=ik2 4891 end do 4892 end do 4893 4894 wmax=0.d0 4895 do ik=1,nmk(c)-1 4896 do ik2=ik+1,nmk(c) 4897 if ( nbp(ik,ik2) == nbm(ik,ik2) ) cycle 4898 W(ik,ik2) = 0.25d0*dble(nbp(ik,ik2)-nbm(ik,ik2))*r(ik,ik2) 4899 if(wmax < dabs(W(ik,ik2)))wmax=dabs(W(ik,ik2)) 4900 end do !ik2 4901 end do !ik 4902 4903 do ik=1,nmk(c) 4904 if(ordref(c,ifem,ik)==12)W(ik,ik)=wmax+1.d0 4905 if(ordref(c,ifem,ik)==21)W(ik,ik)=-wmax-1.d0 4906 ! write(6,*)(W(ik,ik2),ik2=1,nmk(c)) 4907 end do 4908 4909 ok = get_h_from_w(maxval(nmk),W,H,m) 4910 ! 4911 ! 4912 ! stockage des resultats par mere (et non par femelle) 4913 ! et reperage de la phase la plus probable 4914 ! 4915 4916 4917 ! 4918 ! Update the genotypm structure 4919 do ll=1,nmk(c) 4920 if(h(ll) ==1) then 4921 do jm=1,nm 4922 if ( repfem(jm) == ifem ) then 4923 genotypm(c,ll,jm,1)=pheno(c,ll,corref(ifem),1) 4924 genotypm(c,ll,jm,2)=pheno(c,ll,corref(ifem),2) 4925 phm(jm,ll)=1 4926 end if 4927 end do 4928 else 4929 do jm=1,nm 4930 if ( repfem(jm) == ifem ) then 4931 genotypm(c,ll,jm,1)=pheno(c,ll,corref(ifem),2) 4932 genotypm(c,ll,jm,2)=pheno(c,ll,corref(ifem),1) 4933 phm(jm,ll)=2 4934 end if 4935 end do 4936 end if 4937 end do ! ll 4938 end if 4939 end do ! ifem 4940 end do ! chromosome 4941 4942 allocate ( ngend(nchr,nm+1) ) 4943 allocate ( probg(nchr,nm) ) 4944 allocate ( ndesc(nchr,nd) ) 4945 4946 do c=1,nchr 4947 ngenom(c,1)=0 4948 do jm=1,nm 4949 probg(c,jm)=1.d0 4950 ngenom(c,jm+1)=ngenom(c,jm)+1 4951 ngend(c,jm)=ndm(jm) 4952 ngend(c,jm+1)=ndm(jm+1) 4953 do kd=ndm(jm)+1,ndm(jm+1) 4954 ndesc(c,kd)=kd 4955 end do 4956 do ll=1,nmk(c) 4957 if ( phm(jm,ll) == 1 ) then 4958 ptfin(c,ll,ngend(c,jm)+1:ngend(c,jm+1),:)=prot(c,ll,ngend(c,jm)+1:ngend(c,jm+1),:) 4959 else 4960 ptfin(c,ll,ngend(c,jm)+1:ngend(c,jm+1),1)=prot(c,ll,ngend(c,jm)+1:ngend(c,jm+1),2) 4961 ptfin(c,ll,ngend(c,jm)+1:ngend(c,jm+1),2)=prot(c,ll,ngend(c,jm)+1:ngend(c,jm+1),1) 4962 ptfin(c,ll,ngend(c,jm)+1:ngend(c,jm+1),3)=prot(c,ll,ngend(c,jm)+1:ngend(c,jm+1),4) 4963 ptfin(c,ll,ngend(c,jm)+1:ngend(c,jm+1),4)=prot(c,ll,ngend(c,jm)+1:ngend(c,jm+1),3) 4964 end if 4965 end do !! 4966 end do !jm 4967 end do !c 4968 4969 ! ll1=5;ll2=7; 4970 ! do ll=ll1,ll2 4971 ! print *,'MARKER:',ll 4972 ! do ip=1,np 4973 ! do jm=nmp(ip)+1,nmp(ip+1) 4974 ! !print *,trim(mere(jm)),jm 4975 ! if ( trim(mere(jm)) == '634167') then 4976 ! print *,'papa:',get_pheno(1,pheno(1,ll,correp(ip),1)),get_pheno(1,pheno(1,ll,correp(ip),2)) 4977 ! print *,'maman:',get_pheno(1,pheno(1,ll,correm(jm),1)),get_pheno(1,pheno(1,ll,correm(jm),2)) 4978 ! 4979 ! do kd=ndm(jm)+1,ndm(jm+1) 4980 ! print *,get_pheno(1,pheno(1,ll,corred(kd),1)),get_pheno(1,pheno(1,ll,corred(kd),2)) 4981 ! end do 4982 ! end if 4983 ! 4984 ! end do 4985 ! end do 4986 ! end do 4987 ! stop 4988 4989 end subroutine calcul_phases_symmax2sat_dam
calcul_phases_symmax2sat_sire
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
probin
DESCRIPTION
NOTES
SOURCE
4656 subroutine calcul_phases_symmax2sat_sire 4657 4658 integer :: c,nm1,nm2,m1,m2,ip,imark,kkd,nd1,nd2,ik,ik2,ll,jm,kd,m 4659 real(kind=dp) :: p1,p2,wmax 4660 real(kind=dp) , dimension(nd,maxval(nmk)+1) :: trans 4661 integer , dimension(maxval(nmk)) :: H 4662 integer , dimension(maxval(nmk),maxval(nmk)) :: nbp,nbm 4663 real(kind=dp) , dimension(maxval(nmk),maxval(nmk)) :: W 4664 real(kind=dp) , dimension(maxval(nmk),maxval(nmk)) :: r 4665 4666 logical :: ok 4667 integer :: c_11,c_12,c_21,c_22,ll1,ll2 4668 4669 trans=0 4670 ! 4671 ! Etablissement de la phase la plus probable 4672 ! 4673 do c=1,nchr 4674 do ik=1,nmk(c)-1 4675 do ik2=ik+1,nmk(c) 4676 r(ik,ik2) = log( (1.d0-rm(c,ik,ik2))/rm(c,ik,ik2) ) 4677 end do !ik2 4678 end do !ik 4679 4680 do ip=1, np 4681 4682 W=0 4683 ! 4684 ! on reduit le jeu de marqueurs e ceux pour lesquels le pere est heterozygote 4685 ! 4686 ! creation du tableau des transmissions observees 4687 ! 4688 kkd=0 4689 do jm=nmp(ip)+1,nmp(ip+1) 4690 do kd=ndm(jm)+1,ndm(jm+1) 4691 kkd=kkd+1 4692 do ll=1,nmk(c) 4693 trans(kkd,ll)=0 4694 p1=prot(c,ll,kd,1)+prot(c,ll,kd,2) 4695 if(p1.eq.1.d0)trans(kkd,ll)=1 4696 p2=prot(c,ll,kd,3)+prot(c,ll,kd,4) 4697 if(p2.eq.1.d0)trans(kkd,ll)=2 4698 end do !ll 4699 end do !kd 4700 end do !jm 4701 4702 nbp=0 4703 nbm=0 4704 kkd=0 4705 do jm=nmp(ip)+1,nmp(ip+1) 4706 do kd=ndm(jm)+1,ndm(jm+1) 4707 kkd=kkd+1 4708 ik=1 4709 do while (ik<nmk(c)) 4710 do while (ik<=(nmk(c)-1) .and. trans(kkd,ik) == 0 ) 4711 ik=ik+1 4712 end do 4713 if ( ik >= nmk(c)) exit 4714 ik2=ik+1 4715 do while (ik2<=nmk(c) .and. trans(kkd,ik2) == 0 ) 4716 ik2=ik2+1 4717 end do 4718 if ( ik2 > nmk(c)) exit 4719 if ( trans(kkd,ik) == trans(kkd,ik2) ) then 4720 nbp(ik,ik2) = nbp(ik,ik2) + 1 4721 else 4722 nbm(ik,ik2) = nbm(ik,ik2) + 1 4723 end if 4724 4725 ik=ik2 4726 end do !while(ik<nmj(c) 4727 end do !kd 4728 end do !jm 4729 4730 !write(6,*)'ordrep',(ordrep(c,ip,ik),ik=1,nmk(c)) 4731 4732 wmax=0.d0 4733 do ik=1,nmk(c)-1 4734 do ik2=ik+1,nmk(c) 4735 if ( nbp(ik,ik2) == nbm(ik,ik2) ) cycle 4736 W(ik,ik2) = 0.25d0*dble(nbp(ik,ik2)-nbm(ik,ik2))*r(ik,ik2) 4737 if(wmax < dabs(W(ik,ik2)))wmax=dabs(W(ik,ik2)) 4738 end do !ik2 4739 end do !ik 4740 4741 4742 4743 do ik=1,nmk(c) 4744 if(ordrep(c,ip,ik)==12)W(ik,ik)=wmax+1.d0 4745 if(ordrep(c,ip,ik)==21)W(ik,ik)=-wmax-1.d0 4746 ! write(6,*)(W(ik,ik2),ik2=1,nmk(c)) 4747 end do 4748 4749 ok = get_h_from_w(maxval(nmk),W,H,m) 4750 4751 ! 4752 ! stockage du genotype 4753 ! 4754 do ll=1,nmk(c) 4755 if(h(ll) ==1) then 4756 genotyp(c,ll,correp(ip),1)=pheno(c,ll,correp(ip),1) 4757 genotyp(c,ll,correp(ip),2)=pheno(c,ll,correp(ip),2) 4758 else 4759 genotyp(c,ll,correp(ip),1)=pheno(c,ll,correp(ip),2) 4760 genotyp(c,ll,correp(ip),2)=pheno(c,ll,correp(ip),1) 4761 end if 4762 end do !ll 4763 ! 4764 ! Reorganisation du tableau des probabilites de transmission 4765 do jm=nmp(ip)+1,nmp(ip+1) 4766 do kd=ndm(jm)+1,ndm(jm+1) 4767 do ll=1,nmk(c) 4768 if(h(ll).eq.1) then 4769 p1=prot(c,ll,kd,1) 4770 p2=prot(c,ll,kd,2) 4771 prot(c,ll,kd,1)=prot(c,ll,kd,3) 4772 prot(c,ll,kd,2)=prot(c,ll,kd,4) 4773 prot(c,ll,kd,3)=p1 4774 prot(c,ll,kd,4)=p2 4775 end if 4776 end do !ll 4777 end do !kd 4778 end do !jm 4779 end do ! ip 4780 end do ! chromosome 4781 end subroutine calcul_phases_symmax2sat_sire
check_recombination_sire
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
check_recombination_sire
DESCRIPTION
NOTES
verification de la coherence entre la carte et les recombinaisons observees intra pere (meme origine grand parentale) et en opposition. Le taux de recombinaison observe est le rapport nbr en opposition / nbr total>. Connaissant ( d’apres la carte) la distance entre ces deux marqueurs ( donc le taux de recombinaison theorique) , ! on peut calculer la probabilite que ce taux de recombinaison observe soit obtenu ( la loi est binomial de parameters (nombre total de descendants dont on connait la transmission aux deux marqueurs , taux de recombinaison theorique) . Si cette proba est trop faible, on lance un warning
SOURCE
4485 subroutine check_recombination_sire 4486 4487 ! Divers 4488 integer(kind=KIND_PHENO) :: m1,m2 4489 integer imark,ip,jm,kd,kkd,c 4490 integer ll,ll1,ll2 4491 integer trans(nd,maxval(nmk)) 4492 integer oppos,phase,total 4493 integer nd1,nd2,nm1,nm2 4494 real (kind=dp) :: p1,p2,recomb,pp 4495 logical :: hetero(maxval(nmk)) 4496 4497 call log_mess('checking recombination in sires...',INFO_DEF) 4498 ! 4499 4500 do c=1,nchr 4501 do 100 ip=1, np 4502 nm1=nmp(ip)+1 4503 nm2=nmp(ip+1) 4504 ! 4505 ! on cree le tableau hetero qui liste l'heterozygotie des marqueurs 4506 ! 4507 do imark=1,nmk(c) 4508 hetero(imark)=.true. 4509 if(correp(ip).eq.9999) then 4510 m1=nmanque 4511 else 4512 m1=pheno(c,imark,correp(ip),1) 4513 m2=pheno(c,imark,correp(ip),2) 4514 end if 4515 if(m1.eq.nmanque) m2=m1 4516 if(m1.eq.m2) then 4517 hetero(imark)=.false. 4518 end if 4519 end do 4520 ! 4521 ! on reduit le jeu de marqueurs a ceux pour lesquels le pere est heterozygote 4522 ! 4523 ! creation du tableau des transmissions observees 4524 ! 4525 kkd=0 4526 do jm=nm1,nm2 4527 nd1=ndm(jm)+1 4528 nd2=ndm(jm+1) 4529 do kd=nd1,nd2 4530 kkd=kkd+1 4531 do ll=1,nmk(c) 4532 trans(kkd,ll)=0 4533 p1=prot(c,ll,kd,1)+prot(c,ll,kd,2) 4534 if(p1.eq.1.d0)trans(kkd,ll)=1 4535 p2=prot(c,ll,kd,3)+prot(c,ll,kd,4) 4536 if(p2.eq.1.d0)trans(kkd,ll)=2 4537 end do 4538 end do 4539 end do 4540 ! 4541 do ll1=1,nmk(c) 4542 if(hetero(ll1)) go to 10 4543 end do 4544 ! 4545 ! le pere est entierement homozygote 4546 ! 4547 go to 100 4548 ! 4549 10 continue 4550 4551 ! ll1 et ll2 etant les indices de 2 marqueurs heterozygotes successifs 4552 ! (tous les marqueurs intermediaires sont homozygotes), on compare 4553 ! les transmission en ll1 et ll2 4554 ! 4555 ll2=ll1 4556 50 ll2=ll2+1 4557 4558 !DEB OFI : 22/03/2010 : correction bug 4559 if (ll2> nmk(c)) go to 100 4560 !FIN OFI 4561 4562 if(hetero(ll2)) go to 60 4563 if(ll2 < nmk(c)) go to 50 4564 go to 100 4565 60 continue 4566 ! 4567 ! comptage des inversions apparentes des phases locales 4568 ! 4569 phase=0 4570 oppos=0 4571 do kkd=1,ndm(nm2+1)-ndm(nm1) 4572 if(trans(kkd,ll1).ne.0.and.trans(kkd,ll2).ne.0) then 4573 if(trans(kkd,ll1).eq.trans(kkd,ll2))then 4574 phase=phase+1 4575 else 4576 oppos=oppos+1 4577 end if 4578 end if 4579 end do 4580 4581 total=phase+oppos 4582 if(phase < oppos)oppos=phase 4583 4584 recomb=xaldane(posim(c,ll2)-posim(c,ll1)) 4585 call probin(total,recomb,oppos,pp) 4586 if(pp < prob_seuil_recomb) then 4587 write(nficout,*)'Too many recombination events (',oppos,'/',total,') between Marker ['& 4588 //trim(mark(c,ll1))//'] and ['//trim(mark(c,ll2))//& 4589 '] for sire :'//trim(pere(ip))//' ' 4590 ! call log_mess('Too many recombination events between Marker ['& 4591 ! //trim(mark(c,ll1))//'] and ['//trim(mark(c,ll2))//& 4592 ! '] for sire :'//trim(pere(ip))//' ',WARNING_DEF) 4593 end if 4594 4595 ll1=ll2 4596 if(ll2 < nmk(c)) go to 50 4597 4598 100 continue 4599 end do 4600 call log_mess('end check_recombination_sire',INFO_DEF) 4601 return 4602 end subroutine check_recombination_sire
free_internal_struct
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
free_internal_struct
DESCRIPTION
NOTES
SOURCE
157 subroutine free_internal_struct 158 159 deallocate (prot) 160 deallocate (ptfin) 161 deallocate (ordrep) 162 deallocate (ordrem) 163 deallocate (ordre) 164 deallocate (ordref) 165 166 end subroutine free_internal_struct
haplotype_SNP
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
haplotype_SNP
DESCRIPTION
NOTES
SOURCE
248 subroutine haplotype_SNP() 249 250 call set_allocation 251 call ancetre() 252 call set_phasp 253 call set_phasm 254 call gammapf(prot) 255 call pdegp_snp() 256 call pdegm_snp() 257 call setting_alloc_pdd() 258 !call pded_v5() 259 call pded_v5_optim() 260 call check_recombination_sire 261 call free_internal_struct 262 263 end subroutine haplotype_SNP
haplotype_SYMMAX2SAT_SNP
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
haplotype_SNP
DESCRIPTION
NOTES
SOURCE
275 subroutine haplotype_SYMMAX2SAT_SNP 276 277 call set_allocation 278 call ancetre() 279 call set_phasp 280 call set_phasm 281 call gammapf(prot) 282 ! call pdegp_snp() 283 ! call pdegm_snp() 284 285 call calcul_phases_symmax2sat_sire() 286 call calcul_phases_symmax2sat_dam() 287 call setting_alloc_pdd() 288 call pded_v5_optim() 289 call check_recombination_sire 290 call free_internal_struct 291 292 end subroutine haplotype_SYMMAX2SAT_SNP
haplotype_V2
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
haplotype_V2
DESCRIPTION
NOTES
SOURCE
198 subroutine haplotype_V2() 199 200 call set_allocation 201 call ancetre() 202 call set_phasp 203 call set_phasm 204 call gammapf(prot) 205 call pdegp() 206 call pdegm() 207 call setting_alloc_pdd() 208 call pded() 209 call check_recombination_sire 210 call free_internal_struct 211 212 end subroutine haplotype_V2
haplotype_V3
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
haplotype_V3
DESCRIPTION
NOTES
SOURCE
223 subroutine haplotype_V3() 224 call set_allocation 225 call ancetre() 226 call set_phasp 227 call set_phasm 228 call gammapf(prot) 229 call pdegp() 230 call pdegm() 231 call setting_alloc_pdd() 232 !call pded_v5() 233 call pded_v5_optim() 234 call check_recombination_sire 235 call free_internal_struct 236 237 end subroutine haplotype_V3
pded
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pded
DESCRIPTION
NOTES
Calcul des probabilites de transmission le long du chromosome
SOURCE
1262 subroutine pded() 1263 1264 ! 1265 ! Tableaux dimensionnes selon l, le nombre de marqueurs 1266 integer :: cas(maxval(nmk)),trans(maxval(nmk)),infor(maxval(nmk)) 1267 ! dim : nbt*4+nbt,nbt 1268 integer , dimension(:,:),allocatable :: last 1269 real (kind=dp), dimension(:,:),allocatable :: pbt 1270 1271 integer :: index_mark(maxval(nmk)) 1272 logical :: index_mark_valide(maxval(nmk)) 1273 1274 ! 1275 ! Divers 1276 integer(kind=1):: new(11,8) 1277 integer geno,npas,removeCas4,maxcas4,removeCas2,maxCas2 1278 integer ilong,ip,it,jm,jt,rc4,rc2,lit,lkkk,c 1279 integer n,kd,kt,ii 1280 integer mark_selec,im1,im2 1281 integer linf,lk,llk,lmax,lq,lt,maxCas 1282 integer nbt,nd1,nd2,ngeno1,ngeno2,nm1,nm2,nmkinf 1283 real (kind=dp) pos_m,posi_lk,posi_lkinf 1284 real (kind=dp) dg,dd,ddf,ddm,dgf,dgm 1285 real (kind=dp) dx,pbtt,val_min 1286 real (kind=dp) recm,recf,recmlq,recflq,recmqr,recfqr 1287 real (kind=dp) xint,xintm,xintf 1288 !external xaldane 1289 !double precision xaldane 1290 call log_mess('calling...PDED_V4',DEBUG_DEF) 1291 ! 1292 ! tableau des transmission 1293 ! 1294 data new / 1,1,2,2,1,1,1,2,1,1,1, & 1295 1,2,1,2,1,2,1,1,1,2,1, & 1296 0,0,0,0,2,2,1,2,2,2,1, & 1297 0,0,0,0,2,1,2,2,1,2,2, & 1298 0,0,0,0,0,0,0,0,0,0,2, & 1299 0,0,0,0,0,0,0,0,0,0,1, & 1300 0,0,0,0,0,0,0,0,0,0,2, & 1301 0,0,0,0,0,0,0,0,0,0,2/ 1302 1303 1304 ! init pour la suite 1305 allocate (last(2**(((MAX_CAS_ALLOC)/2)+1),2)) 1306 allocate (pbt(4,2**(((MAX_CAS_ALLOC)/2)+1))) 1307 do c=1,nchr 1308 ! 1309 ! Taille du segment explore 1310 ilong=get_ilong(c) 1311 infor = 0 1312 ! 1313 ! Calcul pour chaque descendant de la probabilite de transmission 1314 do 997 ip=1,np 1315 !print *,'***************************PERE:',trim(pere(ip)) 1316 nm1=nmp(ip)+1 1317 nm2=nmp(ip+1) 1318 do 998 jm=nm1,nm2 1319 !print *,'**********************MERE:',trim(mere(jm)) 1320 ngeno1=ngenom(c,jm)+1 1321 ngeno2=ngenom(c,jm+1) 1322 do 999 geno=ngeno1,ngeno2 1323 ! print *,'**********************GENO:',geno 1324 nd1=ngend(c,geno)+1 1325 nd2=ngend(c,geno+1) 1326 !print *,'ND1:',nd1,' ND2:',nd2 1327 do 1000 kd=nd1,nd2 1328 ! 1329 ! le vecteur trans(lk) indique l'information jointe sur la transmission au locus lk 1330 ! le vecteur infor contient la liste des marqeurs doublement informatifs 1331 ! 1332 ! l'origine grand parentale des alleles recu par le descendant k au marqueur lk 1333 ! depend de trans (tableau new) 1334 ! soit 1 pour grand pere et 2 pour grand mere 1335 ! 1336 ! cas 1 cas 2 cas 3 cas 4 1337 ! trans(lk) pere mere pere mere pere mere pere mere 1338 ! 1 1 1 1339 ! 2 1 2 1340 ! 3 2 1 1341 ! 4 2 2 1342 ! 5 1 1 2 2 1343 ! 6 1 2 2 1 1344 ! 7 1 1 1 2 1345 ! 8 2 1 2 2 1346 ! 9 1 1 2 1 1347 ! 10 1 2 2 2 1348 ! 11 1 1 1 2 2 1 2 2 1349 ! 1350 ! 1351 llk=0 1352 do lk=1,nmk(c) 1353 1354 ! 1355 ! cas(lk) est le nombre d'evenements de transmission possible pour le marqueur lk 1356 ! cas(lk) vaut 1 pour trans(lk) = 1 e 4, 1357 ! 2 pour trans(lk) = 5 e 10, 1358 ! 4 pour trans(lk) = 11 1359 ! 1360 cas(lk)=2 1361 ! 1362 ! incertitude totale (donnees manquantes) 1363 ! 1364 trans(lk)=11 1365 ! 1366 ! triplet heterozygotes 1367 ! 1368 if(ptfin(c,lk,kd,2).eq.0.and.ptfin(c,lk,kd,3).eq.0)trans(lk)=5 1369 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=6 1370 ! 1371 ! incertitude sur un parent 1372 ! 1373 if(ptfin(c,lk,kd,3).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=7 1374 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,2).eq.0)trans(lk)=8 1375 if(ptfin(c,lk,kd,2).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=9 1376 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,3).eq.0)trans(lk)=10 1377 ! 1378 ! cas certains 1379 ! 1380 if(ptfin(c,lk,kd,1).eq.1.d0) trans(lk)=1 1381 if(ptfin(c,lk,kd,2).eq.1.d0) trans(lk)=2 1382 if(ptfin(c,lk,kd,3).eq.1.d0) trans(lk)=3 1383 if(ptfin(c,lk,kd,4).eq.1.d0) trans(lk)=4 1384 1385 if(trans(lk).le.4) then 1386 llk=llk+1 1387 infor(llk)=lk 1388 cas(lk)=1 1389 end if 1390 1391 ! 1392 ! reajustement de cas 1393 ! 1394 if(trans(lk).eq.11) cas(lk)=4 1395 end do 1396 1397 nmkinf=llk 1398 1399 ! 1400 ! exploration du groupe de liaison pour la recherche des transmissions possibles 1401 ! le groupe de liaison est separe en segments flanques de MIF 1402 ! avec deux segments externes pouvant ne comporter qu'un seul Marqueur Informatif 1403 ! e droite pour celui de gauche et e gauche pour celui de droite 1404 ! 1405 n=1 1406 linf=1 1407 100 dx=absi(c,n) 1408 1409 ! 1410 ! cas des descendants sans information 1411 ! 1412 if(nmkinf.eq.0) then 1413 linf=0 1414 do ii=1,4 1415 pdd(c,kd,ii,n)=0.25d0 1416 end do 1417 go to 125 1418 end if 1419 1420 1421 if(dx.le.posi(c,infor(linf)).or.linf.eq.nmkinf) go to 105 1422 linf=linf+1 1423 go to 100 1424 1425 105 continue 1426 ! 1427 ! pour des problemes d'arrondis, la position testee pour la qtl 1428 ! peut depasser le bout du groupe de liaison... on l'y ramene 1429 ! 1430 if(dx.ge.posi(c,nmk(c)))dx=posi(c,nmk(c)) 1431 ! 1432 ! cas des positions sur des marqueurs informatifs 1433 ! 1434 if(dx.eq.posi(c,infor(linf))) then 1435 1436 nbt=1 1437 1438 do lq=1,4 1439 pbt(lq,1)=0.d0 1440 end do 1441 1442 pbt(trans(infor(linf)),1)=1.d0 1443 go to 120 1444 1445 end if 1446 ! 1447 ! cas des positions entre des marqueurs informatifs 1448 ! trois situations 1449 ! si dx est avant le premier MIF, on explore les transmissions entre le premier 1450 ! marqueur et le premier MIF 1451 ! si dx est entre deux MIF (cas standard) 1452 ! si dx est apres le dernier MIF, on explore les transmissions entre le dernier MIF et le 1453 ! dernier marqueur 1454 ! 1455 if(dx.lt.posi(c,infor(1))) then 1456 lk=1 1457 lmax=infor(1) 1458 end if 1459 if(dx.gt.posi(c,infor(1)).and.dx.lt.posi(c,infor(nmkinf))) then 1460 lk=infor(linf-1) 1461 lmax=infor(linf) 1462 end if 1463 if(dx.ge.posi(c,infor(nmkinf))) then 1464 lk=infor(nmkinf) 1465 lmax=nmk(c) 1466 end if 1467 1468 ! ---------- AJOUT OFI ALLOCATION STRUCTURES TEMPORAIRES ---------- 1469 !taille des tableaux Max : cas(linf)*cas(linf+1)*..*cas(lmax)=nbt 1470 maxCas = 0 1471 maxCas = MAX_CAS_ALLOC 1472 ! On recherche un echantillonage 1473 removeCas2=0 1474 removeCas4=0 1475 maxcas2=0 1476 maxcas4=0 1477 1478 ! print *,'LK:',lk,'lMAX:',lmax 1479 do lq=lk,lmax 1480 ! print *,lq,':',cas(lq) 1481 if (cas(lq) == 2) maxcas2 = maxcas2 + 1 1482 if (cas(lq) == 4) maxcas4 = maxcas4 + 1 1483 end do 1484 !print *,'MAX CAS:',maxCas, ' nb4:',maxcas4,' nb2',maxcas2 1485 do while ( maxCas >= MAX_CAS_ALLOC) 1486 rc2=removeCas2 1487 rc4=removeCas4 1488 maxCas=0 1489 do lq=lk,lmax 1490 if (cas(lq) > 1) then 1491 if (cas(lq)==4 .and. rc4>0) then 1492 rc4 = rc4-1 1493 cycle 1494 end if 1495 if (cas(lq)==2 .and. rc2>0) then 1496 rc2 = rc2-1 1497 cycle 1498 end if 1499 maxCas = maxCas + cas(lq) 1500 end if 1501 end do 1502 1503 if (maxCas>=MAX_CAS_ALLOC .and. removeCas4<maxcas4) then 1504 removeCas4 = removeCas4 +1 1505 else if (maxCas>=MAX_CAS_ALLOC .and. removeCas2<maxcas2) then 1506 removeCas2 = removeCas2 +1 1507 else if (maxCas>=MAX_CAS_ALLOC) then 1508 call stop_application("-- error dev --") 1509 end if 1510 end do 1511 ! print *,'MAX CAS:',maxCas, ' r4:',removeCas4,' r2',removeCas2 1512 1513 index_mark_valide = .false. 1514 do lq=lk,lmax 1515 index_mark(lq-lk+1)=lq 1516 index_mark_valide(lq-lk+1) = .true. 1517 end do 1518 lq=lk 1519 1520 do while ( removeCas4>0 ) 1521 val_min = 0.d0 1522 mark_selec = -1 1523 ! print *,'lk:',lk 1524 do lit=lk,lmax 1525 if (index_mark_valide(lit-lk+1) .and. cas(lit)==4) then 1526 pos_m=abs(dx-posi(c,lit)) 1527 !print *,'lit ',lit,' pos:',pos_m,' val_min:',val_min 1528 if (pos_m>=val_min) then 1529 ! print *,'min!!! lit:',lit 1530 val_min = pos_m 1531 mark_selec = lit 1532 end if 1533 end if 1534 end do 1535 1536 if (mark_selec == -1) call stop_application("ERROR"); 1537 ! print *,'remov4 :',mark_selec 1538 index_mark_valide(mark_selec-lk+1) = .false. 1539 removeCas4 = removeCas4 -1 1540 end do 1541 1542 do while ( removeCas2>0 ) 1543 val_min = 0.d0 1544 mark_selec = -1 1545 do lit=lk,lmax 1546 if (index_mark_valide(lit-lk+1) .and. cas(lit)==2) then 1547 pos_m=abs(dx-posi(c,lit)) 1548 if (pos_m>=val_min) then 1549 val_min = pos_m 1550 mark_selec = lit 1551 end if 1552 end if 1553 end do 1554 1555 if (mark_selec == -1) call stop_application("ERROR"); 1556 ! print *,'remov2 :',mark_selec 1557 index_mark_valide(mark_selec-lk+1) = .false. 1558 removeCas2 = removeCas2 -1 1559 end do 1560 1561 1562 lkkk=1 1563 do lq=lk,lmax 1564 if ( index_mark_valide(lq-lk+1)) then 1565 index_mark(lkkk)=lq 1566 !print *,lkkk,':',lq,':valide pos:',posi(c,lq) 1567 lkkk=lkkk+1 1568 !else 1569 !print *,lq,':non valide pos:',posi(c,lq) 1570 end if 1571 end do 1572 !stop 1573 1574 lk=1 1575 lmax=lkkk-1 1576 1577 ! 1578 ! nbt est le nombre d'evenements de transmissions possibles entre le premier marqueur 1579 ! et le premier MIF 1580 ! pbt(lq,it) est la probabilite du it emme de ces evenements (PBT), quand la transmission 1581 ! au qtl est lq (1 e 4, correspondant aux 4 transmission grand parentales) 1582 ! 1583 ! on etablit ces donnees iterativement en partant de l'extremite du chromosome 1584 ! 1585 ! Premier marqueur 1586 ! 1587 1588 nbt=cas(index_mark(lk)) 1589 1590 do it=1,nbt 1591 last(it,1)=new(trans(index_mark(lk)),2*(it-1)+1) 1592 last(it,2)=new(trans(index_mark(lk)),2*(it-1)+2) 1593 do lq=1,4 1594 pbt(lq,it)=1.d0/dble(cas(index_mark(lk))) 1595 end do 1596 end do 1597 1598 ! 1599 ! marqueurs suivants 1600 ! 1601 110 lk=lk+1 1602 1603 if(lk.gt.lmax) go to 120 1604 ! 1605 ! le qtl est entre lk-1 et lk, il faut considerer les 4 evenements de tranmissions 1606 ! au qtl 1607 ! 1608 im1 = index_mark(lk-1) 1609 im2 = index_mark(lk) 1610 if(dx.ge.posi(c,im1).and. dx.le.posi(c,im2)) then 1611 1612 xint=posi(c,index_mark(lk))-posi(c,index_mark(lk-1)) 1613 dg=dx-posi(c,index_mark(lk-1)) 1614 dd=posi(c,index_mark(lk))-dx 1615 xintm=posim(c,index_mark(lk))-posim(c,index_mark(lk-1)) 1616 dgm=dg*xintm/xint 1617 ddm=dd*xintm/xint 1618 xintf=posif(c,index_mark(lk))-posif(c,index_mark(lk-1)) 1619 dgf=dg*xintf/xint 1620 ddf=dd*xintf/xint 1621 recmlq=xaldane(dgm) 1622 recflq=xaldane(dgf) 1623 recmqr=xaldane(ddm) 1624 recfqr=xaldane(ddf) 1625 1626 do lq=1,4 1627 do it=1,nbt 1628 if(new(11,1+2*(lq-1)).eq.last(it,1)) then 1629 pbt(lq,it)=pbt(lq,it)*(1.d0-recmlq) 1630 else 1631 pbt(lq,it)=pbt(lq,it)*recmlq 1632 end if 1633 1634 if(new(11,2+2*(lq-1)).eq.last(it,2)) then 1635 pbt(lq,it)=pbt(lq,it)*(1.d0-recflq) 1636 else 1637 pbt(lq,it)=pbt(lq,it)*recflq 1638 end if 1639 end do 1640 end do 1641 1642 do lq=1,4 1643 do it=1,nbt 1644 last(it+nbt*(lq-1),1)=new(11,1+2*(lq-1)) 1645 last(it+nbt*(lq-1),2)=new(11,2+2*(lq-1)) 1646 end do 1647 end do 1648 1649 ! 1650 ! PBTQ est incremente pour les evenements de transmision entre le QTL et le marqueur lk 1651 ! 1652 1653 do lq=1,4 1654 1655 do jt=1,cas(index_mark(lk))-1 1656 do it=1,nbt 1657 pbt(lq,it+nbt*jt)=pbt(lq,it) 1658 end do 1659 end do 1660 1661 do jt=1,cas(index_mark(lk)) 1662 do it=1,nbt 1663 kt=it+nbt*(jt-1) 1664 lt=it+nbt*(lq-1) 1665 if(new(trans(index_mark(lk)),1+2*(jt-1)).eq.last(lt,1)) then 1666 pbt(lq,kt)=pbt(lq,kt)*(1.d0-recmqr) 1667 else 1668 pbt(lq,kt)=pbt(lq,kt)*recmqr 1669 end if 1670 1671 if(new(trans(index_mark(lk)),2+2*(jt-1)).eq.last(lt,2)) then 1672 pbt(lq,kt)=pbt(lq,kt)*(1.d0-recfqr) 1673 else 1674 pbt(lq,kt)=pbt(lq,kt)*recfqr 1675 end if 1676 end do 1677 end do 1678 1679 end do 1680 1681 do jt=1,cas(index_mark(lk)) 1682 do it=1,nbt 1683 kt=it+nbt*(jt-1) 1684 lt=it+nbt*(lq-1) 1685 last(kt,1)=new(trans(index_mark(lk)),1+2*(jt-1)) 1686 last(kt,2)=new(trans(index_mark(lk)),2+2*(jt-1)) 1687 end do 1688 end do 1689 1690 nbt=nbt*cas(index_mark(lk)) 1691 end if 1692 1693 ! 1694 ! le qtl n'est pas entre lk-1 et lk 1695 ! 1696 posi_lkinf = posi(c,index_mark(lk-1)) 1697 posi_lk = posi(c,index_mark(lk)) 1698 if(.not.(dx.ge.posi_lkinf.and.dx.le.posi_lk)) then 1699 ! 1700 ! incrementation de la probabilites des transmission de marqueurs PBT 1701 ! 1702 1703 recm=xaldane(posim(c,index_mark(lk))-posim(c,index_mark(lk-1))) 1704 recf=xaldane(posif(c,index_mark(lk))-posif(c,index_mark(lk-1))) 1705 1706 do lq=1,4 1707 do jt=1,cas(index_mark(lk))-1 1708 do it=1,nbt 1709 pbt(lq,it+nbt*jt)=pbt(lq,it) 1710 end do 1711 end do 1712 1713 do jt=1,cas(index_mark(lk)) 1714 1715 do it=1,nbt 1716 kt=it+nbt*(jt-1) 1717 if(new(trans(index_mark(lk)),1+2*(jt-1)).eq.last(it,1)) then 1718 pbt(lq,kt)=pbt(lq,kt)*(1.d0-recm) 1719 else 1720 pbt(lq,kt)=pbt(lq,kt)*recm 1721 end if 1722 1723 if(new(trans(index_mark(lk)),2+2*(jt-1)).eq.last(it,2)) then 1724 pbt(lq,kt)=pbt(lq,kt)*(1.d0-recf) 1725 else 1726 pbt(lq,kt)=pbt(lq,kt)*recf 1727 end if 1728 end do 1729 1730 end do 1731 end do 1732 1733 do jt=1,cas(index_mark(lk)) 1734 do it=1,nbt 1735 kt=it+nbt*(jt-1) 1736 !lt=it+nbt*(lq-1) 1737 last(kt,1)=new(trans(index_mark(lk)),1+2*(jt-1)) 1738 last(kt,2)=new(trans(index_mark(lk)),2+2*(jt-1)) 1739 end do 1740 end do 1741 1742 nbt=nbt*cas(index_mark(lk)) 1743 end if 1744 1745 go to 110 1746 1747 120 continue 1748 ! 1749 ! on stocke la probabilite de transmission e la position dx 1750 ! 1751 pbtt=0.d0 1752 do lq=1,4 1753 pdd(c,kd,lq,n)=0.d0 1754 do kt=1,nbt 1755 pbtt=pbtt+pbt(lq,kt) 1756 pdd(c,kd,lq,n)=pdd(c,kd,lq,n)+pbt(lq,kt) 1757 end do 1758 end do 1759 if (pbtt /= 0) then 1760 do lq=1,4 1761 pdd(c,kd,lq,n)=pdd(c,kd,lq,n)/pbtt 1762 end do 1763 end if 1764 1765 125 continue 1766 ! 1767 ! on avance d'un pas (en verifiant qu'on ne depasse pas le bout droit) et on recommence 1768 ! 1769 n=n+1 1770 npas=n*pas 1771 if(npas.gt.ilong+1) go to 1000 1772 1773 if(linf.le.nmkinf) go to 100 1774 1775 ! 1776 ! bout du chromosome 1777 ! 1778 1000 continue 1779 999 continue 1780 998 continue 1781 997 continue 1782 end do 1783 1784 if (allocated(last)) deallocate(last) 1785 if (allocated(pbt)) deallocate(pbt) 1786 call log_mess('end PDED_V4',DEBUG_DEF) 1787 1788 end subroutine pded
pded_v5
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pded_v5
DESCRIPTION
NOTES
Calcul des probabilites de transmission le long du chromosome
SOURCE
2377 subroutine pded_v5() 2378 2379 ! 2380 ! Tableaux dimensionnes selon l, le nombre de marqueurs 2381 integer trans(maxval(nmk)) 2382 integer indp(maxval(nmk)),indm(maxval(nmk)),ingp(maxval(nmk)),ingm(maxval(nmk)) 2383 integer indi(maxval(nmk)),ingi(maxval(nmk)) 2384 integer hdp(maxval(nmk)),hdm(maxval(nmk)),hgp(maxval(nmk)),hgm(maxval(nmk)) 2385 integer nbdp(maxval(nmk)),nbdm(maxval(nmk)) 2386 integer nbgp(maxval(nmk)),nbgm(maxval(nmk)) 2387 integer transp(maxval(nmk)),inforp(maxval(nmk)) 2388 2389 double precision trecm(maxval(nmk),maxval(nmk)),trecf(maxval(nmk),maxval(nmk)) 2390 double precision srecm(maxval(nmk),maxval(nmk)),srecf(maxval(nmk),maxval(nmk)) 2391 double precision pbt(4),pbtg(4),pbtd(4),pcum(2),pbante(4,2) 2392 2393 logical infor(maxval(nmk)) 2394 2395 ! Divers 2396 integer new(11,8) 2397 integer geno 2398 integer hdi,hdifin,hdideb,hgi,hgideb,hgifin,hlocp,hlocm 2399 integer idm,igm,idp,igp 2400 integer icoloc,ii,icas,idi,igi,igeno,ilong,ip,iqtl,iqtlp,iqtlm 2401 integer j,jdi,jgi,jj,jk,jm,jn,jnd,jma 2402 integer kcas,kd,kkd 2403 integer lecas,linf,linfp,lk,llk,llkp,lma,lma1,lma2,ln,lq,lqtl 2404 integer maxcas 2405 integer n,nbdi,nbgi 2406 integer nd1,nd2,ndeb,nfin,ngeno1,ngeno2,nm1,nm2,c 2407 integer nmkinf,npas,nposi,nmkinfp 2408 2409 integer kcutopt,temps,tempsmax 2410 integer kcutd,lcutd,ncutd,rcutd,kcutg,lcutg,ncutg,rcutg 2411 2412 double precision apm,app 2413 double precision ddp,ddm,dgp,dgm,dlq,dqr 2414 double precision dx,pbtt,pm,pp,prob 2415 double precision recm,recf,reclq,recqr 2416 double precision recdp,recdm,trecdp,trecdm,srecdp,srecdm 2417 double precision recgp,recgm,trecgp,trecgm,srecgp,srecgm 2418 double precision xint,xintp,xintm 2419 logical coloc,colocdp,colocdm,colocgp,colocgm,finp,finm,possible 2420 ! 2421 ! tableau des transmission 2422 ! 2423 data new / 1,1,2,2,1,1,1,2,1,1,1,& 2424 1,2,1,2,1,2,1,1,1,2,1,& 2425 0,0,0,0,2,2,1,2,2,2,1,& 2426 0,0,0,0,2,1,2,2,1,2,2,& 2427 0,0,0,0,0,0,0,0,0,0,2,& 2428 0,0,0,0,0,0,0,0,0,0,1,& 2429 0,0,0,0,0,0,0,0,0,0,2,& 2430 0,0,0,0,0,0,0,0,0,0,2/ 2431 2432 do c=1,nchr 2433 ! 2434 ! Taille du segment explore 2435 ilong=get_ilong(c) 2436 ! 2437 ! tableau des recombinaison 2438 ! 2439 do lk=1,nmk(c)-1 2440 do jk=lk+1,nmk(c) 2441 recm=xaldane(posim(c,jk)-posim(c,lk)) 2442 recf=xaldane(posif(c,jk)-posif(c,lk)) 2443 trecm(lk,jk)=dlog(recm) 2444 trecf(lk,jk)=dlog(recf) 2445 srecm(lk,jk)=dlog(1.d0-recm) 2446 srecf(lk,jk)=dlog(1.d0-recf) 2447 trecm(jk,lk)=trecm(lk,jk) 2448 trecf(jk,lk)=trecf(lk,jk) 2449 srecm(jk,lk)=srecm(lk,jk) 2450 srecf(jk,lk)=srecf(lk,jk) 2451 end do 2452 end do 2453 ! 2454 ! Calcul pour chaque descendant de la probabilite de transmission 2455 do 1000 ip=1,np 2456 nm1=nmp(ip)+1 2457 nm2=nmp(ip+1) 2458 do 1000 jm=nm1,nm2 2459 ngeno1=ngenom(c,jm)+1 2460 ngeno2=ngenom(c,jm+1) 2461 do 1000 geno=ngeno1,ngeno2 2462 nd1=ngend(c,geno)+1 2463 nd2=ngend(c,geno+1) 2464 do 1000 kd=nd1,nd2 2465 kkd=ndesc(c,kd) 2466 ! 2467 ! le vecteur trans(lk) indique l'information jointe sur la transmission au locus lk 2468 ! le vecteur infor contient la liste des marqeurs doublement informatifs 2469 ! 2470 ! l'origine grand parentale des all�les recu par le descendant k au marqueur lk 2471 ! d�pend de trans (tableau new) 2472 ! soit 1 pour grand p�re et 2 pour grand m�re 2473 ! 2474 ! cas 1 cas 2 cas 3 cas 4 2475 ! trans(lk) pere mere pere mere pere mere pere mere 2476 ! 1 1 1 2477 ! 2 1 2 2478 ! 3 2 1 2479 ! 4 2 2 2480 ! 5 1 1 2 2 2481 ! 6 1 2 2 1 2482 ! 7 1 1 1 2 2483 ! 8 2 1 2 2 2484 ! 9 1 1 2 1 2485 ! 10 1 2 2 2 2486 ! 11 1 1 1 2 2 1 2 2 2487 ! 2488 ! 2489 2490 llkp=0 2491 do lk=1,nmk(c) 2492 2493 ! 2494 ! infor(lk) est � true si le marqueur lk est informatif pour les deux parents 2495 ! 2496 infor(lk)=.false. 2497 ! 2498 ! incertitude totale (donnees manquantes) 2499 ! 2500 trans(lk)=11 2501 ! 2502 ! triplet h�t�rozygotes 2503 ! 2504 if(ptfin(c,lk,kd,2).eq.0.and.ptfin(c,lk,kd,3).eq.0)trans(lk)=5 2505 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=6 2506 ! 2507 ! incertitude sur un parent 2508 ! 2509 if(ptfin(c,lk,kd,3).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=7 2510 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,2).eq.0)trans(lk)=8 2511 if(ptfin(c,lk,kd,2).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=9 2512 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,3).eq.0)trans(lk)=10 2513 ! 2514 ! cas certains 2515 ! 2516 if(ptfin(c,lk,kd,1).eq.1.d0) trans(lk)=1 2517 if(ptfin(c,lk,kd,2).eq.1.d0) trans(lk)=2 2518 if(ptfin(c,lk,kd,3).eq.1.d0) trans(lk)=3 2519 if(ptfin(c,lk,kd,4).eq.1.d0) trans(lk)=4 2520 2521 if(trans(lk).le.4) infor(lk)=.true. 2522 ! 2523 ! 2524 ! cas des familles de demi fr�res 2525 ! 2526 !if(force(jm).or.opt_sib.eq.OPT_SIB_HS) then 2527 if(.not.estfem(repfem(jm)).or.opt_sib.eq.OPT_SIB_HS)then 2528 transp(lk)=0 2529 if(trans(lk).le.2.or.trans(lk).eq.7)transp(lk)=1 2530 if(trans(lk).eq.3.or.trans(lk).eq.4.or.trans(lk).eq.8)transp(lk)=2 2531 if(transp(lk).ne.0) then 2532 llkp=llkp+1 2533 inforp(llkp)=lk 2534 end if 2535 end if 2536 end do 2537 nmkinfp=llkp 2538 ! 2539 ! exploration du groupe de liaison pour la recherche des transmissions possibles 2540 ! 2541 n=1 2542 100 dx=absi(c,n) 2543 ! 2544 ! pour des probl�mes d'arrondis, la position test�e pour la qtl 2545 ! peut d�passer le bout du groupe de liaison... on l'y ramene 2546 ! 2547 if(dx.ge.posi(c,nmk(c)))dx=posi(c,nmk(c)) 2548 ! 2549 ! cas des demi fr�res 2550 ! 2551 ! if(force(jm).or.opt_sib.eq.OPT_SIB_HS) then 2552 if(.not.estfem(repfem(jm)).or.opt_sib.eq.OPT_SIB_HS) then 2553 do ii=1,4 2554 pbt(ii)=0.d0 2555 end do 2556 ! 2557 ! cas des descendants sans information 2558 ! 2559 if(nmkinfp.eq.0) then 2560 do ii=1,3,2 2561 pbt(ii)=0.5d0 2562 end do 2563 go to 130 2564 end if 2565 ! 2566 ! descendants avec information 2567 ! 2568 linfp=1 2569 ! 2570 ! position du qtl avant le premier marqueur informatif 2571 ! 2572 if(dx.lt.posi(c,inforp(1)))then 2573 dqr=posi(c,inforp(1))-dx 2574 if(inforp(1).eq.1) then 2575 if(posi(c,1).ne.0.)dqr=dqr*posim(c,1)/posi(c,1) 2576 else 2577 xint=posi(c,inforp(1))-posi(c,inforp(1)-1) 2578 xintm=posim(c,inforp(1))-posim(c,inforp(1)-1) 2579 if(xint.ne.0.d0)dqr=dqr*xintm/xint 2580 end if 2581 recqr=xaldane(dqr) 2582 2583 if(transp(inforp(1)).eq.1) then 2584 pbt(1)=1.d0-recqr 2585 pbt(3)=recqr 2586 else 2587 pbt(1)=recqr 2588 pbt(3)=1.d0-recqr 2589 end if 2590 2591 go to 130 2592 end if 2593 2594 5 linfp=linfp+1 2595 if(linfp.gt.nmkinfp) then 2596 ! 2597 ! qtl apr�s le dernier marqueur informatif 2598 ! 2599 dlq=dx-posi(c,inforp(nmkinfp)) 2600 if(inforp(nmkinfp).eq.nmk(c)) then 2601 xint=posi(c,inforp(nmkinfp))-posi(c,inforp(nmkinfp)-1) 2602 xintm=posim(c,inforp(nmkinfp))-posim(c,inforp(nmkinfp)-1) 2603 if(xint.ne.0.d0) dlq=dlq*xintm/xint 2604 else 2605 xint=posi(c,inforp(nmkinfp)+1)-posi(c,inforp(nmkinfp)) 2606 xintm=posim(c,inforp(nmkinfp)+1)-posim(c,inforp(nmkinfp)) 2607 if(xint.ne.0.d0) dlq=dlq*xintm/xint 2608 end if 2609 reclq=xaldane(dlq) 2610 if(transp(inforp(nmkinfp)).eq.1) then 2611 pbt(1)=1.d0-reclq 2612 pbt(3)=reclq 2613 else 2614 pbt(1)=reclq 2615 pbt(3)=1.d0-reclq 2616 end if 2617 2618 go to 130 2619 2620 else 2621 2622 if(.not.(dx.ge.posi(c,inforp(linfp-1)) .and. dx.le.posi(c,inforp(linfp)))) go to 5 2623 ! 2624 ! qtl entre deux marqueurs informatifs 2625 ! 2626 2627 dqr=posi(c,inforp(linfp))-dx 2628 dlq=dx-posi(c,inforp(linfp-1)) 2629 xint=posi(c,inforp(linfp))-posi(c,inforp(linfp-1)) 2630 xintm=posim(c,inforp(linfp))-posim(c,inforp(linfp-1)) 2631 if(xint.ne.0.d0)then 2632 dqr=dqr*xintm/xint 2633 recqr=xaldane(dqr) 2634 dlq=dlq*xintm/xint 2635 reclq=xaldane(dlq) 2636 end if 2637 2638 if(transp(inforp(linfp-1)).eq.transp(inforp(linfp)))then 2639 if(transp(inforp(linfp)).eq.1)then 2640 pbt(1)=(1.d0-reclq)*(1.d0-recqr) 2641 pbt(3)=reclq*recqr 2642 else 2643 pbt(1)=reclq*recqr 2644 pbt(3)=(1.d0-reclq)*(1.d0-recqr) 2645 end if 2646 else 2647 if(transp(inforp(linfp)).eq.1)then 2648 pbt(1)=reclq*(1.d0-recqr) 2649 pbt(3)=(1.d0-reclq)*recqr 2650 else 2651 pbt(1)=(1.d0-reclq)*recqr 2652 pbt(3)=reclq*(1.d0-recqr) 2653 end if 2654 end if 2655 2656 go to 130 2657 2658 end if 2659 end if 2660 ! 2661 ! cas des m�lange plein / demi fr�res 2662 ! 2663 ! 2664 ! cas des positions sur des marqueurs informatifs 2665 ! 2666 ! 2667 ! si (coloc=true) le QTL est sur exactement sur le marqueur lma (icoloc =lma) 2668 ! le cas sera trait� sp�cifiquement pour �viter de prendre Log(z�ro) 2669 ! si de surcroit le marqueur est informatif pour les deux parents, on connait 2670 ! directement la probabilit� de transmission 2671 ! 2672 coloc=.false. 2673 do lk=1,nmk(c) 2674 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,lk)))then 2675 coloc=.true. 2676 icoloc=lk 2677 if(infor(lk))then 2678 do lq=1,4 2679 pbt(lq)=0.d0 2680 end do 2681 pbt(trans(lk))=1.d0 2682 go to 130 2683 end if 2684 end if 2685 end do 2686 ! 2687 ! localisation de la position test�e du qtl 2688 ! 2689 do lk=2,nmk(c) 2690 if(dx.ge.posi(c,lk-1).and.dx.lt.posi(c,lk)) lqtl=lk-1 2691 end do 2692 if(dx.eq.posi(c,nmk(c))) lqtl=nmk(c)-1 2693 ! 2694 ! Recherche des zones informatives pour la transmission 2695 ! La zone du p�re (resp. de la m�re) doit se terminer � gauche et � droite 2696 ! par des marqueurs informatifs (dont l'all�le transmis au descendant est 2697 ! connu), sans la pr�sence de marqueurs incertain (trans = 5 ou 6) entre 2698 ! les bornes � gauche (resp � droite ) des deux parents 2699 ! 2700 ! nbdp est le nombre de marqueurs connu pour le pere � droite du qtl 2701 ! nbdm, pour la mere 2702 ! nbdi, nombre de marqueurs incertains 2703 ! 2704 ! hdp est le vecteur des transmision certaines pour le pere � droite 2705 ! hdm, pour la mere 2706 ! 2707 ! indp est le vecteur des indices (dans la lise initiale) des marqueurs 2708 ! connus pour le pere � droite 2709 ! indm, pour la mere 2710 ! indi, le tableau des indices des marqueurs incertains (colonne 1 : le pere, 2 2711 ! : la mere) 2712 ! 2713 nbdi=1 2714 do lk=1,nmk(c) 2715 nbdp(lk)=0 2716 nbdm(lk)=0 2717 end do 2718 jma=lqtl 2719 finp=.false. 2720 finm=.false. 2721 2722 10 jma=jma+1 2723 if(jma.gt.nmk(c)) go to 15 2724 if(trans(jma).le.4) then 2725 nbdp(nbdi)=nbdp(nbdi)+1 2726 indp(nbdp(nbdi))=jma 2727 hdp(nbdp(nbdi))=floor(real(trans(jma)+1)/2.) 2728 nbdm(nbdi)=nbdm(nbdi)+1 2729 indm(nbdm(nbdi))=jma 2730 hdm(nbdm(nbdi))=trans(jma)-2*(hdp(nbdp(nbdi))-1) 2731 finp=.true. 2732 finm=.true. 2733 end if 2734 2735 if(trans(jma).eq.5.or.trans(jma).eq.6) then 2736 nbdp(nbdi)=nbdp(nbdi)+1 2737 indp(nbdp(nbdi))=jma 2738 hdp(nbdp(nbdi))=0 2739 nbdm(nbdi)=nbdm(nbdi)+1 2740 indm(nbdm(nbdi))=jma 2741 hdm(nbdm(nbdi))=0 2742 indi(nbdi)=jma 2743 nbdi=nbdi+1 2744 nbdp(nbdi)=nbdp(nbdi-1) 2745 nbdm(nbdi)=nbdm(nbdi-1) 2746 finp=.false. 2747 finm=.false. 2748 end if 2749 2750 if(trans(jma).eq.7.or.trans(jma).eq.8) then 2751 nbdp(nbdi)=nbdp(nbdi)+1 2752 indp(nbdp(nbdi))=jma 2753 hdp(nbdp(nbdi))=trans(jma)-6 2754 finp=.true. 2755 end if 2756 2757 2758 if(trans(jma).eq.9.or.trans(jma).eq.10) then 2759 nbdm(nbdi)=nbdm(nbdi)+1 2760 indm(nbdm(nbdi))=jma 2761 hdm(nbdm(nbdi))=trans(jma)-8 2762 finm=.true. 2763 end if 2764 2765 if(((.not.finp).or.(.not.finm)).and.jma.lt.nmk(c)) go to 10 2766 2767 ! 2768 ! on cherche les zones � gauche du qtl 2769 ! les elements sont comme � droite , en remplacant d par g 2770 ! 2771 15 nbgi=1 2772 do lk=1,nmk(c) 2773 nbgp(lk)=0 2774 nbgm(lk)=0 2775 end do 2776 jma=lqtl+1 2777 finp=.false. 2778 finm=.false. 2779 2780 20 jma=jma-1 2781 if(jma.eq.0) go to 25 2782 if(trans(jma).le.4) then 2783 nbgp(nbgi)=nbgp(nbgi)+1 2784 ingp(nbgp(nbgi))=jma 2785 hgp(nbgp(nbgi))=floor(real(trans(jma)+1)/2.) 2786 nbgm(nbgi)=nbgm(nbgi)+1 2787 ingm(nbgm(nbgi))=jma 2788 hgm(nbgm(nbgi))=trans(jma)-2*(hgp(nbgp(nbgi))-1) 2789 finp=.true. 2790 finm=.true. 2791 end if 2792 2793 if(trans(jma).eq.5.or.trans(jma).eq.6) then 2794 nbgp(nbgi)=nbgp(nbgi)+1 2795 ingp(nbgp(nbgi))=jma 2796 hgp(nbgp(nbgi))=0 2797 nbgm(nbgi)=nbgm(nbgi)+1 2798 ingm(nbgm(nbgi))=jma 2799 hgm(nbgm(nbgi))=0 2800 ingi(nbgi)=jma 2801 nbgi=nbgi+1 2802 nbgp(nbgi)=nbgp(nbgi-1) 2803 nbgm(nbgi)=nbgm(nbgi-1) 2804 finp=.false. 2805 finm=.false. 2806 end if 2807 2808 if(trans(jma).eq.7.or.trans(jma).eq.8) then 2809 nbgp(nbgi)=nbgp(nbgi)+1 2810 ingp(nbgp(nbgi))=jma 2811 hgp(nbgp(nbgi))=trans(jma)-6 2812 finp=.true. 2813 end if 2814 2815 2816 if(trans(jma).eq.9.or.trans(jma).eq.10) then 2817 nbgm(nbgi)=nbgm(nbgi)+1 2818 ingm(nbgm(nbgi))=jma 2819 hgm(nbgm(nbgi))=trans(jma)-8 2820 finm=.true. 2821 end if 2822 2823 if(((.not.finp).or.(.not.finm)).and.jma.gt.1) go to 20 2824 2825 ! 2826 ! colocgp et colocgm indique (true) si la position du qtl est colocalis�e 2827 ! avec un marqueur (coloc = true) et que ce marqueur est le premier � 2828 ! sa gauche. Dans ce cas, il n'y a pas de recombinaison possible entre 2829 ! le qtl et ce marqueur 2830 ! 2831 ! hlocp et hlocm sont les informations sur la transmission en ce locus (� la 2832 ! fois marqueur et qtl) 2833 ! 2834 25 colocgp=.false. 2835 colocgm=.false. 2836 colocdp=.false. 2837 colocdm=.false. 2838 if(coloc.and.nbgp(1).ne.0)then 2839 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,ingp(1))))colocgp=.true. 2840 end if 2841 if(coloc.and.nbgm(1).ne.0)then 2842 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,ingm(1))))colocgm=.true. 2843 end if 2844 if(coloc.and.nbdp(1).ne.0)then 2845 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,indp(1))))colocdp=.true. 2846 end if 2847 if(coloc.and.nbdm(1).ne.0)then 2848 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,indm(1))))colocdm=.true. 2849 end if 2850 2851 ! 2852 ! les taux de recombinaison avant et apres le QTL sont stock�s dans 2853 ! (t ou s)rec(d ou g)(p ou m) 2854 ! 2855 idp=indp(1) 2856 if(nbdp(1).eq.0)idp=lqtl+1 2857 igp=ingp(1) 2858 if(nbgp(1).eq.0)igp=lqtl 2859 2860 xint=posi(c,lqtl+1)-posi(c,lqtl) 2861 xintp=posim(c,lqtl+1)-posim(c,lqtl) 2862 ! 2863 ! dans le cas de colocalisation, il ne faut pas calculer le taux de 2864 ! recombinaison 2865 ! 2866 dgp=(dx-posi(c,igp))*xintp/xint 2867 recgp=xaldane(dgp) 2868 ddp=(posi(c,idp)-dx)*xintp/xint 2869 recdp=xaldane(ddp) 2870 ! 2871 ! 2872 ! cas des m�res 2873 ! 2874 idm=indm(1) 2875 if(nbdm(1).eq.0)idm=lqtl+1 2876 igm=ingm(1) 2877 if(nbgm(1).eq.0)igm=lqtl 2878 2879 xint=posi(c,lqtl+1)-posi(c,lqtl) 2880 xintm=posif(c,lqtl+1)-posif(c,lqtl) 2881 ! 2882 ! dans le cas de colocalisation, il ne faut pas calculer le taux de 2883 ! recombinaison 2884 ! 2885 dgm=(dx-posi(c,igm))*xintm/xint 2886 recgm=xaldane(dgm) 2887 ddm=(posi(c,idm)-dx)*xintm/xint 2888 recdm=xaldane(ddm) 2889 2890 ! 2891 ! calcul des probabilit�s de transmission pour le descendant kd, position dx 2892 ! 2893 ! 2894 ! on va envisager toutes les possibilit�s pour les marqueurs incertains 2895 ! il y a 2**(nbdi+nbgi) possibilit�s 2896 ! ainsi que les 4 possibilit�s pour la transmission du QTL 2897 ! 2898 ! les hdp, hdm, hgp et hgm sont compl�t�s par des 1 ou 2 pour les marqueurs 2899 ! incertains 2900 ! 2901 ! les boucles en jdi et jgi ont pour fonction de g�n�rer les 2**(nbdi+nbgi) 2902 ! vecteurs d'origines correspondants aux marqueurs incertains 2903 ! 2904 ! initialisation 2905 pbtt=0.d0 2906 do lq=1,4 2907 pbtg(lq)=0.d0 2908 pbtd(lq)=0.d0 2909 end do 2910 2911 ! 2912 ! on va consid�rer les 4 �v�nements de transmission possible au QTL, d�finis 2913 ! pas iqtlp et iqtlm 2914 ! 2915 ! pour chacun des �venements de transmision au qtl on v�rifie la coh�rence 2916 ! avec un �ventuel marqueur colocalis� 2917 ! 2918 iqtlp=0 2919 30 iqtlp=iqtlp+1 2920 if(iqtlp.eq.3)go to 120 2921 2922 iqtlm=0 2923 40 iqtlm=iqtlm+1 2924 if(iqtlm.eq.3)go to 30 2925 ! 2926 iqtl= iqtlm+2*(iqtlp-1) 2927 2928 ! 2929 ! le cas o� il n'y a pas d'incertitude � droite est trait� � part 2930 ! dans ce cas, il faut distinguer les situations o� il n'y a que des 2931 ! transmissions inconnues 2932 ! 2933 if(nbdi.eq.1) then 2934 ! 2935 ! on a nbdp=0 si aucun marqueur � droite n'a trans diff�rent de 9, 10 ou 11 2936 ! on a nbdm=0 si aucun marqueur � droite n'a trans diff�rent de 7, 8 ou 11 2937 ! dans ce cas,les 2 �v�nements iqtlp = 1 ou 2 sont �quiprobables 2938 ! 2939 if(nbdp(1).eq.0.and.nbdm(1).eq.0) then 2940 pbtd(iqtl)=0.25d0 2941 end if 2942 2943 prob=1.d0 2944 if(nbdm(1).ne.0) then 2945 if(colocdm.and.iqtlm.ne.hdm(1))go to 50 2946 if(.not.colocdm)then 2947 if(nbdp(1).eq.0) prob=0.5d0 2948 if(iqtlm.eq.hdm(1))prob=prob*(1.d0-recdm) 2949 if(iqtlm.ne.hdm(1))prob=prob*recdm 2950 end if 2951 end if 2952 2953 if(nbdp(1).ne.0) then 2954 if(colocdp.and.iqtlp.ne.hdp(1))go to 50 2955 if(.not.colocdp)then 2956 if(nbdm(1).eq.0) prob=0.5d0 2957 if(iqtlp.eq.hdp(1))prob=prob*(1.d0-recdp) 2958 if(iqtlp.ne.hdp(1))prob=prob*recdp 2959 end if 2960 end if 2961 2962 pbtd(iqtl)=prob 2963 go to 50 2964 2965 end if 2966 ! 2967 ! cas ou il y a des incertitudes � droite 2968 ! 2969 ! 2970 ! 2971 ! on traite diff�remment les intervalles entre marqueurs incertains : le 1er, les 2972 ! suivants, le dernier 2973 ! 2974 ! cas du premier intervalle 2975 ! nfin est l'indice du marqueur incertain � droite 2976 ! 2977 idi=1 2978 possible=.false. 2979 do hdifin=1,2 2980 ! 2981 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 2982 ! hdp et hdm du premier intervalle 2983 ! 2984 nfin=indi(idi) 2985 hdp(nbdp(1))=new(trans(nfin),2*hdifin-1) 2986 hdm(nbdm(1))=new(trans(nfin),2*hdifin) 2987 2988 if((colocdm.and.iqtlm.ne.hdm(1)).or.(colocdp.and.iqtlp.ne.hdp(1)))then 2989 2990 pbante(iqtl,hdifin)=0.d0 2991 2992 else 2993 2994 possible=.true. 2995 prob=0.d0 2996 2997 if(.not.colocdm)then 2998 if(iqtlm.eq.hdm(1))prob=dlog(1.d0-recdm) 2999 if(recdm.eq.0) print*,'kd,n',kd,n,dx,ndesc(c,kkd),trim(animal(kkd)) 3000 if(iqtlm.ne.hdm(1))prob=dlog(recdm ) 3001 3002 do jma=2,nbdm(1) 3003 if(hdm(jma-1).eq.hdm(jma))then 3004 prob=prob+srecf(indm(jma-1),indm(jma)) 3005 else 3006 prob=prob+trecf(indm(jma-1),indm(jma)) 3007 end if 3008 end do 3009 end if 3010 3011 if(.not.colocdp)then 3012 if(iqtlp.eq.hdp(1))prob=prob+dlog(1.d0-recdp) 3013 if(iqtlp.ne.hdp(1))prob=prob+dlog(recdp) 3014 do jma=2,nbdp(1) 3015 if(hdp(jma-1).eq.hdp(jma))then 3016 prob=prob+srecm(indp(jma-1),indp(jma)) 3017 else 3018 prob=prob+trecm(indp(jma-1),indp(jma)) 3019 end if 3020 end do 3021 end if 3022 3023 pbante(iqtl,hdifin)=dexp(prob) 3024 3025 end if 3026 3027 end do 3028 3029 if(.not.possible) go to 50 3030 ! 3031 ! cas des intervalles suivants 3032 ! 3033 do idi=2,nbdi-1 3034 3035 do hdifin=1,2 3036 ! 3037 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 3038 ! hdp et hdm du idi �me intervalle 3039 ! 3040 nfin=indi(idi) 3041 hdp(nbdp(idi))=new(trans(nfin),2*hdifin-1) 3042 hdm(nbdm(idi))=new(trans(nfin),2*hdifin) 3043 pcum(hdifin)=0.d0 3044 3045 do hdideb=1,2 3046 ! 3047 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 3048 ! de l'intervalle pr�c�dent 3049 ! 3050 ndeb=indi(idi-1) 3051 hdp(nbdp(idi-1))=new(trans(ndeb),2*hdideb-1) 3052 hdm(nbdm(idi-1))=new(trans(ndeb),2*hdideb) 3053 3054 prob=0.d0 3055 do jma=nbdp(idi-1)+1,nbdp(idi) 3056 if(hdp(jma-1).eq.hdp(jma))then 3057 prob=prob+srecm(indp(jma-1),indp(jma)) 3058 else 3059 prob=prob+trecm(indp(jma-1),indp(jma)) 3060 end if 3061 end do 3062 3063 do jma=nbdm(idi-1)+1,nbdm(idi) 3064 if(hdm(jma-1).eq.hdm(jma))then 3065 prob=prob+srecf(indm(jma-1),indm(jma)) 3066 else 3067 prob=prob+trecf(indm(jma-1),indm(jma)) 3068 end if 3069 end do 3070 ! 3071 pcum(hdifin)=pcum(hdifin)+pbante(iqtl,hdideb)*dexp(prob) 3072 3073 end do 3074 end do 3075 3076 do hdifin=1,2 3077 pbante(iqtl,hdifin)=pcum(hdifin) 3078 end do 3079 3080 end do 3081 3082 ! 3083 ! cas du dernier intervalle 3084 ! 3085 idi=nbdi 3086 pcum(1)=0.d0 3087 do hdideb=1,2 3088 ! 3089 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 3090 ! de l'intervalle pr�c�dent 3091 ! 3092 ndeb=indi(idi-1) 3093 hdp(nbdp(idi-1))=new(trans(ndeb),2*hdideb-1) 3094 hdm(nbdm(idi-1))=new(trans(ndeb),2*hdideb) 3095 3096 prob=0.d0 3097 if(nbdp(idi).ne.nbdp(idi-1))then 3098 do jma=nbdp(idi-1)+1,nbdp(idi) 3099 if(hdp(jma-1).eq.hdp(jma))then 3100 prob=prob+srecm(indp(jma-1),indp(jma)) 3101 else 3102 prob=prob+trecm(indp(jma-1),indp(jma)) 3103 end if 3104 end do 3105 end if 3106 3107 if(nbdm(idi).ne.nbdm(idi-1))then 3108 do jma=nbdm(idi-1)+1,nbdm(idi) 3109 if(hdm(jma-1).eq.hdm(jma))then 3110 prob=prob+srecf(indm(jma-1),indm(jma)) 3111 else 3112 prob=prob+trecf(indm(jma-1),indm(jma)) 3113 end if 3114 end do 3115 end if 3116 ! 3117 pcum(1)=pcum(1)+pbante(iqtl,hdideb)*dexp(prob) 3118 3119 end do 3120 3121 pbtd(iqtl)=pcum(1) 3122 3123 50 continue 3124 ! 3125 ! on �num�re les possibilit�s pour les transmissions incertaines � gauche 3126 ! du qtl. 3127 ! 3128 ! 3129 ! le cas o� il n'y a pas d'incertitude � gauche est trait� � part 3130 ! dans ce cas, il faut distinguer les situations o� il n'y a que des 3131 ! transmissions inconnues 3132 ! 3133 3134 if(nbgi.eq.1) then 3135 ! 3136 ! on a nbgp=0 si aucun marqueur � gauche n'a trans diff�rent de 9, 10 ou 11 3137 ! on a nbgm=0 si aucun marqueur � gauche n'a trans diff�rent de 7, 8 ou 11 3138 ! dans ce cas,les 2 �v�nements iqtlp = 1 ou 2 sont �quiprobables 3139 ! 3140 if(nbgp(1).eq.0.and.nbgm(1).eq.0) then 3141 pbtg(iqtl)=0.25d0 3142 end if 3143 3144 prob=1.d0 3145 if(nbgm(1).ne.0) then 3146 if(colocgm.and.iqtlm.ne.hgm(1))go to 60 3147 if(.not.colocgm)then 3148 if(nbgp(1).eq.0) prob=0.5d0 3149 if(iqtlm.eq.hgm(1))prob=prob*(1.d0-recgm) 3150 if(iqtlm.ne.hgm(1))prob=prob*recgm 3151 end if 3152 end if 3153 3154 if(nbgp(1).ne.0) then 3155 if(colocgp.and.iqtlp.ne.hgp(1))go to 60 3156 if(.not.colocgp)then 3157 if(nbgm(1).eq.0) prob=0.5d0 3158 if(iqtlp.eq.hgp(1))prob=prob*(1.d0-recgp) 3159 if(iqtlp.ne.hgp(1))prob=prob*recgp 3160 end if 3161 end if 3162 3163 pbtg(iqtl)=prob 3164 go to 60 3165 3166 end if 3167 ! 3168 ! cas ou il y a des incertitudes � gauche 3169 ! 3170 ! 3171 ! 3172 ! on traite diff�remment les intervalles entre marqueurs incertains : le 1er, les 3173 ! suivants, le dernier 3174 ! 3175 ! cas du premier intervalle 3176 ! nfin est l'indice du marqueur incertain � gauche 3177 ! 3178 igi=1 3179 possible=.false. 3180 do hgifin=1,2 3181 3182 ! 3183 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 3184 ! hgp et hgm du premier intervalle 3185 ! 3186 nfin=ingi(igi) 3187 hgp(nbgp(1))=new(trans(nfin),2*hgifin-1) 3188 hgm(nbgm(1))=new(trans(nfin),2*hgifin) 3189 3190 if( (colocgm.and.iqtlm.ne.hgm(1)).or.(colocgp.and.iqtlp.ne.hgp(1)))then 3191 3192 pbante(iqtl,hgifin)=0.d0 3193 3194 else 3195 3196 possible=.true. 3197 prob=0.d0 3198 3199 if(.not.colocgm)then 3200 if(iqtlm.eq.hgm(1))prob=dlog(1.d0-recgm) 3201 if(iqtlm.ne.hgm(1))prob=dlog(recgm) 3202 3203 do jma=2,nbgm(1) 3204 if(hgm(jma-1).eq.hgm(jma))then 3205 prob=prob+srecf(ingm(jma-1),ingm(jma)) 3206 else 3207 prob=prob+trecf(ingm(jma-1),ingm(jma)) 3208 end if 3209 end do 3210 end if 3211 3212 if(.not.colocgp)then 3213 if(iqtlp.eq.hgp(1))prob=prob+dlog(1.d0-recgp) 3214 if(iqtlp.ne.hgp(1))prob=prob+dlog(recgp) 3215 do jma=2,nbgp(1) 3216 if(hgp(jma-1).eq.hgp(jma))then 3217 prob=prob+srecm(ingp(jma-1),ingp(jma)) 3218 else 3219 prob=prob+trecm(ingp(jma-1),ingp(jma)) 3220 end if 3221 end do 3222 end if 3223 3224 pbante(iqtl,hgifin)=dexp(prob) 3225 3226 end if 3227 3228 end do 3229 3230 if(.not.possible)go to 60 3231 ! 3232 ! cas des intervalles suivants 3233 ! 3234 do igi=2,nbgi-1 3235 3236 do hgifin=1,2 3237 ! 3238 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 3239 ! hgp et hgm du igi �me intervalle 3240 ! 3241 nfin=ingi(igi) 3242 hgp(nbgp(igi))=new(trans(nfin),2*hgifin-1) 3243 hgm(nbgm(igi))=new(trans(nfin),2*hgifin) 3244 pcum(hgifin)=0.d0 3245 3246 do hgideb=1,2 3247 ! 3248 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 3249 ! de l'intervalle pr�c�dent 3250 ! 3251 ndeb=ingi(igi-1) 3252 hgp(nbgp(igi-1))=new(trans(ndeb),2*hgideb-1) 3253 hgm(nbgm(igi-1))=new(trans(ndeb),2*hgideb) 3254 3255 prob=0.d0 3256 do jma=nbgp(igi-1)+1,nbgp(igi) 3257 if(hgp(jma-1).eq.hgp(jma))then 3258 prob=prob+srecm(ingp(jma-1),ingp(jma)) 3259 else 3260 prob=prob+trecm(ingp(jma-1),ingp(jma)) 3261 end if 3262 end do 3263 3264 do jma=nbgm(igi-1)+1,nbgm(igi) 3265 if(hgm(jma-1).eq.hgm(jma))then 3266 prob=prob+srecf(ingm(jma-1),ingm(jma)) 3267 else 3268 prob=prob+trecf(ingm(jma-1),ingm(jma)) 3269 end if 3270 end do 3271 ! 3272 pcum(hgifin)=pcum(hgifin)+pbante(iqtl,hgideb)*dexp(prob) 3273 3274 end do 3275 end do 3276 3277 do hgifin=1,2 3278 pbante(iqtl,hgifin)=pcum(hgifin) 3279 end do 3280 3281 end do 3282 3283 ! 3284 ! cas du dernier intervalle 3285 ! 3286 igi=nbgi 3287 pcum(1)=0.d0 3288 do hgideb=1,2 3289 ! 3290 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 3291 ! de l'intervalle pr�c�dent 3292 ! 3293 ndeb=ingi(igi-1) 3294 hgp(nbgp(igi-1))=new(trans(ndeb),2*hgideb-1) 3295 hgm(nbgm(igi-1))=new(trans(ndeb),2*hgideb) 3296 3297 prob=0.d0 3298 if(nbgp(igi).ne.nbgp(igi-1))then 3299 do jma=nbgp(igi-1)+1,nbgp(igi) 3300 if(hgp(jma-1).eq.hgp(jma))then 3301 prob=prob+srecm(ingp(jma-1),ingp(jma)) 3302 else 3303 prob=prob+trecm(ingp(jma-1),ingp(jma)) 3304 end if 3305 end do 3306 end if 3307 3308 if(nbgm(igi).ne.nbgm(igi-1))then 3309 do jma=nbgm(igi-1)+1,nbgm(igi) 3310 if(hgm(jma-1).eq.hgm(jma))then 3311 prob=prob+srecf(ingm(jma-1),ingm(jma)) 3312 else 3313 prob=prob+trecf(ingm(jma-1),ingm(jma)) 3314 end if 3315 end do 3316 end if 3317 ! 3318 pcum(1)=pcum(1)+pbante(iqtl,hgideb)*dexp(prob) 3319 3320 end do 3321 3322 pbtg(iqtl)= pcum(1) 3323 3324 60 continue 3325 ! 3326 ! on teste si d'autres combinaisons au qtl sont � �valuer 3327 ! 3328 if(iqtlp.eq.2.and.iqtlm.eq.2) go to 120 3329 if(iqtlm.eq.1)go to 40 3330 if(iqtlm.eq.2)go to 30 3331 3332 3333 120 continue 3334 3335 ! 3336 ! on stocke la probabilit� de transmission � la position dx 3337 ! 3338 do lq=1,4 3339 pbt(lq)=pbtg(lq)*pbtd(lq) 3340 end do 3341 3342 130 continue 3343 3344 pbtt=0.d0 3345 do lq=1,4 3346 pbtt=pbtt+pbt(lq) 3347 end do 3348 3349 do lq=1,4 3350 pdd(c,kd,lq,n)=pbt(lq)/pbtt 3351 end do 3352 3353 ! 3354 ! 3355 ! on avance d'un pas (en verifiant qu'on ne depasse pas le bout droit) et on recommence 3356 ! 3357 n=n+1 3358 if ( n > get_npo(c) ) go to 1000 3359 dx=absi(c,n) 3360 if( dx.gt.posi(c,nmk(c)) ) go to 1000 3361 go to 100 3362 1000 continue 3363 3364 end do 3365 call log_mess('END PDED V5....',DEBUG_DEF) 3366 return 3367 end subroutine pded_v5
pded_v5_kd
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pded_v5_kd
DESCRIPTION
NOTES
Calcul des probabilites de transmission le long du chromosome
SOURCE
3458 subroutine pded_v5_kd(c,ip,jm,geno,kd,ilong,recm,recf,trecm,trecf,srecm,srecf) 3459 integer , intent(in) :: c,ip,jm,geno,kd,ilong 3460 real(kind=dp) , intent(in) :: recm,recf 3461 real(kind=dp) , intent(in) :: trecm(maxval(nmk),maxval(nmk)),trecf(maxval(nmk),maxval(nmk)) 3462 real(kind=dp) , intent(in) :: srecm(maxval(nmk),maxval(nmk)),srecf(maxval(nmk),maxval(nmk)) 3463 ! 3464 ! Tableaux dimensionnes selon l, le nombre de marqueurs 3465 integer trans(maxval(nmk)) 3466 integer indp(maxval(nmk)),indm(maxval(nmk)),ingp(maxval(nmk)),ingm(maxval(nmk)) 3467 integer indi(maxval(nmk)),ingi(maxval(nmk)) 3468 integer hdp(maxval(nmk)),hdm(maxval(nmk)),hgp(maxval(nmk)),hgm(maxval(nmk)) 3469 integer nbdp(maxval(nmk)),nbdm(maxval(nmk)) 3470 integer nbgp(maxval(nmk)),nbgm(maxval(nmk)) 3471 integer transp(maxval(nmk)),inforp(maxval(nmk)) 3472 real(kind=dp) :: pbt(4),pbtg(4),pbtd(4),pcum(2),pbante(4,2) 3473 3474 logical infor(maxval(nmk)) 3475 3476 ! Divers 3477 integer new(11,8) 3478 integer hdi,hdifin,hdideb,hgi,hgideb,hgifin,hlocp,hlocm 3479 integer idm,igm,idp,igp 3480 integer icoloc,ii,icas,idi,igi,igeno,iqtl,iqtlp,iqtlm 3481 integer j,jdi,jgi,jj,jk,jn,jnd,jma 3482 integer kcas,kkd 3483 integer lecas,linf,linfp,lk,llk,llkp,lma,lma1,lma2,ln,lq,lqtl 3484 integer maxcas 3485 integer n,nbdi,nbgi 3486 integer ndeb,nfin 3487 integer nmkinf,npas,nposi,nmkinfp 3488 3489 integer kcutopt,temps,tempsmax 3490 integer kcutd,lcutd,ncutd,rcutd,kcutg,lcutg,ncutg,rcutg 3491 3492 real(kind=dp) :: apm,app 3493 real(kind=dp) :: ddp,ddm,dgp,dgm,dlq,dqr 3494 real(kind=dp) :: dx,pbtt,pm,pp,prob 3495 real(kind=dp) :: reclq,recqr 3496 real(kind=dp) :: recdp,recdm,trecdp,trecdm,srecdp,srecdm 3497 real(kind=dp) :: recgp,recgm,trecgp,trecgm,srecgp,srecgm 3498 real(kind=dp) :: xint,xintp,xintm 3499 logical coloc,colocdp,colocdm,colocgp,colocgm,finp,finm,possible 3500 ! 3501 ! tableau des transmission 3502 ! 3503 data new / 1,1,2,2,1,1,1,2,1,1,1,& 3504 1,2,1,2,1,2,1,1,1,2,1,& 3505 0,0,0,0,2,2,1,2,2,2,1,& 3506 0,0,0,0,2,1,2,2,1,2,2,& 3507 0,0,0,0,0,0,0,0,0,0,2,& 3508 0,0,0,0,0,0,0,0,0,0,1,& 3509 0,0,0,0,0,0,0,0,0,0,2,& 3510 0,0,0,0,0,0,0,0,0,0,2/ 3511 3512 3513 kkd=ndesc(c,kd) 3514 ! print *,"new kkd:",kkd 3515 ! 3516 ! le vecteur trans(lk) indique l'information jointe sur la transmission au locus lk 3517 ! le vecteur infor contient la liste des marqeurs doublement informatifs 3518 ! 3519 ! l'origine grand parentale des all�les recu par le descendant k au marqueur lk 3520 ! d�pend de trans (tableau new) 3521 ! soit 1 pour grand p�re et 2 pour grand m�re 3522 ! 3523 ! cas 1 cas 2 cas 3 cas 4 3524 ! trans(lk) pere mere pere mere pere mere pere mere 3525 ! 1 1 1 3526 ! 2 1 2 3527 ! 3 2 1 3528 ! 4 2 2 3529 ! 5 1 1 2 2 3530 ! 6 1 2 2 1 3531 ! 7 1 1 1 2 3532 ! 8 2 1 2 2 3533 ! 9 1 1 2 1 3534 ! 10 1 2 2 2 3535 ! 11 1 1 1 2 2 1 2 2 3536 ! 3537 ! 3538 3539 llkp=0 3540 do lk=1,nmk(c) 3541 3542 ! 3543 ! infor(lk) est � true si le marqueur lk est informatif pour les deux parents 3544 ! 3545 infor(lk)=.false. 3546 ! 3547 ! incertitude totale (donnees manquantes) 3548 ! 3549 trans(lk)=11 3550 ! 3551 ! triplet h�t�rozygotes 3552 ! 3553 if(ptfin(c,lk,kd,2).eq.0.and.ptfin(c,lk,kd,3).eq.0)trans(lk)=5 3554 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=6 3555 ! 3556 ! incertitude sur un parent 3557 ! 3558 if(ptfin(c,lk,kd,3).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=7 3559 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,2).eq.0)trans(lk)=8 3560 if(ptfin(c,lk,kd,2).eq.0.and.ptfin(c,lk,kd,4).eq.0)trans(lk)=9 3561 if(ptfin(c,lk,kd,1).eq.0.and.ptfin(c,lk,kd,3).eq.0)trans(lk)=10 3562 ! 3563 ! cas certains 3564 ! 3565 if(ptfin(c,lk,kd,1).eq.1.d0) trans(lk)=1 3566 if(ptfin(c,lk,kd,2).eq.1.d0) trans(lk)=2 3567 if(ptfin(c,lk,kd,3).eq.1.d0) trans(lk)=3 3568 if(ptfin(c,lk,kd,4).eq.1.d0) trans(lk)=4 3569 3570 if(trans(lk).le.4) infor(lk)=.true. 3571 ! 3572 ! 3573 ! cas des familles de demi fr�res 3574 ! 3575 !if(force(jm).or.opt_sib.eq.OPT_SIB_HS) then 3576 if(.not.estfem(repfem(jm)).or.opt_sib.eq.OPT_SIB_HS)then 3577 transp(lk)=0 3578 if(trans(lk).le.2.or.trans(lk).eq.7)transp(lk)=1 3579 if(trans(lk).eq.3.or.trans(lk).eq.4.or.trans(lk).eq.8)transp(lk)=2 3580 if(transp(lk).ne.0) then 3581 llkp=llkp+1 3582 inforp(llkp)=lk 3583 end if 3584 end if 3585 end do 3586 nmkinfp=llkp 3587 ! 3588 ! exploration du groupe de liaison pour la recherche des transmissions possibles 3589 ! 3590 n=1 3591 100 dx=absi(c,n) 3592 ! 3593 ! pour des probl�mes d'arrondis, la position test�e pour la qtl 3594 ! peut d�passer le bout du groupe de liaison... on l'y ramene 3595 ! 3596 if(dx.ge.posi(c,nmk(c)))dx=posi(c,nmk(c)) 3597 ! 3598 ! cas des demi fr�res 3599 ! 3600 ! if(force(jm).or.opt_sib.eq.OPT_SIB_HS) then 3601 if(.not.estfem(repfem(jm)).or.opt_sib.eq.OPT_SIB_HS) then 3602 do ii=1,4 3603 pbt(ii)=0.d0 3604 end do 3605 ! 3606 ! cas des descendants sans information 3607 ! 3608 if(nmkinfp.eq.0) then 3609 do ii=1,3,2 3610 pbt(ii)=0.5d0 3611 end do 3612 go to 130 ! condition de sortie 3613 end if 3614 ! 3615 ! descendants avec information 3616 ! 3617 linfp=1 3618 ! 3619 ! position du qtl avant le premier marqueur informatif 3620 ! 3621 if(dx.lt.posi(c,inforp(1)))then 3622 dqr=posi(c,inforp(1))-dx 3623 if(inforp(1).eq.1) then 3624 if(posi(c,1).ne.0.)dqr=dqr*posim(c,1)/posi(c,1) 3625 else 3626 xint=posi(c,inforp(1))-posi(c,inforp(1)-1) 3627 xintm=posim(c,inforp(1))-posim(c,inforp(1)-1) 3628 if(xint.ne.0.d0)dqr=dqr*xintm/xint 3629 end if 3630 recqr=xaldane(dqr) 3631 3632 if(transp(inforp(1)).eq.1) then 3633 pbt(1)=1.d0-recqr 3634 pbt(3)=recqr 3635 else 3636 pbt(1)=recqr 3637 pbt(3)=1.d0-recqr 3638 end if 3639 3640 go to 130 3641 end if 3642 3643 5 linfp=linfp+1 3644 if(linfp.gt.nmkinfp) then 3645 ! 3646 ! qtl apr�s le dernier marqueur informatif 3647 ! 3648 dlq=dx-posi(c,inforp(nmkinfp)) 3649 if(inforp(nmkinfp).eq.nmk(c)) then 3650 xint=posi(c,inforp(nmkinfp))-posi(c,inforp(nmkinfp)-1) 3651 xintm=posim(c,inforp(nmkinfp))-posim(c,inforp(nmkinfp)-1) 3652 if(xint.ne.0.d0) dlq=dlq*xintm/xint 3653 else 3654 xint=posi(c,inforp(nmkinfp)+1)-posi(c,inforp(nmkinfp)) 3655 xintm=posim(c,inforp(nmkinfp)+1)-posim(c,inforp(nmkinfp)) 3656 if(xint.ne.0.d0) dlq=dlq*xintm/xint 3657 end if 3658 reclq=xaldane(dlq) 3659 if(transp(inforp(nmkinfp)).eq.1) then 3660 pbt(1)=1.d0-reclq 3661 pbt(3)=reclq 3662 else 3663 pbt(1)=reclq 3664 pbt(3)=1.d0-reclq 3665 end if 3666 3667 go to 130 3668 3669 else 3670 3671 if(.not.(dx.ge.posi(c,inforp(linfp-1)) .and. dx.le.posi(c,inforp(linfp)))) go to 5 3672 ! 3673 ! qtl entre deux marqueurs informatifs 3674 ! 3675 3676 dqr=posi(c,inforp(linfp))-dx 3677 dlq=dx-posi(c,inforp(linfp-1)) 3678 xint=posi(c,inforp(linfp))-posi(c,inforp(linfp-1)) 3679 xintm=posim(c,inforp(linfp))-posim(c,inforp(linfp-1)) 3680 if(xint.ne.0.d0)then 3681 dqr=dqr*xintm/xint 3682 recqr=xaldane(dqr) 3683 dlq=dlq*xintm/xint 3684 reclq=xaldane(dlq) 3685 end if 3686 3687 if(transp(inforp(linfp-1)).eq.transp(inforp(linfp)))then 3688 if(transp(inforp(linfp)).eq.1)then 3689 pbt(1)=(1.d0-reclq)*(1.d0-recqr) 3690 pbt(3)=reclq*recqr 3691 else 3692 pbt(1)=reclq*recqr 3693 pbt(3)=(1.d0-reclq)*(1.d0-recqr) 3694 end if 3695 else 3696 if(transp(inforp(linfp)).eq.1)then 3697 pbt(1)=reclq*(1.d0-recqr) 3698 pbt(3)=(1.d0-reclq)*recqr 3699 else 3700 pbt(1)=(1.d0-reclq)*recqr 3701 pbt(3)=reclq*(1.d0-recqr) 3702 end if 3703 end if 3704 3705 go to 130 3706 3707 end if 3708 end if 3709 ! 3710 ! cas des m�lange plein / demi fr�res 3711 ! 3712 ! 3713 ! cas des positions sur des marqueurs informatifs 3714 ! 3715 ! 3716 ! si (coloc=true) le QTL est sur exactement sur le marqueur lma (icoloc =lma) 3717 ! le cas sera trait� sp�cifiquement pour �viter de prendre Log(z�ro) 3718 ! si de surcroit le marqueur est informatif pour les deux parents, on connait 3719 ! directement la probabilit� de transmission 3720 ! 3721 coloc=.false. 3722 do lk=1,nmk(c) 3723 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,lk)))then 3724 coloc=.true. 3725 icoloc=lk 3726 if(infor(lk))then 3727 do lq=1,4 3728 pbt(lq)=0.d0 3729 end do 3730 pbt(trans(lk))=1.d0 3731 go to 130 3732 end if 3733 end if 3734 end do 3735 ! 3736 ! localisation de la position test�e du qtl 3737 ! 3738 do lk=2,nmk(c) 3739 if(dx.ge.posi(c,lk-1).and.dx.lt.posi(c,lk)) lqtl=lk-1 3740 end do 3741 if(dx.eq.posi(c,nmk(c))) lqtl=nmk(c)-1 3742 ! 3743 ! Recherche des zones informatives pour la transmission 3744 ! La zone du p�re (resp. de la m�re) doit se terminer � gauche et � droite 3745 ! par des marqueurs informatifs (dont l'all�le transmis au descendant est 3746 ! connu), sans la pr�sence de marqueurs incertain (trans = 5 ou 6) entre 3747 ! les bornes � gauche (resp � droite ) des deux parents 3748 ! 3749 ! nbdp est le nombre de marqueurs connu pour le pere � droite du qtl 3750 ! nbdm, pour la mere 3751 ! nbdi, nombre de marqueurs incertains 3752 ! 3753 ! hdp est le vecteur des transmision certaines pour le pere � droite 3754 ! hdm, pour la mere 3755 ! 3756 ! indp est le vecteur des indices (dans la lise initiale) des marqueurs 3757 ! connus pour le pere � droite 3758 ! indm, pour la mere 3759 ! indi, le tableau des indices des marqueurs incertains (colonne 1 : le pere, 2 3760 ! : la mere) 3761 ! 3762 nbdi=1 3763 do lk=1,nmk(c) 3764 nbdp(lk)=0 3765 nbdm(lk)=0 3766 end do 3767 jma=lqtl 3768 finp=.false. 3769 finm=.false. 3770 3771 10 jma=jma+1 3772 if(jma.gt.nmk(c)) go to 15 3773 if(trans(jma).le.4) then 3774 nbdp(nbdi)=nbdp(nbdi)+1 3775 indp(nbdp(nbdi))=jma 3776 hdp(nbdp(nbdi))=floor(real(trans(jma)+1)/2.) 3777 nbdm(nbdi)=nbdm(nbdi)+1 3778 indm(nbdm(nbdi))=jma 3779 hdm(nbdm(nbdi))=trans(jma)-2*(hdp(nbdp(nbdi))-1) 3780 finp=.true. 3781 finm=.true. 3782 end if 3783 3784 if(trans(jma).eq.5.or.trans(jma).eq.6) then 3785 nbdp(nbdi)=nbdp(nbdi)+1 3786 indp(nbdp(nbdi))=jma 3787 hdp(nbdp(nbdi))=0 3788 nbdm(nbdi)=nbdm(nbdi)+1 3789 indm(nbdm(nbdi))=jma 3790 hdm(nbdm(nbdi))=0 3791 indi(nbdi)=jma 3792 nbdi=nbdi+1 3793 nbdp(nbdi)=nbdp(nbdi-1) 3794 nbdm(nbdi)=nbdm(nbdi-1) 3795 finp=.false. 3796 finm=.false. 3797 end if 3798 3799 if(trans(jma).eq.7.or.trans(jma).eq.8) then 3800 nbdp(nbdi)=nbdp(nbdi)+1 3801 indp(nbdp(nbdi))=jma 3802 hdp(nbdp(nbdi))=trans(jma)-6 3803 finp=.true. 3804 end if 3805 3806 3807 if(trans(jma).eq.9.or.trans(jma).eq.10) then 3808 nbdm(nbdi)=nbdm(nbdi)+1 3809 indm(nbdm(nbdi))=jma 3810 hdm(nbdm(nbdi))=trans(jma)-8 3811 finm=.true. 3812 end if 3813 3814 if(((.not.finp).or.(.not.finm)).and.jma.lt.nmk(c)) go to 10 3815 3816 ! 3817 ! on cherche les zones � gauche du qtl 3818 ! les elements sont comme � droite , en remplacant d par g 3819 ! 3820 15 nbgi=1 3821 do lk=1,nmk(c) 3822 nbgp(lk)=0 3823 nbgm(lk)=0 3824 end do 3825 jma=lqtl+1 3826 finp=.false. 3827 finm=.false. 3828 3829 20 jma=jma-1 3830 if(jma.eq.0) go to 25 3831 if(trans(jma).le.4) then 3832 nbgp(nbgi)=nbgp(nbgi)+1 3833 ingp(nbgp(nbgi))=jma 3834 hgp(nbgp(nbgi))=floor(real(trans(jma)+1)/2.) 3835 nbgm(nbgi)=nbgm(nbgi)+1 3836 ingm(nbgm(nbgi))=jma 3837 hgm(nbgm(nbgi))=trans(jma)-2*(hgp(nbgp(nbgi))-1) 3838 finp=.true. 3839 finm=.true. 3840 end if 3841 3842 if(trans(jma).eq.5.or.trans(jma).eq.6) then 3843 nbgp(nbgi)=nbgp(nbgi)+1 3844 ingp(nbgp(nbgi))=jma 3845 hgp(nbgp(nbgi))=0 3846 nbgm(nbgi)=nbgm(nbgi)+1 3847 ingm(nbgm(nbgi))=jma 3848 hgm(nbgm(nbgi))=0 3849 ingi(nbgi)=jma 3850 nbgi=nbgi+1 3851 nbgp(nbgi)=nbgp(nbgi-1) 3852 nbgm(nbgi)=nbgm(nbgi-1) 3853 finp=.false. 3854 finm=.false. 3855 end if 3856 3857 if(trans(jma).eq.7.or.trans(jma).eq.8) then 3858 nbgp(nbgi)=nbgp(nbgi)+1 3859 ingp(nbgp(nbgi))=jma 3860 hgp(nbgp(nbgi))=trans(jma)-6 3861 finp=.true. 3862 end if 3863 3864 3865 if(trans(jma).eq.9.or.trans(jma).eq.10) then 3866 nbgm(nbgi)=nbgm(nbgi)+1 3867 ingm(nbgm(nbgi))=jma 3868 hgm(nbgm(nbgi))=trans(jma)-8 3869 finm=.true. 3870 end if 3871 3872 if(((.not.finp).or.(.not.finm)).and.jma.gt.1) go to 20 3873 3874 ! 3875 ! colocgp et colocgm indique (true) si la position du qtl est colocalis�e 3876 ! avec un marqueur (coloc = true) et que ce marqueur est le premier � 3877 ! sa gauche. Dans ce cas, il n'y a pas de recombinaison possible entre 3878 ! le qtl et ce marqueur 3879 ! 3880 ! hlocp et hlocm sont les informations sur la transmission en ce locus (� la 3881 ! fois marqueur et qtl) 3882 ! 3883 25 colocgp=.false. 3884 colocgm=.false. 3885 colocdp=.false. 3886 colocdm=.false. 3887 if(coloc.and.nbgp(1).ne.0)then 3888 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,ingp(1))))colocgp=.true. 3889 end if 3890 if(coloc.and.nbgm(1).ne.0)then 3891 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,ingm(1))))colocgm=.true. 3892 end if 3893 if(coloc.and.nbdp(1).ne.0)then 3894 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,indp(1))))colocdp=.true. 3895 end if 3896 if(coloc.and.nbdm(1).ne.0)then 3897 if(int(1000.d0*dx).eq.int(1000.d0*posi(c,indm(1))))colocdm=.true. 3898 end if 3899 3900 ! 3901 ! les taux de recombinaison avant et apres le QTL sont stock�s dans 3902 ! (t ou s)rec(d ou g)(p ou m) 3903 ! 3904 idp=indp(1) 3905 if(nbdp(1).eq.0)idp=lqtl+1 3906 igp=ingp(1) 3907 if(nbgp(1).eq.0)igp=lqtl 3908 3909 xint=posi(c,lqtl+1)-posi(c,lqtl) 3910 xintp=posim(c,lqtl+1)-posim(c,lqtl) 3911 ! 3912 ! dans le cas de colocalisation, il ne faut pas calculer le taux de 3913 ! recombinaison 3914 ! 3915 dgp=(dx-posi(c,igp))*xintp/xint 3916 recgp=xaldane(dgp) 3917 ddp=(posi(c,idp)-dx)*xintp/xint 3918 recdp=xaldane(ddp) 3919 ! 3920 ! 3921 ! cas des m�res 3922 ! 3923 idm=indm(1) 3924 if(nbdm(1).eq.0)idm=lqtl+1 3925 igm=ingm(1) 3926 if(nbgm(1).eq.0)igm=lqtl 3927 3928 xint=posi(c,lqtl+1)-posi(c,lqtl) 3929 xintm=posif(c,lqtl+1)-posif(c,lqtl) 3930 ! 3931 ! dans le cas de colocalisation, il ne faut pas calculer le taux de 3932 ! recombinaison 3933 ! 3934 dgm=(dx-posi(c,igm))*xintm/xint 3935 recgm=xaldane(dgm) 3936 ddm=(posi(c,idm)-dx)*xintm/xint 3937 recdm=xaldane(ddm) 3938 3939 ! 3940 ! calcul des probabilit�s de transmission pour le descendant kd, position dx 3941 ! 3942 ! 3943 ! on va envisager toutes les possibilit�s pour les marqueurs incertains 3944 ! il y a 2**(nbdi+nbgi) possibilit�s 3945 ! ainsi que les 4 possibilit�s pour la transmission du QTL 3946 ! 3947 ! les hdp, hdm, hgp et hgm sont compl�t�s par des 1 ou 2 pour les marqueurs 3948 ! incertains 3949 ! 3950 ! les boucles en jdi et jgi ont pour fonction de g�n�rer les 2**(nbdi+nbgi) 3951 ! vecteurs d'origines correspondants aux marqueurs incertains 3952 ! 3953 ! initialisation 3954 pbtt=0.d0 3955 do lq=1,4 3956 pbtg(lq)=0.d0 3957 pbtd(lq)=0.d0 3958 end do 3959 3960 ! 3961 ! on va consid�rer les 4 �v�nements de transmission possible au QTL, d�finis 3962 ! pas iqtlp et iqtlm 3963 ! 3964 ! pour chacun des �venements de transmision au qtl on v�rifie la coh�rence 3965 ! avec un �ventuel marqueur colocalis� 3966 ! 3967 iqtlp=0 3968 30 iqtlp=iqtlp+1 3969 if(iqtlp.eq.3)go to 120 3970 3971 iqtlm=0 3972 40 iqtlm=iqtlm+1 3973 if(iqtlm.eq.3)go to 30 3974 ! 3975 iqtl= iqtlm+2*(iqtlp-1) 3976 3977 ! 3978 ! le cas o� il n'y a pas d'incertitude � droite est trait� � part 3979 ! dans ce cas, il faut distinguer les situations o� il n'y a que des 3980 ! transmissions inconnues 3981 ! 3982 if(nbdi.eq.1) then 3983 ! 3984 ! on a nbdp=0 si aucun marqueur � droite n'a trans diff�rent de 9, 10 ou 11 3985 ! on a nbdm=0 si aucun marqueur � droite n'a trans diff�rent de 7, 8 ou 11 3986 ! dans ce cas,les 2 �v�nements iqtlp = 1 ou 2 sont �quiprobables 3987 ! 3988 if(nbdp(1).eq.0.and.nbdm(1).eq.0) then 3989 pbtd(iqtl)=0.25d0 3990 end if 3991 3992 prob=1.d0 3993 if(nbdm(1).ne.0) then 3994 if(colocdm.and.iqtlm.ne.hdm(1))go to 50 3995 if(.not.colocdm)then 3996 if(nbdp(1).eq.0) prob=0.5d0 3997 if(iqtlm.eq.hdm(1))prob=prob*(1.d0-recdm) 3998 if(iqtlm.ne.hdm(1))prob=prob*recdm 3999 end if 4000 end if 4001 4002 if(nbdp(1).ne.0) then 4003 if(colocdp.and.iqtlp.ne.hdp(1))go to 50 4004 if(.not.colocdp)then 4005 if(nbdm(1).eq.0) prob=0.5d0 4006 if(iqtlp.eq.hdp(1))prob=prob*(1.d0-recdp) 4007 if(iqtlp.ne.hdp(1))prob=prob*recdp 4008 end if 4009 end if 4010 4011 pbtd(iqtl)=prob 4012 go to 50 4013 4014 end if 4015 ! 4016 ! cas ou il y a des incertitudes � droite 4017 ! 4018 ! 4019 ! 4020 ! on traite diff�remment les intervalles entre marqueurs incertains : le 1er, les 4021 ! suivants, le dernier 4022 ! 4023 ! cas du premier intervalle 4024 ! nfin est l'indice du marqueur incertain � droite 4025 ! 4026 idi=1 4027 possible=.false. 4028 do hdifin=1,2 4029 ! 4030 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 4031 ! hdp et hdm du premier intervalle 4032 ! 4033 nfin=indi(idi) 4034 hdp(nbdp(1))=new(trans(nfin),2*hdifin-1) 4035 hdm(nbdm(1))=new(trans(nfin),2*hdifin) 4036 4037 if((colocdm.and.iqtlm.ne.hdm(1)).or.(colocdp.and.iqtlp.ne.hdp(1)))then 4038 4039 pbante(iqtl,hdifin)=0.d0 4040 4041 else 4042 4043 possible=.true. 4044 prob=0.d0 4045 4046 if(.not.colocdm)then 4047 if(iqtlm.eq.hdm(1))prob=dlog(1.d0-recdm) 4048 if(recdm.eq.0) print*,'kd,n',kd,n,dx,ndesc(c,kkd),trim(animal(kkd)) 4049 if(iqtlm.ne.hdm(1))prob=dlog(recdm ) 4050 4051 do jma=2,nbdm(1) 4052 if(hdm(jma-1).eq.hdm(jma))then 4053 prob=prob+srecf(indm(jma-1),indm(jma)) 4054 else 4055 prob=prob+trecf(indm(jma-1),indm(jma)) 4056 end if 4057 end do 4058 end if 4059 4060 if(.not.colocdp)then 4061 if(iqtlp.eq.hdp(1))prob=prob+dlog(1.d0-recdp) 4062 if(iqtlp.ne.hdp(1))prob=prob+dlog(recdp) 4063 do jma=2,nbdp(1) 4064 if(hdp(jma-1).eq.hdp(jma))then 4065 prob=prob+srecm(indp(jma-1),indp(jma)) 4066 else 4067 prob=prob+trecm(indp(jma-1),indp(jma)) 4068 end if 4069 end do 4070 end if 4071 4072 pbante(iqtl,hdifin)=dexp(prob) 4073 4074 end if 4075 4076 end do 4077 4078 if(.not.possible) go to 50 4079 ! 4080 ! cas des intervalles suivants 4081 ! 4082 do idi=2,nbdi-1 4083 4084 do hdifin=1,2 4085 ! 4086 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 4087 ! hdp et hdm du idi �me intervalle 4088 ! 4089 nfin=indi(idi) 4090 hdp(nbdp(idi))=new(trans(nfin),2*hdifin-1) 4091 hdm(nbdm(idi))=new(trans(nfin),2*hdifin) 4092 pcum(hdifin)=0.d0 4093 4094 do hdideb=1,2 4095 ! 4096 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 4097 ! de l'intervalle pr�c�dent 4098 ! 4099 ndeb=indi(idi-1) 4100 hdp(nbdp(idi-1))=new(trans(ndeb),2*hdideb-1) 4101 hdm(nbdm(idi-1))=new(trans(ndeb),2*hdideb) 4102 4103 prob=0.d0 4104 do jma=nbdp(idi-1)+1,nbdp(idi) 4105 if(hdp(jma-1).eq.hdp(jma))then 4106 prob=prob+srecm(indp(jma-1),indp(jma)) 4107 else 4108 prob=prob+trecm(indp(jma-1),indp(jma)) 4109 end if 4110 end do 4111 4112 do jma=nbdm(idi-1)+1,nbdm(idi) 4113 if(hdm(jma-1).eq.hdm(jma))then 4114 prob=prob+srecf(indm(jma-1),indm(jma)) 4115 else 4116 prob=prob+trecf(indm(jma-1),indm(jma)) 4117 end if 4118 end do 4119 ! 4120 pcum(hdifin)=pcum(hdifin)+pbante(iqtl,hdideb)*dexp(prob) 4121 4122 end do 4123 end do 4124 4125 do hdifin=1,2 4126 pbante(iqtl,hdifin)=pcum(hdifin) 4127 end do 4128 4129 end do 4130 4131 ! 4132 ! cas du dernier intervalle 4133 ! 4134 idi=nbdi 4135 pcum(1)=0.d0 4136 do hdideb=1,2 4137 ! 4138 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 4139 ! de l'intervalle pr�c�dent 4140 ! 4141 ndeb=indi(idi-1) 4142 hdp(nbdp(idi-1))=new(trans(ndeb),2*hdideb-1) 4143 hdm(nbdm(idi-1))=new(trans(ndeb),2*hdideb) 4144 4145 prob=0.d0 4146 if(nbdp(idi).ne.nbdp(idi-1))then 4147 do jma=nbdp(idi-1)+1,nbdp(idi) 4148 if(hdp(jma-1).eq.hdp(jma))then 4149 prob=prob+srecm(indp(jma-1),indp(jma)) 4150 else 4151 prob=prob+trecm(indp(jma-1),indp(jma)) 4152 end if 4153 end do 4154 end if 4155 4156 if(nbdm(idi).ne.nbdm(idi-1))then 4157 do jma=nbdm(idi-1)+1,nbdm(idi) 4158 if(hdm(jma-1).eq.hdm(jma))then 4159 prob=prob+srecf(indm(jma-1),indm(jma)) 4160 else 4161 prob=prob+trecf(indm(jma-1),indm(jma)) 4162 end if 4163 end do 4164 end if 4165 ! 4166 pcum(1)=pcum(1)+pbante(iqtl,hdideb)*dexp(prob) 4167 4168 end do 4169 4170 pbtd(iqtl)=pcum(1) 4171 4172 50 continue 4173 ! 4174 ! on �num�re les possibilit�s pour les transmissions incertaines � gauche 4175 ! du qtl. 4176 ! 4177 ! 4178 ! le cas o� il n'y a pas d'incertitude � gauche est trait� � part 4179 ! dans ce cas, il faut distinguer les situations o� il n'y a que des 4180 ! transmissions inconnues 4181 ! 4182 4183 if(nbgi.eq.1) then 4184 ! 4185 ! on a nbgp=0 si aucun marqueur � gauche n'a trans diff�rent de 9, 10 ou 11 4186 ! on a nbgm=0 si aucun marqueur � gauche n'a trans diff�rent de 7, 8 ou 11 4187 ! dans ce cas,les 2 �v�nements iqtlp = 1 ou 2 sont �quiprobables 4188 ! 4189 if(nbgp(1).eq.0.and.nbgm(1).eq.0) then 4190 pbtg(iqtl)=0.25d0 4191 end if 4192 4193 prob=1.d0 4194 if(nbgm(1).ne.0) then 4195 if(colocgm.and.iqtlm.ne.hgm(1))go to 60 4196 if(.not.colocgm)then 4197 if(nbgp(1).eq.0) prob=0.5d0 4198 if(iqtlm.eq.hgm(1))prob=prob*(1.d0-recgm) 4199 if(iqtlm.ne.hgm(1))prob=prob*recgm 4200 end if 4201 end if 4202 4203 if(nbgp(1).ne.0) then 4204 if(colocgp.and.iqtlp.ne.hgp(1))go to 60 4205 if(.not.colocgp)then 4206 if(nbgm(1).eq.0) prob=0.5d0 4207 if(iqtlp.eq.hgp(1))prob=prob*(1.d0-recgp) 4208 if(iqtlp.ne.hgp(1))prob=prob*recgp 4209 end if 4210 end if 4211 4212 pbtg(iqtl)=prob 4213 go to 60 4214 4215 end if 4216 ! 4217 ! cas ou il y a des incertitudes � gauche 4218 ! 4219 ! 4220 ! 4221 ! on traite diff�remment les intervalles entre marqueurs incertains : le 1er, les 4222 ! suivants, le dernier 4223 ! 4224 ! cas du premier intervalle 4225 ! nfin est l'indice du marqueur incertain � gauche 4226 ! 4227 igi=1 4228 possible=.false. 4229 do hgifin=1,2 4230 4231 ! 4232 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 4233 ! hgp et hgm du premier intervalle 4234 ! 4235 nfin=ingi(igi) 4236 hgp(nbgp(1))=new(trans(nfin),2*hgifin-1) 4237 hgm(nbgm(1))=new(trans(nfin),2*hgifin) 4238 4239 if( (colocgm.and.iqtlm.ne.hgm(1)).or.(colocgp.and.iqtlp.ne.hgp(1)))then 4240 4241 pbante(iqtl,hgifin)=0.d0 4242 4243 else 4244 4245 possible=.true. 4246 prob=0.d0 4247 4248 if(.not.colocgm)then 4249 if(iqtlm.eq.hgm(1))prob=dlog(1.d0-recgm) 4250 if(iqtlm.ne.hgm(1))prob=dlog(recgm) 4251 4252 do jma=2,nbgm(1) 4253 if(hgm(jma-1).eq.hgm(jma))then 4254 prob=prob+srecf(ingm(jma-1),ingm(jma)) 4255 else 4256 prob=prob+trecf(ingm(jma-1),ingm(jma)) 4257 end if 4258 end do 4259 end if 4260 4261 if(.not.colocgp)then 4262 if(iqtlp.eq.hgp(1))prob=prob+dlog(1.d0-recgp) 4263 if(iqtlp.ne.hgp(1))prob=prob+dlog(recgp) 4264 do jma=2,nbgp(1) 4265 if(hgp(jma-1).eq.hgp(jma))then 4266 prob=prob+srecm(ingp(jma-1),ingp(jma)) 4267 else 4268 prob=prob+trecm(ingp(jma-1),ingp(jma)) 4269 end if 4270 end do 4271 end if 4272 4273 pbante(iqtl,hgifin)=dexp(prob) 4274 4275 end if 4276 4277 end do 4278 4279 if(.not.possible)go to 60 4280 ! 4281 ! cas des intervalles suivants 4282 ! 4283 do igi=2,nbgi-1 4284 4285 do hgifin=1,2 4286 ! 4287 ! on a deux possibilit�s couvrant l'incertitude, donnant les derniers 4288 ! hgp et hgm du igi �me intervalle 4289 ! 4290 nfin=ingi(igi) 4291 hgp(nbgp(igi))=new(trans(nfin),2*hgifin-1) 4292 hgm(nbgm(igi))=new(trans(nfin),2*hgifin) 4293 pcum(hgifin)=0.d0 4294 4295 do hgideb=1,2 4296 ! 4297 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 4298 ! de l'intervalle pr�c�dent 4299 ! 4300 ndeb=ingi(igi-1) 4301 hgp(nbgp(igi-1))=new(trans(ndeb),2*hgideb-1) 4302 hgm(nbgm(igi-1))=new(trans(ndeb),2*hgideb) 4303 4304 prob=0.d0 4305 do jma=nbgp(igi-1)+1,nbgp(igi) 4306 if(hgp(jma-1).eq.hgp(jma))then 4307 prob=prob+srecm(ingp(jma-1),ingp(jma)) 4308 else 4309 prob=prob+trecm(ingp(jma-1),ingp(jma)) 4310 end if 4311 end do 4312 4313 do jma=nbgm(igi-1)+1,nbgm(igi) 4314 if(hgm(jma-1).eq.hgm(jma))then 4315 prob=prob+srecf(ingm(jma-1),ingm(jma)) 4316 else 4317 prob=prob+trecf(ingm(jma-1),ingm(jma)) 4318 end if 4319 end do 4320 ! 4321 pcum(hgifin)=pcum(hgifin)+pbante(iqtl,hgideb)*dexp(prob) 4322 4323 end do 4324 end do 4325 4326 do hgifin=1,2 4327 pbante(iqtl,hgifin)=pcum(hgifin) 4328 end do 4329 4330 end do 4331 4332 ! 4333 ! cas du dernier intervalle 4334 ! 4335 igi=nbgi 4336 pcum(1)=0.d0 4337 do hgideb=1,2 4338 ! 4339 ! il faut consid�rer les deux possibilit�s couvrant l'incertitude, 4340 ! de l'intervalle pr�c�dent 4341 ! 4342 ndeb=ingi(igi-1) 4343 hgp(nbgp(igi-1))=new(trans(ndeb),2*hgideb-1) 4344 hgm(nbgm(igi-1))=new(trans(ndeb),2*hgideb) 4345 4346 prob=0.d0 4347 if(nbgp(igi).ne.nbgp(igi-1))then 4348 do jma=nbgp(igi-1)+1,nbgp(igi) 4349 if(hgp(jma-1).eq.hgp(jma))then 4350 prob=prob+srecm(ingp(jma-1),ingp(jma)) 4351 else 4352 prob=prob+trecm(ingp(jma-1),ingp(jma)) 4353 end if 4354 end do 4355 end if 4356 4357 if(nbgm(igi).ne.nbgm(igi-1))then 4358 do jma=nbgm(igi-1)+1,nbgm(igi) 4359 if(hgm(jma-1).eq.hgm(jma))then 4360 prob=prob+srecf(ingm(jma-1),ingm(jma)) 4361 else 4362 prob=prob+trecf(ingm(jma-1),ingm(jma)) 4363 end if 4364 end do 4365 end if 4366 ! 4367 pcum(1)=pcum(1)+pbante(iqtl,hgideb)*dexp(prob) 4368 4369 end do 4370 4371 pbtg(iqtl)= pcum(1) 4372 4373 60 continue 4374 ! 4375 ! on teste si d'autres combinaisons au qtl sont � �valuer 4376 ! 4377 if(iqtlp.eq.2.and.iqtlm.eq.2) go to 120 4378 if(iqtlm.eq.1)go to 40 4379 if(iqtlm.eq.2)go to 30 4380 4381 4382 120 continue 4383 4384 ! 4385 ! on stocke la probabilit� de transmission � la position dx 4386 ! 4387 do lq=1,4 4388 pbt(lq)=pbtg(lq)*pbtd(lq) 4389 end do 4390 4391 130 continue 4392 4393 pbtt=0.d0 4394 do lq=1,4 4395 pbtt=pbtt+pbt(lq) 4396 end do 4397 4398 do lq=1,4 4399 pdd(c,kd,lq,n)=pbt(lq)/pbtt 4400 end do 4401 4402 ! 4403 ! 4404 ! on avance d'un pas (en verifiant qu'on ne depasse pas le bout droit) et on recommence 4405 ! 4406 n=n+1 4407 if ( n > get_npo(c) ) return 4408 dx=absi(c,n) 4409 if( dx.gt.posi(c,nmk(c)) ) return 4410 go to 100 4411 4412 return 4413 end subroutine pded_v5_kd
pded_v5_optim
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pded_v5_optim
DESCRIPTION
NOTES
Calcul des probabilites de transmission le long du chromosome
SOURCE
3378 subroutine pded_v5_optim() 3379 real(kind=dp) :: trecm(maxval(nmk),maxval(nmk)),trecf(maxval(nmk),maxval(nmk)) 3380 real(kind=dp) :: srecm(maxval(nmk),maxval(nmk)),srecf(maxval(nmk),maxval(nmk)) 3381 3382 ! Divers 3383 integer :: c,ilong,lk,jk,ip,jm,geno,kd,iITer,nIter,sizeT 3384 real(kind=dp) :: recm,recf 3385 integer ,dimension(:),allocatable :: ipT,jmT,genoT,kdT 3386 3387 sizeT = maxval(ngend(:,(maxval(ngenom(:,nm+1)))+1)) 3388 allocate ( ipT( sizeT ) ) 3389 allocate ( jmT( sizeT ) ) 3390 allocate ( genoT( sizeT ) ) 3391 allocate ( kdT( sizeT ) ) 3392 3393 do c=1,nchr 3394 ! 3395 ! Taille du segment explore 3396 ilong=get_ilong(c) 3397 ! 3398 ! tableau des recombinaison 3399 ! 3400 do lk=1,nmk(c)-1 3401 do jk=lk+1,nmk(c) 3402 recm=xaldane(posim(c,jk)-posim(c,lk)) 3403 recf=xaldane(posif(c,jk)-posif(c,lk)) 3404 trecm(lk,jk)=dlog(recm) 3405 trecf(lk,jk)=dlog(recf) 3406 srecm(lk,jk)=dlog(1.d0-recm) 3407 srecf(lk,jk)=dlog(1.d0-recf) 3408 trecm(jk,lk)=trecm(lk,jk) 3409 trecf(jk,lk)=trecf(lk,jk) 3410 srecm(jk,lk)=srecm(lk,jk) 3411 srecf(jk,lk)=srecf(lk,jk) 3412 end do 3413 end do 3414 ! 3415 ! Calcul pour chaque descendant de la probabilite de transmission 3416 3417 !calcul du nopmbre nombre d iteration 3418 nIter = 0 3419 do ip=1,np 3420 do jm=nmp(ip)+1,nmp(ip+1) 3421 do geno=ngenom(c,jm)+1,ngenom(c,jm+1) 3422 do kd=ngend(c,geno)+1,ngend(c,geno+1) 3423 nIter = nIter + 1 3424 ipT(nIter) = ip 3425 jmT(nIter) = jm 3426 genoT(nIter) = geno 3427 kdT(nIter) = kd 3428 end do 3429 end do 3430 end do 3431 end do 3432 3433 !Calcul + parallelisme 3434 !$OMP PARALLEL DO DEFAULT(SHARED) 3435 do iITer = 1, nIter 3436 call pded_v5_kd(c,ipT(iIter),jmT(iIter),genoT(iIter),kdT(iIter),ilong,recm,recf,trecm,trecf,srecm,srecf) 3437 end do 3438 !$OMP END PARALLEL DO 3439 3440 end do !c 3441 3442 deallocate ( ipT ) 3443 deallocate ( jmT ) 3444 deallocate ( genoT ) 3445 deallocate ( kdT ) 3446 3447 end subroutine pded_v5_optim
pdegm
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pdegm
DESCRIPTION
NOTES
Determination des probabilites des phases pour les meres
SOURCE
788 subroutine pdegm() 789 ! 790 ! Divers 791 integer ,dimension(:,:) ,allocatable :: ngend_t,ndesc_t 792 real (kind=dp),dimension(:,:),allocatable :: probg_t 793 integer geno,t1,t2,t 794 real (kind=dp) :: trans(4,4),s(4),q(4) 795 integer(kind=KIND_PHENO) :: m1,m2 796 integer i,ii,icon,ifem,igeno,imark,j,jj,jm,jnd,kd,c 797 integer ld,ll,ll1,ll2,lmark,lnd,nbcon,nbinc,nd1,nd2 798 integer ngeno,ngeno1,ngeno2,nhomo,nmark,alloc_stat,value_alloc 799 integer nmknmk,nombfem,nphm,ndf(nm+1),repdes(nm,nd),corref(nm) 800 integer , dimension(maxval(nmk)+1) :: hcon,indcon,indinc 801 real (kind=dp) :: pmax,pr,probmax_t,r1,r2,sprob,sprobphm,xlvi 802 logical :: hetero(maxval(nmk)+1) 803 804 integer , dimension(:,:),allocatable ::tphm,h 805 real (kind=dp) , dimension(:),allocatable :: probphm,prob 806 real (kind=dp) , dimension(:,:,:),allocatable ::p 807 real (kind=dp) , dimension(:,:,:,:),allocatable ::ptfin_t 808 809 call log_mess('calling...PDEGM_V4',DEBUG_DEF) 810 allocate (p(maxval(nmk),nd,4),stat=alloc_stat) 811 812 if (alloc_stat/=0) then 813 call stop_application('Not enough memory to compute transmission probabilities. Deacrease number of marker.') 814 end if 815 ! 816 !****************************************************************************** 817 ! 818 ! 819 ! Creation des tableaux 820 ! ndf (nombre total de descendants/femelle) 821 ! repdes (identification des descendants de chaque femelles selon son nuero d'ordre 822 ! en lecture) 823 ! corref (l'equivalent de correm pour les femelles) 824 ! 825 ! cette partie sera plutot e positionner dans lect_genea et lect_typ 826 ! 827 do ifem=1,nfem 828 ndf(ifem)=0 829 end do 830 831 do jm=1,nm 832 nd1=ndm(jm)+1 833 nd2=ndm(jm+1) 834 ifem=repfem(jm) 835 do kd=nd1,nd2 836 ndf(ifem)=ndf(ifem)+1 837 repdes(ifem,ndf(ifem))=kd 838 end do 839 840 corref(ifem)=correm(jm) 841 end do 842 843 ! 844 ! le tableau force indique si une femelle n'a pas assez de descendants 845 ! pour etre estimee 846 ! 847 ! do jm=1,nm 848 ! force(jm)=.false. 849 ! end do 850 ! 851 allocate ( ngend(nchr,1) ) 852 allocate ( probg(nchr,1) ) 853 allocate ( ndesc(nchr,1) ) 854 855 ngenom(:,1)=0 856 ngend(:,1)=0 857 probg=0.d0 858 ! 859 !****************************************************************************** 860 ! 861 ! on traite les femelles (et non les meres) une e une 862 ! 863 ! nombfem est le nombre de femelles dont on peut chercher les phase 864 ! 865 do c=1,nchr 866 nombfem=0 867 do 500 jm=1,nm 868 nd1=ndm(jm)+1 869 nd2=ndm(jm+1) 870 ifem=repfem(jm) 871 ! 872 !****************************************************************************** 873 ! 874 ! Si le nombre de pleins freres est trop faible, le genotype de la mere ne sera pas considere 875 ! Dans ce cas le premier genotype rencontre possible est considere comme le bon 876 877 !if(ndf(ifem).lt.ndmin) then 878 ! force(jm)=.true. 879 if(.not.estfem(ifem))then 880 881 if (allocated(h)) deallocate(h) 882 if (allocated(tphm)) deallocate(tphm) 883 if (allocated(prob)) deallocate(prob) 884 if (allocated(probphm)) deallocate(probphm) 885 allocate (h(nmk(c),1)) 886 allocate (tphm(nmk(c),1)) 887 allocate (prob(1)) 888 allocate (probphm(1)) 889 890 ! on cherche le premier genotype possible 891 ! 892 do nmark=1,nmk(c) 893 tphm(nmark,1)=1 894 if(ordref(c,ifem,nmark).eq.21)tphm(nmark,1)=2 895 end do 896 nphm=1 897 898 probphm(1)=1.d0 899 sprobphm=probphm(1) 900 901 go to 400 902 903 end if 904 ! 905 !****************************************************************************** 906 ! 907 ! 908 ! on cree le tableau hetero qui liste l'heterozygotie des marqueurs 909 ! 910 911 nhomo=0 912 do nmark=1,nmk(c) 913 hetero(nmark)=.true. 914 if(corref(ifem).eq.9999) then 915 m1=nmanque 916 else 917 m1=pheno(c,nmark,corref(ifem),1) 918 m2=pheno(c,nmark,corref(ifem),2) 919 end if 920 if(m1.eq.nmanque) m2=m1 921 if(m1.eq.m2) then 922 nhomo=nhomo+1 923 hetero(nmark)=.false. 924 end if 925 end do 926 ! 927 ! creation du tableau h (anciennement dans combine) 928 ! en se restreignant aux phases utiles 929 ! 930 ! on commence par lister les marqueurs dont l'origine est connue 931 ! soit par l'ascendance, soit parcequ'ils sont homozygotes 932 ! (dans ce cas on donne l'ordrep=12) 933 ! 934 935 nbcon=0 936 nbinc=0 937 do lmark=1,nmk(c) 938 939 icon=0 940 if(ordref(c,ifem,lmark).eq.12)then 941 nbcon=nbcon+1 942 hcon(nbcon)=1 943 indcon(nbcon)=lmark 944 icon=1 945 end if 946 947 if(ordref(c,ifem,lmark).eq.21)then 948 nbcon=nbcon+1 949 hcon(nbcon)=2 950 indcon(nbcon)=lmark 951 icon=1 952 end if 953 954 if(.not.hetero(lmark))then 955 nbcon=nbcon+1 956 hcon(nbcon)=1 957 indcon(nbcon)=lmark 958 icon=1 959 end if 960 961 if(icon.eq.0) then 962 nbinc=nbinc+1 963 indinc(nbinc)=lmark 964 end if 965 966 end do 967 968 !--------------------- ALLOC DYNAMIQUE --------------------- 969 ! ---- A FAIRE VALIDER PAR JM ---------------------------- 970 if (allocated(h)) deallocate(h) 971 if (allocated(tphm)) deallocate(tphm) 972 if (allocated(prob)) deallocate(prob) 973 if (allocated(probphm)) deallocate(probphm) 974 975 if (nbinc>0) then 976 value_alloc= (2**(nbinc-1))*2 977 allocate (h(nmk(c),value_alloc),STAT=alloc_stat) 978 allocate (tphm(nmk(c),value_alloc),STAT=alloc_stat) 979 allocate (prob(value_alloc),STAT=alloc_stat) 980 allocate (probphm(value_alloc),STAT=alloc_stat) 981 else 982 allocate (h(nmk(c),1),STAT=alloc_stat) 983 allocate (tphm(nmk(c),1),STAT=alloc_stat) 984 allocate (prob(1),STAT=alloc_stat) 985 allocate (probphm(1),STAT=alloc_stat) 986 end if 987 988 if (alloc_stat/=0) then 989 call stop_application('Not enough memory to compute transmission probabilities. Deacrease number of marker.') 990 end if 991 992 !--------------------------- END OFI --------------------- 993 ! 994 ! s'il reste des inconnues 995 ! on remplit le tableau h, en completant les elements connus issus 996 ! de indcon et hcon par les phases possibles des elements inconnus 997 ! 998 if(nbinc.ne.0) then 999 nmknmk=2**(nbinc-1) 1000 ngeno=nmknmk 1001 do geno=1,nmknmk 1002 1003 do lmark=1,nbcon 1004 h(indcon(lmark),geno)=hcon(lmark) 1005 end do 1006 1007 do lmark=1,nbinc 1008 lnd=2**(lmark-1) 1009 jnd=floor(float(lnd-1+geno)/float(lnd)) 1010 h(indinc(nbinc+1-lmark),geno)=1+mod(jnd+1,2) 1011 end do 1012 1013 end do 1014 ! 1015 ! quand le genotype est oriente (phasm(c,jm) e true), il faut aussi calculer les 1016 ! probabilites des phases miroirs 1017 ! 1018 1019 if(phasm(c,jm)) then 1020 ngeno=ngeno+nmknmk 1021 do geno=1,nmknmk 1022 1023 do lmark=1,nbcon 1024 h(indcon(lmark),nmknmk+geno)=hcon(lmark) 1025 end do 1026 1027 do lmark=1,nbinc 1028 lnd=2**(lmark-1) 1029 jnd=floor(float(lnd-1+geno)/float(lnd)) 1030 h(indinc(nbinc+1-lmark),nmknmk+geno)=2-mod(jnd+1,2) 1031 end do 1032 1033 end do 1034 end if 1035 ! 1036 ! cas ou tout est connu 1037 ! 1038 else 1039 ngeno=1 1040 do lmark=1,nbcon 1041 h(indcon(lmark),ngeno)=hcon(lmark) 1042 end do 1043 end if 1044 1045 ! 1046 ! on va calculer, pour les ngeno configurations possibles 1047 ! la probabilite de la phase d'apres la descendance 1048 ! 1049 probmax_t=-1.d6 1050 ! 1051 do 200 geno=1,ngeno 1052 ! 1053 xlvi=0.d0 1054 do 100 ld=1,ndf(ifem) 1055 kd=repdes(ifem,ld) 1056 ! 1057 ! Reorganisation du tableau des probabilites de transmission 1058 1059 do ll=1,nmk(c) 1060 do i=1,4 1061 p(ll,kd,i)=prot(c,ll,kd,i) 1062 end do 1063 if(h(ll,geno).eq.2) then 1064 p(ll,kd,1)=prot(c,ll,kd,2) 1065 p(ll,kd,2)=prot(c,ll,kd,1) 1066 p(ll,kd,3)=prot(c,ll,kd,4) 1067 p(ll,kd,4)=prot(c,ll,kd,3) 1068 end if 1069 end do 1070 ! 1071 ! Initialisation 1072 do i=1,4 1073 q(i)=p(1,kd,i) 1074 end do 1075 ! 1076 ! Calcul de la probabilite du phenotype du descendant 1077 do ll2=2,nmk(c) 1078 ll1=ll2-1 1079 r1=rm(c,ll1,ll2) 1080 r2=rf(c,ll1,ll2) 1081 call transmi(r1,r2,trans) 1082 do t2=1,4 1083 s(t2)=0.d0 1084 do t1=1,4 1085 s(t2)=s(t2)+trans(t1,t2)*q(t1) 1086 end do 1087 end do 1088 do i=1,4 1089 q(i)=p(ll2,kd,i)*s(i) 1090 end do 1091 end do 1092 pr=0.d0 1093 do t=1,4 1094 pr=pr+q(t) 1095 end do 1096 xlvi=xlvi+log(pr) 1097 1098 100 continue 1099 prob(geno)=xlvi 1100 if(probmax_t.le.xlvi) probmax_t=xlvi 1101 1102 200 continue 1103 ! 1104 !****************************************************************************** 1105 1106 ! 1107 ! on rescale les prob(geno) en retranchant pmax e toutes les valeurs 1108 ! ce qui fait que la plus grande vaut 0 et les autres sont augmentees 1109 ! de -pmax. On prend l'exponentielle de ce resultat 1110 ! 1111 sprob=0.d0 1112 do geno=1,ngeno 1113 if (prob(geno).ne.0.d0) then 1114 prob(geno)=dexp(prob(geno)-probmax_t) 1115 sprob=sprob+prob(geno) 1116 end if 1117 end do 1118 1119 do geno=1,ngeno 1120 prob(geno)=prob(geno)/sprob 1121 end do 1122 ! 1123 ! elimination des phases improbables 1124 ! 1125 nphm=0 1126 sprobphm=0.d0 1127 do geno=1,ngeno 1128 if(prob(geno).ge.PRSEUIL)then 1129 nphm=nphm+1 1130 do imark=1,nmk(c) 1131 tphm(imark,nphm)=h(imark,geno) 1132 end do 1133 probphm(nphm)=prob(geno) 1134 sprobphm=sprobphm+prob(geno) 1135 end if 1136 end do 1137 1138 do geno=1,nphm 1139 probphm(geno)=probphm(geno)/sprobphm 1140 end do 1141 300 continue 1142 1143 !****************************************************************************** 1144 ! 1145 ! Impression des resultats 1146 ! 1147 400 continue 1148 1149 if(nphm.eq.0)then 1150 call stop_application('Dam '//trim(femelle(ifem))//' has no possible haplotype '// & 1151 '(none has probability greater than '// trim(str(PRSEUIL))//') : '// & 1152 'check the genotype compatibility in the all sire family') 1153 end if 1154 1155 ! if (corref(ifem).ne.9999) then 1156 ! if (ndf(ifem).ge.ndmin) then 1157 ! nombfem=nombfem+1 1158 ! end if 1159 ! end if 1160 ! 1161 ! stockage des resultats par mere (et non par femelle) 1162 ! et reperage de la phase la plus probable 1163 ! 1164 1165 pmax=0.d0 1166 allocate ( probg_t(nchr,size(probg,2)) ) 1167 probg_t(:,1:size(probg,2)) = probg 1168 1169 deallocate ( probg ) 1170 allocate ( probg(nchr,size(probg_t,2)+nphm) ) 1171 probg=0.d0 1172 probg(:,1:size(probg_t,2)) = probg_t 1173 1174 deallocate(probg_t) 1175 1176 ngenom(c,jm+1)=ngenom(c,jm) 1177 do geno=1,nphm 1178 ngenom(c,jm+1)=ngenom(c,jm+1)+1 1179 ! 1180 ! ajout de genotypm 1181 do ii=1,nmk(c) 1182 if (correm(jm) <= size(genotyp,3).and. corref(ifem) <= size(pheno,3)) then 1183 genotypm(c,ii,ngenom(c,jm+1),1)=pheno(c,ii,corref(ifem),tphm(ii,geno)) 1184 genotypm(c,ii,ngenom(c,jm+1),2)=pheno(c,ii,corref(ifem),3-tphm(ii,geno)) 1185 end if 1186 end do 1187 ! 1188 ! 1189 probg(c,ngenom(c,jm+1))=probphm(geno) 1190 end do 1191 1192 ngeno1=ngenom(c,jm)+1 1193 ngeno2=ngenom(c,jm+1) 1194 igeno=0 1195 1196 allocate ( ngend_t(nchr,size(ngend,2)) ) 1197 ngend_t(:,1:size(ngend,2)) = ngend 1198 deallocate ( ngend ) 1199 allocate ( ngend(nchr,size(ngend_t,2)+(ngeno2-ngeno1)+1) ) 1200 ngend=0.d0 1201 ngend(:,1:size(ngend_t,2))=ngend_t 1202 deallocate(ngend_t) 1203 1204 if ( size(ptfin,3) < maxval(ngend)+(ngeno2-ngeno1+1)*(nd2-nd1+1) ) then 1205 allocate (ptfin_t(nchr,maxval(nmk),size(ptfin,3),4)) 1206 ptfin_t=ptfin 1207 deallocate (ptfin) 1208 !pour eviter de reallouer sans cesse et faire des copies de tableau, on alloue nd valeurs supplementaires 1209 !normalement ca devrait etre : maxval(ngend)+(ngeno2-ngeno1+1)*(nd2-nd1+1) 1210 allocate (ptfin(nchr,maxval(nmk),maxval(ngend)+nd,4)) 1211 ptfin(:,:,1:size(ptfin_t,3),:)=ptfin_t 1212 deallocate (ptfin_t) 1213 end if 1214 1215 do geno=ngeno1,ngeno2 1216 igeno=igeno+1 1217 ngend(c,geno+1)=ngend(c,geno) 1218 allocate ( ndesc_t(nchr,size(ndesc,2)) ) 1219 ndesc_t(:,1:size(ndesc,2)) = ndesc 1220 deallocate ( ndesc ) 1221 allocate ( ndesc(nchr,size(ndesc_t,2)+(nd2-nd1)+1) ) 1222 ndesc=0 1223 ndesc(:,1:size(ndesc_t,2))=ndesc_t 1224 deallocate( ndesc_t ) 1225 do kd=nd1,nd2 1226 ngend(c,geno+1)=ngend(c,geno+1)+1 1227 ndesc(c,ngend(c,geno+1))=kd 1228 do ll=1,nmk(c) 1229 do i=1,4 1230 ptfin(c,ll,ngend(c,geno+1),i)=prot(c,ll,kd,i) 1231 end do 1232 if(tphm(ll,igeno).eq.2) then 1233 ptfin(c,ll,ngend(c,geno+1),1)=prot(c,ll,kd,2) 1234 ptfin(c,ll,ngend(c,geno+1),2)=prot(c,ll,kd,1) 1235 ptfin(c,ll,ngend(c,geno+1),3)=prot(c,ll,kd,4) 1236 ptfin(c,ll,ngend(c,geno+1),4)=prot(c,ll,kd,3) 1237 end if 1238 1239 end do 1240 end do 1241 end do 1242 500 continue 1243 end do 1244 1245 deallocate (p) 1246 if (allocated(h)) deallocate(h) 1247 if (allocated(tphm)) deallocate(tphm) 1248 if (allocated(prob)) deallocate(prob) 1249 if (allocated(probphm)) deallocate(probphm) 1250 call log_mess('end PDEGM_V4',DEBUG_DEF) 1251 end subroutine pdegm
pdegm_snp
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pdegm_snp
DESCRIPTION
NOTES
Determination des probabilites des phases pour les meres
SOURCE
2036 subroutine pdegm_snp() 2037 2038 ! 2039 ! Divers 2040 integer :: geno,igeno 2041 integer(kind=KIND_PHENO) :: m1,m2 2042 integer i,ii,icon,ifem,indmin,jm,kd,kkd,c 2043 integer ll,ll1,ll2,lmark,nbcon,nd1,nd2 2044 integer ngeno1,ngeno2,nmark 2045 integer nombfem,ndf(nm+1),repdes(nm,nd),corref(nm) 2046 integer oppos,phase,phm(maxval(nmk)) 2047 integer trans(nd,maxval(nmk)),hcon(maxval(nmk)),indcon(maxval(nmk)) 2048 logical :: hetero(maxval(nmk)), mafor(maxval(nmk)) 2049 real (kind=dp) :: p1,p2,ratio,ratmin 2050 2051 call log_mess('calling...PDEGM_FAST',DEBUG_DEF) 2052 ! 2053 !****************************************************************************** 2054 ! 2055 ! 2056 ! Creation des tableaux 2057 ! ndf (nombre total de descendants/femelle) 2058 ! repdes (identification des descendants de chaque femelles selon son nuero d'ordre 2059 ! en lecture) 2060 ! corref (l'equivalent de correm pour les femelles) 2061 ! 2062 ! cette partie sera plutot e positionner dans lect_genea et lect_typ 2063 ! 2064 2065 !init OFI : 2066 trans = 0 2067 2068 do ifem=1,nfem 2069 ndf(ifem)=0 2070 end do 2071 2072 do jm=1,nm 2073 nd1=ndm(jm)+1 2074 nd2=ndm(jm+1) 2075 ifem=repfem(jm) 2076 do kd=nd1,nd2 2077 ndf(ifem)=ndf(ifem)+1 2078 repdes(ifem,ndf(ifem))=kd 2079 end do 2080 2081 corref(ifem)=correm(jm) 2082 end do 2083 2084 ! 2085 ! le tableau force indique si une femelle n'a pas assez de descendants 2086 ! pour etre estimee 2087 ! 2088 ! do jm=1,nm 2089 ! force(jm)=.false. 2090 ! end do 2091 ! 2092 allocate ( ngend(nchr,nm+1) ) 2093 allocate ( probg(nchr,nm) ) 2094 allocate ( ndesc(nchr,nd) ) 2095 2096 ngenom(:,1)=0 2097 ngend=0 2098 probg=0.d0 2099 ! 2100 !****************************************************************************** 2101 ! 2102 ! on traite les femelles (et non les meres) une a une 2103 ! 2104 ! nombfem est le nombre de femelles dont on peut chercher les phases 2105 ! 2106 do c=1,nchr 2107 nombfem=0 2108 do 500 jm=1,nm 2109 nd1=ndm(jm)+1 2110 nd2=ndm(jm+1) 2111 ifem=repfem(jm) 2112 ! 2113 !****************************************************************************** 2114 ! 2115 ! Si le nombre de pleins freres est trop faible, le genotype de la mere ne sera pas considere 2116 ! Dans ce cas le premier genotype rencontre possible est considere comme le bon 2117 2118 !if(ndf(ifem).lt.ndmin) then 2119 ! force(jm)=.true. 2120 if(estfem(ifem))then 2121 ! on cherche le premier genotype possible 2122 ! 2123 do nmark=1,nmk(c) 2124 phm(nmark)=1 2125 if(ordref(c,ifem,nmark).eq.21)phm(nmark)=2 2126 end do 2127 2128 go to 400 2129 2130 end if 2131 ! 2132 !****************************************************************************** 2133 ! 2134 ! 2135 ! on cree le tableau hetero qui liste l'heterozygotie des marqueurs 2136 ! 2137 do nmark=1,nmk(c) 2138 hetero(nmark)=.true. 2139 if(corref(ifem).eq.9999) then 2140 m1=nmanque 2141 else 2142 m1=pheno(c,nmark,corref(ifem),1) 2143 m2=pheno(c,nmark,corref(ifem),2) 2144 end if 2145 if(m1.eq.nmanque) m2=m1 2146 if(m1.eq.m2) then 2147 hetero(nmark)=.false. 2148 end if 2149 end do 2150 ! 2151 ! creation du tableau des transmissions observees 2152 ! 2153 kkd=0 2154 do kd=nd1,nd2 2155 kkd=kkd+1 2156 do ll=1,nmk(c) 2157 trans(kkd,ll)=0 2158 p1=prot(c,ll,kd,1)+prot(c,ll,kd,3) 2159 if(p1.eq.1.d0)trans(kkd,ll)=1 2160 p2=prot(c,ll,kd,2)+prot(c,ll,kd,4) 2161 if(p2.eq.1.d0)trans(kkd,ll)=2 2162 end do 2163 end do 2164 ! 2165 ! on commence par lister les marqueurs dont l'origine est connue 2166 ! soit par l'ascendance, soit parcequ'ils sont homozygotes 2167 ! (dans ce cas on donne l'ordrem=12) 2168 ! 2169 2170 nbcon=0 2171 do lmark=1,nmk(c) 2172 2173 mafor(lmark)=.false. 2174 hcon(lmark)=0 2175 phm(lmark)=0 2176 2177 icon=0 2178 if(ordref(c,ifem,lmark).eq.12)then 2179 nbcon=nbcon+1 2180 hcon(lmark)=1 2181 indcon(nbcon)=lmark 2182 mafor(lmark)=.true. 2183 end if 2184 2185 if(ordref(c,ifem,lmark).eq.21)then 2186 nbcon=nbcon+1 2187 hcon(lmark)=2 2188 indcon(nbcon)=lmark 2189 mafor(lmark)=.true. 2190 end if 2191 2192 end do 2193 ! 2194 ! reconstruction de la phase la plus probable par approximations successives 2195 ! en partant du premier marqueur heterozygote chez le pere 2196 ! 2197 do ll1=1,nmk(c) 2198 if(hetero(ll1)) go to 10 2199 end do 2200 ! 2201 ! le pere est entierement homozygote 2202 ! 2203 do lmark=1,nmk(c) 2204 phm(lmark)=1 2205 end do 2206 2207 go to 100 2208 2209 10 phm(ll1)=1 2210 if(nbcon.ne.0)phm(ll1)=hcon(indcon(1)) 2211 ! 2212 ! si le premier marqueur informatif n'est pas le premier, on pose que 2213 ! les phases locales e sa gauche sont identiques 2214 ! 2215 if(ll1.gt.1) then 2216 do ll=ll1-1,1,-1 2217 phm(ll)=phm(ll1) 2218 end do 2219 end if 2220 ! 2221 ! on part du premier marqueur informatif (s'il existe) et on construit 2222 ! les phases locales en glissant vers la droite 2223 ! ll1 et ll2 etant les indices de 2 marqueurs heterozygotes successifs 2224 ! (tous les marqueurs intermediaires sont homozygotes), on compare 2225 ! les transmission en ll1 et ll2, et on en deduit la phase locale 2226 ! 2227 ! en cas d'incoherence entre cette reconstruction et la phase donnee par 2228 ! l'ascendance au prochain marqueur informatif, il faudra inverser le resultat 2229 ! e partir du marqueur (indmin) le moins clair (rapport d'effectifs en phase 2230 ! et en opposition ratmin) le plus proche de 0.5 2231 ! 2232 40 continue 2233 2234 ratmin=0.d0 2235 indmin=ll1 2236 ll2=ll1 2237 50 ll2=ll2+1 2238 2239 if ( ll2 <= nmk(c) ) then 2240 if(hetero(ll2)) go to 60 2241 end if 2242 ! 2243 ! cas du dernier marqueur heterozygote 2244 ! 2245 if(ll2.ge.nmk(c)) then 2246 2247 do ll=ll1+1,nmk(c) 2248 phm(ll)=phm(ll1) 2249 end do 2250 go to 100 2251 end if 2252 2253 go to 50 2254 60 continue 2255 ! 2256 ! comptage des inversions apparentes des phases locales 2257 ! 2258 phase=0 2259 oppos=0 2260 do kkd=1,ndf(ifem) 2261 if(trans(kkd,ll1).ne.0.and.trans(kkd,ll2).ne.0) then 2262 if(trans(kkd,ll1).eq.trans(kkd,ll2))then 2263 phase=phase+1 2264 else 2265 oppos=oppos+1 2266 end if 2267 end if 2268 end do 2269 if((phase+oppos).eq.0)go to 50 2270 ! 2271 ! etablissement des phases locales 2272 ! 2273 if(phase.gt.oppos)then 2274 do ll=ll1+1,ll2 2275 phm(ll)=phm(ll1) 2276 end do 2277 ratio=dble(oppos)/dble(phase) 2278 else 2279 do ll=ll1+1,ll2 2280 phm(ll)=3-phm(ll1) 2281 end do 2282 ratio=dble(phase)/dble(oppos) 2283 end if 2284 2285 if(ratio.ge.ratmin) then 2286 ratmin=ratio 2287 indmin=ll2 2288 end if 2289 2290 2291 ll1=ll2 2292 if(ll2.eq.nmk(c)) go to 100 2293 if(.not.mafor(ll2)) go to 50 2294 ! 2295 ! test de la coherence entre phase reconstruite d'apres la descendance 2296 ! et information de l'ascendance 2297 ! en cas d'incoherence on inverse e partir de indmin 2298 ! 2299 if(phm(ll2).ne.hcon(ll2)) then 2300 do ll = indmin,ll2 2301 phm(ll)=3-phm(ll) 2302 end do 2303 go to 40 2304 end if 2305 go to 50 2306 ! 2307 100 continue 2308 2309 400 continue 2310 2311 ! 2312 ! stockage des resultats par mere (et non par femelle) 2313 ! et reperage de la phase la plus probable 2314 ! 2315 ngenom(c,jm+1)=ngenom(c,jm)+1 2316 ! 2317 ! 2318 probg(c,ngenom(c,jm+1))=1.d0 2319 2320 do ii=1,nmk(c) 2321 if (correm(jm) <= size(genotypm,3).and. corref(ifem) <= size(pheno,3)) then 2322 genotypm(c,ii,ngenom(c,jm+1),1)=pheno(c,ii,corref(ifem),phm(ii)) 2323 genotypm(c,ii,ngenom(c,jm+1),2)=pheno(c,ii,corref(ifem),3-phm(ii)) 2324 2325 if ( correm(jm) > size(genotypm,3) ) then 2326 call log_mess('dam:'//trim(mere(jm))//' have no genotype.',WARNING_DEF) 2327 end if 2328 if ( corref(ifem) > size(pheno,3) ) then 2329 call log_mess('dam:'//trim(mere(jm))//' have no genotype.',WARNING_DEF) 2330 2331 exit 2332 end if 2333 2334 end if 2335 end do 2336 2337 2338 ngeno1=ngenom(c,jm)+1 2339 ngeno2=ngenom(c,jm+1) 2340 igeno=0 2341 do geno=ngeno1,ngeno2 2342 igeno=igeno+1 2343 ngend(c,geno+1)=ngend(c,geno) 2344 do kd=nd1,nd2 2345 ngend(c,geno+1)=ngend(c,geno+1)+1 2346 ndesc(c,ngend(c,geno+1))=kd 2347 do ll=1,nmk(c) 2348 do i=1,4 2349 ptfin(c,ll,ngend(c,geno+1),i)=prot(c,ll,kd,i) 2350 end do 2351 if(phm(ll).eq.2) then 2352 ptfin(c,ll,ngend(c,geno+1),1)=prot(c,ll,kd,2) 2353 ptfin(c,ll,ngend(c,geno+1),2)=prot(c,ll,kd,1) 2354 ptfin(c,ll,ngend(c,geno+1),3)=prot(c,ll,kd,4) 2355 ptfin(c,ll,ngend(c,geno+1),4)=prot(c,ll,kd,3) 2356 end if 2357 2358 end do 2359 end do 2360 end do 2361 2362 500 continue 2363 end do 2364 call log_mess('end PDEGM_FAST',DEBUG_DEF) 2365 return 2366 end subroutine pdegm_snp
pdegp
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pdegp
DESCRIPTION
NOTES
Determination de la phase la plus probable pour les peres
Tableaux Arbitraire des 4 cas de transmission
Male \ Femelle Chr1 Chr2
Chr1 1 2
Chr2 3 4
Utilisation des structures de donnees
prot(c,desc,ll,cas) : probabilite de transmission du marqueur ll pour le descenant desc
du Cas 'cas' (1,2,3 ou 4)
SOURCE
399 subroutine pdegp() 400 ! Divers 401 integer geno,t,tder 402 integer(kind=KIND_PHENO) :: m1,m2 403 404 integer i,ii,imark,ip,j,jm,jmax1,jnd,kd,kkd 405 integer linf,ll,ll1,ll2,lmark,lnd,alloc_stat,c 406 integer trans(nd,maxval(nmk)),hcon(maxval(nmk)),indcon(maxval(nmk)),indinc(maxval(nmk)) 407 integer comptrans(maxval(nmk)) 408 integer nbcon,nbinc,nd1,nd2,ngeno,nhomo,nlvi,nm1,nm2 409 integer nmark,nmkinf,nmknmk,nphp,value_alloc 410 real (kind=dp) :: p1,p2,pmax,pr,probmax_t,sprob,xlvi 411 logical :: hetero(maxval(nmk)) 412 !Les valeurs sont [1,2] un octet suffit 413 integer(kind=1), dimension (:,:),allocatable :: h,tphp 414 real (kind=dp) , dimension (:),allocatable ::prob,probphp 415 416 call log_mess('calling...PDEGP_V4',DEBUG_DEF) 417 ! 418 ! Etablissement de la phase la plus probable 419 ! 420 do c=1,nchr 421 do 300 ip=1, np 422 nm1=nmp(ip)+1 423 nm2=nmp(ip+1) 424 ! 425 ! on cree le tableau hetero qui liste l'heterozygotie des marqueurs 426 ! 427 nmark=nmk(c) 428 nhomo=0 429 430 do imark=1,nmk(c) 431 hetero(imark)=.true. 432 if(correp(ip).eq.9999) then 433 m1=nmanque 434 else 435 m1=pheno(c,imark,correp(ip),1) 436 m2=pheno(c,imark,correp(ip),2) 437 end if 438 if(m1.eq.nmanque) m2=m1 439 if(m1.eq.m2) then 440 nhomo=nhomo+1 ! le pere est homozygote au marqueur imark 441 hetero(imark)=.false. 442 end if 443 end do 444 ! 445 ! tableau de comptage des transmissions connues 446 ! 447 do ll=1,nmk(c) 448 comptrans(ll)=0 449 end do 450 ! 451 ! creation du tableau des transmissions observees 452 ! 453 kkd=0 454 do jm=nm1,nm2 455 nd1=ndm(jm)+1 456 nd2=ndm(jm+1) 457 do kd=nd1,nd2 458 kkd=kkd+1 459 do ll=1,nmk(c) 460 trans(kkd,ll)=0 461 p1=prot(c,ll,kd,1)+prot(c,ll,kd,2) 462 if(p1.eq.1.d0)trans(kkd,ll)=1 ! recu 1er Chr du pere 463 p2=prot(c,ll,kd,3)+prot(c,ll,kd,4) 464 if(p2.eq.1.d0)trans(kkd,ll)=2 ! recu 2eme Chr du pere 465 if(trans(kkd,ll).ne.0)comptrans(ll)=comptrans(ll)+1 466 end do 467 end do 468 end do 469 ! 470 ! comptage des marqueurs informatifs 471 ! 472 nmkinf=0 473 do ll=1,nmk(c) 474 if(comptrans(ll).ne.0)nmkinf=nmkinf+1 475 end do 476 ! 477 ! creation du tableau h (anciennement dans combine) 478 ! en se restreignant aux phases utiles 479 ! 480 ! on commence par lister les marqueurs dont l'origine est connue 481 ! soit par l'ascendance, soit parcequ'ils sont homozygotes 482 ! (dans ce cas on donne l'ordrep=12) 483 ! 484 485 nbcon=0 486 nbinc=0 487 do lmark=1,nmk(c) 488 489 !1er allele paternel / 2eme allele maternel 490 if(ordrep(c,ip,lmark).eq.12)then 491 nbcon=nbcon+1 492 hcon(nbcon)=1 493 indcon(nbcon)=lmark 494 !1er allele maternel / 2eme allele paternel 495 else if(ordrep(c,ip,lmark).eq.21)then 496 nbcon=nbcon+1 497 hcon(nbcon)=2 498 indcon(nbcon)=lmark 499 !Le pere est homozygote sur ce marqueur 500 else if(.not.hetero(lmark))then 501 nbcon=nbcon+1 502 hcon(nbcon)=1 503 indcon(nbcon)=lmark 504 else 505 nbinc=nbinc+1 506 indinc(nbinc)=lmark 507 end if 508 509 end do 510 511 !--------------------- ALLOC DYNAMIQUE --------------------- 512 ! ---- A FAIRE VALIDER PAR JM ---------------------------- 513 if (allocated(h)) deallocate(h) 514 if (allocated(tphp)) deallocate(tphp) 515 if (allocated(prob)) deallocate(prob) 516 if (allocated(probphp)) deallocate(probphp) 517 518 if (nbinc>0) then 519 if ( (nbinc*2) > MAX_MARKER) then 520 call stop_application("You can not used this version of haplotype calculation. Please use --snp option") 521 end if 522 523 value_alloc= (2**(nbinc-1))*2 524 allocate (h(nmk(c),value_alloc),STAT=alloc_stat) 525 allocate (tphp(nmk(c),value_alloc),STAT=alloc_stat) 526 allocate (prob(value_alloc),STAT=alloc_stat) 527 allocate (probphp(value_alloc),STAT=alloc_stat) 528 else 529 allocate (h(nmk(c),1),STAT=alloc_stat) 530 allocate (tphp(nmk(c),1),STAT=alloc_stat) 531 allocate (prob(1),STAT=alloc_stat) 532 allocate (probphp(1),STAT=alloc_stat) 533 end if 534 535 if (alloc_stat/=0) then 536 call stop_application('Not enough memory to compute '// & 537 'transmission probabilities. Deacrease number of marker.') 538 end if 539 540 !--------------------------- END OFI --------------------- 541 542 !***************** WARNING 543 !***************** WARNING 544 !***************** WARNING 545 !***************** WARNING 546 h = 1 ! AJOUT OFI, sinon il y a des cas ou h n est pas initialise et ca fait planter l application.... 547 !***************** WARNING 548 !***************** WARNING 549 !***************** WARNING 550 !***************** WARNING 551 ! 552 ! traitement des cas limites (nmk(c) = 0 ou 1 marqueur) 553 ! 554 if(nmk(c).eq.0.or.nmkinf.eq.0)then 555 call stop_application(' Sire '//trim(pere(ip))// & 556 ' had no informative progeny '//& 557 ' for any of the markers'//' It should not be'//& 558 ' included in the analysis') 559 end if 560 561 562 if(nmk(c).eq.1.or.nmkinf.eq.1) then 563 ngeno=1 564 prob(1)=1.d0 565 if(nbcon.eq.1)then 566 h(1,1)=hcon(1) 567 else 568 h(1,1)=1 569 end if 570 go to 200 571 end if 572 573 574 ! 575 ! s'il reste des inconnues 576 ! on remplit le tableau h, en completant les elements connus issus 577 ! de indcon et hcon par les phases possibles des elements inconnus 578 ! 579 if(nbinc.ne.0) then 580 nmknmk=2**(nbinc-1) 581 ngeno=nmknmk 582 583 do geno=1,nmknmk 584 do lmark=1,nbcon 585 h(indcon(lmark),geno)=hcon(lmark) 586 end do 587 588 do lmark=1,nbinc 589 lnd=2**(lmark-1) 590 jnd=floor(float(lnd-1+geno)/float(lnd)) 591 h(indinc(nbinc+1-lmark),geno)=1+mod(jnd+1,2) 592 end do 593 594 end do 595 596 597 ! 598 ! quand le genotype est oriente (phasp(ip) e true), il faut aussi calculer les 599 ! probabilites des phases miroirs 600 ! 601 602 if(phasp(c,ip)) then 603 ngeno=ngeno+nmknmk 604 do geno=1,nmknmk 605 606 do lmark=1,nbcon 607 h(indcon(lmark),nmknmk+geno)=hcon(lmark) 608 end do 609 610 do lmark=1,nbinc 611 lnd=2**(lmark-1) 612 jnd=floor(float(lnd-1+geno)/float(lnd)) 613 h(indinc(nbinc+1-lmark),nmknmk+geno)=2-mod(jnd+1,2) 614 end do 615 end do 616 end if 617 ! 618 ! cas ou tout est connu 619 ! 620 else 621 ngeno=1 622 do lmark=1,nbcon 623 h(indcon(lmark),ngeno)=hcon(lmark) 624 end do 625 end if 626 ! 627 ! traitement des cas limites (0 ou 1 marqueur heterozygote) 628 ! 629 630 if(nmk(c)-nhomo.le.1)then 631 ngeno=1 632 prob(1)=1.d0 633 go to 200 634 end if 635 probmax_t=-1.d6 636 ! 637 ! 638 do 100 geno=1,ngeno 639 ! 640 ! calcul de la probabilite de la phase d'apres la descendance 641 ! 642 xlvi=0.d0 643 kkd=0 644 do jm=nm1,nm2 645 nd1=ndm(jm)+1 646 nd2=ndm(jm+1) 647 nlvi=0 648 do kd=nd1,nd2 649 kkd=kkd+1 650 ! 651 ! Calcul de la probabilite du phenotype du descendant 652 ll1=1 653 do while (trans(kkd,ll1).eq.0.and.ll1.lt.nmk(c)) 654 ll1=ll1+1 655 end do 656 tder=trans(kkd,ll1) 657 if (tder.ne.0) then 658 if (h(ll1,geno).eq.2) tder=3-tder 659 linf=ll1+1 660 do ll2=linf,nmk(c) 661 t=trans(kkd,ll2) 662 pr=0.d0 663 if (t.ne.0) then 664 if (h(ll2,geno).eq.2) t=3-t 665 if (t.eq.tder) then 666 pr=1.d0-rm(c,ll1,ll2) 667 else 668 pr=rm(c,ll1,ll2) 669 end if 670 tder=t 671 ll1=ll2 672 xlvi=xlvi+dlog(pr) 673 end if 674 end do 675 end if 676 677 end do 678 end do 679 680 prob(geno)=xlvi 681 if(probmax_t.le.xlvi) probmax_t=xlvi 682 683 100 continue 684 685 ! 686 ! Standardisation 687 ! 688 sprob=0.d0 689 do geno=1,ngeno 690 if (prob(geno).ne.0.d0) then 691 prob(geno)=dexp(prob(geno)-probmax_t) 692 sprob=sprob+prob(geno) 693 end if 694 end do 695 696 !AJOUT OFI : sprob peut valoir 0 697 if ( sprob /= 0 ) then 698 do geno=1,ngeno 699 prob(geno)=prob(geno)/sprob 700 end do 701 end if 702 ! 703 ! elimination des phases improbables 704 ! 705 200 continue 706 nphp=0 707 do geno=1,ngeno 708 if(prob(geno).ge.PHPSEUIL)then 709 nphp=nphp+1 710 do imark=1,nmark 711 tphp(imark,nphp)=h(imark,geno) 712 end do 713 probphp(nphp)=prob(geno) 714 end if 715 end do 716 717 718 if(nphp.eq.0) then 719 call stop_application(' Sire '//trim(pere(ip))// & 720 ' had no phase with a probability above the threshold '//trim(str(PHPSEUIL))// & 721 ' The model is not appropriate to your data We suggest the use of the linealised model') 722 end if 723 724 ! 725 ! cas des inidvidus totalement homozygotes 726 ! 727 728 if(.not.phasp(c,ip).and.nhomo.eq.nmk(c))then 729 prob(1)=1.d0 730 do i=2,ngeno 731 prob(i)=0.d0 732 end do 733 end if 734 735 ! 736 ! Recherche de la phase la plus probable 737 pmax=0.d0 738 do geno=1,nphp 739 if(pmax.lt.probphp(geno)) then 740 pmax=probphp(geno) 741 jmax1=geno 742 end if 743 end do 744 ! 745 ! stockage du genotype 746 ! 747 do i=1,nmk(c) 748 genotyp(c,i,correp(ip),1)=pheno(c,i,correp(ip),tphp(i,jmax1)) 749 genotyp(c,i,correp(ip),2)=pheno(c,i,correp(ip),3-tphp(i,jmax1)) 750 end do 751 ! 752 ! Reorganisation du tableau des probabilites de transmission 753 do jm=nm1,nm2 754 nd1=ndm(jm)+1 755 nd2=ndm(jm+1) 756 do kd=nd1,nd2 757 do ll=1,nmk(c) 758 if(tphp(ll,jmax1).eq.2) then 759 p1=prot(c,ll,kd,1) 760 p2=prot(c,ll,kd,2) 761 prot(c,ll,kd,1)=prot(c,ll,kd,3) 762 prot(c,ll,kd,2)=prot(c,ll,kd,4) 763 prot(c,ll,kd,3)=p1 764 prot(c,ll,kd,4)=p2 765 end if 766 end do 767 end do 768 end do 769 300 continue 770 end do 771 772 if (allocated(h)) deallocate (h) 773 if (allocated(tphp)) deallocate(tphp) 774 if (allocated(prob)) deallocate (prob) 775 if (allocated(probphp)) deallocate (probphp) 776 call log_mess('end PDEGP_V4',DEBUG_DEF) 777 end subroutine pdegp
pdegp_snp
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
pdegp_snp
DESCRIPTION
NOTES
Determination de la phase la plus probable pour les peres
SOURCE
1799 subroutine pdegp_snp() 1800 ! Divers 1801 integer(kind=KIND_PHENO) :: m1,m2 1802 1803 integer ii,imark,ip,indmin,jm,kd,kkd,c 1804 integer ll,ll1,ll2,lmark 1805 integer trans(nd,maxval(nmk)),hcon(maxval(nmk)),indcon(maxval(nmk)) 1806 integer oppos,phase,php(maxval(nmk)) 1807 integer nbcon,nd1,nd2,nm1,nm2 1808 real (kind=dp) :: p1,p2,ratio,ratmin 1809 logical :: hetero(maxval(nmk)), mafor(maxval(nmk)) 1810 1811 call log_mess('calling...PDEGP_FAST',DEBUG_DEF) 1812 1813 ! 1814 ! Etablissement de la phase la plus probable 1815 ! 1816 do c=1,nchr 1817 do 300 ip=1, np 1818 nm1=nmp(ip)+1 1819 nm2=nmp(ip+1) 1820 ! 1821 ! on cree le tableau hetero qui liste l'heterozygotie des marqueurs 1822 ! 1823 do imark=1,nmk(c) 1824 hetero(imark)=.true. 1825 if(correp(ip).eq.9999) then 1826 m1=nmanque 1827 else 1828 m1=pheno(c,imark,correp(ip),1) 1829 m2=pheno(c,imark,correp(ip),2) 1830 end if 1831 if(m1.eq.nmanque) m2=m1 1832 if(m1.eq.m2) then 1833 hetero(imark)=.false. 1834 end if 1835 end do 1836 ! 1837 ! on reduit le jeu de marqueurs e ceux pour lesquels le pere est heterozygote 1838 ! 1839 ! creation du tableau des transmissions observees 1840 ! 1841 kkd=0 1842 do jm=nm1,nm2 1843 nd1=ndm(jm)+1 1844 nd2=ndm(jm+1) 1845 do kd=nd1,nd2 1846 kkd=kkd+1 1847 do ll=1,nmk(c) 1848 trans(kkd,ll)=0 1849 p1=prot(c,ll,kd,1)+prot(c,ll,kd,2) 1850 if(p1.eq.1.d0)trans(kkd,ll)=1 1851 p2=prot(c,ll,kd,3)+prot(c,ll,kd,4) 1852 if(p2.eq.1.d0)trans(kkd,ll)=2 1853 end do 1854 end do 1855 end do 1856 ! 1857 ! Liste des marqueurs dont l'origine est connue par l'ascendance 1858 ! 1859 nbcon=0 1860 do lmark=1,nmk(c) 1861 1862 mafor(lmark)=.false. 1863 hcon(lmark)=0 1864 php(lmark)=0 1865 1866 if(ordrep(c,ip,lmark).eq.12)then 1867 nbcon=nbcon+1 1868 hcon(lmark)=1 1869 indcon(nbcon)=lmark 1870 mafor(lmark)=.true. 1871 end if 1872 1873 if(ordrep(c,ip,lmark).eq.21)then 1874 nbcon=nbcon+1 1875 hcon(lmark)=2 1876 indcon(nbcon)=lmark 1877 mafor(lmark)=.true. 1878 end if 1879 1880 end do 1881 ! 1882 ! reconstruction de la phase la plus probable par approximations successives 1883 ! en partant du premier marqueur heterozygote chez le pere 1884 ! 1885 do ll1=1,nmk(c) 1886 if(hetero(ll1)) go to 10 1887 end do 1888 ! 1889 ! le pere est entierement homozygote 1890 ! 1891 do lmark=1,nmk(c) 1892 php(lmark)=1 1893 end do 1894 1895 go to 100 1896 ! 1897 10 php(ll1)=1 1898 if(nbcon.ne.0)php(ll1)=hcon(indcon(1)) 1899 ! 1900 ! si le premier marqueur informatif n'est pas le premier, on pose que 1901 ! les phases locales e sa gauche sont identiques 1902 ! 1903 if(ll1.gt.1) then 1904 do ll=ll1-1,1,-1 1905 php(ll)=php(ll1) 1906 end do 1907 end if 1908 ! 1909 ! on part du premier marqueur informatif (s'il existe) et on construit 1910 ! les phases locales en glissant vers la droite 1911 ! ll1 et ll2 etant les indices de 2 marqueurs heterozygotes successifs 1912 ! (tous les marqueurs intermediaires sont homozygotes), on compare 1913 ! les transmission en ll1 et ll2, et on en deduit la phase locale 1914 ! 1915 ! en cas d'incoherence entre cette reconstruction et la phase donnee par 1916 ! l'ascendance au prochain marqueur informatif, il faudra inverser le resultat 1917 ! e partir du marqueur (indmin) le moins clair (rapport d'effectifs en phase 1918 ! et en opposition ratmin) le plus proche de 0.5 1919 ! 1920 40 continue 1921 1922 ratmin=0.d0 1923 indmin=ll1 1924 ll2=ll1 1925 50 ll2=ll2+1 1926 if ( ll2 <= nmk(c) ) then 1927 if(hetero(ll2)) go to 60 1928 end if 1929 ! 1930 ! cas du dernier marqueur heterozygote 1931 ! 1932 if(ll2.ge.nmk(c)) then 1933 1934 do ll=ll1+1,nmk(c) 1935 php(ll)=php(ll1) 1936 end do 1937 go to 100 1938 end if 1939 1940 go to 50 1941 60 continue 1942 ! 1943 ! comptage des inversions apparentes des phases locales 1944 ! 1945 phase=0 1946 oppos=0 1947 do kkd=1,ndm(nm2+1)-ndm(nm1) 1948 if(trans(kkd,ll1).ne.0.and.trans(kkd,ll2).ne.0) then 1949 if(trans(kkd,ll1).eq.trans(kkd,ll2))then 1950 phase=phase+1 1951 else 1952 oppos=oppos+1 1953 end if 1954 end if 1955 end do 1956 if((phase+oppos).eq.0)go to 50 1957 ! 1958 ! etablissement des phases locales 1959 ! 1960 if(phase.gt.oppos)then 1961 do ll=ll1+1,ll2 1962 php(ll)=php(ll1) 1963 end do 1964 ratio=dble(oppos)/dble(phase) 1965 else 1966 do ll=ll1+1,ll2 1967 php(ll)=3-php(ll1) 1968 end do 1969 ratio=dble(phase)/dble(oppos) 1970 end if 1971 1972 if(ratio.ge.ratmin) then 1973 ratmin=ratio 1974 indmin=ll2 1975 end if 1976 1977 ll1=ll2 1978 if(ll2.eq.nmk(c)) go to 100 1979 if(.not.mafor(ll2)) go to 50 1980 ! 1981 ! test de la coherence entre phase reconstruite d'apres la descendance 1982 ! et information de l'ascendance 1983 ! en cas d'incoherence on inverse e partir de indmin 1984 ! 1985 if(php(ll2).ne.hcon(ll2)) then 1986 do ll = indmin,ll2 1987 php(ll)=3-php(ll) 1988 end do 1989 go to 40 1990 end if 1991 go to 50 1992 ! 1993 100 continue 1994 1995 ! 1996 ! stockage du genotype 1997 ! 1998 do ll=1,nmk(c) 1999 genotyp(c,ll,correp(ip),1)=pheno(c,ll,correp(ip),php(ll)) 2000 genotyp(c,ll,correp(ip),2)=pheno(c,ll,correp(ip),3-php(ll)) 2001 end do 2002 2003 ! 2004 ! Reorganisation du tableau des probabilites de transmission 2005 do jm=nm1,nm2 2006 nd1=ndm(jm)+1 2007 nd2=ndm(jm+1) 2008 do kd=nd1,nd2 2009 do ll=1,nmk(c) 2010 if(php(ll).eq.2) then 2011 p1=prot(c,ll,kd,1) 2012 p2=prot(c,ll,kd,2) 2013 prot(c,ll,kd,1)=prot(c,ll,kd,3) 2014 prot(c,ll,kd,2)=prot(c,ll,kd,4) 2015 prot(c,ll,kd,3)=p1 2016 prot(c,ll,kd,4)=p2 2017 end if 2018 end do 2019 end do 2020 end do 2021 300 continue 2022 end do 2023 call log_mess('end PDEGP_FAST',DEBUG_DEF) 2024 return 2025 end subroutine pdegp_snp
probin
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
probin
DESCRIPTION
NOTES
probin calcule la fonction de repartition d une binomial de parametres n et r pp est la probabilite de trouver au moins m individus recombinants parmi n pour un tqux de recombinaison r
SOURCE
4615 subroutine probin(n,r,m,pp) 4616 double precision r,pp,c,do 4617 integer n,m,l 4618 4619 if(r**n == 0.d0) then 4620 if(m == 0) then 4621 pp = 1.d0 4622 else 4623 pp = 0.d0 4624 end if 4625 return 4626 end if 4627 4628 if(r == 1.d0) then 4629 if(m == n) then 4630 pp = 1.d0 4631 else 4632 pp = 0.d0 4633 end if 4634 return 4635 end if 4636 4637 c=1.d0 4638 do=r**n 4639 pp=c*do 4640 do l=n,m,-1 4641 c=c * dble(l) / dble(n-l+1) 4642 do=do*(1.d0-r)/ r 4643 pp=pp+c*do 4644 end do 4645 end subroutine probin
set_allocation
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
set_allocation
DESCRIPTION
NOTES
SOURCE
114 subroutine set_allocation 115 116 integer :: stat 117 118 allocate (ordrep(nchr,np,maxval(nmk)),STAT=stat) 119 call check_allocate(stat,'ordrep [m_qtlmap_haplotype_V2]') 120 allocate (ordrem(nchr,nm,maxval(nmk)),STAT=stat) 121 call check_allocate(stat,'ordrem [m_qtlmap_haplotype_V2]') 122 allocate (ordre(nchr,nr,maxval(nmk)),STAT=stat) 123 call check_allocate(stat,'ordre [m_qtlmap_haplotype_V2]') 124 allocate (ordref(nchr,nm,maxval(nmk)),STAT=stat) 125 call check_allocate(stat,'ordref [m_qtlmap_haplotype_V2]') 126 allocate (prot(nchr,maxval(nmk),nd,4),STAT=stat) 127 call check_allocate(stat,'prot [m_qtlmap_haplotype_V2]') 128 allocate (ptfin(nchr,maxval(nmk),nd,4),STAT=stat) 129 call check_allocate(stat,'ptfin [m_qtlmap_haplotype_V2]') 130 131 ! a voir allocation genotyp et genotypm...... : la 3eme dim est approximatif... 132 133 allocate (genotyp(nchr,maxval(nmk),(size(numero)),2),STAT=stat) 134 call check_allocate(stat,'genotyp [m_qtlmap_haplotype_V2]') 135 allocate (genotypm(nchr,maxval(nmk),(size(numero)),2),STAT=stat) 136 call check_allocate(stat,'genotypm [m_qtlmap_haplotype_V2]') 137 138 genotyp=nmanque 139 allocate (ngenom(nchr,nm+1),STAT=stat) 140 call check_allocate(stat,'ngenom [m_qtlmap_haplotype_V2]') 141 allocate (phasp(nchr,np),STAT=stat) 142 call check_allocate(stat,'phasp [m_qtlmap_haplotype_V2]') 143 allocate (phasm(nchr,nm),STAT=stat) 144 call check_allocate(stat,'phasm [m_qtlmap_haplotype_V2]') 145 146 end subroutine set_allocation
set_phasm
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
set_phasm
DESCRIPTION
NOTES
SOURCE
4452 subroutine set_phasm() 4453 integer :: ll,im,c 4454 4455 if (size(phasm,2) /= size(ordrem,2)) then 4456 call stop_application("Devel error: using set_phasm with dim(phasm) :dim(ordrem,2)") 4457 end if 4458 phasm=.false. 4459 do c=1,nchr 4460 do ll=1,nmk(c) 4461 do im=1,nm 4462 if(ordrem(c,im,ll) /= 0)phasm(c,im)=.true. 4463 end do 4464 end do 4465 end do 4466 end subroutine set_phasm
set_phasp
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
set_phasp
DESCRIPTION
NOTES
SOURCE
4424 subroutine set_phasp() 4425 4426 integer :: ll,ip,c 4427 4428 if (size(phasp,2) /= size(ordrep,2)) then 4429 call stop_application("Devel error: using set_phasp with dim(phasp) :dim(ordrep,2)") 4430 end if 4431 phasp=.false. 4432 do c=1,nchr 4433 do ll=1,nmk(c) 4434 do ip=1,np 4435 if(ordrep(c,ip,ll) /= 0 ) then 4436 phasp(c,ip)=.true. 4437 end if 4438 end do 4439 end do 4440 end do 4441 end subroutine set_phasp
setting_alloc_pdd
[ Top ] [ m_qtlmap_haplotype_V2 ] [ Subroutines ]
NAME
setting_alloc_pdd
DESCRIPTION
NOTES
SOURCE
177 subroutine setting_alloc_pdd 178 integer :: stat,c 179 integer ,dimension(nchr) :: valnpo 180 181 call log_mess('Second dim of pdd:'//trim(str(maxval(ngend))),DEBUG_DEF) 182 do c=1,nchr 183 valnpo(c)=get_npo(c) 184 end do 185 allocate( pdd(nchr,maxval(ngend),4,maxval(valnpo)) ) 186 pdd=0.d0 187 end subroutine setting_alloc_pdd