m_qtlmap_genotype
NAME
m_qtlmap_genotype -- Genotype routines module
SYNOPSIS
DESCRIPTION
NOTES
SEE ALSO
unit_genotype
[ Top ] [ m_qtlmap_genotype ] [ Constants ]
NAME
unit_genotype
DESCRIPTION
index_position_marker
[ Top ] [ m_qtlmap_genotype ] [ Variables ]
NAME
index_position_marker
DESCRIPTION
give the position (column index) according a marker index
NOTES
each i element contains the position in the file of for the marker mark(i)
nmarker_in_genotype_file
[ Top ] [ m_qtlmap_genotype ] [ Variables ]
NAME
nmarker_in_genotype_file
DESCRIPTION
number of marker in the genotype file
NOTES
allocate_vector_
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
allocate_vector_
DESCRIPTION
Allocation of main arrays coming from m_qtlmap_data : corregp,corregm,correr,correp,correm,corred,presentg
SOURCE
320 subroutine allocate_vector_() 321 integer :: ios 322 323 call log_mess('SUBROUTINE : allocate_vector_correspond',DEBUG_DEF) 324 allocate (corregp(size(gpere)), stat=ios) 325 call check_allocate(ios,'corregp') 326 327 allocate (corregm(size(gmere)), stat=ios) 328 call check_allocate(ios,'corregm') 329 330 allocate (correr(size(repro)), stat=ios) 331 call check_allocate(ios,'correr') 332 333 allocate (correp(size(pere)), stat=ios) 334 call check_allocate(ios,'correp') 335 336 allocate (correm(size(mere)), stat=ios) 337 call check_allocate(ios,'correm') 338 339 allocate (corred(size(animal)), stat=ios) 340 call check_allocate(ios,'corred') 341 342 allocate (presentg(nchr,size(animal)), stat=ios) 343 call check_allocate(ios,'presentg') 344 345 call log_mess('END SUBROUTINE : allocate_vector_correspond',DEBUG_DEF) 346 end subroutine allocate_vector_
check_HWE
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
check_HWE
DESCRIPTION
check_Hardy Weinberg : check the equilibrium of marker transmission within each family
INPUTS
ndmin : minimum number of progenies to build full sib family
SOURCE
1502 subroutine check_HWE 1503 1504 integer :: ll,ip,jm,kd,nb_hete_prog,nb_known_prog 1505 integer :: c,borne(2) 1506 1507 call log_mess('checking HWE in heterozygous sires...',INFO_DEF) 1508 1509 write (nficout,*) 1510 write (nficout,*) "** Check the equilibrium of marker transmission within each family **" 1511 write (nficout,*) 1512 1513 do c=1,nchr 1514 !for each marker 1515 do ll=1,nmk(c) 1516 if ( nall(c,ll) > 2 ) cycle 1517 do ip=1,np ! for each sire 1518 1519 if( pheno(c,ll,correp(ip),1) /= pheno(c,ll,correp(ip),2)) then 1520 nb_known_prog=0 1521 nb_hete_prog=0 1522 1523 do jm=nmp(ip)+1,nmp(ip+1) 1524 do kd=ndm(jm)+1,ndm(jm+1) 1525 if ( corred(kd) /= INT_NOT_DEFINED ) then 1526 if ( pheno(c,ll,corred(kd),1) /= nmanque .and. pheno(c,ll,corred(kd),2) /= nmanque) then 1527 nb_known_prog=nb_known_prog+1 1528 if (pheno(c,ll,corred(kd),1) /= pheno(c,ll,corred(kd),2)) nb_hete_prog=nb_hete_prog+1 1529 end if 1530 end if 1531 end do ! end progeny 1532 end do ! end dam 1533 1534 call confbin(nb_known_prog,pseuilHWE,borne) 1535 1536 if(nb_hete_prog < borne(1) .or. nb_hete_prog > borne(2)) then 1537 write(nficout,*)'Marker ['//trim(mark(c,ll))//'] for sire :'//trim(pere(ip))//& 1538 ' not in HWE : ',nb_hete_prog, ' heterozygous progeny amongst ',nb_known_prog 1539 ! call log_mess('Marker ['//trim(mark(c,ll))//'] for sire :'//trim(pere(ip))//& 1540 ! ' not in HWE ... should be excluded of the analysis' ,WARNING_DEF) 1541 1542 end if 1543 end if 1544 end do ! end sire 1545 end do ! end marker 1546 end do ! end chr 1547 1548 call log_mess('check_HWE ok',INFO_DEF) 1549 end subroutine
check_marker_name
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
check_marker_name
DESCRIPTION
read and check the header genotypic file. fill internal arrays index_chr,index_position_marker,
INPUTS
file : genotypic file
SOURCE
110 subroutine check_marker_name(file) 111 character(len=LENGTH_MAX_FILE),intent(in) :: file 112 character(len=LEN_DEF) :: token 113 integer :: ios,i,j,k,c,long,s1 114 logical :: is_ok 115 integer, parameter :: LENBUF=LEN_LINE*10 116 character(len=LENBUF) ,pointer :: line_read 117 118 character(len=LEN_DEF) , dimension (:),allocatable ::mark00 119 120 call log_mess('SUBROUTINE : check_marker_name',DEBUG_DEF) 121 122 rewind (unit_genotype) 123 124 ios = 0 125 allocate (line_read) 126 line_read = '' 127 128 do while ( trim(line_read) == '' .and. (ios == 0) ) 129 call GET(unit_genotype,line_read,IOSTAT=ios) 130 end do 131 132 if (ios /= 0) then 133 call stop_application('No marker name are found in the genotype file. [file:'//trim(file)//']') 134 end if 135 136 allocate (mark00(50000)) 137 138 line_read=trim(line_read) 139 long = len(trim(line_read)) 140 nmarker_in_genotype_file=1 141 is_ok=.false. 142 s1=1 143 do i=1,long 144 if ( line_read(i:i) == ' ') then 145 if (is_ok) then 146 mark00(nmarker_in_genotype_file)=adjustl(line_read(s1:i)) 147 s1=i 148 nmarker_in_genotype_file=nmarker_in_genotype_file+1 149 is_ok=.false. 150 end if 151 else if ( .not. is_ok ) then 152 is_ok=.true. 153 end if 154 end do 155 mark00(nmarker_in_genotype_file)=adjustl(line_read(s1:long)) 156 !stop 157 call log_mess('number of marker in the header genotype file:'// trim(str(nmarker_in_genotype_file)) , VERBOSE_DEF) 158 159 allocate (index_position_marker(nmarker_in_genotype_file),stat=ios) 160 call check_allocate(ios,'index_position_marker') 161 162 allocate (index_chr(nmarker_in_genotype_file),stat=ios) 163 call check_allocate(ios,'index_chr') 164 165 index_position_marker = INT_NOT_DEFINED 166 index_chr = 0 167 168 do c=1,nchr 169 do i=1,nmk(c) 170 ! we search the marker associated 171 j = 1 172 do while ( j <= size(mark00) .and. mark(c,i) /= mark00(j) ) 173 !print *,mark(c,i), mark00(j) 174 j = j+1 175 if ( j > size(mark00) ) exit 176 end do 177 178 if ( j > size(mark00) ) then 179 call log_mess('bad construction of genotype file. a marker from map file' & 180 //' is not finding on the genotype file ['//trim(mark(c,i))//']',ERROR_DEF) 181 stop 1 182 end if 183 call log_mess( 'index in genotype file:'//trim(str(j))//& 184 ' --> index in analyse:'//trim(str(i)),DEBUG_DEF) 185 index_position_marker(j) = i 186 index_chr(j) = c 187 call log_mess('marker from map file ['//trim(mark(c,i))//'] is finding in position:'// trim(str(j)),VERBOSE_DEF) 188 end do 189 end do 190 191 deallocate (mark00) 192 deallocate (line_read) 193 194 call log_mess('END SUBROUTINE : check_marker_name',DEBUG_DEF) 195 end subroutine check_marker_name
check_typage
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
check_typage
DESCRIPTION
check_typage
check phenotype marker of the progeny
p1,p2 phenotype marker of the kd progeny
sp1,sp2 phenotype marker of kd's sire
dp1,dp2 phenotype marker of kd's dam
if p1==sp1 .or. sp2 then p2==dp1 .or. dp2
if p1==dp1 .or. dp2 then p2==sp1 .or. sp2
SOURCE
622 subroutine check_typage 623 integer :: ll,ip,jm,kd,igp,jgm,kr 624 logical :: typp,typm,typgp,typgm 625 integer :: i,c,nberr,typerr 626 integer ,parameter :: MAX_ERR=100 627 integer (kind=KIND_PHENO) :: error1(MAX_ERR,2),error2(MAX_ERR,2),error3(MAX_ERR,2) 628 integer :: kderror(MAX_ERR),iperror(MAX_ERR),jmerror(MAX_ERR),chrerror(MAX_ERR) 629 character(len=LEN_DEF) :: merror(MAX_ERR) 630 631 call log_mess('checking phenotype marker of progenies...',INFO_DEF) 632 633 nberr=0 634 do c=1,nchr 635 !for each marker 636 do ll=1,nmk(c) 637 do ip=1,np ! for each sire 638 typp=.false. 639 if(correp(ip) /= INT_NOT_DEFINED)then 640 if( pheno(c,ll,correp(ip),1) /= nmanque ) then 641 typp=.true. 642 endif 643 end if 644 645 do jm=nmp(ip)+1,nmp(ip+1) 646 typm=.false. 647 if(correm(jm) /= INT_NOT_DEFINED)then 648 if(pheno(c,ll,correm(jm),1) /= nmanque) then 649 typm=.true. 650 endif 651 end if 652 653 ! MODIF OFI 21 janv 2010 : (.not. typp) .or. (.not. typm) => (.not. typp) .and. (.not. typm) 654 !no information about parents 655 if ( (.not. typp) .and. (.not. typm) ) then 656 call log_mess('none information find about sire ['//trim(pere(ip))//& 657 '] dam ['//trim(mere(jm))//'] for marker ['//trim(mark(c,ll))//']',DEBUG_DEF) 658 cycle 659 end if 660 661 !checking progenies 662 do kd=ndm(jm)+1,ndm(jm+1) 663 if ( corred(kd) == INT_NOT_DEFINED ) then 664 cycle 665 end if 666 667 if ( pheno(c,ll,corred(kd),1) == nmanque .or. pheno(c,ll,corred(kd),2) == nmanque) then 668 cycle 669 end if 670 671 if (.not. typp .or. pheno(c,ll,corred(kd),1) == pheno(c,ll,correp(ip),1) .or. & 672 pheno(c,ll,corred(kd),1) == pheno(c,ll,correp(ip),2)) then 673 !if ( correm(jm) == INT_NOT_DEFINED ) cycle 674 if (.not. typm .or. (pheno(c,ll,corred(kd),2) == pheno(c,ll,correm(jm),1)) .or. & 675 (pheno(c,ll,corred(kd),2) == pheno(c,ll,correm(jm),2)) ) then 676 cycle 677 end if 678 end if 679 680 if ( .not. typm .or. & 681 pheno(c,ll,corred(kd),1) == pheno(c,ll,correm(jm),1) .or. & 682 pheno(c,ll,corred(kd),1) == pheno(c,ll,correm(jm),2)) then 683 684 if ( .not. typp .or. & 685 (pheno(c,ll,corred(kd),2) == pheno(c,ll,correp(ip),1)) .or. & 686 (pheno(c,ll,corred(kd),2) == pheno(c,ll,correp(ip),2)) ) then 687 cycle 688 end if 689 end if 690 691 if (nberr >= MAX_ERR) exit 692 nberr=nberr+1 693 error1(nberr,1)=pheno(c,ll,correp(ip),1) 694 error2(nberr,1)=pheno(c,ll,correm(jm),1) 695 error3(nberr,1)=pheno(c,ll,corred(kd),1) 696 error1(nberr,2)=pheno(c,ll,correp(ip),2) 697 error2(nberr,2)=pheno(c,ll,correm(jm),2) 698 error3(nberr,2)=pheno(c,ll,corred(kd),2) 699 kderror(nberr)=kd 700 iperror(nberr)=ip 701 jmerror(nberr)=jm 702 merror(nberr)=mark(c,ll) 703 chrerror(nberr)=c 704 705 end do ! end progeny 706 end do ! end dam 707 end do ! end sire 708 end do ! end marker 709 end do ! end chr 710 typerr=nberr 711 if ( nberr < MAX_ERR ) then 712 !! CHECKING Reproducter.......... 713 do c=1,nchr 714 do ll=1,nmk(c) 715 do igp=1,ngp 716 typgp=.false. 717 if(corregp(igp) /= INT_NOT_DEFINED) then 718 if(pheno(c,ll,corregp(igp),1) /= nmanque) then 719 typgp=.true. 720 end if 721 end if 722 723 do jgm=ngmgp(igp)+1,ngmgp(igp+1) 724 typgm = .false. 725 if(corregm(jgm) /= INT_NOT_DEFINED) then 726 if(pheno(c,ll,corregm(jgm),1) /= nmanque) then 727 typgm=.true. 728 end if 729 end if 730 !no information about parents 731 if ( (.not. typgp) .or. (.not. typgm) ) then 732 call log_mess('none information find about grand sire ['//trim(gpere(igp))//& 733 '] grand dam ['//trim(gmere(jgm))//'] for marker ['//trim(mark(c,ll))//']',DEBUG_DEF) 734 cycle 735 end if 736 737 do kr=nrgm(jgm)+1,nrgm(jgm+1) 738 if ( correr(kr) == INT_NOT_DEFINED ) then 739 cycle 740 end if 741 742 if ( pheno(c,ll,correr(kr),1) == nmanque .or. pheno(c,ll,correr(kr),2) == nmanque) then 743 cycle 744 end if 745 if ( .not. typgp .or. & 746 pheno(c,ll,correr(kr),1) == pheno(c,ll,corregp(igp),1) .or. & 747 pheno(c,ll,correr(kr),1) == pheno(c,ll,corregp(igp),2)) then 748 if ( .not. typgm .or. & 749 (pheno(c,ll,correr(kr),2) == pheno(c,ll,corregm(jgm),1)) .or. & 750 (pheno(c,ll,correr(kr),2) == pheno(c,ll,corregm(jgm),2)) ) then 751 cycle 752 end if 753 end if 754 755 if ( .not. typgm .or. & 756 pheno(c,ll,correr(kr),1) == pheno(c,ll,corregm(jgm),1) .or. & 757 pheno(c,ll,correr(kr),1) == pheno(c,ll,corregm(jgm),2)) then 758 759 if ( .not. typgp .or. & 760 (pheno(c,ll,correr(kr),2) == pheno(c,ll,corregp(igp),1)) .or. & 761 (pheno(c,ll,correr(kr),2) == pheno(c,ll,corregp(igp),2)) ) then 762 cycle 763 end if 764 end if 765 if (nberr >= MAX_ERR) exit 766 nberr=nberr+1 767 error1(nberr,1)=pheno(c,ll,corregp(igp),1) 768 error2(nberr,1)=pheno(c,ll,corregm(jgm),1) 769 error3(nberr,1)=pheno(c,ll,correr(kr),1) 770 error1(nberr,2)=pheno(c,ll,corregp(igp),2) 771 error2(nberr,2)=pheno(c,ll,corregm(jgm),2) 772 error3(nberr,2)=pheno(c,ll,correr(kr),2) 773 kderror(nberr)=kr 774 iperror(nberr)=igp 775 jmerror(nberr)=jgm 776 merror(nberr)=mark(c,ll) 777 chrerror(nberr)=c 778 end do ! end repro 779 end do ! end grand dam 780 end do ! end grand sire 781 end do ! end marker 782 end do ! end chr 783 end if 784 785 if ( nberr > 0 ) then 786 do i=1,typerr 787 c=chrerror(i) 788 call log_mess(' ---------- ',ERROR_DEF) 789 call log_mess('Phenotype sire :'//trim(get_pheno(c,error1(i,1)))//' '//& 790 trim(get_pheno(c,error1(i,2))),ERROR_DEF) 791 call log_mess('Phenotype dam :'//trim(get_pheno(c,error2(i,1)))//' '//& 792 trim(get_pheno(c,error2(i,2))),ERROR_DEF) 793 call log_mess('Phenotype progeny :'//trim(get_pheno(c,error3(i,1)))//' '//& 794 trim(get_pheno(c,error3(i,2))),ERROR_DEF) 795 call log_mess('progeny ['//trim(animal(kderror(i)))//'] from sire ['//trim(pere(iperror(i)))//& 796 '] dam ['//trim(mere(jmerror(i)))//'] for marker ['//trim(merror(i))//']',ERROR_DEF) 797 798 end do 799 800 do i=typerr+1,nberr 801 c=chrerror(i) 802 call log_mess(' ---------- ',ERROR_DEF) 803 call log_mess('Phenotype grand sire :'//trim(get_pheno(c,error1(i,1)))//' '//& 804 trim(get_pheno(c,error1(i,2))),ERROR_DEF) 805 call log_mess('Phenotype grand dam :'//trim(get_pheno(c,error2(i,1)))//' '//& 806 trim(get_pheno(c,error2(i,2))),ERROR_DEF) 807 call log_mess('Phenotype reproductor :'//trim(get_pheno(c,error3(i,1)))//' '//& 808 trim(get_pheno(c,error3(i,2))),ERROR_DEF) 809 call log_mess('reproductor ['//trim(repro(kderror(i)))//'] from grand sire ['//trim(gpere(iperror(i)))//& 810 '] grand dam ['//trim(gmere(jmerror(i)))//'] for marker ['//trim(merror(i))//']',ERROR_DEF) 811 end do 812 call stop_application(' ** NB ERROR :'// trim(str(nberr))//'** ') 813 else 814 call log_mess('ok',INFO_DEF) 815 end if 816 end subroutine check_typage
confbin
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
confbin
DESCRIPTION
INPUTS
n : p : borne :
SOURCE
1563 subroutine confbin(n,p,borne) 1564 real(kind=dp) :: p,c,seuil,tc 1565 integer :: n,m,l,borne(2) 1566 1567 seuil=log(p/2.d0)+dble(n) *log(2.d0) 1568 tc=1.d0 1569 c=1.d0 1570 m=n 1571 l=1 1572 1573 do while (log(tc) < seuil) 1574 c=c *dble(m) / dble(l) 1575 tc=tc+c 1576 l=l+1.d0 1577 m=m-1.d0 1578 end do 1579 1580 borne(1)=n-m 1581 borne(2)=m 1582 1583 end subroutine confbin
internal_read_file
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
internal_read_file
DESCRIPTION
read the genotypic file. This file contains the animals phenotypes at the markers. The first line gives the marker names, the markers must belong to the marker map file. For each animal, a line gives : * its ID (as decribed in the pedigree file) * the markers phenotypes, ranked following in the first line order . Each phenotype is made of 2 alleles, unordered. When an animal has no phenotype for a marker, both alleles must be given the missing value code as given in the parametrisation of the analysis.
INPUTS
file : genotypic file
SOURCE
214 subroutine internal_read_file(file) 215 character(len=*),intent(in) :: file 216 character(len=LEN_DEF) :: buffer_char,name 217 character(len=LEN_DEF) :: val_nmanque,w 218 integer :: eof,ios,l,i_ind,i,j,k 219 logical :: is_ok 220 character(len=LEN_DEF), dimension(:,:,:) , allocatable :: temp_marker_list 221 character(len=LEN_DEF), dimension(:,:,:,:) , allocatable :: pheno_t 222 223 call log_mess('SUBROUTINE : internal_read_file',DEBUG_DEF) 224 rewind (unit_genotype) 225 226 eof = 0 227 nmes = 0 228 229 do while ( eof == 0 ) 230 read(unit_genotype,*,iostat=eof) buffer_char 231 if ( trim(buffer_char) /= '' .and. eof == 0 ) then 232 nmes = nmes+1 233 end if 234 end do 235 236 ! the first line is the header (marker name) 237 nmes = nmes - 1 238 call log_mess('find ['//trim(str(nmes))//'] animals in the genotype file.',INFO_DEF) 239 240 ! fix bug : parents and grand parent can be add later.... 241 allocate ( numero(nmes+ngp+ngm+np+nm) , stat=ios) 242 call check_allocate(ios,'numero') 243 numero=STRING_NOT_DEFINED 244 245 allocate ( pheno(nchr,maxval(nmk),size(numero),2) , stat=ios) 246 call check_allocate(ios,'pheno') 247 pheno=nmanque 248 249 allocate (temp_marker_list(nmes,nmarker_in_genotype_file,2), stat=ios) 250 call check_allocate(ios,'temp_marker_list') 251 252 rewind (unit_genotype) 253 read(unit_genotype,*,iostat=eof) buffer_char 254 !Recupere les nom de tous les individus pour l affichage d erreur 255 i=1 256 do while ( eof == 0 .and. i <= nmes ) 257 read (unit_genotype,*,iostat=eof) numero(i),((temp_marker_list(i,k,j),j=1,2),k=1,nmarker_in_genotype_file) 258 if ( trim(numero(i)) /= '' .and. eof == 0 ) then 259 i = i+1 260 end if 261 end do 262 263 allocate (pheno_t(nchr,maxval(nmk),size(numero),2)) 264 !first encode, unknown char 265 val_nmanque = get_pheno(1,nmanque) 266 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i) 267 !$OMP DO 268 do i_ind=1,nmes 269 !****************************************************** 270 ! 05/07/2010 : Test si l individu a deja ete lu (plusieurs definition de typage existe) 271 ! test beaucoup trop long 272 ! do i=1,(i_ind-1) 273 ! if ( numero(i_ind) == numero(i)) then 274 ! call stop_on_error (1,file,l,'animal ['//trim(numero(i_ind))//& 275 ! '] have two lines definitions of phenotype marker !!') 276 ! end if 277 ! end do 278 279 ! save now in the pheno structure.... 280 do i=1,nmarker_in_genotype_file ! for all markers in the genotype file 281 if ( index_position_marker(i) == INT_NOT_DEFINED) cycle 282 if ( trim(temp_marker_list(i_ind,i,1)) == val_nmanque .or. trim(temp_marker_list(i_ind,i,2)) == val_nmanque) then 283 if ( .not. (trim(temp_marker_list(i_ind,i,1)) == val_nmanque .and. trim(temp_marker_list(i_ind,i,2)) == val_nmanque) ) then 284 call log_mess("*** Only one phenotype marker for [ animal:"//trim(numero(i_ind))//& 285 "] [marker:"// trim(mark(index_chr(i),index_position_marker(i)))//& 286 "] is not allowed *** [setting -->" & 287 //trim(val_nmanque)//" "//trim(val_nmanque)//"]",WARNING_DEF) 288 end if 289 pheno(index_chr(i),index_position_marker(i),i_ind,1) = nmanque 290 pheno(index_chr(i),index_position_marker(i),i_ind,2) = nmanque 291 else 292 !$OMP CRITICAL 293 !w=trim(temp_marker_list(i_ind,i,1)) 294 pheno(index_chr(i),index_position_marker(i),i_ind,1) = set_pheno(index_chr(i),temp_marker_list(i_ind,i,1)) 295 !w=trim(temp_marker_list(i_ind,i,2)) 296 pheno(index_chr(i),index_position_marker(i),i_ind,2) = set_pheno(index_chr(i),temp_marker_list(i_ind,i,2)) 297 !w='' 298 !$OMP END CRITICAL 299 end if 300 end do 301 end do 302 !$OMP END DO 303 !$OMP END PARALLEL 304 305 deallocate (temp_marker_list) 306 deallocate (index_position_marker) 307 deallocate (index_chr) 308 call log_mess('number of animals defined in genotype file:'// trim(str(nmes)),VERBOSE_DEF) 309 call log_mess('END SUBROUTINE : internal_read_file',DEBUG_DEF) 310 end subroutine internal_read_file
read_genotype
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
read_genotype
DESCRIPTION
Read the user genotypic file
INPUTS
ndmin : minimum number of progenies to build full sib family
NOTES
ndmin used by set_estfem
SOURCE
76 subroutine read_genotype(ndmin) 77 integer ,intent(in) :: ndmin 78 integer :: ios 79 call log_mess('START SUBROUTINE : read_genotype',DEBUG_DEF) 80 81 open (unit_genotype,FILE=in_typage,IOSTAT=ios,recl=2**20) 82 83 if ( ios /= 0 ) then 84 call log_mess('Cannot open genotype file ['//trim(in_typage)//']',ERROR_DEF) 85 stop 1 86 endif 87 88 call log_mess('reading genotype file...',INFO_DEF) 89 call allocate_vector_() 90 call check_marker_name(in_typage) 91 call internal_read_file(in_typage) 92 close(unit_genotype) 93 call set_corresponding_vector() 94 call set_allele_info_vector() 95 call set_estfem(ndmin) 96 97 call log_mess('END SUBROUTINE : read_genotype',INFO_DEF) 98 end subroutine read_genotype
recup
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
recup
DESCRIPTION
NOTES
Deduction des phenotypes des reproducteurs non types d apres leur descendance
SOURCE
1076 subroutine recup 1077 1078 ! Divers 1079 logical typp,typm,typgp,typgm 1080 integer :: ll,igp,nagp,ngm1,ngm2,jgm,nagm,nr1,nr2,kr,il 1081 integer :: jmes,imes,ip,nap,nm1,nm2,jm,nam,nd1,nd2,kd,c 1082 integer(kind=KIND_PHENO) :: m1,m2,mc1,mc2,md1,md2,mx 1083 1084 ! 1085 !****************************************************************************** 1086 ! Reconstitution des phenotypes des grands parents 1087 !****************************************************************************** 1088 ! 1089 1090 call log_mess('finding phenotype marker...',INFO_DEF) 1091 ! Recherche des phenotypes des grands parents 1092 do c=1,nchr 1093 do ll=1,nmk(c) 1094 do 10 igp=1,ngp 1095 if(corregp(igp).eq.9999) then 1096 typgp=.false. 1097 nagp=0 1098 else if(pheno(c,ll,corregp(igp),1).eq.nmanque) then 1099 typgp=.false. 1100 nagp=0 1101 else 1102 typgp=.true. 1103 nagp=2 1104 m1=pheno(c,ll,corregp(igp),1) 1105 m2=pheno(c,ll,corregp(igp),2) 1106 end if 1107 ngm1=ngmgp(igp)+1 1108 ngm2=ngmgp(igp+1) 1109 do 11 jgm=ngm1,ngm2 1110 if(corregm(jgm).eq.9999) then 1111 typgm=.false. 1112 nagm=0 1113 else if(pheno(c,ll,corregm(jgm),1).eq.nmanque) then 1114 typgm=.false. 1115 nagm=0 1116 else 1117 typgm=.true. 1118 nagm=2 1119 mc1=pheno(c,ll,corregm(jgm),1) 1120 mc2=pheno(c,ll,corregm(jgm),2) 1121 end if 1122 if(typgp.and.typgm)go to 11 1123 1124 ! if ( corregp(igp) /= 9999 ) then 1125 ! call log_mess('Grand sire ['//trim(gpere(igp))//'] -->'//trim(get_pheno(c,pheno(c,ll,corregp(igp),1)))//& 1126 ! ' '//trim(get_pheno(c,pheno(c,ll,corregp(igp),2)))//' marker:'//trim(mark(c,ll)),DEBUG_DEF) 1127 ! end if 1128 ! 1129 ! if ( corregm(jgm) /= 9999 ) then 1130 ! call log_mess('Grand dam ['//trim(gmere(jgm))//'] -->'//trim(get_pheno(c,pheno(c,ll,corregm(jgm),1)))//' '//& 1131 ! trim(get_pheno(c,pheno(c,ll,corregm(jgm),2)))//' marker:'//trim(mark(c,ll)),DEBUG_DEF) 1132 ! end if 1133 ! 1134 ! Passage en revue des parents 1135 nr1=nrgm(jgm)+1 1136 nr2=nrgm(jgm+1) 1137 do 12 kr=nr1,nr2 1138 if(correr(kr).eq.9999) go to 12 1139 if(pheno(c,ll,correr(kr),1).eq.nmanque) go to 12 1140 md1=pheno(c,ll,correr(kr),1) 1141 md2=pheno(c,ll,correr(kr),2) 1142 mx=nmanque 1143 1144 ! if ( correr(kr) /= 9999 ) then 1145 ! call log_mess('Repro ['//trim(repro(kr))//'] -->'//trim(get_pheno(c,pheno(c,ll,correr(kr),1)))//' '//& 1146 ! trim(get_pheno(c,pheno(c,ll,correr(kr),2)))//' marker:'//trim(mark(c,ll)),DEBUG_DEF) 1147 ! end if 1148 ! 1149 ! Parent homozygote 1150 ! Grand pere inconnu 1151 if(md1.eq.md2)then 1152 ! call log_mess('Reproductor homo ['//trim(get_pheno(c,md1))//' '//trim(get_pheno(c,md2))//']',DEBUG_DEF) 1153 if(nagp.ne.2) then 1154 if(nagp.eq.0) then 1155 nagp=1 1156 m1=md1 1157 ! call log_mess('finding one allele for grand sire ['//trim(get_pheno(c,m1))//']',DEBUG_DEF) 1158 else 1159 if(md1.ne.m1) then 1160 nagp=2 1161 m2=md1 1162 ! call log_mess('finding one allele for grand sire ['//trim(get_pheno(c,m2))//']',DEBUG_DEF) 1163 end if 1164 end if 1165 end if 1166 ! 1167 ! Grand mere inconnue 1168 if(nagm.ne.2) then 1169 if(nagm.eq.0) then 1170 nagm=1 1171 mc1=md1 1172 ! call log_mess('finding one allele for grand dam ['//trim(get_pheno(c,mc1))//']',DEBUG_DEF) 1173 else 1174 if(md1.ne.mc1) then 1175 nagm=2 1176 mc2=md1 1177 ! call log_mess('finding one allele for grand dam ['//trim(get_pheno(c,mc2))//']',DEBUG_DEF) 1178 end if 1179 end if 1180 end if 1181 ! 1182 ! Parent heterozygote 1183 ! Grand mere connue et grand pere inconnu 1184 else 1185 ! call log_mess('Reproductor hetero ['//trim(get_pheno(c,md1))//','//trim(get_pheno(c,md2))//']',DEBUG_DEF) 1186 if(nagm.eq.2.and.nagp.ne.2) then 1187 if((md1.eq.mc1.and.md2.ne.mc2).or.(md1.eq.mc2.and.md2.ne.mc1)) mx=md2 1188 if((md2.eq.mc1.and.md1.ne.mc2).or.(md2.eq.mc2.and.md1.ne.mc1)) mx=md1 1189 if(mx.ne.nmanque)then 1190 if(nagp.eq.0) then 1191 nagp=1 1192 m1=mx 1193 ! call log_mess('finding one allele for grand sire ['//trim(get_pheno(c,m1))//']',DEBUG_DEF) 1194 else 1195 if(mx.ne.m1) then 1196 nagp=2 1197 m2=mx 1198 ! call log_mess('finding one allele for grand sire ['//trim(get_pheno(c,m2))//']',DEBUG_DEF) 1199 end if 1200 end if 1201 end if 1202 end if 1203 ! 1204 ! Grand pere connu et grand mere inconnue 1205 if(nagp.eq.2.and.nagm.ne.2) then 1206 if((md1.eq.m1.and.md2.ne.m2).or.(md1.eq.m2.and.md2.ne.m1)) mx=md2 1207 if((md2.eq.m1.and.md1.ne.m2).or.(md2.eq.m2.and.md1.ne.m1)) mx=md1 1208 if(mx.ne.nmanque)then 1209 if(nagm.eq.0) then 1210 nagm=1 1211 mc1=mx 1212 ! call log_mess('finding one allele for grand dam ['//trim(get_pheno(c,mc1))//']',DEBUG_DEF) 1213 else 1214 if(mx.ne.mc1) then 1215 nagm=2 1216 mc2=mx 1217 ! call log_mess('finding one allele for grand dam ['//trim(get_pheno(c,mc2))//']',DEBUG_DEF) 1218 end if 1219 end if 1220 end if 1221 end if 1222 end if 1223 12 continue 1224 ! 1225 ! Stockage du phenotype reconstitue de la grand mere 1226 if(typgm.or.nagm.lt.2) go to 11 1227 if(corregm(jgm).eq.9999) then 1228 nmes=nmes+1 1229 numero(nmes)=gmere(jgm) 1230 corregm(jgm)=nmes 1231 pheno(c,:,corregm(jgm),1)=nmanque 1232 pheno(c,:,corregm(jgm),2)=nmanque 1233 end if 1234 jmes=corregm(jgm) 1235 pheno(c,ll,jmes,1)=mc1 1236 pheno(c,ll,jmes,2)=mc2 1237 1238 call log_mess('setting phenotype marker for grand dam ['//trim(gmere(jgm))//'] :'//& 1239 trim(get_pheno(c,pheno(c,ll,corregm(jgm),1)))//' '//trim(get_pheno(c,pheno(c,ll,corregm(jgm),2)))& 1240 //' marker ['//trim(mark(c,ll))//']',DEBUG_DEF) 1241 1242 11 continue 1243 ! 1244 ! Stockage du phenotype reconstitue du grand pere 1245 if(typgp.or.nagp.lt.2) go to 10 1246 if(corregp(igp).eq.9999) then 1247 nmes=nmes+1 1248 numero(nmes)=gpere(igp) 1249 corregp(igp)=nmes 1250 pheno(c,:,corregp(igp),1)=nmanque 1251 pheno(c,:,corregp(igp),2)=nmanque 1252 end if 1253 imes=corregp(igp) 1254 pheno(c,ll,imes,1)=m1 1255 pheno(c,ll,imes,2)=m2 1256 ! call log_mess('setting phenotype marker for grand sire ['//trim(gpere(igp))//'] :'//& 1257 ! trim(get_pheno(c,pheno(c,ll,corregp(igp),1)))//' '//trim(get_pheno(c,pheno(c,ll,corregp(igp),2)))//& 1258 ! ' marker ['//trim(mark(c,ll))//']',DEBUG_DEF) 1259 10 continue 1260 end do 1261 ! 1262 !****************************************************************************** 1263 ! Reconstitution des phenotypes des parents 1264 !****************************************************************************** 1265 ! 1266 ! Recherche des phenotypes des parents 1267 do ll=1,nmk(c) 1268 do 20 ip=1,np 1269 if(correp(ip).eq.9999) then 1270 typp=.false. 1271 nap=0 1272 else if (pheno(c,ll,correp(ip),1).eq.nmanque) then 1273 typp=.false. 1274 nap=0 1275 else 1276 typp=.true. 1277 nap=2 1278 m1=pheno(c,ll,correp(ip),1) 1279 m2=pheno(c,ll,correp(ip),2) 1280 end if 1281 nm1=nmp(ip)+1 1282 nm2=nmp(ip+1) 1283 do 21 jm=nm1,nm2 1284 if(correm(jm).eq.9999) then 1285 typm=.false. 1286 nam=0 1287 else if(pheno(c,ll,correm(jm),1).eq.nmanque) then 1288 typm=.false. 1289 nam=0 1290 else 1291 typm=.true. 1292 nam=2 1293 mc1=pheno(c,ll,correm(jm),1) 1294 mc2=pheno(c,ll,correm(jm),2) 1295 end if 1296 if(typp.and.typm)go to 21 1297 1298 ! if ( correp(ip) /= 9999 ) then 1299 ! call log_mess('Sire ['//trim(pere(ip))//'] -->'//trim(get_pheno(c,pheno(c,ll,correp(ip),1)))//& 1300 ! ' '//trim(get_pheno(c,pheno(c,ll,correp(ip),2)))//' marker:'//trim(mark(c,ll)),DEBUG_DEF) 1301 ! end if 1302 ! 1303 ! if ( correm(jm) /= 9999 ) then 1304 ! call log_mess('Dam ['//trim(mere(jm))//'] -->'//trim(get_pheno(c,pheno(c,ll,correm(jm),1)))//' '& 1305 ! //trim(get_pheno(c,pheno(c,ll,correm(jm),2)))//' marker:'//trim(mark(c,ll)),DEBUG_DEF) 1306 ! end if 1307 ! Passage en revue des descendants 1308 nd1=ndm(jm)+1 1309 nd2=ndm(jm+1) 1310 do 22 kd=nd1,nd2 1311 if(corred(kd).eq.9999)go to 22 1312 if(pheno(c,ll,corred(kd),1).eq.nmanque)go to 22 1313 md1=pheno(c,ll,corred(kd),1) 1314 md2=pheno(c,ll,corred(kd),2) 1315 mx=nmanque 1316 1317 ! if ( corred(kd) /= 9999 ) then 1318 ! call log_mess('Progeny ['//trim(animal(kd))//'] -->'//trim(get_pheno(c,pheno(c,ll,corred(kd),1)))//& 1319 ! ' '//trim(get_pheno(c,pheno(c,ll,corred(kd),2)))//' marker:'//trim(mark(c,ll)),DEBUG_DEF) 1320 ! end if 1321 ! 1322 ! Descendant homozygote 1323 ! Pere inconnu 1324 if(md1.eq.md2)then 1325 ! call log_mess('Progeny homo',DEBUG_DEF) 1326 if(nap.ne.2) then 1327 if(nap.eq.0) then 1328 nap=1 1329 m1=md1 1330 ! call log_mess('finding one allele for sire ['//trim(get_pheno(c,m1))//']',DEBUG_DEF) 1331 else 1332 if(md1.ne.m1) then 1333 nap=2 1334 m2=md1 1335 ! call log_mess('finding one allele for sire ['//trim(get_pheno(c,m2))//']',DEBUG_DEF) 1336 end if 1337 end if 1338 end if 1339 ! 1340 ! Mere inconnue 1341 if(nam.ne.2) then 1342 if(nam.eq.0) then 1343 nam=1 1344 mc1=md1 1345 ! call log_mess('finding one allele for dam ['//trim(get_pheno(c,mc1))//']',DEBUG_DEF) 1346 else 1347 if(md1.ne.mc1) then 1348 nam=2 1349 mc2=md1 1350 ! call log_mess('finding one allele for dam ['//trim(get_pheno(c,mc2))//']',DEBUG_DEF) 1351 end if 1352 end if 1353 end if 1354 ! 1355 ! Descendant heterozygote 1356 ! Mere connue et pere inconnu 1357 else 1358 ! call log_mess('Progeny hetero',DEBUG_DEF) 1359 if(nam.eq.2.and.nap.ne.2) then 1360 if((md1.eq.mc1.and.md2.ne.mc2).or. & 1361 (md1.eq.mc2.and.md2.ne.mc1)) mx=md2 1362 if((md2.eq.mc1.and.md1.ne.mc2).or. & 1363 (md2.eq.mc2.and.md1.ne.mc1)) mx=md1 1364 if(mx.ne.nmanque) then 1365 if(nap.eq.0) then 1366 nap=1 1367 m1=mx 1368 ! call log_mess('finding one allele for sire ['//trim(get_pheno(c,m1))//']',DEBUG_DEF) 1369 else 1370 if(mx.ne.m1) then 1371 nap=2 1372 m2=mx 1373 ! call log_mess('finding one allele for sire ['//trim(get_pheno(c,m2))//']',DEBUG_DEF) 1374 end if 1375 end if 1376 end if 1377 end if 1378 ! 1379 ! Pere connu et mere inconnue 1380 if(nap.eq.2.and.nam.ne.2) then 1381 if((md1.eq.m1.and.md2.ne.m2).or. & 1382 (md1.eq.m2.and.md2.ne.m1)) mx=md2 1383 if((md2.eq.m1.and.md1.ne.m2).or. & 1384 (md2.eq.m2.and.md1.ne.m1)) mx=md1 1385 if(mx.ne.nmanque) then 1386 if(nam.eq.0) then 1387 nam=1 1388 mc1=mx 1389 ! call log_mess('finding one allele for dam ['//trim(get_pheno(c,mc1))//']',DEBUG_DEF) 1390 else 1391 if(mx.ne.mc1) then 1392 nam=2 1393 mc2=mx 1394 ! call log_mess('finding one allele for dam ['//trim(get_pheno(c,mc2))//']',DEBUG_DEF) 1395 end if 1396 end if 1397 end if 1398 end if 1399 end if 1400 22 continue 1401 ! 1402 ! Stockage du phenotype reconstitue de la mere 1403 if(typm.or.nam.lt.2) go to 21 1404 if(correm(jm).eq.9999)then 1405 nmes=nmes+1 1406 numero(nmes)= mere(jm) 1407 correm(jm)=nmes 1408 pheno(c,:,correm(jm),1)=nmanque 1409 pheno(c,:,correm(jm),2)=nmanque 1410 end if 1411 1412 1413 1414 jmes=correm(jm) 1415 pheno(c,ll,jmes,1)=mc1 1416 pheno(c,ll,jmes,2)=mc2 1417 ! call log_mess('setting phenotype marker for dam ['//trim(mere(jm))//'] :'//& 1418 ! trim(get_pheno(c,pheno(c,ll,correm(jm),1)))//' '//trim(get_pheno(c,pheno(c,ll,correm(jm),2)))//& 1419 ! ' marker ['//trim(mark(c,ll))//']',DEBUG_DEF) 1420 21 continue 1421 ! 1422 ! Stockage du phenotype reconstitue du pere 1423 if(typp.or.nap.lt.2) go to 20 1424 if(correp(ip).eq.9999) then 1425 nmes=nmes+1 1426 numero(nmes)=pere(ip) 1427 correp(ip)=nmes 1428 pheno(c,:,correp(ip),1)=nmanque 1429 pheno(c,:,correp(ip),2)=nmanque 1430 end if 1431 imes=correp(ip) 1432 pheno(c,ll,imes,1)=m1 1433 pheno(c,ll,imes,2)=m2 1434 ! call log_mess('setting phenotype marker for sire ['//trim(pere(ip))//'] :'//& 1435 ! trim(get_pheno(c,pheno(c,ll,correp(ip),1)))//' '//trim(get_pheno(c,pheno(c,ll,correp(ip),2)))//& 1436 ! ' marker ['//trim(mark(c,ll))//']',DEBUG_DEF) 1437 20 continue 1438 1439 end do 1440 end do ! end nchr 1441 end subroutine recup
set_allele_info_vector
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
set_allele_info_vector
DESCRIPTION
Set basic statistic and information about genotypic information (frequency, number of allele by marker,...)
SOURCE
493 subroutine set_allele_info_vector() 494 integer :: i,j,kd,imes,k,nb,kmes,ios,c 495 logical , dimension(:) , allocatable :: adress 496 logical :: tri_termine 497 character(len=LEN_DEF) :: temp 498 ! dim : nb marker,2 * nmes 499 integer(kind=KIND_PHENO) , dimension (:,:),allocatable :: mm 500 501 call log_mess('START SUBROUTINE : set_allele_info_vector',DEBUG_DEF) 502 503 allocate (mm(maxval(nmk),size(numero)*2), stat=ios) 504 call check_allocate(ios,'mm (m_qtlmap_genotype.f95:set_allele_info_vector)') 505 506 allocate (alleles(nchr,maxval(nmk),size(numero)*2), stat=ios) 507 call check_allocate(ios,'alleles (m_qtlmap_genotype.f95:set_allele_info_vector)') 508 509 allocate (pc_all(nchr,maxval(nmk),size(numero)*2), stat=ios) 510 call check_allocate(ios,'pc_all (m_qtlmap_genotype.f95:set_allele_info_vector)') 511 512 allocate (nall(nchr,maxval(nmk)), stat=ios) 513 call check_allocate(ios,'nall (m_qtlmap_genotype.f95:set_allele_info_vector)') 514 515 do c=1,nchr 516 mm = nmanque 517 call log_mess('['//trim(chromo(c))//'] count number of allele by marker,sort and computing frequencies...',VERBOSE_DEF) 518 do i=1,nmk(c) 519 do j=1,nmes 520 mm(i,j) = pheno(c,i,j,1) 521 mm(i,nmes+j) = pheno(c,i,j,2) 522 end do 523 end do 524 !if none information are found for a 525 do i=1,nmes 526 do j=1,nmk(c) 527 if ( pheno(c,j,i,1) /= nmanque .and. pheno(c,j,i,2) /= nmanque) then 528 exit ! this animal is ok!! 529 end if 530 end do 531 532 if ( j > nmk(c) ) then 533 do kd=1,size(animal) 534 if (animal(kd) == numero(i)) then 535 presentg(c,kd)=.false. 536 exit 537 end if 538 end do 539 end if 540 end do 541 542 543 ! COMPTAGE DES ALLELES PAR MARQUEUR ET CREATION D'UN VECTEUR D'ALLELES par MARQUEUR 544 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(adress,imes,kmes,tri_termine,temp,k,nb) 545 allocate (adress(size(numero)*2), stat=ios) 546 call check_allocate(ios,'adress (m_qtlmap_genotype.f95:set_allele_info_vector)') 547 548 !$OMP DO 549 do i=1,nmk(c) 550 nall(c,i)=0 551 adress= .true. 552 do imes=1,nmes*2 553 if (mm(i,imes)==nmanque) then 554 adress(imes)=.false. 555 else 556 do kmes=1,imes-1 557 ! 558 ! debut modif JME le 11/03/2010 559 560 if (mm(i,imes)==mm(i,kmes)) then 561 adress(imes)=.false. 562 exit 563 end if 564 ! fin modif JME le 11/03/2010 565 ! if (mm(i,imes)==mm(i,kmes)) adress(imes)=.false. 566 enddo 567 endif 568 enddo 569 do imes=1,nmes*2 570 if (adress(imes))then 571 nall(c,i)=nall(c,i)+1 572 alleles(c,i,nall(c,i))=get_pheno(c,mm(i,imes)) 573 endif 574 ! TRIE DES VECTEUR D'ALLELES PAR MARQUEUR 575 tri_termine = .true. 576 do while ( tri_termine ) 577 tri_termine = .false. 578 do k=2,nall(c,i) 579 if ( alleles(c,i,k).LT.alleles(c,i,k-1) ) then 580 tri_termine = .true. 581 temp = alleles(c,i,k-1) 582 alleles(c,i,k-1) = alleles(c,i,k) 583 alleles(c,i,k) = temp 584 endif 585 enddo 586 enddo 587 enddo 588 589 ! CALCUL DE LA FREQUENCE DE CHAQUE ALLELE 590 do k=1,nall(c,i) 591 nb=0 592 do imes=1,nmes*2 593 if (get_pheno(c,mm(i,imes))==alleles(c,i,k)) nb=nb+1 594 enddo 595 pc_all(c,i,k)=dble(nb)*100.d0/(dble(nmes)*2.d0) 596 enddo 597 enddo 598 !$OMP END DO 599 deallocate ( adress ) 600 !$OMP END PARALLEL 601 end do ! end chr 602 603 deallocate ( mm ) 604 call log_mess('END SUBROUTINE : set_allele_info_vector',DEBUG_DEF) 605 end subroutine set_allele_info_vector
set_corresponding_vector
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
set_corresponding_vector
DESCRIPTION
* Build the correspondence vector between genealogy array and the array numero (name of genotyped animal) : corregp,corregm,correp,correm,correr,corred * initialize presentg array
SOURCE
357 subroutine set_corresponding_vector() 358 integer :: igp,jgm,kr,ip,jm,kd,mes,ich 359 integer :: origin_nmes 360 call log_mess('SUBROUTINE : set_corresponding_vector',DEBUG_DEF) 361 362 origin_nmes = nmes 363 ! Creation d'un vecteur de correspondance pour les grands peres (corregp) 364 do igp=1,size(gpere) 365 corregp(igp)=INT_NOT_DEFINED 366 do mes=1,nmes 367 if(gpere(igp) == numero(mes)) then 368 corregp(igp)=mes 369 exit 370 endif 371 end do 372 !Ajout OFI Juin 2009 --> Add unknown genotype to avoid error with array acces 373 if ( corregp(igp) == INT_NOT_DEFINED ) then 374 call log_mess('The grand sire '// trim(gpere(igp))// ' has no genotype',WARNING_DEF) 375 nmes = nmes + 1 376 numero(nmes) = gpere(igp) 377 corregp(igp)=nmes 378 pheno(:,:,corregp(igp),:)=nmanque 379 end if 380 381 end do 382 ! if (t_imp) write(16,4000)gpere(igp) 383 !4000 format(1x,'The grand sire ',i6,' has no genotype') 384 ! 42 continue 385 386 ! Creation d'un vecteur de correspondance pour les grands meres (corregm) 387 do jgm=1,size(gmere) 388 corregm(jgm)=INT_NOT_DEFINED 389 do mes=1,nmes 390 if(gmere(jgm) == numero(mes)) then 391 corregm(jgm)=mes 392 exit 393 endif 394 end do 395 396 !Ajout OFI Juin 2009 --> Add unknown genotype to avoid error with array acces 397 if ( corregm(jgm) == INT_NOT_DEFINED ) then 398 call log_mess('The grand dam '// trim(gmere(jgm))// ' has no genotype',WARNING_DEF) 399 nmes = nmes + 1 400 numero(nmes) = gmere(jgm) 401 corregm(jgm)=nmes 402 pheno(:,:,corregm(jgm),:)=nmanque 403 end if 404 end do 405 406 ! Creation d'un vecteur de correspondance pour les peres (correp) 407 do ip=1,size(pere) 408 correp(ip)=INT_NOT_DEFINED 409 do mes=1,nmes 410 if(pere(ip) == numero(mes)) then 411 correp(ip)=mes 412 exit 413 endif 414 end do 415 416 !Ajout OFI Juin 2009 --> Add unknown genotype to avoid error with array acces 417 if ( correp(ip) == INT_NOT_DEFINED ) then 418 call log_mess('The sire '// trim(pere(ip))// ' has no genotype',WARNING_DEF) 419 nmes = nmes + 1 420 numero(nmes) = pere(ip) 421 correp(ip)=nmes 422 pheno(:,:,correp(ip),:)=nmanque 423 end if 424 425 end do 426 ! if (t_imp) write(16,4002)pere(ip) 427 !4002 format(1x,'The sire ',i6,' has no genotype') 428 ! 44 continue 429 430 !Creation d'un vecteur de correspondance pour les meres (correm) 431 do jm=1,size(mere) 432 correm(jm)=INT_NOT_DEFINED 433 do mes=1,nmes 434 if(mere(jm) == numero(mes)) then 435 correm(jm)=mes 436 exit 437 endif 438 end do 439 440 !Ajout OFI Juin 2009 --> Add unknown genotype to avoid error with array acces 441 if ( correm(jm) == INT_NOT_DEFINED ) then 442 call log_mess('The dam '// trim(mere(jm))// ' has no genotype',WARNING_DEF) 443 nmes = nmes + 1 444 numero(nmes) = mere(jm) 445 correm(jm)=nmes 446 pheno(:,:,correm(jm),:)=nmanque 447 end if 448 ! if (t_imp) write(16,4003)mere(jm) 449 ! 4003 format(1x,'The dam ',i6,' has no genotype') 450 ! 45 continue 451 end do 452 453 ! Creation d'un vecteur de correspondance pour les parents (correr) 454 do kr=1,size(repro) 455 correr(kr)=INT_NOT_DEFINED 456 do mes=1,nmes 457 if(repro(kr) == numero(mes)) then 458 correr(kr)=mes 459 exit 460 endif 461 end do 462 463 end do 464 465 ! Creation d'un vecteur de correspondance pour les descendants (corred) 466 do kd=1,size(animal) 467 presentg(:,kd)=.false. 468 corred(kd)=INT_NOT_DEFINED 469 do mes=1,nmes 470 if(animal(kd) == numero(mes)) then 471 corred(kd)=mes 472 do ich=1,nchr 473 presentg(ich,kd)=.not. (allocated(pheno(ich,:,mes,:)==nmanque)) 474 end do 475 exit 476 endif 477 end do 478 end do 479 480 nmes = origin_nmes 481 482 call log_mess('END SUBROUTINE : set_corresponding_vector',DEBUG_DEF) 483 end subroutine set_corresponding_vector
set_estfem
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
set_estfem
DESCRIPTION
initialize the array estfem.
INPUTS
ndmin : minimum number of progenies to build full sib family
SOURCE
1478 subroutine set_estfem(ndmin) 1479 integer , intent(in) :: ndmin 1480 integer :: ip,jm 1481 if ( allocated (estfem)) deallocate ( estfem ) 1482 allocate (estfem(size(mere)) ) 1483 estfem=.false. 1484 1485 do ip=1,size(pere) 1486 do jm=nmp(ip)+1,nmp(ip+1) 1487 if((ndm(jm+1)-ndm(jm)).ge.ndmin) estfem(repfem(jm))=.true. 1488 1489 end do 1490 end do 1491 end subroutine set_estfem
sim_typ
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
sim_typ
DESCRIPTION
NOTES
Simulation des typages sur trois generations , reperage des alleles au qtl des descendants.
SOURCE
827 subroutine sim_typ(ndmin) 828 ! 829 integer , intent(in) :: ndmin 830 ! Divers 831 logical :: typ2 832 integer :: c,p 833 double precision liminf,limsup 834 integer gpare,k,pare,i,j,iq,jgm,ngm1,ngm2,igp 835 integer l1,l2,nr1,nr2,ir,ip,nm1,nm2,jm,id,ic,mes,nd1,nd2 836 real (kind=dp) :: r 837 real :: x_rand 838 839 real, external :: ranf 840 ! 841 ! 842 !******************************************************************* 843 ! Lecture parametres de simulation du QTL 844 !******************************************************************* 845 ! 846 call log_mess('simulation of genotype...',VERBOSE_DEF) 847 pheno=nmanque 848 ! 849 !*********************************************************************** 850 ! Simulation des typages des grands-parents 851 ! 852 !*********************************************************************** 853 ! Alleles des Grands peres 854 ! 855 if (.not. allocated (pc_all) ) then 856 call stop_application('Devel error: pc_all array is not allocated. Can not simulate type.') 857 end if 858 859 if (.not. allocated (alleles) ) then 860 call stop_application('Devel error: all array is not allocated. Can not simulate type.') 861 end if 862 863 864 mes=0 865 do igp=1,ngp 866 mes=mes+1 867 numero(mes)=gpere(igp) 868 corregp(igp)=mes 869 ! Alleles aux marqueurs 870 do c=1,nchr 871 do j=1,2 872 do i=1,nmk(c) 873 x_rand=ranf() 874 limsup=0.d0 875 do k=1,nall(c,i) 876 liminf=limsup 877 limsup=liminf+pc_all(c,i,k)/1.d2 878 if(dble(x_rand).gt.liminf.and.dble(x_rand).le.limsup) then 879 pheno(c,i,mes,j)=set_pheno(c,alleles(c,i,k)) 880 end if 881 end do 882 end do 883 end do 884 end do 885 end do 886 ! 887 ! Alleles des Grands meres 888 do jgm=1,ngm 889 mes=mes+1 890 numero(mes)=gmere(jgm) 891 corregm(jgm)=mes 892 ! 893 ! Alleles aux marqueurs 894 do c=1,nchr 895 do j=1,2 896 do i=1,nmk(c) 897 x_rand=ranf() 898 limsup=0.d0 899 do k=1,nall(c,i) 900 liminf=limsup 901 limsup=liminf+pc_all(c,i,k)/1.d2 902 if(dble(x_rand).gt.liminf.and.dble(x_rand).le.limsup) then 903 pheno(c,i,mes,j)=set_pheno(c,alleles(c,i,k)) 904 end if 905 end do 906 end do 907 end do 908 end do 909 end do !! jgm 910 ! 911 ! 912 !********************************************************************** 913 ! Simulation de la transmission des phases aux parents 914 ! Les parents qui sont aussi gparents doivent etre enregistres les premiers 915 ! 916 !********************************************************************** 917 ! 918 do igp=1,ngp 919 ngm1=ngmgp(igp)+1 920 ngm2=ngmgp(igp+1) 921 do jgm=ngm1,ngm2 922 nr1=nrgm(jgm)+1 923 nr2=nrgm(jgm+1) 924 do ir=nr1,nr2 925 mes=mes+1 926 numero(mes)=repro(ir) 927 correr(ir)=mes 928 ! 929 ! Deux chromosomes a la suite : p=1 gd pere ; p=2 gd mere 930 do p=1,2 931 if (p.eq.1) gpare=corregp(igp) 932 if (p.eq.2) gpare=corregm(jgm) 933 l1=1 934 k=0 935 iq=1 936 ! Premier marqueur 937 ! Choix de la phase heritee au 1er marqueur : k 938 x_rand=ranf() 939 if (dble(x_rand).le.0.5d0) then 940 pheno(:,l1,mes,p)=pheno(:,l1,gpare,1) 941 k=1 942 else 943 pheno(:,l1,mes,p)=pheno(:,l1,gpare,2) 944 k=2 945 end if 946 ! 947 ! Marqueurs suivants 948 do c=1,nchr 949 do l2=2,nmk(c) 950 if (p.eq.1) r=(rm(c,l1,l2)) 951 if (p.eq.2) r=(rf(c,l1,l2)) 952 x_rand=ranf() 953 typ2=.true. 954 ! 955 ! Il y a recombinaison 956 if (dble(x_rand).lt.r) then 957 ! Phase 1 heritee -> echange 958 if (k .eq. 1 .and. typ2) then 959 pheno(c,l2,mes,p)=pheno(c,l2,gpare,3-k) 960 ! 961 k=2 962 typ2=.false. 963 end if 964 ! 965 ! Phase 2 heritee -> echange si pas deja fait de phase 1 a 2 966 if (k .eq. 2 .and. typ2) then 967 pheno(c,l2,mes,p)=pheno(c,l2,gpare,3-k) 968 k=1 969 ! 970 end if 971 end if 972 973 ! 974 ! Pas de recombinaison 975 if (dble(x_rand).ge.r)pheno(c,l2,mes,p)=pheno(c,l2,gpare,k) 976 ! 977 l1=l1+1 978 end do !! fin chr 979 end do 980 end do !!fin deux gparents 981 ! 982 end do !! repro 983 end do !! gmere 984 end do !! gpere 985 ! 986 ! 987 !********************************************************************** 988 ! Simulation de la transmission des phases aux descendants 989 ! 990 !********************************************************************** 991 ! 992 do ip=1,np 993 correp(ip)=correr(reppere(ip)) 994 nm1=nmp(ip)+1 995 nm2=nmp(ip+1) 996 do jm=nm1,nm2 997 correm(jm)=correr(repmere(jm)) 998 nd1=ndm(jm)+1 999 nd2=ndm(jm+1) 1000 do id=nd1,nd2 1001 mes=mes+1 1002 numero(mes)=animal(id) 1003 corred(id)=mes 1004 do p=1,2 1005 if (p.eq.1) pare=correp(ip) 1006 if (p.eq.2) pare=correm(jm) 1007 l1=1 1008 k=0 1009 iq=1 1010 ! Premier marqueur 1011 x_rand=ranf() 1012 if (dble(x_rand).le.0.5d0) then 1013 pheno(:,l1,mes,p)=pheno(:,l1,pare,1) 1014 k=1 1015 else 1016 pheno(:,l1,mes,p)=pheno(:,l1,pare,2) 1017 k=2 1018 end if 1019 ! Marqueurs suivants 1020 do c=1,nchr 1021 do l2=2,nmk(c) 1022 if (p.eq.1) r=(rm(c,l1,l2)) 1023 if (p.eq.2) r=(rf(c,l1,l2)) 1024 x_rand=ranf() 1025 typ2=.true. 1026 ! 1027 if (dble(x_rand).lt.r) then 1028 if (k .eq. 1 .and. typ2) then 1029 pheno(c,l2,mes,p)=pheno(c,l2,pare,2) 1030 k=2 1031 typ2=.false. 1032 goto 10 1033 end if 1034 ! 1035 if (k .eq. 2 .and. typ2) then 1036 pheno(c,l2,mes,p)=pheno(c,l2,pare,1) 1037 k=1 1038 goto 10 1039 end if 1040 end if 1041 ! 1042 if (dble(x_rand).ge.r) then 1043 if (k.eq.1) pheno(c,l2,mes,p)=pheno(c,l2,pare,1) 1044 if (k.eq.2) pheno(c,l2,mes,p)=pheno(c,l2,pare,2) 1045 goto 10 1046 end if 1047 ! 1048 10 continue 1049 l1=l1+1 1050 end do !! fin chr 1051 end do 1052 end do !!fin deux gparents 1053 1054 end do !! animal 1055 end do !! mere 1056 end do !! pere 1057 ! 1058 !********************************************************************** 1059 ! Ecriture des genotypes dans le fichier typ_sim 1060 !********************************************************************** 1061 ! 1062 nmes=mes 1063 1064 call set_estfem(ndmin) 1065 end subroutine sim_typ
write_typ
[ Top ] [ m_qtlmap_genotype ] [ Subroutines ]
NAME
write_typ
DESCRIPTION
INPUTS
file_name : name of the ouput file
SOURCE
1452 subroutine write_typ(file_name) 1453 character(len=*),intent(in) :: file_name 1454 integer :: mes,i,j,c,lk 1455 1456 1457 open(1111,file=file_name) 1458 1459 1460 write (1111,*) ((trim(mark(c,lk))//' ',lk=1,nmk(c)),c=1,nchr) 1461 1462 do mes=1,nmes 1463 write (1111,*) trim(numero(mes))//' ',(((trim(get_pheno(c,pheno(c,i,mes,j)))//' ',j=1,2),i=1,nmk(c)),c=1,nchr) 1464 end do 1465 1466 close (1111) 1467 end subroutine write_typ