m_qtlmap_haplotype_ldla
NAME
m_qtlmap_haplotype_ldla
DESCRIPTION
NOTES
SEE ALSO
comp_haplo
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
comp_haplo
DESCRIPTION
DIMENSIONS
1: 2:
NOTES
SOURCE
182 logical , dimension (:,:) , allocatable :: comp_haplo
count_haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
count_haplo_complet
DESCRIPTION
DIMENSIONS
1:
NOTES
SOURCE
131 double precision , dimension (:) , allocatable :: count_haplo_complet
gamete
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
gamete
DESCRIPTION
DIMENSIONS
1: 1,2 2: marker (1 <= <= longhap) 3: 2**longhap
NOTES
gamete(1,lk,j_gam) l allele paternel au lk eme marqueur pour le j_gam eme gamete gamete(2,lk,j_gam) l allele maternel au lk eme marqueur pour le j_gam eme gamete
SOURCE
94 integer(kind=KIND_PHENO), dimension (:,:,:) , allocatable :: gamete
haplo
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
haplo
DESCRIPTION
list of the possible haplotype (with unknwon allele) at the current position
DIMENSIONS
1: index of the haplotype 2: index of the allele (1 < index allele < longap)
NOTES
SOURCE
40 integer(kind=KIND_PHENO), dimension (:,:) , allocatable :: haplo
haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
haplo_complet
DESCRIPTION
list of the possible haplotype (without unknwon allele) at the current position
DIMENSIONS
1: index of the complete haplotype 2: index of the allele (1 < index allele < longhap)
NOTES
see liste_haplo_complet
SOURCE
66 integer(kind=KIND_PHENO), dimension (:,:) , allocatable :: haplo_complet
loc_haplo
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
loc_haplo
DESCRIPTION
!-----!
A A A T T T A A A A T T T T A T A T A T A T T T T
| | | |
| | | ---------------
| | ------------- |
| ------- | |
loc_haplo(1) | | |
loc_haplo(2) | |
loc_haplo(3) |
loc_haplo(4)
DIMENSIONS
1: index of the marker (size:longhap)
NOTES
loc_haplo, le vecteur de dimension longhap qui contient les numeros des marqueurs pour l haplotype entourant la position testee dx, et lk le premier marqueur flanquant a droite
SOURCE
119 integer , dimension (:) , allocatable :: loc_haplo
nb_haplo_possible
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
nb_haplo_possible
DESCRIPTION
DIMENSIONS
1:
NOTES
SOURCE
170 integer , dimension (:) , allocatable :: nb_haplo_possible
pb_haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
pb_haplo_complet
DESCRIPTION
DIMENSIONS
1: 2: 3:
NOTES
SOURCE
145 double precision , dimension (:) , allocatable :: pb_haplo_complet
race_h
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
race_h
DESCRIPTION
breed origin of haplotype haplo
DIMENSIONS
1: index of the haplotype
NOTES
SOURCE
52 character(len=LEN_DEF), dimension (:) , allocatable :: race_h
race_haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
race_haplo_complet
DESCRIPTION
breed origin of haplotype haplo_complet
DIMENSIONS
1: index of the complete haplotype
NOTES
SOURCE
78 character(len=LEN_DEF), dimension (:) , allocatable :: race_haplo_complet
tab_haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Variables ]
NAME
tab_haplo_complet
DESCRIPTION
DIMENSIONS
1: 2:
NOTES
SOURCE
158 integer(kind=KIND_PHENO), dimension (:,:) , allocatable :: tab_haplo_complet
bilan_shared_haplo
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
bilan_shared_haplo
DESCRIPTION
NOTES
bilan_shared_haplo informe sur le nombre et les carateristiques des haplotypes paternels retrouves chez les meres Utilise nmk,np,posi,correp,nm de m_qtlmap_data Utilise genotyp de m_qtlmap_haplotype_data
SOURCE
2191 subroutine bilan_shared_haplo(c,dx,n,ip) 2192 2193 ! 2194 integer , intent(in) :: c,ip 2195 real(kind=dp) , intent(in) :: dx 2196 integer, intent (in) :: n 2197 integer itab,jtab,ktab,stat,imax1,imax2 2198 integer , dimension (:,:) , allocatable :: tab 2199 ! allocation 2200 ! 2201 allocate (tab(2,maxval(nmk)),STAT=stat) 2202 call check_allocate(stat,'tab [m_qtlmap_haplotype_LD]') 2203 2204 if(n ==1 .and. ip ==1)then 2205 write(nficout,20) ip 2206 20 format(' Dam carried Sire ',i6,' haplotypes length distribution (grouped by 20 SNP)') 2207 write(nficout,30) 2208 30 format(' position * sire * 1st haplotype / 2nd haplotype'//72('*')//& 2209 ' * * 10 20 30 40 50 60 70 80 90 100 110 120',& 2210 ' 130 140 150 160 170 180 190 200.....'// 72('*')) 2211 end if 2212 2213 tab=0 2214 jtab=0 2215 ktab=1 2216 do itab =1,nmk(c) 2217 jtab=jtab+1 2218 tab(1,ktab)=tab(1,ktab)+ tab_shared_haplo(ip,1,1+itab) 2219 tab(2,ktab)=tab(2,ktab)+ tab_shared_haplo(ip,2,1+itab) 2220 if(jtab == 10) then 2221 jtab=0 2222 ktab=ktab+1 2223 end if 2224 end do 2225 do imax1 =ktab,1,-1 2226 if(tab(1,imax1) /=0) exit 2227 end do 2228 do imax2 =ktab,1,-1 2229 if(tab(2,imax2) /=0) exit 2230 end do 2231 2232 write(nficout,31) dx,ip,(tab(1,itab),itab=1,imax1) 2233 31 format(3x,f5.3,' *',i8,'*',30(1x,i3)) 2234 write(nficout,32) (tab(2,itab),itab=1,imax2) 2235 32 format(3x,5x,11x,'*',30(1x,i3)) 2236 deallocate (tab) 2237 2238 end subroutine bilan_shared_haplo
check_change_maker_windows
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
liste_haplo
DESCRIPTION
NOTES
SOURCE
575 function check_change_maker_windows(c,dx) result(change) 576 integer , intent(in) :: c 577 real(kind=dp) , intent(in) :: dx 578 579 logical :: change 580 integer :: lkmin,lkmax,lk 581 582 change = .true. 583 if (.not. allocated(loc_haplo) ) then 584 ! print *,'not allocated...' 585 return 586 end if 587 ! lkmax est la limite pour le marqueur flanquant l haplotype sur sa droite 588 ! 589 lkmin=longhap/2 590 if ( lkmin == 0 ) lkmin = 1 591 lkmax=nmk(c)-lkmin+1 592 ! la position testee doit etre compatible avec la longueur de l'haplotype 593 ! 594 595 if((dx < posi(c,lkmin)) .or. (dx > posi(c,lkmax))) then 596 call stop_application("liste_haplo : incompatibilty between haplotype length and tested position.") 597 end if 598 599 lk=lkmin +1 600 601 do while (dx > posi(c,lk) .and. lk < lkmax) 602 lk=lk+1 603 end do 604 605 if ( loc_haplo(1) == lk-lkmin ) change=.false. 606 607 end function check_change_maker_windows
free_internal_struct_ldla
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
free_internal_struct_ldla
DESCRIPTION
NOTES
SOURCE
285 subroutine free_internal_struct_ldla 286 if (allocated(haplo)) deallocate (haplo) 287 if (allocated(haplo_complet)) deallocate (haplo_complet) 288 if (allocated(race_h)) deallocate (race_h) 289 if (allocated(loc_haplo)) deallocate (loc_haplo) 290 if (allocated(race_haplo_complet)) deallocate (race_haplo_complet) 291 if (allocated(gamete)) deallocate (gamete) 292 if (allocated(count_haplo_complet)) deallocate (count_haplo_complet) 293 if (allocated(nb_gam)) deallocate (nb_gam) 294 if (allocated(prob_gam)) deallocate (prob_gam) 295 if (allocated(num_haplo_pere)) deallocate (num_haplo_pere) 296 if (allocated(num_haplo_mere)) deallocate (num_haplo_mere) 297 if (allocated(num_haplo_desc)) deallocate (num_haplo_desc) 298 if (allocated(nb_gam_pere)) deallocate (nb_gam_pere) 299 if (allocated(nb_gam_mere)) deallocate (nb_gam_mere) 300 if (allocated(nb_gam_desc)) deallocate (nb_gam_desc) 301 ! if (allocated(prob_gam_pere)) deallocate (prob_gam_pere) 302 ! if (allocated(prob_gam_mere)) deallocate (prob_gam_mere) 303 ! if (allocated(prob_gam_desc)) deallocate (prob_gam_desc) 304 if (allocated(tab_IBS)) deallocate (tab_IBS) 305 if (allocated(name_haplo)) deallocate (name_haplo) 306 if (allocated(name_haplo_reduit)) deallocate (name_haplo_reduit) 307 if (allocated(race_haplo_reduit)) deallocate (race_haplo_reduit) 308 if (allocated(name_haplo_complet)) deallocate (name_haplo_complet) 309 if (allocated(count_haplo)) deallocate (count_haplo) 310 if (allocated(tab_shared_haplo)) deallocate( tab_shared_haplo) 311 if (allocated(comp_haplo)) deallocate (comp_haplo) 312 if (allocated(comp_haplo_reduit)) deallocate (comp_haplo_reduit) 313 if (allocated(liste_haplo_reduit)) deallocate (liste_haplo_reduit) 314 if (allocated(pb_haplo_reduit)) deallocate (pb_haplo_reduit) 315 if (allocated(pb_haplo_complet)) deallocate (pb_haplo_complet) 316 if (allocated(tab_haplo_complet)) deallocate (tab_haplo_complet) 317 if (allocated(nb_haplo_possible)) deallocate (nb_haplo_possible) 318 319 end subroutine free_internal_struct_ldla
global_gamete
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
global_gamete
DESCRIPTION
NOTES
global_gamete trouve les gametes parentaux transmis par ses parents (ip et jm) pour une phase geno de la mere a un individu particulier (kd) sur l'ensemble du chromosome (c) les doubles recombinaisons entre marqueurs informatifs flanquants sont supposees absentes en sortie, global_gamete fournit tab_IBS(kd,lk,i) pour lk=1,nmk(c) et i=1 (gamete paternel) et 2 gamete maternel) Utilise nmk,nmanque,correp,corred,correm,pheno de m_qtlmap_data Utilise genotyp de m_qtlmap_haplotype_data HYSTORY correction BUG : leak memory 05/11/2010
SOURCE
1996 subroutine global_gamete(ip,jm,geno,kd,c) 1997 integer , intent(in) ::ip,jm,geno,kd,c 1998 1999 ! divers 2000 integer lk,llk 2001 integer (kind=KIND_PHENO) orig_left,orig_right 2002 logical found 2003 integer (kind=KIND_PHENO) :: allele(2) 2004 integer :: orig(2) 2005 integer :: stat 2006 integer (kind=KIND_PHENO) :: unknown=0 2007 ! dimension allocatable 2008 integer(kind=KIND_PHENO) , dimension (:) , allocatable :: origine 2009 ! allocation 2010 ! 2011 allocate (origine(maxval(nmk)),STAT=stat) 2012 call check_allocate(stat,'origine [m_qtlmap_haplotype_LD]') 2013 2014 tab_IBS(kd,:,:)=0 2015 2016 do lk = 1,nmk(c) 2017 ! 2018 ! pour chaque marqueur de l'intervalle, on cherche (avec point_gamete) le nom des alleles recus 2019 ! de chacun des parents (allele), ainsi que leur origine grand parentale 2020 ! orig vaut 1 si origine = grand pere; 2 = grand mere ; 2021 ! 0 = indetermine (duo ou triplet heterozygote identique); 3= parent homozygote, 4= manquant 2022 ! 2023 call point_gamete(ip,jm,geno,kd,c,lk,allele,orig) 2024 origine(lk)=orig(1) 2025 if(origine(lk) ==1 .or. origine(lk) ==2 .or. origine(lk) == 3) tab_IBS(kd,lk,1)=allele(1) 2026 ! 2027 ! si l'allele paternel est indermine on regarde ce qui s'est passe a gauche et a droite 2028 if(origine(lk) == 0) then 2029 2030 orig_left=0 2031 found=.false. 2032 llk=lk 2033 do while((.not. found ).and. llk > 1) 2034 llk=llk-1 2035 if(origine(llk) ==1 .or. origine(llk) ==2 ) then 2036 found=.true. 2037 orig_left=origine(llk) 2038 end if 2039 end do ! found 2040 orig_right=0 2041 found=.false. 2042 llk=lk 2043 do while((.not. found ) .and. llk < nmk(c)) 2044 llk=llk+1 2045 call point_gamete(ip,jm,geno,kd,c,llk,allele,orig) 2046 if(orig(1) ==1 .or. orig(1) ==2 ) then 2047 found=.true. 2048 orig_right=orig(1) 2049 end if 2050 end do ! found 2051 2052 if(orig_left == orig_right) tab_IBS(kd,lk,1)=genotyp(c,lk,correp(ip),orig_left) 2053 !if(orig_left == 0) tab_IBS(kd,lk,1)=genotyp(c,lk,correp(ip),orig_right) !!ATTENTION orig_right peut valoir =0 2054 !if(orig_right == 0) tab_IBS(kd,lk,1)=genotyp(c,lk,correp(ip),orig_left) 2055 if( ((orig_left == 0).or.(orig_left == 4) ) .and. ( (orig_right /= 0) .and. (orig_right /= 4)) ) then 2056 if ( orig_right /= 3 ) then 2057 tab_IBS(kd,lk,1)=genotyp(c,lk,correp(ip),orig_right) 2058 else 2059 tab_IBS(kd,lk,1)=genotyp(c,lk,correp(ip),1) 2060 end if 2061 end if 2062 2063 if( ((orig_right == 0).or.(orig_right == 4) ) .and. ( (orig_left /= 0) .and. (orig_left /= 4)) ) then 2064 if ( orig_left /= 3 ) then 2065 tab_IBS(kd,lk,1)=genotyp(c,lk,correp(ip),orig_left) 2066 else 2067 tab_IBS(kd,lk,1)=genotyp(c,lk,correp(ip),1) 2068 end if 2069 end if
liste_gamete
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
liste_gamete
DESCRIPTION
NOTES
liste_gamete trouve les gametes du chromosome c ayant forme chacun des descendant Utilise nmk,np,posi,correp,nm,repfem,correm,estfem de m_qtlmap_data Utilise genotyp,ngenom,genotypm de m_qtlmap_haplotype_data
SOURCE
1946 subroutine liste_gamete(c) 1947 1948 ! 1949 integer , intent(in) :: c 1950 1951 integer ip,jm,kd,geno,i,lk!,indtab(10000) 1952 1953 geno=0 1954 ! 1955 do ip=1,np 1956 do jm=nmp(ip)+1,nmp(ip+1) 1957 do kd=ndm(jm)+1,ndm(jm+1) 1958 if (presentg(c,kd)) call global_gamete(ip,jm,geno,kd,c) 1959 1960 ! do i=1,2 1961 ! do lk=1,nmk(c) 1962 ! indtab(lk)=0 1963 ! if(tab_IBS(kd,lk,i) /= 0)indtab(lk)=tab_IBS(kd,lk,i)- VAL_MIN_INDEX_PHENO 1964 ! end do ! lk 1965 ! end do ! i 1966 end do !kd 1967 end do !jm 1968 end do !ip 1969 ! 1970 1971 end subroutine liste_gamete
liste_haplo
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
liste_haplo
DESCRIPTION
NOTES
liste_haplo cherche les haplotypes de longeur longhap possibles autour de la position dx du chromosome c Utilise nmk,np,posi,correp,nm,repfem,correm,estfem de m_qtlmap_data Utilise genotyp,ngenom,genotypm de m_qtlmap_haplotype_data
SOURCE
622 subroutine liste_haplo(c,dx,n,hsire,hdam) 623 integer , intent(in) :: c 624 real(kind=dp) , intent(in) :: dx 625 integer, intent (in) :: n 626 logical , intent(in) :: hsire,hdam 627 628 integer lkmax,lkmin,i,ip,j,jm,geno,j_gam,kd,lk 629 630 nb_gam=0 631 632 ! lkmax est la limite pour le marqueur flanquant l haplotype sur sa droite 633 ! 634 lkmin=longhap/2 635 if ( lkmin == 0 ) lkmin = 1 636 lkmax=nmk(c)-lkmin+1 637 ! la position testee doit etre compatible avec la longueur de l'haplotype 638 ! 639 640 if((dx < posi(c,lkmin)) .or. (dx > posi(c,lkmax))) then 641 call stop_application("liste_haplo : incompatibilty between haplotype length and tested position.") 642 end if 643 644 ! 645 ! initialisation de loc_haplo, le vecteur de dimension longhap qui contient les numeros des marqueurs 646 ! pour l haplotype entourant la position testee dx, et lk le premier marqueur flanquant a droite 647 do i=1,longhap 648 loc_haplo(i)=i 649 end do !i 650 651 lk=loc_haplo(lkmin) +1 652 653 do while (dx > posi(c,lk) .and. lk < lkmax) 654 lk=lk+1 655 end do 656 657 ! 658 ! lk est le premier marqueur a droite de dx 659 do i=1,longhap 660 loc_haplo(i)=lk-lkmin-1+i 661 end do !i 662 663 !print *,loc_haplo 664 ! 665 ! le premier haplotype dans la liste sera celui des donnees manqantes, il y en aura un par race quand il y a plusieurs races 666 ! 667 nb_haplo=0 668 race_h='UNKNOWN' 669 do j=1,NB_RACES 670 nb_haplo=nb_haplo+1 671 do i=1,longhap 672 haplo(nb_haplo,i)=nmanque 673 end do !i 674 race_h(nb_haplo)=nom_race(j) 675 enddo 676 ! 677 ! haplotypes trouves chez les peres 678 ! 679 if(hsire)then 680 do ip=1,np 681 do j=1,2 682 nb_haplo=nb_haplo+1 683 do i=1,longhap 684 haplo(nb_haplo,i)=genotyp(c,loc_haplo(i),correp(ip),j) 685 end do !i 686 !! Si les parents du père sont connus (reppere.ne.INT_NOT_DEFINED) et sont 687 !! soit de meme race soit génotypés sur au moins un marqueur 688 if (raceFileDefined) then 689 if (reppere(ip).eq.INT_NOT_DEFINED) then 690 call stop_application('Both parents of sires and dams have to be defined in the ped file '//& 691 'to use breed origin of haplotypes. The parents of sire '//trim(pere(ip))//& 692 ' are not in the ped file! ') 693 else 694 if (( racep(reppere(ip))==racem(reppere(ip)) ) & 695 .or. (phasp(c,ip) )) then 696 if (j==1) race_h(nb_haplo)= racep(reppere(ip)) 697 if (j==2) race_h(nb_haplo)= racem(reppere(ip)) 698 endif 699 endif 700 endif 701 end do !j 702 703 end do !ip 704 end if ! hsire 705 706 if(hdam)then 707 708 ! 709 ! haplotypes trouves chez les meres phasees 710 ! 711 712 do ip=1,np 713 do jm=nmp(ip)+1,nmp(ip+1) 714 if(estfem(repfem(jm)))then 715 ! 716 ! cas des femelles dont on estime la phase 717 ! 718 do geno=ngenom(c,jm)+1,ngenom(c,jm+1) 719 do j=1,2 720 nb_haplo=nb_haplo+1 721 do i=1,longhap 722 haplo(nb_haplo,i)=genotypm(c,loc_haplo(i),geno,j) 723 end do !i 724 !! Si les parents de la mère sont connus (repmere.ne.INT_NOT_DEFINED) et sont 725 !! soit de meme race soitgénotypés sur au moins un marqueur 726 if (raceFileDefined) then 727 if (repmere(jm).eq.INT_NOT_DEFINED) then 728 call stop_application ("Both parents of sires and dams have to be defined in the ped file "//& 729 "to use breed origin of haplotypes. the parents of dam "//trim(mere(jm))//& 730 " are not in the ped file!") 731 else 732 if ( ( racep(repmere(jm))==racem(repmere(jm)) ) & 733 .or. (phasm(c,jm) )) then 734 if (j==1) race_h(nb_haplo)= racep(repmere(jm)) 735 if (j==2) race_h(nb_haplo)= racem(repmere(jm)) 736 end if 737 end if 738 end if 739 end do !j 740 end do !geno 741 ! 742 ! cas des femelles non phasees (issues d'une seule race lorsqu'elle est renseignée) 743 ! 744 else 745 geno=0 746 do kd=ndm(jm)+1,ndm(jm+1) 747 if(presentg(c,kd)) then 748 call local_gamete(ip,jm,geno,kd,c,dx) 749 do j_gam=1,nb_gam(kd) 750 nb_haplo=nb_haplo+1 751 do i=1,longhap 752 haplo(nb_haplo,i)=gamete(2,i,j_gam) 753 end do !i 754 if (raceFileDefined) then 755 if (repmere(jm).eq.INT_NOT_DEFINED) then 756 call stop_application ("Both parents of sires and dams have to be defined in the ped file "//& 757 "to use breed origin of haplotypes. the parents of dam "//trim(mere(jm))//& 758 " are not in the ped file!") 759 STOP 760 else 761 if (racep(repmere(jm))==racem(repmere(jm))) race_h(nb_haplo)= racep(repmere(jm)) 762 endif 763 endif 764 end do !j_gam 765 end if !presentg 766 end do !kd 767 end if 768 769 end do !jm 770 end do !ip 771 end if ! hdam 772 !ATTENTION: lorsqu'un haplotype n'a pas une race connue -->le programme ne continue pas 773 ! normalement cette erreur a été detectée avant 774 775 do j=1,nb_haplo 776 if (raceFileDefined.and.race_h(j)=='UNKNOWN') then 777 call stop_application("An haplotype has an unknown breed origin, PLEASE check if you give"//& 778 " the breeds of all funder parents!") 779 endif 780 !print *, 'HAPLO_pop:',j, race_h(j), (haplo(j,i),i=1,longhap) 781 enddo 782 end subroutine liste_haplo
liste_haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
liste_haplo_complet
DESCRIPTION
NOTES
liste_haplo_complet etablit la liste des haplotypes complets possible en entree, il lui faut loc_haplo , longhap et loc_haplo créés par la routine liste_haplo du module m_qtlmap_haplotype_ldla nall et alleles créés par la routine set_allele_info_vector du module m_qtlmap_genotype en sortie il fournit tab_haplo_complet et nb_haplo_complet
SOURCE
799 subroutine liste_haplo_complet(c) 800 integer , intent(in) :: c 801 integer i1,i2,j,j1,j2 802 integer nb_row,nb_col 803 integer :: stat 804 logical :: allocated 805 806 ! allocation 807 808 nb_haplo_complet = 1 809 do j=1,longhap 810 nb_haplo_complet = nb_haplo_complet * nall(c,loc_haplo(j)) 811 end do !j 812 813 allocated = allocated(tab_haplo_complet) 814 815 if ( .not. allocated .or. (size(tab_haplo_complet,1)<nb_haplo_complet) ) then 816 if (allocated) deallocate (tab_haplo_complet) 817 allocate (tab_haplo_complet(nb_haplo_complet*NB_RACES ,longhap),STAT=stat) 818 call check_allocate(stat,'tab_haplo_complet [m_qtlmap_haplotype_ldla:liste_haplo_complet]') 819 end if 820 821 822 823 tab_haplo_complet = nmanque 824 825 826 do i1 = 1,nall(c,loc_haplo(1)) 827 tab_haplo_complet(i1,1) = set_pheno(c,alleles(c,loc_haplo(1),i1)) 828 end do !i1 829 830 831 nb_row=nall(c,loc_haplo(1)) 832 nb_col=1 833 834 do j=2,longhap 835 836 do i1 = 1,nb_row 837 tab_haplo_complet(i1,nb_col+1) = set_pheno(c,alleles(c,loc_haplo(j),1)) 838 end do !i1 839 840 do j1= 2,nall(c,loc_haplo(j)) 841 842 do i2 = 1,nb_row 843 do j2= 1,nb_col 844 tab_haplo_complet(i2+(j1-1)*nb_row,j2) = tab_haplo_complet(i2,j2) 845 end do !j2 846 tab_haplo_complet(i2+(j1-1)*nb_row,nb_col+1) = set_pheno(c,alleles(c,loc_haplo(j),j1)) 847 end do !i2 848 849 end do !j1 850 nb_row=nb_row * nall(c,loc_haplo(j)) 851 nb_col= nb_col+1 852 end do !j 853 854 855 856 end subroutine liste_haplo_complet
liste_shared_haplo
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
liste_shared_haplo
DESCRIPTION
NOTES
liste_shared_haplo determine si un descendant a recu de sa mere , autour de la positon dx, un allele ibs a un de ceux de son pere, et sur quelle longueur l'effet haplotypique sera ajoute au modle de description des performances si cette longueur est superieure a LONG_MIN_IBS Utilise nmk,np,posi,correp,nm de m_qtlmap_data Utilise genotyp de m_qtlmap_haplotype_data
SOURCE
2113 subroutine liste_shared_haplo(c,dx,n) 2114 2115 ! 2116 integer , intent(in) :: c 2117 real(kind=dp) , intent(in) :: dx 2118 integer, intent (in) :: n 2119 2120 integer lkleft,ip,jm,kd,jp,iall,icount,icount_max,lk,stat 2121 2122 if (allocated(shared_haplo)) deallocate(shared_haplo) 2123 allocate (shared_haplo(nd),STAT=stat) 2124 call check_allocate(stat,'shared_haplo [m_qtlmap_haplotype_LD]') 2125 2126 2127 ! la position testee doit etre entre les bornes du chromsome 2128 ! 2129 if((dx < posi(c,1)) .or. (dx > posi(c,nmk(c)))) then 2130 call stop_application("liste_shared_area : incompatibilty between haplotype length and tested position.") 2131 end if 2132 2133 tab_shared_haplo=0 2134 ! 2135 ! recherche du marqueur a gauche de dx 2136 do lkleft=1,nmk(c)-1 2137 if(dx <= posi(c,lkleft)) exit 2138 end do 2139 ! 2140 do ip=1,np 2141 do jm=nmp(ip)+1,nmp(ip+1) 2142 do kd=ndm(jm)+1,ndm(jm+1) 2143 icount_max=0 2144 shared_haplo(kd)=0 2145 ! 2146 ! icount est la longueur du segment de l'haplotype maternel ibs au chromosome iall du pere jp 2147 do jp=1,np 2148 do iall=1,2 2149 icount=0 2150 do lk=lkleft,1,-1 2151 if(tab_IBS(kd,lk,2) /= genotyp(c,lk,correp(jp),iall)) exit 2152 if(genotyp(c,lk,correp(jp),iall)/= nmanque) icount=icount+1 2153 end do 2154 do lk=lkleft+1,nmk(c) 2155 if(tab_IBS(kd,lk,2) /= genotyp(c,lk,correp(jp),iall)) exit 2156 if(genotyp(c,lk,correp(jp),iall)/= nmanque) icount=icount+1 2157 end do 2158 2159 if(icount > icount_max .and. icount >= LONG_MIN_IBS) then 2160 icount_max=icount 2161 shared_haplo(kd)=iall+2*(jp-1) 2162 end if 2163 end do ! iall 2164 end do ! jp 2165 2166 if( shared_haplo(kd) /= 0) then 2167 iall=2-modulo(shared_haplo(kd),2) 2168 jp=(shared_haplo(kd)+modulo(shared_haplo(kd),2))/2 2169 tab_shared_haplo(jp,iall,1+icount_max)=tab_shared_haplo(jp,iall,1+icount_max)+1 !nombre d'haplotyped e longueur count_max 2170 tab_shared_haplo(jp,iall,1)=tab_shared_haplo(jp,iall,1)+1 ! nombre d haplotype partage 2171 end if 2172 end do !kd 2173 end do !jm 2174 end do !ip 2175 ! 2176 end subroutine liste_shared_haplo
local_gamete
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
local_gamete
DESCRIPTION
NOTES
local_gamete trouve les gametes parentaux transmis par ses parents (ip et jm) pour une phase geno de la mere a un individu particulier (kd) autour d'un point precis defini par le chromosome (c) et la position (dx) ainsi que le nombre de marqueurs l'entourant (longhap) en sortie, local_gamete fournit nb_gam = le nombre de couples de gametes parentaux (2 ** nb indeterminations) gamete(1,lk,j_gam) l allele paternel au lk eme marqueur pour le j_gam eme gamete gamete(2,lk,j_gam) l allele maternel au lk eme marqueur pour le j_gam eme gamete prob_gam(j_gam) la probabilite de cet allele Utilise posi,nmk,nmanque,correp,corred,correm,pheno de m_qtlmap_data Utilise genotyp,genotypm de m_qtlmap_haplotype_data
SOURCE
1338 subroutine local_gamete(ip,jm,geno,kd,c,dx) 1339 integer , intent(in) ::ip,jm,geno,kd,c 1340 real (kind=dp), intent(in) :: dx 1341 1342 ! divers 1343 integer lk,lkleft,lkright,l,llk 1344 integer left_known_sire,left_known_dam,right_known_sire,right_known_dam 1345 integer lks,lkd,rks,rkd 1346 integer isex,igam 1347 integer i,j,k,idup 1348 integer ngamkept,nb_inc 1349 integer value_alloc 1350 logical known_sire_left,known_dam_left,known_sire_right,known_dam_right,manquant 1351 logical ksr,ksl,kdr,kdl 1352 integer (kind=KIND_PHENO) :: allele(2),unknown 1353 integer :: orig(2) 1354 real(kind=dp) recm,recf,probtot 1355 integer :: stat 1356 1357 ! dimension allocatable 1358 integer , dimension (:,:) , allocatable :: base 1359 integer(kind=KIND_PHENO) , dimension (:,:,:) , allocatable :: origine 1360 integer , dimension (:,:) , allocatable :: tabdup 1361 ! allocation 1362 ! 1363 allocate (base(2,nmk(c)),STAT=stat) 1364 call check_allocate(stat,'base [m_qtlmap_haplotype_LD]') 1365 1366 value_alloc= 2**(2*longhap) 1367 1368 allocate (origine(2,nmk(c),value_alloc),STAT=stat) 1369 call check_allocate(stat,'origine [m_qtlmap_haplotype_LD]') 1370 if ( allocated(tabdup)) deallocate (tabdup) 1371 allocate (tabdup(value_alloc,2*longhap),STAT=stat) 1372 call check_allocate(stat,'tabdup[m_qtlmap_haplotype_LD]') 1373 unknown=0 1374 1375 gamete=nmanque 1376 1377 ! la position testee doit etre compatible avec la longueur de l'haplotype 1378 ! 1379 1380 if ( longhap/2 > 0 ) then 1381 if((dx < posi(c,longhap/2)) .or. (dx > posi(c,nmk(c)-longhap/2+1)))then 1382 call stop_application("local_gamete : incompatibilty between haplotype length and tested position.") 1383 end if 1384 end if 1385 ! 1386 ! on cherche le marqueur flanquant a gauche la position dx 1387 ! 1388 lk=1 1389 do while (dx > posi(c,lk)) 1390 lk=lk+1 1391 end do 1392 lk=lk-1 1393 ! 1394 ! lkleft et lkright seront les numeros des marqueurs bornant l haplotype a determiner 1395 ! 1396 lkleft =lk-longhap/2+1 1397 if(lk <= longhap/2) lkleft=1 1398 lkright =lk+longhap/2 1399 if(lk >= nmk(c)-longhap/2) lkright=nmk(c) 1400 ! 1401 ! du fait des incertitudes 1402 ! (pere et descendant heterozygotes identiques si la mere n est pas marquee 1403 ! deux parents et descendants heterozygotes identiques si la mere est marquee) 1404 ! plusieurs gametes peuvent avoir forme un descendant 1405 ! 1406 ! nb_gam = 2**(nb d'incertitudes) est le nombre de gametes possibles 1407 ! gamete(1,lk,j_gam) est le nom de l'allele au i eme marqueur transmis par le pere 1408 ! (par la mere pour gamete(2,lk,j_gam)) 1409 ! pour la j_gam eme configuration gametique 1410 ! origine(1,lk,j_gam) en est l origine grand paternel 1411 ! (0 pour inconnue, 1 si grand pere, 2 si grand mere, 3 le pere est homozygote, 4 si manquant) 1412 ! (grand maternelle pour origine(2,lk)) 1413 ! 1414 ! construction des gametes correspondant a l haplotype 1415 ! pour chaque marqueur de l'intervalle, on cherche (avec point_gamete) le nom des alleles recus 1416 ! de chacun des parents, ainsi que leur origine grand parentale (mis dans base) 1417 ! base vaut 1 si origine = grand pere; 2 = grand mere ; 1418 ! 0 = indetermine (duo ou triplet heterozygote identique); 3= parent homozygote, 4= manquant 1419 ! 1420 origine =0 1421 base=0 1422 1423 manquant=.false. 1424 do lk = lkleft,lkright 1425 call point_gamete(ip,jm,geno,kd,c,lk,allele,base(:,lk)) 1426 if(allele(1) == nmanque .or. allele(2) == nmanque) manquant=.true. 1427 origine(1,lk,1)=base(1,lk) 1428 origine(2,lk,1)=base(2,lk) 1429 end do !lk 1430 1431 ! 1432 ! je mets a manque total les haplotypes compremnant une donnee manquante 1433 ! 1434 if(manquant) then 1435 nb_gam(kd)=1 1436 prob_gam(kd,1)=1.d0 1437 1438 llk=0 1439 do lk = lkleft,lkright 1440 llk=llk+1 1441 do isex=1,2 1442 gamete(isex,llk,1)=nmanque 1443 end do ! isex 1444 end do !lk 1445 1446 return 1447 1448 end if 1449 ! 1450 ! il y a 2**(nombre d'origines inconnues) gametes possibles dont on doit calculer les probabilites 1451 ! 1452 ! quand les marqueurs flanquants sont d'origine inconnue, il faut rechercher une information 1453 ! a l exterieur de l'haplotype 1454 ! 1455 ! vers la gauche 1456 ! 1457 left_known_sire=0 1458 left_known_dam=0 1459 known_sire_left=.false. 1460 if(base(1,lkleft) ==1 .or.base(1,lkleft) == 2) then 1461 known_sire_left=.true. 1462 left_known_sire=lkleft 1463 end if 1464 known_dam_left=.false. 1465 if(base(2,lkleft) ==1 .or.base(2,lkleft) == 2) then 1466 known_dam_left=.true. 1467 left_known_dam=lkleft 1468 end if 1469 1470 if(.not.manquant)then 1471 if(estfem(repfem(jm))) then 1472 1473 if(.not. (known_sire_left .and. known_dam_left)) then 1474 do lk=lkleft-1,1,-1 1475 call point_gamete(ip,jm,geno,kd,c,lk,allele,base(:,lk)) 1476 1477 if(.not. known_sire_left .and. (base(1,lk) == 1 .or. base(1,lk) == 2)) then 1478 known_sire_left=.true. 1479 left_known_sire=lk 1480 origine(1,lk,1)=base(1,lk) 1481 end if 1482 if(.not. known_dam_left .and. (base(2,lk) == 1 .or. base(2,lk) == 2)) then 1483 known_dam_left=.true. 1484 left_known_dam=lk 1485 origine(2,lk,1)=base(2,lk) 1486 end if 1487 1488 if(known_sire_left .and. known_dam_left) exit 1489 end do !lk 1490 end if 1491 1492 else 1493 1494 if(.not.known_sire_left) then 1495 do lk=lkleft-1,1,-1 1496 call point_gamete(ip,jm,geno,kd,c,lk,allele,base(:,lk)) 1497 1498 if(.not. known_sire_left .and. (base(1,lk) == 1 .or. base(1,lk) == 2)) then 1499 known_sire_left=.true. 1500 left_known_sire=lk 1501 origine(1,lk,1)=base(1,lk) 1502 end if 1503 if(known_sire_left) exit 1504 end do !lk 1505 end if 1506 1507 end if 1508 end if 1509 ! 1510 ! vers la droite 1511 ! 1512 right_known_sire=0 1513 right_known_dam=0 1514 known_sire_right=.false. 1515 if(base(1,lkright) ==1 .or.base(1,lkright) == 2) then 1516 known_sire_right=.true. 1517 right_known_sire=lkright 1518 end if 1519 known_dam_right=.false. 1520 if(base(2,lkright) ==1 .or.base(2,lkright) == 2) then 1521 known_dam_right=.true. 1522 right_known_dam=lkright 1523 end if 1524 1525 if(.not.manquant) then 1526 if(estfem(repfem(jm))) then 1527 1528 if(.not. (known_sire_right .and. known_dam_right)) then 1529 1530 do lk=lkright+1,nmk(c) 1531 1532 call point_gamete(ip,jm,geno,kd,c,lk,allele,base(:,lk)) 1533 1534 1535 if(.not. known_sire_right .and. (base(1,lk) == 1 .or. base(1,lk) == 2)) then 1536 known_sire_right=.true. 1537 right_known_sire=lk 1538 origine(1,lk,1)=base(1,lk) 1539 end if 1540 if(.not. known_dam_right .and. (base(2,lk) == 1 .or. base(2,lk) == 2)) then 1541 known_dam_right=.true. 1542 right_known_dam=lk 1543 origine(2,lk,1)=base(2,lk) 1544 end if 1545 1546 if(known_sire_right .and. known_dam_right) exit 1547 end do !lk 1548 end if 1549 else 1550 1551 if(.not.known_sire_right) then 1552 do lk=lkright+1,nmk(c) 1553 call point_gamete(ip,jm,geno,kd,c,lk,allele, base(:,lk)) 1554 1555 if(.not. known_sire_right .and. ( base(1,lk) == 1 .or. base(1,lk) == 2)) then 1556 known_sire_right=.true. 1557 right_known_sire=lk 1558 origine(1,lk,1)= base(1,lk) 1559 end if 1560 if(known_sire_right) exit 1561 end do !lk 1562 end if 1563 end if 1564 end if 1565 ! 1566 ! Construction de la liste des gametes possibles 1567 ! en parcourant les marqueurs le long de l haplotype et en dedoublant la liste 1568 ! a chaque fois qu'une incertitude (base =0) est rencontree 1569 ! nb_inc est le nombre d'incertitudes 1570 ! nb_gam(kd) le nombre de gametes possibles pour kd (nb_gam(kd)=2**nb_inc dans un premier temps) 1571 ! 1572 nb_inc=0 1573 1574 do lk = lkleft,lkright 1575 do isex=1,2 1576 if(base(isex,lk) == 0) nb_inc=nb_inc + 1 1577 end do !isex 1578 end do !lk 1579 nb_gam(kd)=2**nb_inc 1580 ! 1581 !tabdup contient la liste des vecteurs d'origines possibles pour les nb_gam(kd) incertitudes 1582 ! 1583 if(nb_inc > 0) then 1584 tabdup(1,nb_inc)=1 1585 tabdup(2,nb_inc)=2 1586 do i=2,nb_inc 1587 do j=1,2**(i-1) 1588 tabdup(j,nb_inc-i+1)=1 1589 tabdup(2**(i-1)+j,nb_inc-i+1)=2 1590 do k=1,i-1 1591 tabdup(2**(i-1)+j,nb_inc-k+1)=tabdup(j,nb_inc-k+1) 1592 end do !k 1593 end do !j 1594 end do !i 1595 1596 do igam=1,nb_gam(kd) 1597 idup=0 1598 do lk = lkleft,lkright 1599 do isex=1,2 1600 origine(isex,lk,igam)=origine(isex,lk,1) 1601 if(base(isex,lk) == 0) then 1602 idup=idup+1 1603 origine(isex,lk,igam)=tabdup(igam,idup) 1604 end if 1605 end do ! isex 1606 end do ! lk 1607 if(known_sire_left) origine(1,left_known_sire,igam) =origine(1,left_known_sire,1) 1608 if(known_sire_right) origine(1,right_known_sire,igam)=origine(1,right_known_sire,1) 1609 if(known_dam_left) origine(2,left_known_dam,igam) = origine(2,left_known_dam,1) 1610 if(known_dam_right) origine(2,right_known_dam,igam)= origine(2,right_known_dam,1) 1611 end do ! igam 1612 1613 end if !nb_gam 1614 1615 1616 ksr=known_sire_right 1617 ksl=known_sire_left 1618 kdr=known_dam_right 1619 kdl=known_dam_left 1620 lks=left_known_sire 1621 rks=right_known_sire 1622 lkd=left_known_dam 1623 rkd=right_known_dam 1624 1625 ! 1626 ! calcul des probabilites des gametes 1627 ! 1628 do igam=1,nb_gam(kd) 1629 prob_gam(kd,igam)=1.d0 1630 known_sire_right=ksr 1631 known_sire_left=ksl 1632 known_dam_right=kdr 1633 known_dam_left=kdl 1634 left_known_sire=lks 1635 right_known_sire=rks 1636 left_known_dam=lkd 1637 right_known_dam=rkd 1638 1639 do lk=lkleft,lkright 1640 1641 if(origine(1,lk,igam) < 3) then 1642 if(known_sire_left) then 1643 recm=xaldane(posim(c,lk)-posim(c,left_known_sire)) 1644 if(origine(1,lk,igam) == origine(1,left_known_sire,igam)) then 1645 prob_gam(kd,igam)=prob_gam(kd,igam) * (1.d0 -recm) 1646 else 1647 prob_gam(kd,igam)=prob_gam(kd,igam) * recm 1648 end if 1649 else 1650 known_sire_left=.true. 1651 end if 1652 left_known_sire=lk 1653 end if 1654 1655 if(geno /=0 .and. origine(2,lk,igam) < 3) then 1656 if(known_dam_left) then 1657 recf=xaldane(posif(c,lk)-posif(c,left_known_dam)) 1658 if(origine(2,lk,igam) == origine(2,left_known_dam,igam)) then 1659 prob_gam(kd,igam)=prob_gam(kd,igam) * (1.d0 -recf) 1660 else 1661 prob_gam(kd,igam)=prob_gam(kd,igam) * recf 1662 end if 1663 else 1664 known_dam_left=.true. 1665 end if 1666 left_known_dam=lk 1667 end if 1668 1669 end do ! lk 1670 ! 1671 ! bout du segment 1672 if(known_sire_right .and. right_known_sire > lkright) then 1673 recm=xaldane(posim(c,right_known_sire)-posim(c,lkright)) 1674 if(origine(1,lkright,igam) == origine(1,right_known_sire,igam)) then 1675 prob_gam(kd,igam)=prob_gam(kd,igam) * (1.d0 -recm) 1676 else 1677 prob_gam(kd,igam)=prob_gam(kd,igam) * recm 1678 end if 1679 end if 1680 if(geno /= 0 .and. known_dam_right .and. right_known_dam > lkright) then 1681 recf=xaldane(posif(c,right_known_dam)-posif(c,lkright)) 1682 if(origine(2,lkright,igam) == origine(2,right_known_dam,igam)) then 1683 prob_gam(kd,igam)=prob_gam(kd,igam) * (1.d0 -recf) 1684 else 1685 prob_gam(kd,igam)=prob_gam(kd,igam) * recf 1686 end if 1687 end if 1688 1689 end do !igam 1690 1691 ! au cas ou, standardisation des probabilites 1692 ! 1693 probtot=0.d0 1694 do igam=1,nb_gam(kd) 1695 probtot=probtot+prob_gam(kd,igam) 1696 end do 1697 do igam=1,nb_gam(kd) 1698 prob_gam(kd,igam)=prob_gam(kd,igam)/probtot 1699 end do 1700 ! 1701 ! on enrichir les informations sur les gametes pour chacune des configurations possibles 1702 ! 1703 gamete=nmanque 1704 do igam=1,nb_gam(kd) 1705 llk=0 1706 do lk=lkleft,lkright 1707 llk=llk+1 1708 if(origine(1,lk,igam) == 4) gamete(1,llk,igam) = nmanque 1709 if(origine(1,lk,igam) == 3) gamete(1,llk,igam) = genotyp(c,lk,correp(ip),1) 1710 if(origine(1,lk,igam) == 2) gamete(1,llk,igam) = genotyp(c,lk,correp(ip),2) 1711 if(origine(1,lk,igam) == 1) gamete(1,llk,igam) = genotyp(c,lk,correp(ip),1) 1712 1713 1714 if(origine(1,lk,igam) <= 3) then 1715 1716 if(gamete(1,llk,igam) == pheno(c,lk,corred(kd),1)) then 1717 gamete(2,llk,igam) = pheno(c,lk,corred(kd),2) 1718 else 1719 gamete(2,llk,igam) = pheno(c,lk,corred(kd),1) 1720 end if 1721 end if 1722 end do ! lk 1723 end do !igam 1724 ! 1725 ! 1726 ! elimination des cas improbables 1727 ! 1728 ngamkept=0 1729 do igam=1,nb_gam(kd) 1730 if(prob_gam(kd,igam) >= prob_seuil_gam)then 1731 ngamkept=ngamkept+1 1732 prob_gam(kd,ngamkept)=prob_gam(kd,igam) 1733 do llk=1,longhap 1734 gamete(1,llk,ngamkept)=gamete(1,llk,igam) 1735 gamete(2,llk,ngamkept)=gamete(2,llk,igam) 1736 end do !i 1737 end if 1738 end do !igam 1739 ! 1740 ! standardisation 1741 ! 1742 probtot=0.d0 1743 do igam=1,ngamkept 1744 probtot=probtot+prob_gam(kd,igam) 1745 end do 1746 do igam=1,ngamkept 1747 prob_gam(kd,igam)=prob_gam(kd,igam)/probtot 1748 end do 1749 1750 nb_gam(kd)=ngamkept 1751 1752 ! correctiopn BUG : leak memory 05/11/2010 1753 deallocate (base) 1754 deallocate (origine) 1755 end subroutine local_gamete
point_gamete
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
point_gamete
DESCRIPTION
NOTES
point_gamete determine le nom des alleles (allele) au marqueur lk du chromosome c recus par les descendant kd du pere ip et de la mere jm affectee de la phase geno ainsi que leur origine grand parentale (orig) allele(1) est l'allele recu du pere, allele(2) de la mere allele(1)=allele(2)=0 si le duo ou le triplet est heterozygte identique allele(1)=allele(2)=nmanque si le phenotype du descednant et ou du pere sont manquants orig(1) vaut 1 si origine = grand pere paternel ; 2 = grand mere paternelle; 0 = indetermine (duo ou triplet heterozygote identique); 3= pere homozygote orig(2) est le symetrique pour l'origne maternelle Utilise nmanque,correp,pheno,corred,correm de m_qtlmap_data Utilise genotyp,genotypm de m_qtlmap_genotype_data
SOURCE
1778 subroutine point_gamete(ip,jm,geno,kd,c,lk,allele,orig) 1779 1780 1781 integer(kind=KIND_PHENO) :: allele(2),unknown 1782 1783 integer :: orig(2) 1784 integer , intent(in) :: ip,jm,geno,kd,c,lk 1785 logical found 1786 1787 allele(1)=nmanque 1788 allele(2)=nmanque 1789 unknown=0 1790 found=.false. 1791 ! le genotype du pere et le phenotype du descendant ne doivent pas etre manquants 1792 ! 1793 if(genotyp(c,lk,correp(ip),1) /= nmanque .and. pheno(c,lk,corred(kd),1) /= nmanque) then !1 1794 ! 1795 ! cas du pere homozygote 1796 ! 1797 if(genotyp(c,lk,correp(ip),1) == genotyp(c,lk,correp(ip),2)) then !2 1798 1799 allele(1)=genotyp(c,lk,correp(ip),1) 1800 if(pheno(c,lk,corred(kd),1) == allele(1)) then 1801 allele(2)=pheno(c,lk,corred(kd),2) 1802 else 1803 allele(2)=pheno(c,lk,corred(kd),1) 1804 end if 1805 found=.true. 1806 ! 1807 ! pere heterozygote 1808 1809 else 1810 ! descendant homozygote 1811 ! 1812 if(pheno(c,lk,corred(kd),1) == pheno(c,lk,corred(kd),2)) then !3 1813 allele(1)=pheno(c,lk,corred(kd),1) 1814 allele(2)=pheno(c,lk,corred(kd),2) 1815 found=.true. 1816 end if 1817 1818 ! descendant heterozygote 1819 ! son premier allele ne pouvant pas venir du pere 1820 ! 1821 if(pheno(c,lk,corred(kd),1) /= genotyp(c,lk,correp(ip),1) .and. & 1822 pheno(c,lk,corred(kd),1) /= genotyp(c,lk,correp(ip),2)) then !4 1823 allele(1)=pheno(c,lk,corred(kd),2) 1824 allele(2)=pheno(c,lk,corred(kd),1) 1825 found=.true. 1826 end if 1827 ! 1828 ! descendant heterozygote 1829 ! son deuxieme allele ne pouvant pas venir du pere 1830 ! 1831 if(pheno(c,lk,corred(kd),2) /= genotyp(c,lk,correp(ip),1) .and. & 1832 pheno(c,lk,corred(kd),2) /= genotyp(c,lk,correp(ip),2)) then !5 1833 allele(1)=pheno(c,lk,corred(kd),1) 1834 allele(2)=pheno(c,lk,corred(kd),2) 1835 found=.true. 1836 end if 1837 ! 1838 ! incertitude pere - descendant, la mere etant connue 1839 ! 1840 if (.not. found ) then !6 1841 allele(1)=unknown 1842 allele(2)=unknown 1843 if ( pheno(c,lk,correm(jm),1) /= nmanque ) then !6 1844 ! 1845 ! cas de la mere homozygote 1846 ! 1847 if(pheno(c,lk,correm(jm),1) == pheno(c,lk,correm(jm),2)) then !7 1848 1849 allele(2)=pheno(c,lk,correm(jm),1) 1850 if(pheno(c,lk,corred(kd),1) == allele(2)) then 1851 allele(1)=pheno(c,lk,corred(kd),2) 1852 else 1853 allele(1)=pheno(c,lk,corred(kd),1) 1854 end if 1855 found=.true. 1856 end if 1857 ! 1858 ! mere heterozygote 1859 ! descendant heterozygote 1860 ! son premier allele ne pouvant pas venir de la mere 1861 ! 1862 if(pheno(c,lk,corred(kd),1) /= pheno(c,lk,correm(jm),1) .and. & 1863 pheno(c,lk,corred(kd),1) /= pheno(c,lk,correm(jm),2)) then !8 1864 allele(1)=pheno(c,lk,corred(kd),1) 1865 allele(2)=pheno(c,lk,corred(kd),2) 1866 found=.true. 1867 end if 1868 ! 1869 ! descendant heterozygote 1870 ! son deuxieme allele ne pouvant pas venir de la mere 1871 ! 1872 if(pheno(c,lk,corred(kd),2) /= pheno(c,lk,correm(jm),1) .and. & 1873 pheno(c,lk,corred(kd),2) /= pheno(c,lk,correm(jm),2)) then !9 1874 allele(1)=pheno(c,lk,corred(kd),2) 1875 allele(2)=pheno(c,lk,corred(kd),1) 1876 found=.true. 1877 end if 1878 1879 end if 1880 end if 1881 end if 1882 end if 1883 1884 ! 1885 ! on peut alors determiner les origines grands parentales des gametes recus par kd 1886 ! 1887 ! 0 = indetermine (le descendant est heterozygote comme ses parents) 1888 ! 3 = indeterminable (les parents sont homozygotes ) 1889 ! 4 = donnees manquantes) 1890 ! 1891 orig(1)=0 1892 if(genotyp(c,lk,correp(ip),1) == nmanque .or. pheno(c,lk,corred(kd),1) == nmanque) orig(1) =4 1893 if(orig(1) == 0 .and. (genotyp(c,lk,correp(ip),1) == genotyp(c,lk,correp(ip),2))) orig(1)=3 1894 if(orig(1) == 0 .and. allele(1) /= unknown) then 1895 if(allele(1) == genotyp(c,lk,correp(ip),1)) then 1896 orig(1)=1 1897 else 1898 orig(1)=2 1899 end if 1900 end if 1901 1902 orig(2)=0 1903 if(geno /= 0) then 1904 ! 1905 ! cas ou la mere est phasee (orig correspond aux origines grand parentales) 1906 ! 1907 if(genotypm(c,lk,geno,1) == nmanque .or. pheno(c,lk,corred(kd),1) == nmanque) orig(2)= 4 1908 if(orig(2) == 0 .and. (genotypm(c,lk,geno,1) == genotypm(c,lk,geno,2))) orig(2)=3 1909 if(orig(2) == 0 .and. allele(2) /= unknown) then 1910 if(allele(2) == genotypm(c,lk,geno,1)) then 1911 orig(2)=1 1912 else 1913 orig(2)=2 1914 end if 1915 end if 1916 else 1917 ! 1918 ! cas ou la mere n'est pas phasee (orig correspond a l'ordre de lecture des phenotypes) 1919 ! 1920 if(pheno(c,lk,correm(jm),1) == nmanque .or. pheno(c,lk,corred(kd),1) == nmanque) orig(2)= 4 1921 if(orig(2) == 0 .and. (pheno(c,lk,correm(jm),1) == pheno(c,lk,correm(jm),2))) orig(2)=3 1922 if(orig(2) == 0 .and. allele(2) /= unknown) then 1923 if(allele(2) == pheno(c,lk,correm(jm),1)) then 1924 orig(2)=1 1925 else 1926 orig(2)=2 1927 end if 1928 end if 1929 end if 1930 end subroutine point_gamete
proba_haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
proba_haplo_complet
DESCRIPTION
NOTES
proba_haplo_complet etablit les probabilites des haplotypes complets en entree, il lui faut nb_haplo_complet et tab_haplo_complet créés par la routine liste_haplo_complet du module m_qtlmap_haplotype_ldla haplo et nb_haplo créés par la routine set_allele_info_vector du module m_qtlmap_genotype en sortie il fournit pb_haplo_complet
SOURCE
875 subroutine proba_haplo_complet(c,dx) 876 integer , intent(in) :: c 877 real(kind=dp) , intent(in) :: dx 878 double precision erreur 879 integer iht,jhc,i,iter,ri 880 logical compatible 881 integer :: stat 882 real(kind=dp) :: S 883 ! dim allocatable 884 885 real(kind=dp) , dimension (:) , allocatable :: phc 886 real(kind=dp) , dimension (:) , allocatable :: pht 887 888 ! allocation 889 if ( allocated(pb_haplo_complet)) deallocate (pb_haplo_complet) 890 allocate (pb_haplo_complet(nb_haplo_complet*NB_RACES),STAT=stat) 891 call check_allocate(stat,'pb_haplo_complet [m_qtlmap_haplotype_ldla:proba_haplo_complet]') 892 allocate (comp_haplo(nb_haplo,nb_haplo_complet*NB_RACES),STAT=stat) 893 call check_allocate(stat,'comp_haplo [m_qtlmap_haplotype_ldla]') 894 895 896 allocate (phc(nb_haplo_complet*NB_RACES),STAT=stat) 897 call check_allocate(stat,'phc [m_qtlmap_haplotype_ldla:proba_haplo_complet]') 898 allocate (pht(nb_haplo),STAT=stat) 899 call check_allocate(stat,'pht [m_qtlmap_haplotype_ldla:proba_haplo_complet]') 900 901 902 pb_haplo_complet = 0.d0 903 comp_haplo=.false. 904 race_haplo_complet='UNKNOWN' 905 do ri=1,NB_RACES 906 do jhc = 1,nb_haplo_complet 907 do i=1,longhap 908 tab_haplo_complet(jhc+nb_haplo_complet*(ri-1),i)=tab_haplo_complet(jhc,i) 909 race_haplo_complet(jhc+nb_haplo_complet*(ri-1))=nom_race(ri) 910 enddo 911 enddo 912 enddo 913 nb_haplo_complet= nb_haplo_complet*NB_RACES 914 915 ! 916 ! creation de la table de compatibilites entre les haplotypes observes et les haplotypes complets possibles 917 ! 918 do iht=1,nb_haplo 919 do jhc = 1,nb_haplo_complet 920 compatible =.true. 921 !haplotypes de races différentes 922 if (race_haplo_complet(jhc) /= race_h(iht)) then 923 compatible = .false. 924 !exit 925 else 926 !haplotypes différents de meme race 927 do i=1,longhap 928 if( (haplo(iht,i) /= tab_haplo_complet(jhc,i)) .and. (haplo(iht,i) /= nmanque) )then 929 compatible = .false. 930 !exit 931 end if 932 end do !i 933 endif 934 if(compatible) comp_haplo(iht,jhc)=.true. 935 end do !jhc 936 end do !iht 937 938 939 ! 940 ! Algorithme EM pour l'estimation des probabilites des haplotypes complets 941 ! 942 pb_haplo_complet =1.d0 / dble(nb_haplo_complet) 943 944 erreur=1.d0 945 iter=0 946 947 do while( erreur > EPS_EM) 948 949 pht=0.d0 950 do iht = 1, nb_haplo 951 do jhc=1,nb_haplo_complet 952 if(comp_haplo(iht,jhc)) pht(iht)=pht(iht)+pb_haplo_complet(jhc) 953 end do !jhc 954 end do !iht 955 956 957 erreur=0.d0 958 phc=0.d0 959 do jhc=1,nb_haplo_complet 960 do iht = 1, nb_haplo 961 if(comp_haplo(iht,jhc)) phc(jhc)=phc(jhc)+pb_haplo_complet(jhc)/pht(iht) 962 end do !iht 963 phc(jhc)=phc(jhc)/dble(nb_haplo) 964 erreur=erreur+dabs(phc(jhc)-pb_haplo_complet(jhc)) 965 end do !jhc 966 967 iter=iter+1 968 if ( iter > ITER_EM_MAX) then 969 call stop_application( ' the EM algorithm for the estimation of haplotypes frequency does not converge ') 970 end if 971 do jhc=1,nb_haplo_complet 972 pb_haplo_complet(jhc)=phc(jhc) 973 end do !jhc 974 end do !while 975 S=0.D0 976 ! print *,'RACE PB_HAPLO' 977 do jhc=1,nb_haplo_complet 978 S=S+pb_haplo_complet(jhc) 979 !print *,trim(race_haplo_complet(jhc)), pb_haplo_complet(jhc) 980 end do !jhc 981 ! print *, 'Somme PROB_haplo_complet', S, nb_haplo_complet 982 983 deallocate (phc) 984 deallocate (pht) 985 986 end subroutine proba_haplo_complet
set_allocation
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
set_allocation
DESCRIPTION
NOTES
SOURCE
218 subroutine set_allocation_ldla 219 integer :: stat,value_alloc_1,value_alloc_2,c 220 221 222 allocate (loc_haplo(longhap),STAT=stat) 223 call check_allocate(stat,'loc_haplo [m_qtlmap_haplotype_ldla]') 224 allocate (nb_gam(nd),STAT=stat) 225 call check_allocate(stat,'nb_gam [m_qtlmap_haplotype_ldla]') 226 allocate (nb_gam_pere(np,2),STAT=stat) 227 call check_allocate(stat,'nb_gam_pere [m_qtlmap_haplotype_ldla]') 228 allocate (nb_gam_mere(maxval(ngenom(:,nm+1)),2),STAT=stat) 229 call check_allocate(stat,'nb_gam_mere [m_qtlmap_haplotype_ldla]') 230 allocate (tab_shared_haplo(np,2,1+maxval(nmk)),STAT=stat) 231 call check_allocate(stat,'tab_shared_haplo [m_qtlmap_haplotype_ldla]') 232 233 value_alloc_1= (2*(np+nm)+nd*2**longhap)*NB_RACES 234 value_alloc_2=maxval(nall)**(2*longhap)+1 235 236 allocate (haplo(value_alloc_1,longhap),STAT=stat) 237 call check_allocate(stat,'haplo [m_qtlmap_haplotype_ldla]') 238 ! allocate (haplo_complet(value_alloc_1,value_alloc_2),STAT=stat) 239 ! call check_allocate(stat,'haplo_complet [m_qtlmap_haplotype_ldla]') 240 allocate (nb_haplo_possible(value_alloc_1),STAT=stat) 241 call check_allocate(stat,'nb_haplo_possible [m_qtlmap_haplotype_ldla]') 242 !add race info 243 allocate (race_h(value_alloc_1),STAT=stat) 244 call check_allocate(stat,'race_h [m_qtlmap_haplotype_ldla]') 245 allocate (race_haplo_complet(value_alloc_1),STAT=stat) 246 call check_allocate(stat,'race_haplo_complet [m_qtlmap_haplotype_ldla]') 247 248 249 allocate (count_haplo_complet(value_alloc_1),STAT=stat) 250 call check_allocate(stat,'count_haplo_complet [m_qtlmap_haplotype_ldla]') 251 252 allocate (prob_gam(nd,value_alloc_2),STAT=stat) 253 call check_allocate(stat,'prob_gam [m_qtlmap_haplotype_ldla]') 254 allocate (nb_gam_desc(nd),STAT=stat) 255 call check_allocate(stat,'nb_gam_desc [m_qtlmap_haplotype_ldla]') 256 allocate (gamete(2,longhap,value_alloc_2),STAT=stat) 257 call check_allocate(stat,'gamete [m_qtlmap_haplotype_ldla]') 258 259 allocate (tab_IBS(nd,maxval(nmk),2),STAT=stat) 260 call check_allocate(stat,'tab_IBS [m_qtlmap_haplotype_ldla]') 261 262 allocate (name_haplo(value_alloc_2),STAT=stat) 263 call check_allocate(stat,'name_haplo [m_qtlmap_haplotype_ldla]') 264 allocate (name_haplo_reduit(value_alloc_2),STAT=stat) 265 call check_allocate(stat,'name_haplo_reduit [m_qtlmap_haplotype_ldla]') 266 allocate (race_haplo_reduit(value_alloc_2),STAT=stat) 267 call check_allocate(stat,'race_haplo_reduit [m_qtlmap_haplotype_ldla]') 268 allocate (name_haplo_complet(value_alloc_2),STAT=stat) 269 call check_allocate(stat,'name_haplo_complet [m_qtlmap_haplotype_ldla]') 270 allocate (count_haplo(value_alloc_2),STAT=stat) 271 call check_allocate(stat,'count_haplo [m_qtlmap_haplotype_ldla]') 272 273 274 end subroutine set_allocation_ldla
set_haplo_final
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
set_haplo_final
DESCRIPTION
NOTES
set_haplo_final met au propre les listes de correspondances des haplotypes paternels et maternels num_haplo_pere(ip,j,l) pour le pere ip, chromosome j et haplotype complet de la liste reduite possible l (avec = 1,... nb_gam_pere(ip,j)) num_haplo_mere(geno,j,l) pour la combinaison mere x phase geno, chromosome j et haplotype complet de la liste reduite possible l (avec = 1,... nb_gam_mere(geno,j)) num_haplo_desc(kd,j,l) pour le descendant kd, chromosome j et haplotype complet de la liste reduite possible l (avec = 1,... nb_gam_desc(kd)) Utilise np,nm,ndm,estfem,repfem de m_qtlmap_data utilise nb_haplo_reduit et comp_haplo_reduit de m_qtlmap_haplotype_ldla (tri_haplo_complet) utilise ngenom de m_qtlmap_haplotpe_data utilise presentg utilise nb_gam
SOURCE
1190 subroutine set_haplo_final(c,n,hsire,hdam) 1191 1192 logical, intent (in) :: hsire,hdam 1193 integer, intent (in) :: n,c 1194 integer :: stat 1195 ! divers 1196 integer ip,j,j_haplo,k_haplo,jm,geno,kd,j_gam 1197 integer :: nb_gam_t,i_gam,value_alloc,value_alloc_2 1198 real (kind=dp) :: su 1199 1200 if ( allocated(num_haplo_pere)) deallocate (num_haplo_pere) 1201 allocate (num_haplo_pere(np,2,nb_haplo_complet),STAT=stat) 1202 call check_allocate(stat,'num_haplo_pere [m_qtlmap_haplotype_ldla]') 1203 num_haplo_pere=0 1204 if ( allocated(num_haplo_mere)) deallocate (num_haplo_mere) 1205 allocate (num_haplo_mere(maxval(ngenom(:,nm+1)),2,nb_haplo_complet),STAT=stat) 1206 call check_allocate(stat,'num_haplo_mere [m_qtlmap_haplotype_ldla]') 1207 num_haplo_mere = 0 1208 value_alloc=nb_haplo_reduit*maxval(nb_gam) 1209 value_alloc_2=maxval(nall)**(2*longhap)+1 1210 if ( allocated(num_haplo_desc)) deallocate (num_haplo_desc) 1211 ! allocate (num_haplo_desc(nd,nb_haplo,nb_haplo_complet),STAT=stat) 1212 allocate (num_haplo_desc(nd,value_alloc),STAT=stat) 1213 call check_allocate(stat,'num_haplo_desc [m_qtlmap_haplotype_ldla]') 1214 num_haplo_desc = 0 1215 if ( allocated(pb_haplo_desc)) deallocate (pb_haplo_desc) 1216 ! allocate (pb_haplo_desc(nd,nb_haplo_reduit),STAT=stat) 1217 allocate (pb_haplo_desc(nd,value_alloc_2),STAT=stat) 1218 call check_allocate(stat,'pb_haplo_desc [m_qtlmap_haplotype_ldla]') 1219 pb_haplo_desc = 0.d0 1220 1221 j_haplo=1 1222 1223 ! 1224 ! pour les peres 1225 ! 1226 nb_gam_pere=0 1227 if(hsire)then 1228 1229 do ip=1,np 1230 do j=1,2 1231 j_haplo=j_haplo+1 1232 do k_haplo=1,nb_haplo_reduit 1233 if(comp_haplo_reduit(j_haplo,k_haplo)) then 1234 nb_gam_pere(ip,j)=nb_gam_pere(ip,j)+1 1235 num_haplo_pere(ip,j,nb_gam_pere(ip,j))= k_haplo 1236 end if 1237 end do !k_haplo 1238 end do !j 1239 end do !ip 1240 1241 end if ! opt_model(hsire) 1242 ! 1243 ! pour les meres 1244 ! 1245 nb_gam_mere=0 1246 nb_gam_desc=0 1247 if(hdam)then 1248 1249 do jm=1,nm 1250 if(estfem(repfem(jm)))then 1251 ! 1252 ! meres phasees 1253 ! 1254 do geno=ngenom(c,jm)+1,ngenom(c,jm+1) 1255 do j=1,2 1256 j_haplo=j_haplo+1 1257 do k_haplo=1,nb_haplo_reduit 1258 if(comp_haplo_reduit(j_haplo,k_haplo)) then 1259 nb_gam_mere(geno,j)=nb_gam_mere(geno,j)+1 1260 num_haplo_mere(geno,j,nb_gam_mere(geno,j))= k_haplo 1261 end if 1262 end do !k_haplo 1263 end do !j 1264 end do !geno 1265 else 1266 ! 1267 ! non phasees 1268 ! do kd=ndm(jm)+1,ndm(jm+1) 1269 ! if(presentg(c,kd))then 1270 ! do j_gam=1,nb_gam(kd) 1271 ! j_haplo=j_haplo+1 1272 ! do k_haplo=1,nb_haplo_reduit 1273 ! if(comp_haplo_reduit(j_haplo,k_haplo)) then 1274 ! nb_gam_desc(kd,j_gam)=nb_gam_desc(kd,j_gam)+1 1275 ! num_haplo_desc(kd,j_gam,nb_gam_desc(kd,j_gam))= k_haplo 1276 ! end if 1277 ! end do !k_haplo 1278 ! end do !j_gam 1279 ! end if 1280 ! end do !kd 1281 1282 do kd=ndm(jm)+1,ndm(jm+1) 1283 if(presentg(c,kd))then 1284 do j_gam=1,nb_gam(kd) 1285 j_haplo=j_haplo+1 1286 nb_gam_t=nb_gam_desc(kd)+1 1287 do k_haplo=1,nb_haplo_reduit 1288 if(comp_haplo_reduit(j_haplo,k_haplo)) then 1289 nb_gam_desc(kd)=nb_gam_desc(kd)+1 1290 num_haplo_desc(kd,nb_gam_desc(kd))=k_haplo 1291 end if 1292 end do !k 1293 su = dble(sum(pb_haplo_reduit(num_haplo_desc(kd,nb_gam_t:nb_gam_desc(kd))))) 1294 if ( su == 0 ) cycle 1295 do i_gam=nb_gam_t,nb_gam_desc(kd) 1296 ! pb_haplo_desc(kd,i_gam) = pb_haplo_desc(kd,i_gam) & 1297 ! +prob_gam(kd,j_gam) * (pb_haplo_reduit(num_haplo_desc(kd,i_gam)) /su) 1298 pb_haplo_desc(kd,num_haplo_desc(kd,i_gam)) = & 1299 pb_haplo_desc(kd,num_haplo_desc(kd,i_gam)) & 1300 +prob_gam(kd,j_gam) * (pb_haplo_reduit(num_haplo_desc(kd,i_gam)) /su) 1301 end do !i_gam 1302 end do !j_gam 1303 end if !presentg(c,kd) 1304 end do !kd 1305 1306 1307 end if ! estfem 1308 end do !jm 1309 end if ! opt_model (hdam) 1310 1311 end subroutine set_haplo_final
set_haplo_for_ldla
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
set_haplo_for_ldla
DESCRIPTION
NOTES
set_haplo construit la liste des haplotypes de longueur longhap Utilise nchr,posi,pas,nmk de m_qtlmap_data utilise longhap de m_qtlmap_data
SOURCE
333 subroutine set_haplo_for_ldla(c,dx,n,hsire,hdam) 334 integer c,i,j,k 335 double precision dx 336 integer, intent (in) :: n 337 logical,intent(in) :: hsire,hdam 338 339 if (.not. check_change_maker_windows(c,dx)) return 340 ! 341 ! longhap est la longueur de l haplotype dont on estime l effet 342 ! longhap est defini en nombre de marqueurs ( multiple de 2) 343 ! PROB_HAPLO_MIN=0.0d0 344 call free_internal_struct_ldla 345 call set_allocation_ldla 346 call liste_haplo(c,dx,n,hsire,hdam) 347 call liste_haplo_complet(c) 348 call proba_haplo_complet(c,dx) 349 call tri_haplo_complet(c,dx) 350 call set_haplo_final(c,n,hsire,hdam) 351 end subroutine set_haplo_for_ldla
set_tab_ibs
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
set_tab_ibs
DESCRIPTION
NOTES
set_tab_ibs construit la liste des haplotypes correspondant au chromosome entier
SOURCE
362 subroutine set_tab_ibs(c) 363 ! 364 ! set_tab_ibs construit la liste des haplotypes correspondant au chromosome entier 365 ! 366 ! Utilise nchr,posi,pas,nmk de m_qtlmap_data 367 ! utilise longhap de m_qtlmap_data 368 ! 369 integer, intent (in) :: c 370 ! 371 ! longhap est la longueur de l haplotype dont on estime l effet 372 ! longhap est defini en nombre de marqueurs ( multiple de 2) 373 call free_internal_struct_ldla 374 call set_allocation_ldla 375 call liste_gamete(c) 376 return 377 378 end subroutine set_tab_ibs
sort_haplo
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
sort_haplo
DESCRIPTION
NOTES
sort_haplo trie les haplotypes dans la liste haplo et les recodifie chaque reproducteur est affecte de ses haplotype recodifies la liste des noms d'haplotype (name_haplo) est cree Utilise np,nm,ndm,estfem,repfem de m_qtlmap_data utilise ngenom de m_qtlmap_haplotpe_data
SOURCE
394 subroutine sort_haplo(c,n,hsire,hdam) 395 logical :: hsire,hdam 396 ! dim allocatable 397 integer(kind=KIND_PHENO) , dimension (:,:) , allocatable :: kept_haplo 398 integer , dimension (:) , allocatable :: corr_haplo 399 integer , dimension (:) , allocatable :: corr_corr_haplo 400 integer , dimension (:) , allocatable :: inv_corr_haplo 401 ! divers 402 integer c 403 integer i,i_haplo,i_kept,ip,j,jm,j_gam,kd,geno,previous,j_haplo 404 integer nb_kept,value_alloc 405 integer :: stat,lstr1,lstr2 406 integer, intent (in) :: n 407 character(len=LEN_DEF) :: buf 408 double precision dencount 409 logical inconnu, new 410 411 ! allocation 412 value_alloc=nb_haplo 413 allocate (kept_haplo(value_alloc,longhap),STAT=stat) 414 call check_allocate(stat,'kept_haplo [m_qtlmap_haplotype_LD:sort_haplo]') 415 416 allocate (corr_haplo(value_alloc),STAT=stat) 417 call check_allocate(stat,'corr_haplo [m_qtlmap_haplotype_LD:sort_haplo]') 418 corr_haplo=0 419 420 allocate (corr_corr_haplo(value_alloc),STAT=stat) 421 call check_allocate(stat,'corr_corr_haplo [m_qtlmap_haplotype_LD:sort_haplo]') 422 423 allocate (inv_corr_haplo(value_alloc),STAT=stat) 424 call check_allocate(stat,'inv_corr_haplo [m_qtlmap_haplotype_LD:sort_haplo]') 425 426 427 428 ! on trie les haplotypes trouves en eliminant les redondances 429 ! 430 ! la liste reduite sera dans kept_haplo 431 ! de longueur nb_kept 432 ! le i eme ancien haplotype est a la corr_haplo(i) eme place dans kept_haplo 433 ! 434 name_haplo='' 435 kept_haplo(1,1)=haplo(1,1) 436 lstr1=len(trim(get_pheno(c,haplo(1,1)))) 437 name_haplo(1)(1:lstr1)=trim(get_pheno(c,haplo(1,1))) 438 439 do i=2,longhap 440 kept_haplo(1,i)=haplo(1,i) 441 lstr2=len(trim(get_pheno(c,haplo(1,i)))) 442 name_haplo(1)(lstr1+1:lstr1+lstr2)=trim(get_pheno(c,haplo(1,i))) 443 lstr1=lstr1+lstr2 444 end do !i 445 corr_haplo(1)=1 446 nb_kept=1 447 448 do i_haplo=2, nb_haplo 449 inconnu=.true. 450 do i_kept =1,nb_kept 451 new=.false. 452 do i=1,longhap 453 if(kept_haplo(i_kept,i) /= haplo(i_haplo,i))new=.true. 454 if(new)exit 455 end do 456 if(.not.new) previous=i_kept 457 if(.not.(inconnu.and.new))inconnu=.false. 458 end do ! i_kept 459 if(inconnu)then 460 nb_kept=nb_kept+1 461 kept_haplo(nb_kept,1)=haplo(i_haplo,1) 462 lstr1=len(trim(get_pheno(c,haplo(i_haplo,1)))) 463 name_haplo(nb_kept)(1:lstr1)=trim(get_pheno(c,haplo(i_haplo,1))) 464 do i=2,longhap 465 kept_haplo(nb_kept,i)=haplo(i_haplo,i) 466 lstr2=len(trim(get_pheno(c,haplo(i_haplo,i)))) 467 name_haplo(nb_kept)(lstr1+1:lstr1+lstr2)=trim(get_pheno(c,haplo(i_haplo,i))) 468 lstr1=lstr1+lstr2 469 buf=get_pheno(c,haplo(i_haplo,i)) 470 end do !i 471 corr_haplo(i_haplo)=nb_kept 472 else 473 corr_haplo(i_haplo)=previous 474 end if 475 end do !i_haplo 476 ! 477 ! Compteur d'haplotypes 478 ! 479 dencount=1.d0/dble(nb_haplo) 480 count_haplo(:)=0.d0 481 do i_haplo=1,nb_haplo 482 count_haplo(corr_haplo(i_haplo))=count_haplo(corr_haplo(i_haplo))+dencount 483 end do 484 485 nb_max_haplo=nb_kept 486 487 if(nb_max_haplo > NB_HAPLO_PRIOR)then 488 call stop_application('more haplotypes found than expected ['//trim(str(nb_haplo))//& 489 '] opt_nb_haplo_prior must be increased') 490 end if 491 492 ! 493 ! regroupement des haplotypes rares dans le paquet des donnees manquantes 494 ! 495 nb_kept=1 496 corr_corr_haplo(1)=1 497 inv_corr_haplo(1)=1 498 do i_haplo=2,nb_max_haplo 499 if(count_haplo(i_haplo) < PROB_HAPLO_MIN) then 500 corr_corr_haplo(i_haplo)=1 501 inv_corr_haplo(1)=1!peut etre que c est i_haplo.... 502 else 503 nb_kept=nb_kept+1 504 corr_corr_haplo(i_haplo)=nb_kept 505 inv_corr_haplo(nb_kept)=i_haplo 506 end if 507 end do ! i_haplo 508 count_haplo(:)=0.d0 509 do i_haplo=1,nb_haplo 510 j_haplo=corr_haplo(i_haplo) 511 count_haplo(corr_corr_haplo(j_haplo))=count_haplo(corr_corr_haplo(j_haplo))+dencount 512 end do 513 nb_max_haplo=nb_kept 514 515 ! creation d un indice d haplotype par reproducteur 516 ! 517 ! pour les peres 518 j_haplo=1 519 if(hsire)then 520 do ip=1,np 521 do j=1,2 522 j_haplo=j_haplo+1 523 num_haplo_pere(ip,j,1)=corr_corr_haplo(corr_haplo(j_haplo)) 524 end do !j 525 end do !ip 526 end if ! opt_model 527 528 if(hdam)then 529 ! 530 ! pour les meres 531 ! 532 do jm=1,nm 533 if(estfem(repfem(jm)))then 534 ! 535 ! meres phasees 536 ! 537 do geno=ngenom(c,jm)+1,ngenom(c,jm+1) 538 do j=1,2 539 j_haplo=j_haplo+1 540 num_haplo_mere(geno,j,1)=corr_corr_haplo(corr_haplo(j_haplo)) 541 end do !j 542 end do !geno 543 else 544 ! 545 ! non phasees 546 do kd=ndm(jm)+1,ndm(jm+1) 547 if(presentg(c,kd))then 548 do j_gam=1,nb_gam(kd) 549 j_haplo=j_haplo+1 550 ! num_haplo_desc(kd,j_gam,1)=corr_corr_haplo(corr_haplo(j_haplo)) 551 num_haplo_desc(kd,j_gam)=corr_corr_haplo(corr_haplo(j_haplo)) 552 end do !j_gam 553 end if 554 end do !kd 555 end if 556 end do !jm 557 end if ! opt_model 558 559 deallocate(kept_haplo) 560 deallocate(corr_haplo) 561 deallocate(corr_corr_haplo) 562 deallocate(inv_corr_haplo) 563 564 end subroutine sort_haplo
tri_haplo_complet
[ Top ] [ m_qtlmap_haplotype_ldla ] [ Subroutines ]
NAME
tri_haplo_complet
DESCRIPTION
NOTES
tri_haplo_complet elimine de la liste les haplotypes complets peu probables en entree, il lui faut nb_haplo_complet et tab_haplo_complet créés par la routine liste_haplo_complet du module m_qtlmap_haplotype_ldla haplo et nb_haplo créés par la routine set_allele_info_vector du module m_qtlmap_genotype en sortie il fournit pb_haplo_complet recalcule en, eliminant les haplotypes peu probables name_haplo_complet le nom compacte des haplotypes complets
SOURCE
1006 subroutine tri_haplo_complet(c,dx) 1007 integer , intent(in) :: c 1008 real(kind=dp) , intent(in) :: dx 1009 integer :: stat,lstr1,lstr2 1010 integer jhc,i,idx,ipb,khr,iht,ir 1011 integer nb_possible 1012 double precision tphc 1013 logical une_race_rare 1014 character(kind=1) :: nm_race_rare 1015 real(kind=dp) :: spb_haplo_reduit 1016 ! dim allocatable 1017 double precision , dimension (:) , allocatable :: phc 1018 ! allocation 1019 allocate (phc(nb_haplo_complet*NB_RACES),STAT=stat) 1020 call check_allocate(stat,'phc [m_qtlmap_haplotype_ldla:tri_haplo_complet]') 1021 allocate (liste_haplo_reduit(nb_haplo_complet*NB_RACES),STAT=stat) 1022 call check_allocate(stat,'liste_haplo_reduit [m_qtlmap_haplotype_ldla]') 1023 1024 1025 !correction BUG OFI : 14/09/2010 1026 if ( size (name_haplo_complet)<nb_haplo_complet ) then 1027 deallocate(name_haplo_complet) 1028 allocate (name_haplo_complet(nb_haplo_complet*NB_RACES)) 1029 end if 1030 1031 ! 1032 ! nom compacte des haplotypes mis dans name_haplo_complet 1033 ! 1034 name_haplo_complet='' 1035 1036 do jhc=1,nb_haplo_complet 1037 lstr1=len(trim(get_pheno(c,tab_haplo_complet(jhc,1)))) 1038 name_haplo_complet(jhc)(1:lstr1)=trim(get_pheno(c,tab_haplo_complet(jhc,1))) 1039 1040 do i=2,longhap 1041 lstr2=len(trim(get_pheno(c,tab_haplo_complet(jhc,i)))) 1042 name_haplo_complet(jhc)(lstr1+1:lstr1+lstr2)=trim(get_pheno(c,tab_haplo_complet(jhc,i))) 1043 lstr1=lstr1+lstr2 1044 end do !i 1045 1046 end do ! jhc 1047 ! 1048 ! comptage des haplotypes complets poossibles par haplotype trouvé 1049 ! 1050 do iht =1,nb_haplo 1051 nb_haplo_possible(iht)=0 1052 do jhc=1,nb_haplo_complet 1053 if(comp_haplo(iht,jhc)) nb_haplo_possible(iht)=nb_haplo_possible(iht)+1 1054 end do ! jhc 1055 end do !iht 1056 1057 ! regroupement des haplotypes rares dans le paquet des donnees manquantes 1058 ! on prepare le calcul des probabilites des haplotypes non rares (phc(jhc)) 1059 ! et on les compte (nb_haplo_reduit) 1060 ! 1061 phc=0.d0 1062 tphc=0.d0 1063 i=0 1064 nb_haplo_reduit=0 1065 liste_haplo_reduit=0 1066 spb_haplo_reduit=0.D0 1067 une_race_rare=.true. 1068 nm_race_rare='UNKNOWN' 1069 !print *,PROB_HAPLO_MIN 1070 do jhc=1,nb_haplo_complet 1071 if(pb_haplo_complet(jhc) > PROB_HAPLO_MIN) then 1072 phc(jhc)=pb_haplo_complet(jhc) 1073 tphc=tphc+phc(jhc) 1074 nb_haplo_reduit= nb_haplo_reduit+1 1075 liste_haplo_reduit(nb_haplo_reduit)=jhc 1076 spb_haplo_reduit=spb_haplo_reduit+phc(jhc) 1077 else 1078 1079 if (race_haplo_complet(jhc).ne.'UNKNOWN') i=i+1 1080 if (i==1) nm_race_rare=race_haplo_complet(jhc) 1081 if(nm_race_rare.ne.'UNKNOWN'.and.race_haplo_complet(jhc).ne.nm_race_rare) & 1082 une_race_rare=.false. 1083 end if 1084 end do ! jhc 1085 !print *,'somme des p_haplo=', spb_haplo_reduit , nb_haplo_complet 1086 ! 1087 ! arret si aucun haplotype possible 1088 ! 1089 if ( tphc == 0.d0) then 1090 return 1091 idx=int(100.*dx) 1092 ipb=int(100.*PROB_HAPLO_MIN) 1093 call stop_application(' No possible haplotypes found at location '//trim(str(idx))//& 1094 ' cM with a threshold '//trim(str(ipb))//' %') 1095 end if 1096 ! 1097 ! allocations selon le nobre d'haplotypes non rare 1098 ! 1099 if ( allocated(comp_haplo_reduit)) deallocate (comp_haplo_reduit) 1100 allocate (comp_haplo_reduit(nb_haplo,nb_haplo_reduit+nb_races),STAT=stat) 1101 call check_allocate(stat,'comp_haplo_reduit [m_qtlmap_haplotype_ldla]') 1102 1103 if ( allocated(pb_haplo_reduit)) deallocate (pb_haplo_reduit) 1104 allocate (pb_haplo_reduit(nb_haplo_reduit+nb_races),STAT=stat) 1105 call check_allocate(stat,'pb_haplo_reduit [m_qtlmap_haplotype_ldla]') 1106 1107 if ( allocated(name_haplo_reduit)) deallocate (name_haplo_reduit) 1108 allocate (name_haplo_reduit(nb_haplo_reduit+nb_races),STAT=stat) 1109 call check_allocate(stat,'name_haplo_reduit [m_qtlmap_haplotype_ldla]') 1110 ! 1111 ! creation des probabilites et correspondances 1112 ! 1113 pb_haplo_reduit =0.d0 1114 race_haplo_reduit='UNKNOWN' 1115 1116 do khr=1,nb_haplo_reduit 1117 pb_haplo_reduit(khr)=phc(liste_haplo_reduit(khr)) 1118 ! pb_haplo_reduit(khr)=phc(liste_haplo_reduit(khr))/tphc 1119 name_haplo_reduit(khr)=name_haplo_complet(liste_haplo_reduit(khr)) 1120 race_haplo_reduit(khr)=race_haplo_complet(liste_haplo_reduit(khr)) 1121 end do ! khr 1122 1123 1124 do iht =1,nb_haplo 1125 nb_possible=0 1126 do khr=1,nb_haplo_reduit 1127 comp_haplo_reduit(iht,khr)=comp_haplo(iht,liste_haplo_reduit(khr)) 1128 if(comp_haplo_reduit(iht,khr)) nb_possible=nb_possible+1 1129 end do ! khr 1130 if(nb_possible < nb_haplo_possible(iht))comp_haplo_reduit(iht,nb_haplo_reduit+1)=.true. 1131 end do !iht 1132 1133 if(nb_haplo_reduit < nb_haplo_complet) then 1134 if (raceFileDefined) then 1135 do ir=1, nb_races 1136 spb_haplo_reduit=0.d0 1137 nb_haplo_reduit= nb_haplo_reduit+1 1138 do jhc=1,nb_haplo_complet 1139 if ((race_haplo_complet(jhc)==nom_race(ir)) & 1140 .and.(pb_haplo_complet(jhc).LT.PROB_HAPLO_MIN)) & 1141 spb_haplo_reduit=spb_haplo_reduit+pb_haplo_complet(jhc) 1142 end do !jhc 1143 pb_haplo_reduit(nb_haplo_reduit)= spb_haplo_reduit 1144 name_haplo_reduit(nb_haplo_reduit)='RARE' 1145 race_haplo_reduit(nb_haplo_reduit)=nom_race(ir) 1146 1147 enddo 1148 else 1149 nb_haplo_reduit= nb_haplo_reduit+1 1150 pb_haplo_reduit(nb_haplo_reduit)=1.d0-tphc 1151 name_haplo_reduit(nb_haplo_reduit)='RARE' 1152 if (une_race_rare) race_haplo_reduit(nb_haplo_reduit)=nm_race_rare 1153 end if 1154 end if 1155 1156 nb_max_haplo=nb_haplo_reduit 1157 1158 ! write(nficout,*) 'Num_HAPLO RACE HAPLO pb_HAPLO' 1159 ! do i= 1, nb_haplo_reduit 1160 ! write(nficout,*) i,trim(race_haplo_reduit(i)),' ',& 1161 ! trim(name_haplo_reduit(i)), & 1162 ! pb_haplo_reduit(i) 1163 ! enddo 1164 1165 ! correctiopn BUG : leak memory 05/11/2010 1166 deallocate (phc) 1167 end subroutine tri_haplo_complet