m_qtlmap_simulation
NAME
m_qtlmap_simulation -- Simulation routines.
SYNOPSIS
DESCRIPTION
NOTES
ios_simul_file
[ Top ] [ m_qtlmap_simulation ] [ Constants ]
NAME
ios_simul_file
DESCRIPTION
Unit fortran record to read the simulation file
current_line
[ Top ] [ m_qtlmap_simulation ] [ Variables ]
NAME
ios_simul_file
DESCRIPTION
number of the current line while the parse
init_genotype_simul
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
init_genotype_simul
DESCRIPTION
initialize the array of marker coming from m_qtlmap_data
NOTES
SOURCE
436 subroutine init_genotype_simul 437 integer :: ios 438 439 nmes = ngp+ngm+np+nm+nd 440 allocate ( numero(nmes) , stat=ios) 441 call check_allocate(ios,'numero') 442 443 allocate ( pheno(size(nmk),maxval(nmk),nmes,2) , stat=ios) 444 call check_allocate(ios,'pheno') 445 446 allocate (corregp(size(gpere)), stat=ios) 447 call check_allocate(ios,'corregp') 448 449 allocate (corregm(size(gmere)), stat=ios) 450 call check_allocate(ios,'corregm') 451 452 allocate (correr(size(repro)), stat=ios) 453 call check_allocate(ios,'correr') 454 455 allocate (correp(size(pere)), stat=ios) 456 call check_allocate(ios,'correp') 457 458 allocate (correm(size(mere)), stat=ios) 459 call check_allocate(ios,'correm') 460 461 allocate (corred(size(animal)), stat=ios) 462 call check_allocate(ios,'corred') 463 464 allocate (presentg(nchr,size(animal)), stat=ios) 465 call check_allocate(ios,'presentg') 466 467 end subroutine init_genotype_simul
init_permutation
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
init_permutation
DESCRIPTION
initialise a permutation (none simulation user file is readed)
INPUTS
is_transcriptome : data are transcriptomic ?
NOTES
SOURCE
208 subroutine init_permutation(is_transcriptome) 209 logical , intent(in) :: is_transcriptome 210 211 simulMap=.false. 212 simulGenea=.false. 213 call read_genealogy 214 call read_model(in_param_ef) 215 call set_simulation_traits() !! initilialize filter_car private integer array 216 call fixe_structure_model(filter_car) 217 call read_traits(is_transcriptome,filter_car) 218 219 end subroutine init_permutation
init_simul_genealogy
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
init_simul_genealogy
DESCRIPTION
initialize the array of marker coming from m_qtlmap_data
NOTES
SOURCE
392 subroutine init_simul_genealogy(ngp,ngm,np,nm,nd,nr) 393 integer , intent(in) :: ngp,ngm,np,nm,nd,nr 394 integer :: stat 395 396 allocate (gmere(ngm), STAT = stat) 397 call check_allocate(stat,'gmere') 398 allocate (gpere(ngp), STAT = stat) 399 call check_allocate(stat,'gpere') 400 allocate (repro(nr), STAT = stat) 401 call check_allocate(stat,'repro') 402 allocate (animal(nd), STAT = stat) 403 call check_allocate(stat,'animal') 404 allocate (pere(np), STAT = stat) 405 call check_allocate(stat,'pere') 406 allocate (mere(nm), STAT = stat) 407 call check_allocate(stat,'mere') 408 allocate (femelle(nm), STAT = stat) 409 call check_allocate(stat,'femelle') 410 nfem = 0 411 412 allocate (ngmgp(ngp+1), STAT = stat) 413 call check_allocate(stat,'ngmgp') 414 allocate (nrgm(ngm+1), STAT = stat) 415 call check_allocate(stat,'nrgm') 416 allocate (nmp(np+1), STAT = stat) 417 call check_allocate(stat,'nmp') 418 allocate (ndm(nm+1), STAT = stat) 419 call check_allocate(stat,'ndm') 420 allocate (reppere(np), STAT = stat) 421 call check_allocate(stat,'reppere') 422 allocate (repmere(nm), STAT = stat) 423 call check_allocate(stat,'repmere') 424 allocate (repfem(nm), STAT = stat) 425 call check_allocate(stat,'repfem') 426 end subroutine init_simul_genealogy
init_simul_marker
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
init_simul_marker
DESCRIPTION
initialize the array of marker coming from m_qtlmap_data
NOTES
SOURCE
292 subroutine init_simul_marker() 293 integer :: alloc_stat,nballelepossible,i 294 295 allocate (mark(size(nmk,1),maxval(nmk)), stat = alloc_stat) 296 call check_allocate(alloc_stat,'mark') 297 298 do i=1,nmk(1) 299 mark(1,i)='gen_mark_'//str(i) 300 end do 301 302 allocate (posi(size(nmk,1),maxval(nmk)), stat = alloc_stat) 303 call check_allocate(alloc_stat,'posi') 304 305 allocate (posif(size(nmk,1),maxval(nmk)), stat = alloc_stat) 306 call check_allocate(alloc_stat,'posif') 307 308 allocate (posim(size(nmk,1),maxval(nmk)), stat = alloc_stat) 309 call check_allocate(alloc_stat,'posim') 310 311 allocate (rm(size(nmk,1),maxval(nmk),maxval(nmk)), stat = alloc_stat) 312 call check_allocate(alloc_stat,'rm') 313 314 allocate (rf(size(nmk,1),maxval(nmk),maxval(nmk)), stat = alloc_stat) 315 call check_allocate(alloc_stat,'rf') 316 317 allocate ( nall(size(nmk,1),maxval(nmk)) , stat = alloc_stat) 318 call check_allocate(alloc_stat,'nall') 319 320 nballelepossible = int(VAL_MAX_INDEX_PHENO) - int(VAL_MIN_INDEX_PHENO) + 1 321 322 allocate ( alleles(size(nmk,1),maxval(nmk),nballelepossible) , stat = alloc_stat) 323 call check_allocate(alloc_stat,'alleles') 324 325 allocate ( pc_all(size(nmk,1),maxval(nmk),nballelepossible) , stat = alloc_stat) 326 call check_allocate(alloc_stat,'pc_all') 327 328 end subroutine init_simul_marker
init_traits_simul
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
init_traits_simul
DESCRIPTION
initialize the array of marker coming from m_qtlmap_data
NOTES
SOURCE
477 subroutine init_traits_simul 478 integer :: alloc_stat,ic 479 480 allocate (carac(ncar) , stat=alloc_stat) 481 call check_allocate(alloc_stat,'carac') 482 483 allocate (corperf(size(animal)), stat = alloc_stat) 484 call check_allocate(alloc_stat,'corperf') 485 486 allocate (y(ncar,size(animal)), stat = alloc_stat) 487 call check_allocate(alloc_stat,'y') 488 489 allocate (presentc(ncar,size(animal)), stat = alloc_stat) 490 call check_allocate(alloc_stat,'presentc') 491 492 allocate(natureY(ncar), stat = alloc_stat) 493 call check_allocate(alloc_stat,'natureY') 494 495 allocate (cd(ncar,size(animal)), stat = alloc_stat) 496 call check_allocate(alloc_stat,'cd') 497 cd=1.d0 498 presentc = .true. 499 500 allocate (ndelta(ncar,size(animal)), stat = alloc_stat) 501 call check_allocate(alloc_stat,'ndelta') 502 503 allocate (covar(nd,ncov), stat = alloc_stat) 504 call check_allocate(alloc_stat,'covar') 505 506 allocate (namefix(nfix), stat = alloc_stat) 507 call check_allocate(alloc_stat,'namefix') 508 509 allocate (namecov(ncov), stat = alloc_stat) 510 call check_allocate(alloc_stat,'namecov') 511 512 allocate (listelev(nfix,nd), stat = alloc_stat) 513 call check_allocate(alloc_stat,'listelev') 514 515 allocate (modele(ncar,3+(2*nfix)+ncov), stat = alloc_stat) 516 call check_allocate(alloc_stat,'modele') 517 modele=0 518 519 allocate (listModelTrait(ncar)) 520 do ic=1,ncar 521 allocate ( listModelTrait(ic)%indexFixedEffect(0)) 522 listModelTrait(ic)%nbfe=0 523 allocate ( listModelTrait(ic)%indexCovariate(0)) 524 listModelTrait(ic)%nbco=0 525 allocate (listModelTrait(ic)%nbint(MAX_QTL)) 526 listModelTrait(ic)%nbint=0 527 end do 528 529 allocate (nlev(nfix), stat = alloc_stat) 530 call check_allocate(alloc_stat,'nlev') 531 532 allocate (nivx(nd,nfix), stat = alloc_stat) 533 call check_allocate(alloc_stat,'nivx') 534 535 end subroutine init_traits_simul
read_attributes_discrete_trait
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
read_attributes_discrete_trait
DESCRIPTION
NOTES
SOURCE
1018 subroutine read_attributes_discrete_trait(nc,ic,line_read) 1019 integer , intent(in) :: nc,ic 1020 character(len=LEN_DEF) ,intent(inout) :: line_read 1021 character(len=LEN_DEF) :: word 1022 integer :: nbmod,i 1023 logical :: is_ok 1024 real :: summod 1025 real (kind=dp),dimension(:,:) , allocatable :: freqmodsimultrait_t 1026 1027 if ( .not. allocated(nbmodsimultrait)) then 1028 allocate(nbmodsimultrait(nc)) 1029 nbmodsimultrait=0 1030 end if 1031 !number of modality 1032 word = trim(next_word(line_read,is_ok)) 1033 1034 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,'number of modality ['//trim(str(ic))//'] is not defined.') 1035 nbmod=get_int(word,is_ok) 1036 1037 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,& 1038 'modality of traits ['//trim(str(ic))//'] have to be an integer.') 1039 nbmodsimultrait(ic)=nbmod 1040 1041 if ( .not. allocated(freqmodsimultrait)) then 1042 allocate(freqmodsimultrait(nc,nbmod)) 1043 freqmodsimultrait=0.d0 1044 else 1045 if (size(freqmodsimultrait,2)<nbmod) then 1046 allocate(freqmodsimultrait_t(size(freqmodsimultrait,1),size(freqmodsimultrait,2))) 1047 freqmodsimultrait_t=freqmodsimultrait 1048 deallocate(freqmodsimultrait) 1049 allocate(freqmodsimultrait(size(freqmodsimultrait_t,1),nbmod)) 1050 freqmodsimultrait=0.0 1051 freqmodsimultrait(:,1:size(freqmodsimultrait_t,2))=freqmodsimultrait_t 1052 deallocate(freqmodsimultrait_t) 1053 end if 1054 end if 1055 1056 summod=0.d0 1057 do i=1,nbmod 1058 word = trim(next_word(line_read,is_ok)) 1059 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,'freq [trait:'//& 1060 trim(str(ic))//' mod:"//trim(str(i))//"] is not defined.') 1061 freqmodsimultrait(ic,i)=get_real(word,is_ok) 1062 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,& 1063 'frequency of modality :'//trim(str(i))//'for trait ['//trim(str(ic))//'] have to be a real.') 1064 if ( (freqmodsimultrait(ic,i) > 1.d0) .or. (freqmodsimultrait(ic,i) < 0.d0) ) then 1065 call stop_on_error (1,in_parsim,current_line,& 1066 'frequency of modality :'//trim(str(i))//'for trait ['//trim(str(ic))//'] have to be 0<=mod<=1.') 1067 end if 1068 summod = summod + freqmodsimultrait(ic,i) 1069 end do 1070 1071 if ( summod /= 1.d0 ) then 1072 call stop_on_error (1,in_parsim,current_line,& 1073 'sum of frequency have to be equal to 1. Trait ['//trim(str(ic))//' : sum:'//trim(str(summod))//']') 1074 end if 1075 1076 end subroutine read_attributes_discrete_trait
read_correlation_matrix
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
read_correlation_matrix
DESCRIPTION
NOTES
SOURCE
731 subroutine read_correlation_matrix(nc) 732 integer ,intent(in) :: nc 733 character(len=LEN_DEF) :: word 734 character(len=LEN_LINE) :: line_read 735 character(len=1000) :: buf_line 736 logical :: is_ok 737 integer :: j,ic,ios,jc 738 739 740 if ( nc <= 1 ) return 741 !correlation matrix definition 742 line_read='' 743 744 do while ( trim(line_read)=='') 745 current_line = current_line + 1 746 call GET(ios_simul_file,line_read,IOSTAT=ios) 747 if ( (ios /=0) ) then 748 !read(ios_simul_file,*,iostat=ios) (rho(ic,jc),jc=1,ic-1) 749 call stop_on_error(ios,in_parsim, current_line,"The correlation matrix is not defined") 750 end if 751 end do 752 753 !FORMAT [ [corr(1)] [ corr(2,1) corr(2,2)] ... ] 754 755 !!keyword correlation 756 word = trim(next_word(line_read,is_ok)) 757 if (word /= 'correlation') then 758 call stop_application("'correlation' keyword expected! **"//trim(word)//"**") 759 end if 760 !! first "[" 761 buf_line = trim(line_read) 762 763 j=index(buf_line,"[") 764 if ( j == 0 ) then 765 call stop_on_error(ios,in_parsim, current_line,"Expecting '[' :**"//trim(line_read)//"**") 766 end if 767 buf_line=buf_line(j+1:) 768 769 do ic=2,nc 770 call parse_real_array(buf_line,ic-1,RhoP(ic,:),is_ok) 771 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,'Correlation matrix definition for ['//& 772 trim(str(ic))//'] unknown.') 773 do jc=1,ic-1 774 775 !Add test : a correlation is between -1 and 1 776 if ( RhoP(ic,jc) > 1 .or. RhoP(ic,jc) < -1 ) then 777 call stop_application("Correlation have to be defined betwwen -1 and 1 bad value :["//trim(str(RhoP(ic,jc)))//"]") 778 end if 779 call log_mess ('Correlation matrix definition ['//trim(str(ic))//']:'//trim(str(RhoP(ic,jc))),VERBOSE_DEF) 780 end do 781 end do 782 783 j=index(buf_line,"]") 784 if ( j == 0 ) then 785 call stop_on_error(1,in_parsim, current_line,"Expecting ']' to close the matrix definition") 786 end if 787 buf_line=buf_line(:j-1) 788 if ( trim(buf_line) /= '' ) then 789 call stop_on_error(1,in_parsim, current_line,"Too many information in correlation matrix :"//trim(buf_line)//']') 790 end if 791 792 !fill the matrix 793 do ic=2,nc 794 do jc=1,ic-1 795 RhoP(jc,ic)=RhoP(ic,jc) 796 end do 797 end do 798 799 end subroutine read_correlation_matrix
read_qtl_effect_on_trait
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
read_qtl_effect_on_trait
DESCRIPTION
NOTES
SOURCE
809 subroutine read_qtl_effect_on_trait(nc) 810 integer, intent(in) :: nc 811 character(len=LEN_DEF) :: word 812 character(len=LEN_LINE) :: line_read 813 logical :: is_ok 814 integer :: ic,iq,ios 815 816 allocate (ue(nc,nqtl)) 817 ue=0.d0 818 if (nqtl > 0) then 819 !qtl effect 820 821 current_line = current_line + 1 822 line_read='' 823 do while ( trim(line_read) == '' ) 824 call GET(ios_simul_file,line_read,IOSTAT=ios) 825 if ( (ios /= 0) ) then 826 call stop_on_error(1,in_parsim, current_line,"QTL effects on traits") 827 end if 828 end do 829 word = trim(next_word(line_read,is_ok)) 830 if (word /= 'qtleffect') then 831 call stop_application("'qtleffect' keyword expected! **"//trim(word)//"**") 832 end if 833 834 835 do ic=1,nc 836 do iq=1,nqtl 837 word = next_word(line_read,is_ok) 838 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,& 839 'QTL ['//trim(str(iq))//'] effects on traits ['//trim(str(ic))//'] unknown.') 840 ue(ic,iq) = get_real(word,is_ok) 841 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,& 842 'QTL ['//trim(str(iq))//'] effects on traits ['//trim(str(ic))//'] have to be a real.') 843 call log_mess ( 'QTL ['//trim(str(iq))//'] effects on traits ['//trim(str(ic))//']:'& 844 //trim(str(ue(ic,iq))),VERBOSE_DEF) 845 !read(ios_simul_file,*,iostat=ios) ((ue(ic,iq),iq=1,nqtl),ic=1,ncar) 846 !call stop_on_error(ios,in_parsim, current_line,"QTL effects on traits") 847 end do 848 end do 849 850 else 851 call log_mess('NQTL=0 --> No qtl effect to define',VERBOSE_DEF) 852 end if 853 end subroutine read_qtl_effect_on_trait
read_simulation_file
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
read_simulation_file
DESCRIPTION
read simulation user file and initialize internal array
INPUTS
parsim_file : the simulation user file is_transcriptome : data are transcriptomic ? ndmin : minimum of progenies for consider full sib family
OUTPUTS
croisement : OUTBRED_KEYWORD, F2_KEYWORD, BC_KEYWORD nqtlsimul : number of qtl to simulate simulQtl : at least one qtl is simulated simMap : simulation map case ?
NOTES
SOURCE
91 subroutine read_simulation_file(is_transcriptome,simulQtl,& 92 nqtlsimul,simMap,ndmin,croisement) 93 logical , intent(in) :: is_transcriptome 94 logical,intent(inout) :: simulQtl 95 logical , intent(inout) :: simMap 96 integer , intent(in) :: ndmin 97 integer , intent(out) :: nqtlsimul 98 character(len=LEN_BUFFER_WORD),intent(out):: croisement 99 character(len=LEN_BUFFER_WORD) :: token 100 integer :: ios 101 logical :: label_traits_present,label_simultraits_present 102 103 call log_mess('SUBROUTINE : read_simulation_file',DEBUG_DEF) 104 call log_mess('reading simulation file...',INFO_DEF) 105 !in_parsim = parsim_file 106 simulQtl = .false. 107 simMap = .false. 108 label_traits_present = .false. 109 label_simultraits_present = .false. 110 nqtlsimul=0 111 112 call file_exist(in_parsim) 113 114 115 open(unit=ios_simul_file,file=in_parsim,action="read",form="formatted") 116 ios = 0 117 do while ( ios == 0 ) 118 current_line = current_line + 1 119 read(ios_simul_file,*,iostat=ios) token 120 if ( ios == 0 ) then 121 select case (trim(token)) 122 case (LABEL_MARKERS) 123 call log_mess('setting simulation of markers...',INFO_DEF) 124 call set_simulation_marker(simMap) 125 case (LABEL_GENEALOGY) 126 call log_mess('setting simulation of genealogy...',INFO_DEF) 127 call set_simulation_genealogy(croisement) 128 call init_simul_genealogy(ngp,ngm,np,nm,nd,nr) 129 call sim_genea(inmp,indm,croisement) 130 case (LABEL_QTL) 131 call log_mess('setting simulation of qtl...',INFO_DEF) 132 call set_simulation_qtl(simulQtl) 133 nqtlsimul=nqtl 134 case (LABEL_TRAITS) ! trait are described in the model and unknown data are taken 135 if ( .not. traitsFileDefined ) then 136 call stop_application("None traits file are defined but a LABEL "//trim(LABEL_TRAITS)& 137 //" is detected. Please use a LABEL "//trim(LABEL_SIMULTRAITS)//" to not use real traits.") 138 end if 139 if ( label_simultraits_present ) then 140 call stop_application(trim(LABEL_SIMULTRAITS)//" have been defined. You cannot define "//& 141 trim(LABEL_TRAITS)//"." ) 142 end if 143 call log_mess('setting simulation of traits...',INFO_DEF) 144 if (simulGenea) then 145 call stop_application("You cannot use real data traits with a simulated genealogy.") 146 end if 147 call log_mess('setting simul with real traits...',INFO_DEF) 148 croisement=F2_KEYWORD 149 150 if ( .not. geneaFileDefined) then 151 call stop_application("Genealogy file is not defined and none genealogy are simulated.") 152 end if 153 154 call read_genealogy 155 call read_model(in_param_ef) 156 call set_simulation_traits() !! initilialize filter_car private integer array 157 call fixe_structure_model(filter_car) 158 call read_traits(is_transcriptome,filter_car) 159 label_traits_present=.true. 160 161 case (LABEL_SIMULTRAITS) ! traits are not real (not defined in the model file) 162 if ( label_traits_present ) then 163 call stop_application(LABEL_TRAITS//" have been defined. You cannot define "//LABEL_SIMULTRAITS//"." ) 164 end if 165 if (.not. simulGenea) then 166 if ( .not. geneaFileDefined) then 167 call stop_application("Genealogy file is not defined and none genealogy are simulated.") 168 end if 169 croisement=F2_KEYWORD 170 call read_genealogy 171 end if 172 call log_mess('setting simul traits...',INFO_DEF) 173 call set_simulation_simultraits() 174 label_simultraits_present=.true. 175 case default 176 call stop_on_error(-1,in_parsim, current_line,"label unknown:"//trim(token)) 177 end select 178 end if 179 end do 180 181 ! checking at least simulation traits 182 if ( .not. simulTraits ) then 183 call stop_application('Check simulation file. '//trim(LABEL_SIMULTRAITS)//' Label is not defined.') 184 end if 185 186 if ( .not. simulMap .and. .not. genotypeFileDefined) then 187 call stop_application("Genotype file is not defined and none genotype are simulated.") 188 end if 189 if ( .not. simulMap .and. .not. mapFileDefined) then 190 call stop_application("Map file is not defined and none map are simulated.") 191 end if 192 ! 193 close(ios_simul_file) 194 195 call log_mess('END SUBROUTINE : read_simulation_file',DEBUG_DEF) 196 end subroutine read_simulation_file
set_simulation_genealogy
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
set_simulation_genealogy
DESCRIPTION
read the GENEALOGY part
OUTPUTS
croisement : OUTBRED_KEYWORD, F2_KEYWORD, BC_KEYWORD
NOTES
SOURCE
341 subroutine set_simulation_genealogy(croisement) 342 character(len=LEN_DEF) :: message 343 integer :: ios 344 character(len=LEN_BUFFER_WORD),intent(out):: croisement 345 346 simulGenea = .true. 347 current_line = current_line + 1 348 read(ios_simul_file,*,iostat=ios) croisement 349 message = "None genealogy description is finding." 350 call stop_on_error(ios,in_parsim, current_line,trim(message)) 351 352 if ( croisement /= F2_KEYWORD .and. & 353 croisement /= BC_KEYWORD .and. & 354 croisement /= OUTBRED_KEYWORD ) then 355 message = in_parsim // "-Problem detected at the genealogy description:"//& 356 ". SIMULATION FOR NEW PEDIGREE IS ONLY AVAILABLE FOR F2 TYPE CROSSES : value ["//trim(croisement)//"]" 357 call stop_application(trim(message)) 358 endif 359 360 361 current_line = current_line + 1 362 read(ios_simul_file,*,iostat=ios) nbpere, inmp, indm 363 call stop_on_error(ios,in_parsim, current_line,"number of sires , number of dams per sire" // & 364 ", number of progeny per dam") 365 366 if (croisement /= OUTBRED_KEYWORD ) then 367 np = nbpere 368 ngp= np 369 ngm= np 370 nm=np*inmp 371 nd=nm*indm 372 nr=np+nm 373 else 374 np = nbpere 375 ngp = nbpere*inmp + nbpere 376 ngm = nbpere*inmp + nbpere 377 nm = np*inmp 378 nd= nm*indm 379 nr=np+nm 380 end if 381 382 end subroutine set_simulation_genealogy
set_simulation_marker
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
set_simulation_marker
DESCRIPTION
read the MARKER part
OUTPUTS
simMap : simulation map ?
NOTES
SOURCE
231 subroutine set_simulation_marker(simMap) 232 logical , intent(inout) :: simMap 233 character(len=LEN_DEF) :: message 234 integer :: ios,nb_marker_used,alloc_stat,iq,ic 235 real :: temp 236 simMap = .true. 237 simulMap = .true. 238 current_line = current_line + 1 239 read(ios_simul_file,*,iostat=ios) dens,nalle,taille 240 241 if ( dens <= 0.d0 ) then 242 call stop_application("Density of marker can not be equal or less than 0.") 243 end if 244 245 if ( taille <= 0.d0 ) then 246 call stop_application("Size of map can not be equal or less than 0.") 247 end if 248 249 if ( nalle <= 0.d0 ) then 250 call stop_application("Number allele per marker can not be equal or less than 0.") 251 end if 252 253 if ( dens > taille ) then 254 call stop_application("Density can not be greater than the size map.") 255 end if 256 257 258 message = "Problem detected at the markers description:" // & 259 " Allowed: density(float) allele number(int) chromosome size(float)" // & 260 " unity(int) unknwon char(char)" 261 262 call stop_on_error(ios,in_parsim, current_line,trim(message)) 263 ! size marker... 264 !allocate (nmk(1),stat=ios) 265 !call check_allocate(ios,'nmk (m_qtlmap_simulation)') 266 !print *,'HOLA' 267 temp = taille 268 nmk(1) = 0 269 do while ( temp >= 0 ) 270 nmk(1) = nmk(1)+1 271 temp = temp - dens 272 end do 273 274 nb_marker_used = nmk(1) 275 276 do iq=1,nqtl 277 do ic=1,ncar 278 call log_mess('Trait :'//trim(str(ic))//' Qtl:'//trim(str(iq))//' val:'//trim(str(ue(ic,iq)))//' M') 279 end do !! fin ic 280 end do !! fin iq 281 282 end subroutine set_simulation_marker
set_simulation_qtl
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
set_simulation_qtl
DESCRIPTION
NOTES
SOURCE
545 subroutine set_simulation_qtl(simQtl) 546 logical,intent(inout) :: simQtl 547 integer :: ios,i,iq,jq,chr 548 real :: l 549 character(len=LEN_DEF) :: word 550 character(len=LEN_LINE) :: line_read 551 logical :: is_ok 552 553 simtyp = .true. 554 555 current_line = current_line + 1 556 read(ios_simul_file,*,iostat=ios) nqtl 557 558 call stop_on_error(ios,in_parsim, current_line,"QTL number") 559 560 if (nqtl.gt.0) then 561 simQtl=.true. ! More than one qtl are simulated 562 end if 563 564 allocate(posiqtl(nqtl)) 565 allocate(chrqtl(nqtl)) 566 allocate (xlim(nqtl)) 567 568 569 !*** POSITION *** 570 !***************** 571 current_line = current_line + 1 572 call GET(ios_simul_file,line_read,IOSTAT=ios) 573 call stop_on_error(ios,in_parsim, current_line,"QTL : waiting for keyword 'position' and values' ") 574 575 word = trim(next_word(line_read,is_ok)) 576 if (word /= 'position') then 577 call stop_application("'position' keyword expected! **"//trim(word)//"**") 578 end if 579 580 do i=1,nqtl 581 word = trim(next_word(line_read,is_ok)) 582 if (.not. is_ok ) then 583 call stop_application("qtl "// trim(str(i))//" position expected!") 584 end if 585 posiqtl(i)=get_real(word,is_ok) 586 if ( .not. is_ok ) then 587 call stop_application("qtl "// trim(str(i))//" position expected! **"//trim(word)//"**") 588 end if 589 end do 590 591 !*** CHROMOSOME *** 592 !******************* 593 current_line = current_line + 1 594 call GET(ios_simul_file,line_read,IOSTAT=ios) 595 call stop_on_error(ios,in_parsim, current_line,"QTL : waiting for keyword 'chromosome' and values' ") 596 597 word = trim(next_word(line_read,is_ok)) 598 if (word /= 'chromosome') then 599 call stop_application("'chromosome' keyword expected! **"//trim(word)//"**") 600 end if 601 602 do i=1,nqtl 603 chrqtl(i) = trim(next_word(line_read,is_ok)) 604 if (.not. is_ok ) then 605 call stop_application("qtl "// trim(str(i))//" chromosome expected!") 606 end if 607 end do 608 609 !*** FREQUENCIES *** 610 !******************** 611 current_line = current_line + 1 612 call GET(ios_simul_file,line_read,IOSTAT=ios) 613 call stop_on_error(ios,in_parsim, current_line,"QTL : waiting for keyword 'frequency' and values' ") 614 615 word = trim(next_word(line_read,is_ok)) 616 if (word /= 'frequency') then 617 call stop_application("'frequency' keyword expected! **"//trim(word)//"**") 618 end if 619 620 do i=1,nqtl 621 word = trim(next_word(line_read,is_ok)) 622 if (.not. is_ok ) then 623 call stop_application("qtl "// trim(str(i))//" frequency expected!") 624 end if 625 xlim(i)=get_real(word,is_ok) 626 if ( .not. is_ok ) then 627 call stop_application("qtl "// trim(str(i))//" frequency expected! **"//trim(word)//"**") 628 end if 629 end do 630 631 632 ! check : * position have to be ordered and not at the same position 633 do iq=1,nqtl-1 634 do jq=iq+1,nqtl 635 if (posiqtl(iq) == posiqtl(jq)) then 636 call stop_application("QTLs"//trim(str(iq))//" and "//trim(str(jq))//"are simulated on the same position") 637 end if 638 639 if (posiqtl(iq) > posiqtl(jq)) then 640 call stop_application("QTLs have to be ordered : qtl "//trim(str(iq))//"> qtl "//trim(str(jq))) 641 end if 642 end do 643 end do 644 645 call stop_on_error(ios,in_parsim, current_line,"QTL frequencies") 646 647 end subroutine set_simulation_qtl
set_simulation_simultraits
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
set_simulation_simultraits
DESCRIPTION
NOTES
SOURCE
1086 subroutine set_simulation_simultraits 1087 integer :: ic,jc,iq,j 1088 integer :: ios 1089 character(len=LEN_DEF) :: word 1090 character(len=LEN_LINE) :: line_read 1091 logical :: is_ok 1092 1093 simulTraits = .true. 1094 current_line = current_line + 1 1095 ios=0 1096 word="" 1097 1098 do while (ios == 0 .and. word == "" ) 1099 read (ios_simul_file,*,iostat=ios) word 1100 ncar = get_int(word,is_ok) 1101 if ( .not. is_ok ) then 1102 call stop_on_error (1,in_parsim,current_line,'bad definition of number trait: ['//trim(word)//'].') 1103 end if 1104 end do 1105 1106 call stop_on_error(ios,in_parsim, current_line,"Number of traits") 1107 1108 call init_traits_simul 1109 allocate(h2(ncar)) 1110 h2=0.d0 1111 allocate(RhoP(ncar,ncar)) 1112 RhoP=0.d0 1113 allocate(RhoG(ncar,ncar)) 1114 RhoG=0.d0 1115 1116 ! FORMAT TRAITNAME NATURE HERITABILITY N=NUMBER_OF_MOD FREQ1 FREQ2..FREQN 1117 ! Constraint : FREQ1+FREQ2+..+FREQN=1 1118 do ic=1,ncar 1119 current_line = current_line + 1 1120 call GET(ios_simul_file,line_read,IOSTAT=ios) 1121 if (trim(line_read)=='') cycle ! empty line 1122 1123 if ( (ios /=0) ) then 1124 call stop_on_error (1,in_parsim,current_line,'simul traits ['//trim(str(ic))//'] is not defined.') 1125 end if 1126 1127 word = trim(next_word(line_read,is_ok)) 1128 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,'traits ['//trim(str(ic))//'] is not defined.') 1129 carac(ic)=word 1130 1131 word = trim(next_word(line_read,is_ok)) 1132 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,'nature of traits ['//trim(str(ic))//'] is not defined.') 1133 1134 if ( word /= 'i' .and. word /= 'r') then 1135 call stop_on_error (1,in_parsim,current_line,& 1136 "Availbale value for nature of traits is 'i' or 'r' **"//trim(word)//'** .') 1137 end if 1138 1139 natureY(ic)=word 1140 1141 word = trim(next_word(line_read,is_ok)) 1142 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,'Heritability of traits ['//& 1143 trim(carac(j))//'] is not defined.') 1144 h2(ic) = get_real(word,is_ok) 1145 1146 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,& 1147 'Heritability of traits ['//trim(carac(j))//'] have to be a real.') 1148 call log_mess ('Heritability of traits ['//trim(str(ic))//']:'//trim(str(h2(ic))),VERBOSE_DEF) 1149 1150 if (natureY(ic) == 'i') call read_attributes_discrete_trait(ncar,ic,line_read) 1151 1152 end do 1153 call set_soglia(ncar,natureY) 1154 !checking heritability 1155 do ic=1,ncar 1156 if (h2(ic)>1) then 1157 call stop_application("Heritability for the "//trim(str(ic))//" is greater than 1") 1158 end if 1159 if (h2(ic)<0) then 1160 call stop_application("Heritability for the "//trim(str(ic))//" is less than 0") 1161 end if 1162 end do 1163 1164 1165 call read_correlation_matrix(ncar) 1166 call read_qtl_effect_on_trait(ncar) 1167 1168 end subroutine set_simulation_simultraits
set_simulation_traits
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
set_simulation_traits
DESCRIPTION
set the array filter_car. This filter contains index trait referenced in the model.
NOTES
SOURCE
657 subroutine set_simulation_traits() 658 integer :: ic,jc,iq,j 659 integer :: ios 660 character(len=LEN_DEF) :: word 661 character(len=LEN_LINE) :: line_read 662 logical :: is_ok 663 integer :: nc 664 665 simulTraits = .true. 666 current_line = current_line + 1 667 read (ios_simul_file,*,iostat=ios) nc 668 669 if ( nc > ncar ) then 670 call log_mess(" Number of trait in the model file:"//trim(str(ncar)),ERROR_DEF) 671 call log_mess("Number of trait in the simulation parameter file:"//trim(str(nc)),ERROR_DEF) 672 call stop_application("Your model file do not match with the simulation parameter file."//& 673 " Too many traits defined in the parameter file.") 674 end if 675 676 call stop_on_error(ios,in_parsim, current_line,"Number of traits") 677 678 ! allocation of filter 679 allocate(filter_car(nc)) 680 filter_car=0 681 682 if ( (.not. allocated(carac_t)) .or. size(carac_t) <= 0 ) then 683 call stop_application("DEVEL ERROR- set_simulation_traits :: carac_t is empty!") 684 end if 685 ic=1 686 ! FORMAT: TRAIT_NAME H2 687 do while (ic<=nc) 688 current_line = current_line + 1 689 call GET(ios_simul_file,line_read,IOSTAT=ios) 690 if ( (ios /=0) ) then 691 call stop_on_error (1,in_parsim,current_line,'traits ['//trim(str(ic))//'] is not defined.') 692 end if 693 if (trim(line_read)=='' ) cycle ! empty line 694 695 word = trim(next_word(line_read,is_ok)) 696 if ( .not. is_ok ) call stop_on_error (1,in_parsim,current_line,'traits ['//trim(str(ic))//'] is not defined.') 697 698 ! carac_t is a buffer array read from the model file 699 do j=1,size(carac_t) 700 ! first case : real name, second case : name generated by the keyword all 701 if ( carac_t(j) == word .or. (LABEL_NAME_TRAIT//trim(str(j)) == carac_t(j)) ) exit 702 end do 703 704 if ( j>size(carac_t) ) then 705 call stop_application("Trait["//trim(word)//"] is not defined in the model!") 706 end if 707 filter_car(ic)=j 708 709 if (natureY_t(ic) == 'i') call read_attributes_discrete_trait(nc,ic,line_read) 710 ic=ic+1 711 end do 712 713 call set_soglia(nc,natureY_t) 714 715 if (nc == 1) call log_mess('Only one trait --> no correlation matrix to define',VERBOSE_DEF) 716 717 !call read_correlation_matrix(nc) 718 719 call read_qtl_effect_on_trait(nc) 720 721 end subroutine set_simulation_traits
set_soglia
[ Top ] [ m_qtlmap_simulation ] [ Subroutines ]
NAME
set_soglia
DESCRIPTION
NOTES
SOURCE
863 subroutine set_soglia(nc,nature) 864 integer ,intent(in):: nc 865 character(len=*),dimension(nc),intent(in) :: nature 866 double precision :: ppf,freqtot 867 integer :: ic,j 868 869 if ( .not. allocated(freqmodsimultrait) ) return 870 871 allocate (soglia(nc,size(freqmodsimultrait,2))) 872 873 do ic=1,nc 874 875 if ( nature(ic) /='i') cycle 876 877 freqtot=0.d0 878 do j=1,nbmodsimultrait(ic)-1 879 freqtot=freqtot+freqmodsimultrait(ic,j) 880 call NORPPF(freqtot,PPF) 881 soglia(ic,j)=PPF 882 end do 883 end do 884 885 end subroutine set_soglia