INPUT
NAME
INPUT
DESCRIPTION
Package input :
- Initilize Data user routines
m_qtlmap_traits
NAME
m_qtlmap_traits -- Performances routines.
SYNOPSIS
DESCRIPTION
NOTES
SEE ALSO
ALL_LABEL_MODEL
[ Top ] [ m_qtlmap_traits ] [ Constants ]
NAME
ALL_LABEL_MODEL
DESCRIPTION
Keyword to apply a generic model for all files
NOTES
MAX_QTL
[ Top ] [ m_qtlmap_traits ] [ Constants ]
NAME
MAX_QTL
DESCRIPTION
Maximum number of interaction
NOTES
nombre maximum possible de qtl en interaction definit dans le fichier model
NB_DES_MIN
[ Top ] [ m_qtlmap_traits ] [ Constants ]
NAME
NB_DES_MIN
DESCRIPTION
The minimum number of progenies used by the permutation method
NOTES
nbre minimum de descendant a considerer dans une famille (pere, pere-mere) pour faire des permutations
unit_mod
[ Top ] [ m_qtlmap_traits ] [ Constants ]
NAME
unit_mod
DESCRIPTION
the unit fortran assigned to the model file
NOTES
unit_perf
[ Top ] [ m_qtlmap_traits ] [ Constants ]
NAME
unit_perf
DESCRIPTION
the unit fortran assigned to the phenotypic file
NOTES
carac_t
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
carac_t
DESCRIPTION
buffer array for carac
DIMENSIONS
ncar
NOTES
nombre de qtl en interaction definit pour le caractere
cdt
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
cdt
DESCRIPTION
value of censured data for a trait ic/progeny kd
DIMENSIONS
nc,nbete
NOTES
cov
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
cov
DESCRIPTION
value of the covariate cov for a progeny kd
DIMENSIONS
nd,ncov
NOTES
h2_t
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
h2_t
DESCRIPTION
buffer array for h2
DIMENSIONS
ncar
NOTES
int_qtl
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
int_qtl
DESCRIPTION
indicate if a qtl is in interaction with a fixed effect in the model trait (1 otherwise 0)
DIMENSIONS
nc,nqtl,nfix
NOTES
int_qtl_t
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
int_qtl_t
DESCRIPTION
buffer array for int_qtl
DIMENSIONS
nc,nqtl,nfix
NOTES
natureY_t
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
natureY_t
DESCRIPTION
buffer array for natureY
DIMENSIONS
ncar
NOTES
nb_qtl_def
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
nb_qtl_def
DESCRIPTION
number of qtl - fixed effect interaction defined for each trait
DIMENSIONS
ncar
NOTES
nombre de qtl en interaction definit pour le caractere
ndelt
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
ndelt
DESCRIPTION
validity of a phenotypic value for a progeny kd
DIMENSIONS
nc,nbete
NOTES
niv
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
niv
DESCRIPTION
level of the fixed effect fe for a progeny kd
DIMENSIONS
nd,nfix
NOTES
nuis
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
nuis
DESCRIPTION
indicate if fixed effect / covariate are take in care (1 otherwise 0) for a specific trait ic
DIMENSIONS
nc,nfix+ncov
NOTES
nuis_t
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
nuis_t
DESCRIPTION
buffer array for nuis
DIMENSIONS
nc,nfix+ncov
NOTES
RhoG_t
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
RhoG_t
DESCRIPTION
buffer array for RhoG
DIMENSIONS
ncar,ncar
NOTES
RhoP_t
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
RhoP_t
DESCRIPTION
buffer array for RhoP
DIMENSIONS
ncar,ncar
NOTES
valeur
[ Top ] [ m_qtlmap_traits ] [ Variables ]
NAME
valeur
DESCRIPTION
phenotypic value for a progeny kd
DIMENSIONS
nc,nbete
NOTES
check_cd
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
check_cd
DESCRIPTION
set cd to 1.0 if the trait is not 'average' and cd /= 0 for each animal
NOTES
SOURCE
2344 subroutine check_cd 2345 integer :: i,j 2346 2347 call log_mess("check cd...",INFO_DEF) 2348 do i=1,ncar 2349 if ( natureY(i) /= 'a' ) then 2350 do j=1,size(cd,2) 2351 if (cd(i,j)/=0.d0) cd(i,j)=1.d0 2352 end do 2353 end if 2354 end do ! i 2355 2356 end subroutine check_cd
check_traits_and_fathers
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
check_traits_and_fathers
DESCRIPTION
check if number of progenies is greater than 1 otheriwse (stop the program with an error message)
NOTES
SOURCE
1420 subroutine check_traits_and_fathers() 1421 integer :: i,ic,nm1,nm2,nd1,nd2,jm,kd,effp 1422 call log_mess("check traits of sire family...",INFO_DEF) 1423 1424 if ( .not. allocated(nmp) ) then 1425 call stop_application("Dev error : try to read traits file before genealogy initialization...") 1426 end if 1427 1428 do ic=1,size(carac) 1429 do i=1,size(pere) 1430 nm1=nmp(i)+1 1431 nm2=nmp(i+1) 1432 effp=0 1433 do jm=nm1,nm2 1434 nd1=ndm(jm)+1 1435 nd2=ndm(jm+1) 1436 effp = effp+count(presentc(ic,nd1:nd2)) 1437 end do 1438 if (effp == 0) then 1439 call stop_application('Father ['//trim(pere(i))// & 1440 '] have no child with the trait ['//trim(carac(ic))//']') 1441 end if 1442 if (effp == 1) then 1443 call stop_application('Father ['//trim(pere(i))// & 1444 '] have only one child with the trait ['//trim(carac(ic))//']') 1445 end if 1446 end do 1447 end do 1448 call log_mess("end check traits...",INFO_DEF) 1449 end subroutine check_traits_and_fathers
fixe_structure_model
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
fixe_structure_model
DESCRIPTION
Initialize the general array : carac,nuis,int_qtl,natureY,h2,RhoP,RhoG and release buffer array associated
INPUTS
filter : vector of index to select a subset of trait defined in the model file
NOTES
SOURCE
712 subroutine fixe_structure_model(filter) 713 integer, dimension(:),intent(in) :: filter ! index filter to initialize model structure 714 integer :: i,nc,j 715 nc = size(filter) 716 if ( nc /= 0 ) then 717 allocate (carac(nc)) 718 allocate (nuis(nc,nfix+ncov)) 719 allocate (int_qtl(nc,MAX_QTL,nfix)) 720 allocate (natureY(nc)) 721 722 allocate (h2(nc)) 723 allocate (RhoP(nc,nc)) 724 allocate (RhoG(nc,nc)) 725 726 do i=1,nc 727 carac(i)=carac_t(filter(i)) 728 nuis(i,:)=nuis_t(filter(i),:) 729 int_qtl(i,:,:)=int_qtl_t(filter(i),:,:) 730 natureY(i)=natureY_t(filter(i)) 731 h2(i) = h2_t(filter(i)) 732 733 do j=i+1,nc 734 RhoP(i,j) = RhoP_t(filter(i),filter(j)) 735 RhoP(j,i) = RhoP(i,j) 736 RhoG(i,j) = RhoG_t(filter(i),filter(j)) 737 RhoG(j,i) = RhoG(i,j) 738 end do 739 end do 740 741 else 742 allocate (carac(ncar)) 743 allocate (nuis(ncar,nfix+ncov)) 744 allocate (int_qtl(ncar,MAX_QTL,nfix)) 745 allocate (natureY(ncar)) 746 allocate (h2(ncar)) 747 allocate (RhoP(ncar,ncar)) 748 allocate (RhoG(ncar,ncar)) 749 750 carac=carac_t 751 nuis=nuis_t 752 int_qtl=int_qtl_t 753 natureY=natureY_t 754 h2=h2_t 755 RhoP=RhoP_t 756 RhoG=RhoG_t 757 end if 758 759 deallocate(carac_t) 760 deallocate(nuis_t) 761 deallocate(int_qtl_t) 762 deallocate(natureY_t) 763 deallocate (h2_t) 764 deallocate (RhoP_t) 765 deallocate (RhoG_t) 766 767 end subroutine fixe_structure_model
init_model_struct
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
init_model_struct
DESCRIPTION
fill the array model and listModelTrait variable
NOTES
* the array model is used by the m_qtlmap_modlin routines and derived * the listModelTrait variable is used by the m_qtlmap_incidence routines and derived
SOURCE
780 subroutine init_model_struct() 781 integer :: alloc_stat 782 integer :: i,j,ic,k,nf,nc,ifix,ico,nqf,nl,il,kc,iqtl 783 784 call log_mess('init_model_struct',INFO_DEF) 785 786 allocate (modele(ncar,3+(3*nfix)+ncov), stat = alloc_stat) 787 call check_allocate(alloc_stat,'model') 788 allocate (listModelTrait(ncar)) 789 790 do ic=1,ncar 791 allocate ( listModelTrait(ic)%indexFixedEffect(nfix)) 792 nf=0 793 ! Pour les effets fixes et covariables 794 do ifix=1,nfix 795 if(nuis(ic,ifix) == 1) then 796 nf=nf+1 797 modele(ic,3+nf)=ifix 798 listModelTrait(ic)%indexFixedEffect(nf)=ifix 799 end if 800 end do 801 listModelTrait(ic)%nbfe = nf 802 modele(ic,1)=nf 803 804 allocate ( listModelTrait(ic)%indexCovariate(ncov)) 805 nc=0 806 do ico=1,ncov 807 if(nuis(ic,nfix+ico) == 1) then 808 nc=nc+1 809 modele(ic,3+nf+nc)=ico 810 listModelTrait(ic)%indexCovariate(nc)=ico 811 end if 812 end do 813 modele(ic,2)=nc 814 listModelTrait(ic)%nbco = nc 815 ! Pour les interactions effets fixes, QTL 816 nqf=0 817 818 if ( nb_qtl_def(ic) >= 1 ) then 819 do ifix=1,nfix 820 if(int_qtl(ic,1,ifix) == 1) then 821 nqf=nqf+1 822 modele(ic,3+nf+nc+nqf)=ifix 823 end if 824 end do 825 end if 826 modele(ic,3)=nqf 827 828 allocate (listModelTrait(ic)%nbint(MAX_QTL)) 829 listModelTrait(ic)%nbint=0 830 allocate ( listModelTrait(ic)%indexFixedEffectWithInteraction(nb_qtl_def(ic),nfix)) 831 do iqtl=1,nb_qtl_def(ic) 832 nqf = 0 833 do ifix=1,nfix 834 if(int_qtl(ic,iqtl,ifix) == 1) then 835 nqf=nqf+1 836 listModelTrait(ic)%indexFixedEffectWithInteraction(iqtl,nqf)=ifix 837 end if 838 end do 839 listModelTrait(ic)%nbint(iqtl)=nqf 840 end do 841 842 end do 843 844 deallocate (nb_qtl_def) 845 846 end subroutine init_model_struct
init_random
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
init_random
DESCRIPTION
NOTES
faire un module random
SOURCE
2718 subroutine init_random() 2719 integer :: n,clock,i 2720 integer, dimension(:), allocatable :: seed 2721 2722 !initialisation pour les routines intrasec de fortran 2723 call random_seed(size=n) 2724 allocate(seed(n)) 2725 call system_clock(count=clock) 2726 seed = clock + 37 * (/ (i - 1, i = 1, n) /) 2727 call random_seed(put = seed) 2728 deallocate(seed) 2729 2730 end subroutine init_random
initialise_struct_internal
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
initialise_struct_internal
DESCRIPTION
allocation of main array of dataset user : corperf,niveau,covar,y,ycategorial,presentc,ndelta,cd
NOTES
SOURCE
858 subroutine initialise_struct_internal() 859 integer :: alloc_stat 860 integer :: i,j,ic,k,nf,nc,ifix,ico,nqf,nl,il,icCategorial,kc,iqtl 861 !dim : nd 862 integer , dimension (:) ,allocatable :: nbvx 863 !dim: ndx,nfix 864 character(len=LEN_DEF) , dimension (:,:),allocatable :: niveau 865 866 call log_mess('initialise_struct_internal model and performance',INFO_DEF) 867 868 allocate (corperf(nd), stat = alloc_stat) 869 call check_allocate(alloc_stat,'corperf') 870 871 allocate (niveau(nd,nfix), stat = alloc_stat) 872 call check_allocate(alloc_stat,'niveau') 873 874 allocate (covar(nd,ncov), stat = alloc_stat) 875 call check_allocate(alloc_stat,'covar') 876 877 allocate (y(ncar,nd), stat = alloc_stat) 878 call check_allocate(alloc_stat,'y') 879 880 if ( ncarcat /=0 ) then 881 allocate (ycategorial(ncarcat,nd), stat = alloc_stat) 882 call check_allocate(alloc_stat,'ycategorial') 883 end if 884 885 allocate (presentc(ncar,nd), stat = alloc_stat) 886 call check_allocate(alloc_stat,'presentc') 887 888 allocate (ndelta(ncar,nd), stat = alloc_stat) 889 call check_allocate(alloc_stat,'ndelta') 890 891 allocate (cd(ncar,nd), stat = alloc_stat) 892 call check_allocate(alloc_stat,'cdt') 893 894 y=REAL_NOT_DEFINED 895 presentc=.false. 896 covar=REAL_NOT_DEFINED 897 niveau=STRING_NOT_DEFINED 898 corperf=INT_NOT_DEFINED 899 cd=0.d0 900 !---------------------- corperf et presentc 901 !for each descendant 902 do i=1,nd 903 !find the animal associated 904 do j=1,size(bete) 905 if ( animal(i) == bete(j) ) then 906 exit ! exit from this do 907 endif 908 end do 909 910 if ( j >= size(bete)+1 ) then !we does not find the animal 911 cycle ! next descendant 912 else ! animal find.... 913 corperf(i)=j 914 icCategorial = 0 915 do ic=1,ncar 916 if (natureY(ic) == 'c') icCategorial = icCategorial + 1 917 cd(ic,i)=cdt(ic,j) 918 ndelta(ic,i)=ndelt(ic,j) 919 if(cdt(ic,j) /= 0.d0) then 920 if (natureY(ic) == 'c') then 921 ycategorial(icCategorial,i)=str(valeur(ic,j)) 922 else 923 y(ic,i)=valeur(ic,j) 924 end if 925 presentc(ic,i)=.true. 926 endif 927 end do 928 929 do ifix=1,nfix 930 niveau(i,ifix)=niv(j,ifix) 931 932 ! un phenotype manquant est considere lorsque l'effet fixe, 933 ! considere dans le modele, est non renseigne 934 935 if (niveau(i,ifix) == STRING_NOT_DEFINED) then 936 do ic=1, ncar 937 do k=1, modele(ic,1) 938 if(modele(ic,3+k) == ifix) then 939 presentc(ic,i)=.false. 940 end if 941 enddo 942 do k=1, modele(ic,3) 943 if(modele(ic,3+modele(ic,1)+modele(ic,2)+k)==ifix) then 944 presentc(ic,i)=.false. 945 end if 946 enddo 947 enddo 948 endif 949 end do 950 951 do ico=1,ncov 952 covar(i,ico)=cov(j,ico) 953 954 ! un phenotype manquant est considere lorsque la covariable 955 ! est non renseignee 956 957 if (covar(i,ico) == REAL_NOT_DEFINED) then 958 do ic=1, ncar 959 do k=1, modele(ic,2) 960 if (modele(ic,3+modele(ic,1)+k) == ico) then 961 presentc(ic,i)=.false. 962 endif 963 enddo 964 enddo 965 endif 966 end do 967 end if 968 end do 969 !----------------------------------------------------------------------------- 970 ! DET DU NOMBRE DE NIVEAUX DES EFFETS FIXES 971 !----------------------------------------------------------------------------- 972 allocate (nlev(nfix), stat = alloc_stat) 973 call check_allocate(alloc_stat,'nlev') 974 975 allocate (nbvx(nd), stat = alloc_stat) 976 call check_allocate(alloc_stat,'nbvx') 977 978 allocate (nivx(nd,nfix), stat = alloc_stat) 979 call check_allocate(alloc_stat,'nivx') 980 981 !call log_mess("DEVEL TODO : listelev can be optimized",WARNING_DEF) 982 allocate (listelev(nfix,nd), stat = alloc_stat) 983 call check_allocate(alloc_stat,'listelev') 984 985 do ifix=1,nfix 986 nl=0 987 nlev(ifix)=0 988 do i=1,nd 989 nbvx(i)=0 990 enddo 991 !Comptage du nombre de fois qu'un niveau est rencontre pour le premier i et mise a 992 ! 0 pour les autres 993 do i=1,nd 994 if (niveau(i,ifix) /= STRING_NOT_DEFINED) then 995 if (nbvx(i) == 2) then 996 nbvx(i)=0 997 cycle 998 else 999 nbvx(i)=1 1000 do j=i+1,nd 1001 if (niveau(i,ifix) == niveau(j,ifix)) then 1002 nbvx(i)=nbvx(i)+1 1003 nbvx(j)=2 1004 end if 1005 end do 1006 end if 1007 end if 1008 end do 1009 1010 do i=1,nd 1011 if ( nbvx(i) /= 0 ) then 1012 nl=nl+1 1013 listelev(ifix,nl)= niveau(i,ifix) 1014 endif 1015 enddo 1016 1017 nlev(ifix)=nl 1018 1019 ! affectation des differents niveaux d'effets fixes renumerotes de 1 a nl 1020 do i=1,nd 1021 if (niveau(i,ifix) == STRING_NOT_DEFINED) nivx(i,ifix)=0 1022 do il=1,nlev(ifix) 1023 if (listelev(ifix,il) == niveau(i,ifix)) nivx(i,ifix)=il 1024 enddo 1025 enddo 1026 enddo 1027 1028 deallocate(niveau) 1029 deallocate(nbvx) 1030 call log_mess('END SUBROUTINE initialise_struct_internal',DEBUG_DEF) 1031 end subroutine initialise_struct_internal
log_infoqtldefinedbyuser
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
log_infoqtldefinedbyuser
DESCRIPTION
print in the console information about localisation of qtl simulated.
NOTES
SOURCE
2680 subroutine log_infoqtldefinedbyuser(nqtl) 2681 integer , intent(in) :: nqtl 2682 integer :: i,chr 2683 real(kind=dp) :: l 2684 logical :: ok 2685 2686 if ( nqtl <=0 ) return 2687 ! LOG***** 2688 call log_mess("--------------- QTL DEFINED --------------------",INFO_DEF) 2689 2690 do i=1,nqtl 2691 chr=0 2692 ok=.false. 2693 do while ( .not. ok ) 2694 chr = chr + 1 2695 if ( chr > nchr ) call stop_application("none chromosome "//trim(chrqtl(i))//" are defined !"); 2696 ok = chrqtl(i) == chromo(chr) 2697 end do 2698 2699 call log_mess(trim(str(i))//":Position defined by the user :"//trim(str(posiqtl(i))),INFO_DEF) 2700 l=(get_pos((posiqtl(i)-posi(chr,1)))*((posi(chr,nmk(chr))-posi(chr,1))/get_npo(chr)))+posi(chr,1) 2701 call log_mess(trim(str(i))//":Position according the sampling :"//trim(str(l)),INFO_DEF) 2702 end do 2703 2704 call log_mess("--------------------------------------------------------",INFO_DEF) 2705 end subroutine log_infoqtldefinedbyuser
manage_data
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
manage_data
FUNCTION
Normalize data for continue data and initilialized data structure sigt, xmut Count discrete data for discrete data and compute Proportions (prop structure) and Threshold (seuil structure)
SYNOPSIS
Mange process according to the type of the data (continue, discrete or categorial)
INPUTS
is_transcriptom -- boolean to set for no printing information whith transcriptomic data
is_simul -- boolean to set for not printing information in simulation case
RESULT
none.
EXAMPLE
NOTES
SEE ALSO
is_transcriptom,
set_count_discrete,
set_proportion_discrete
SOURCE
2648 subroutine manage_data(is_transcriptom,is_simul,normalize) 2649 2650 logical , intent(in) :: is_transcriptom 2651 logical , intent(in) :: is_simul 2652 logical , intent(in) :: normalize 2653 2654 if ( .not. allocated(xmut) )allocate(xmut(ncar)) 2655 if ( .not. allocated(sigt) )allocate(sigt(ncar)) 2656 2657 if ( normalize ) then 2658 call normalize_data(is_transcriptom,is_simul) 2659 else 2660 xmut=0.d0 2661 sigt=1.d0 2662 end if 2663 2664 call set_count_discrete() 2665 call set_proportion_discrete() 2666 2667 end subroutine manage_data
manage_simulator_traits
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
manage_simulator_traits
DESCRIPTION
Manage the type of simulator (depends opt_cal/nature of traits)
4 cases :
1) permutation are switch on : -> permuted data according to the type analysis
a) MULTITRAIT ANALYSIS -> PERMUTATION BY LINE
b) UNITRAIT ANALYSIS -> PERMUTATION on each trait
2) permutation are switch off and analysis is a multitrait
--> simulation of traits with heritability and correlation matrix (according to the nature)
3) permutation are switch off and analysis is unitrait
--> simulation of traits with heritability according to the nature of the trait
NOTES
Simulation des typages sur trois generations , reperage des alleles au qtl des descendants.
SOURCE
1603 subroutine manage_simulator_traits(analyse_is_multi_traits,permute_mode,step,simqtl,nqtl,qtl,croisement) 1604 logical ,intent(in) :: analyse_is_multi_traits ! generaly the simulation 1605 logical ,intent(in) :: permute_mode 1606 integer,intent(in) :: nqtl,step 1607 1608 integer ,dimension(np+nfem+nd,nqtl,2) ,intent(inout) :: qtl 1609 character(len=LEN_BUFFER_WORD) , intent(in) :: croisement 1610 logical ,intent(in) :: simqtl 1611 1612 real (kind=dp),dimension (1) :: h2_t 1613 1614 real (kind=dp),dimension (1,size(animal)) :: ys 1615 1616 real (kind=dp),dimension (1,1) :: bidonrho 1617 integer :: i,j,ic 1618 character(len=1) :: nat 1619 1620 if (ncar <= 0 ) then 1621 call stop_application("DEVEL ERROR : NCAR is not initialized") 1622 end if 1623 1624 if ( .not. allocated(natureY) ) then 1625 call stop_application("DEVEL ERROR : NATUREY is not initialized") 1626 endif 1627 if ( size(natureY) <= 0 ) then 1628 call stop_application("DEVEL ERROR : NATUREY is empty") 1629 end if 1630 1631 !the permutation mode is active, none nimulator are called 1632 if ( permute_mode ) then 1633 call set_estime() 1634 if ( .not. analyse_is_multi_traits ) then 1635 call log_mess("** Permutation for unitrait analysis **",INFO_DEF) 1636 call sim_perf_shuffling(.false.) 1637 else !traits line are permuted between animals 1638 call log_mess("** Permutation for multitraits analysis **",INFO_DEF) 1639 call sim_perf_shuffling(.true.) 1640 end if 1641 return 1642 end if 1643 1644 1645 if ( analyse_is_multi_traits .or. allocated(natureY=='r')) then 1646 !check the nature of all traits (have to bo the same) 1647 nat=natureY(1) 1648 do i=1,size(natureY)-1 1649 if (natureY(i) /= natureY(i+1)) then 1650 call stop_application("You have to defined a model value with the same nature of trait"//& 1651 " for a multitrait analysis - [trait "//trim(str(i))//" nature:"//natureY(i)//"] [trait "//& 1652 trim(str(i+1))//" nature:"//natureY(i+1)//"]") 1653 end if 1654 end do 1655 1656 if ( nat == 'r' ) then 1657 call log_mess("** Simulation for real multitrait analysis ** ",INFO_DEF) 1658 call sim_perf_tirage(ncar,RhoP,h2,y) 1659 if (simqtl) then 1660 do ic=1,ncar 1661 call sim_QTL(ic,nqtl,step,qtl,croisement) 1662 end do 1663 end if 1664 end if 1665 1666 if ( nat == 'i' ) then 1667 call log_mess(" ** Simulation for discrete multitrait analysis ** ",INFO_DEF) 1668 call sim_perf_tirage(ncar,RhoP,h2,y) 1669 if (simqtl) then 1670 do ic=1,ncar 1671 call sim_QTL(ic,nqtl,step,qtl,croisement) 1672 end do 1673 end if 1674 do i=1,ncar 1675 call sim_transform_discret(i) 1676 end do 1677 end if 1678 1679 if ( nat == 'c' ) then 1680 call stop_application("None simulator are defined for categorial data") 1681 end if 1682 1683 if ( nat == 'a' ) then 1684 call log_mess("** Simulation for real multitrait analysis ** ",INFO_DEF) 1685 call sim_perf_tirage(ncar,RhoP,h2,y) 1686 do i=1,ncar 1687 do j=1,size(y,2) 1688 if (cd(i,j)/=0) y(i,j)=ys(1,j)/sqrt(cd(i,j)) 1689 end do 1690 end do 1691 1692 if (simqtl) then 1693 do ic=1,ncar 1694 call sim_QTL(ic,nqtl,step,qtl,croisement) 1695 end do 1696 end if 1697 end if 1698 1699 else !otherwise we call for each trait the simulator corresponding to the nature of the traits 1700 bidonrho=0.d0 1701 do i=1,ncar 1702 h2_t(1) = h2(i) 1703 call log_mess("** Simulator for TRAIT ["//trim(str(i))//"]->Nature:"//trim(natureY(i)//"] **"),INFO_DEF) 1704 select case (natureY(i)) 1705 case ('r') !! real / continue value 1706 call sim_perf_tirage(1,bidonrho,h2_t,ys) 1707 y(i,:)=ys(1,:) 1708 if (simqtl) then 1709 call sim_QTL(i,nqtl,step,qtl,croisement) 1710 end if 1711 case ('i') !! discrete value 1712 call sim_perf_tirage(1,bidonrho,h2_t,ys) 1713 y(i,:)=ys(1,:) 1714 if (simqtl) then 1715 call sim_QTL(i,nqtl,step,qtl,croisement) 1716 end if 1717 call sim_transform_discret(i) 1718 case ('c') 1719 call stop_application("None simulator are defined for categorial data") 1720 case ('a') 1721 call sim_perf_tirage(1,bidonrho,h2_t,ys) 1722 1723 do j=1,size(y,2) 1724 if (cd(i,j)/=0) y(i,j)=ys(1,j)/sqrt(cd(i,j)) 1725 end do 1726 1727 if (simqtl) then 1728 call sim_QTL(i,nqtl,step,qtl,croisement) 1729 end if 1730 1731 !call stop_application("None simulator are defined for average data") 1732 case default 1733 call stop_application("Unknown nature of data ["//natureY(i)//"]") 1734 end select 1735 end do 1736 end if 1737 1738 ! 1739 call set_estime() 1740 1741 end subroutine manage_simulator_traits
normal
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
normal
DESCRIPTION
NOTES
faire un module random
SOURCE
2743 function normal(mean,sigma) !returns a normal distribution 2744 real(kind=dp) , intent(in) :: mean,sigma 2745 real(kind=dp) normal,tmp 2746 integer flag 2747 real(kind=dp) fac,gsave,rsq,r1,r2,u 2748 save flag,gsave 2749 data flag /0/ 2750 !$omp threadprivate (flag) 2751 if (flag.eq.0) then 2752 rsq=2.0 2753 do while(rsq.ge.1.0.or.rsq.eq.0.0) ! new from for do 2754 call random_number(u) 2755 r1=2.0*u-1.0 2756 call random_number(u) 2757 r2=2.0*u-1.0 2758 rsq=r1*r1+r2*r2 2759 enddo 2760 fac=sqrt(-2.0*log(rsq)/rsq) 2761 gsave=r1*fac 2762 tmp=r2*fac 2763 flag=1 2764 else 2765 tmp=gsave 2766 flag=0 2767 endif 2768 normal=tmp*sigma+mean 2769 return 2770 end function normal
normalize_data
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
normalize_data
FUNCTION
Compute xmut (Mean) and sigt (Standart deviation) and normalize Y array.
SYNOPSIS
Normalizing data for continue traits
INPUTS
* is_transcriptom -- boolean to set for no printing information whith transcriptomic data
* is_simul -- boolean to set for not printing information in simulation case
RESULT
none.
EXAMPLE
NOTES
If none traits are found found for one founder, the analyse stop
SEE ALSO
xmut sigt natureY MATH_QTLMAP_G01FAF
SOURCE
2384 subroutine normalize_data(is_transcriptom,is_simul) 2385 2386 logical , intent(in) :: is_transcriptom 2387 logical , intent(in) :: is_simul 2388 integer :: kd,i 2389 real (kind=dp) :: sy,sy2,eff 2390 2391 call log_mess('SUBROUTINE normalize_data',DEBUG_DEF) 2392 2393 xmut = 0.d0 2394 sigt = 1.d0 2395 2396 do i=1,ncar 2397 if (.not.(natureY(i) == 'r' .or. natureY(i) == 'a')) cycle 2398 sy=0.d0 2399 sy2=0.d0 2400 eff=0.d0 2401 xmut(i)=0.d0 2402 do kd=1,nd 2403 if(presentc(i,kd)) then 2404 sy2=sy2+y(i,kd)*y(i,kd) 2405 sy=sy+y(i,kd) 2406 eff=eff+1.d0 2407 endif 2408 enddo 2409 if (eff == 0.d0) then 2410 call stop_application("Can not find animal with performance for trait :"//trim(carac(i))) 2411 end if 2412 2413 xmut(i)=sy/eff 2414 sigt(i)=dble(sqrt((sy2/eff)-xmut(i)*xmut(i))) 2415 2416 ! normalisation des donn�es 2417 if ( .not. is_transcriptom .and. .not. is_simul ) then 2418 call log_mess('Trait ['//trim(carac(i))//'] qtlmap normalize data with means:'//& 2419 str(xmut(i))//'+-'//str(sigt(i)),INFO_DEF) 2420 end if 2421 2422 do kd=1,nd 2423 if(presentc(i,kd)) then 2424 y(i,kd)=(y(i,kd)-xmut(i))/sigt(i) 2425 end if 2426 end do 2427 end do 2428 2429 call log_mess('END SUBROUTINE normalize_data',DEBUG_DEF) 2430 end subroutine normalize_data
read_model
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
read_model
DESCRIPTION
Read the model file
INPUTS
file : model file name filter_car : vector of index to select a subset of trait defined in the model file
NOTES
SOURCE
341 subroutine read_model(file,filter_car) 342 character(len=LENGTH_MAX_FILE),intent(in) :: file 343 integer , dimension(:), pointer ,optional :: filter_car 344 integer :: ios,alloc_stat 345 character(len=LEN_LINE) ::line_read,line_orig,saveLineRead 346 character(len=LEN_DEF) :: word,word2 347 integer :: i,j,l,k,ii,ncar2,iqtl,ic 348 logical :: is_ok,all_mode,interQtl 349 character(len=LEN_DEF),dimension(:),allocatable :: lcar,carac_t2,natureY_t2 350 351 call log_mess('SUBROUTINE read_model',DEBUG_DEF) 352 call log_mess('reading model file...',INFO_DEF) 353 open(unit_mod,FILE=file,IOSTAT=ios) 354 355 if ( ios /= 0 ) then 356 call stop_application('Cannot open model file ['//trim(file)//']') 357 endif 358 359 call log_mess('Model description....',VERBOSE_DEF) 360 361 read(unit_mod,*,IOSTAT=ios) ncar 362 if (ios /= 0) then 363 call stop_application('Cannot read number of trait ['//trim(file)//']' & 364 //' line:1') 365 end if 366 367 call log_mess('Number of traits :'//trim(str(ncar)),VERBOSE_DEF) 368 369 read(unit_mod,*,IOSTAT=ios) nfix,ncov 370 if (ios /= 0) then 371 call stop_application('Cannot read number of fixed effect and number of covariate ['//trim(file)//']' & 372 //' line:2') 373 end if 374 call log_mess('Number of fixed effects :'//trim(str(nfix)),VERBOSE_DEF) 375 call log_mess('Number of covariates :'//trim(str(ncov)),VERBOSE_DEF) 376 377 !Lecture des cov et effets fixe 378 CALL GET(unit_mod,line_read,IOSTAT=ios) 379 if ( ios /= 0) then 380 call stop_application('give name of fixed effect and covariate in the file ['//trim(file)//']' & 381 //' line:3'//' by default write "none"') 382 endif 383 384 allocate (namefix(nfix)) 385 allocate (namecov(ncov)) 386 387 l = 3 388 do i=1,nfix 389 namefix(i) = trim(next_word(line_read,is_ok)) 390 if ( .not. is_ok ) call stop_on_error (1,file,l,'fixed effect name ['//trim(str(i))//']') 391 call log_mess('Fixed effect :'//trim(trim(namefix(i))),VERBOSE_DEF) 392 end do 393 do i=1,ncov 394 namecov(i) = next_word(line_read,is_ok) 395 if ( .not. is_ok ) call stop_on_error (1,file,l,'covariate name ['//trim(str(i))//']') 396 call log_mess('Covariate :'//trim(trim(namecov(i))),VERBOSE_DEF) 397 end do 398 399 !lecture des caracteres 400 !---------------------- 401 ALLOCATE (carac_t(ncar), stat = alloc_stat) 402 CALL check_allocate(alloc_stat,'carac') 403 ALLOCATE (nuis_t(ncar,nfix+ncov), stat = alloc_stat) 404 CALL check_allocate(alloc_stat,'nuis') 405 ALLOCATE (int_qtl_t(ncar,MAX_QTL,nfix), stat = alloc_stat) 406 CALL check_allocate(alloc_stat,'int_qtl') 407 int_qtl_t=0 408 409 ALLOCATE (nb_qtl_def(ncar)) 410 nb_qtl_def=0 411 allocate(natureY_t(ncar), stat = alloc_stat) 412 call check_allocate(alloc_stat,'natureY') 413 414 all_mode = .false. 415 ncarcat = 0 416 do i=1,ncar 417 !for each trait we read the model description 418 l = l+1 419 call GET(unit_mod,line_read,IOSTAT=ios) 420 421 if ( ios /= 0) then 422 call stop_application('none description trait find for the'// trim(str(i))//' trait ' & 423 //' line:'//trim(str(l))//'model file:['//trim(file)//']') 424 end if 425 word = next_word(line_read,is_ok) 426 if ( .not. is_ok ) call stop_on_error (1,file,l,'trait name') 427 428 !! ALL is a generic name for apply model on all traits 429 if (word == ALL_LABEL_MODEL) then 430 call log_mess('Keywords ['//ALL_LABEL_MODEL//'] detected. apply the description for each trait',INFO_DEF) 431 432 do j=i,ncar 433 word2=LABEL_NAME_TRAIT//trim(str(i)) 434 carac_t(i) = trim(word2) 435 end do 436 all_mode = .true. 437 end if 438 if (.not. all_mode) then 439 carac_t(i) = trim(word) 440 else 441 do j=i,ncar 442 carac_t(j) = LABEL_NAME_TRAIT//trim(str(j)) 443 end do 444 end if 445 446 word = trim(next_word(line_read,is_ok)) 447 448 if ( .not. is_ok ) call stop_on_error (1,file,l,'Trait nature') 449 450 if ( word /= 'r' .and. word /= 'i' .and. word /= 'c' .and.word /= 'a') then 451 call stop_application('Trait nature have to be r(real), i(integer), c(categorial) or a(average)'//& 452 ' token:'//trim(word)) 453 end if 454 455 if ( .not. all_mode ) then 456 natureY_t(i)=trim(word) 457 else 458 do j=i,ncar 459 natureY_t(j)=trim(word) 460 end do 461 end if 462 463 if ( .not. all_mode ) then 464 if ( natureY_t(i) == 'c') ncarcat = ncarcat +1 465 else 466 467 do j=i,ncar 468 if ( natureY_t(j) == 'c') ncarcat = ncarcat +1 469 end do 470 end if 471 472 do j=1,(nfix+ncov) 473 word = trim(next_word(line_read,is_ok)) 474 if ( .not. is_ok ) call stop_on_error (1,file,l,'model for trait ['//trim(carac_t(i))//']') 475 476 nuis_t (i,j) = get_int(word,is_ok) 477 478 if ( .not. is_ok ) call stop_on_error (1,file,l,'possible value for model : 0,1 [val:'//trim(word)//']') 479 !ALL mode 480 if (all_mode) then 481 do k=i+1,ncar 482 nuis_t (k,j) = nuis_t(i,j) 483 end do 484 end if 485 end do 486 487 iqtl=0 488 interQtl=(nfix/=0) 489 490 do while (interQtl) 491 saveLineRead=line_read 492 word = trim(next_word(saveLineRead,is_ok)) 493 494 if ( is_ok ) then 495 iqtl=iqtl+1 496 if ( MAX_QTL < iqtl ) call stop_on_error(1,file,l,& 497 "You can not define more than ["//trim(str(MAX_QTL))//"] qtl interaction.") 498 do j=1,nfix 499 word = trim(next_word(line_read,is_ok)) 500 if ( .not. is_ok ) call stop_on_error (1,file,l,'possible value for model : 0,1 [val:'//trim(word)//'] numQTL '//& 501 '['// trim(str(iqtl))//']') 502 503 int_qtl_t(i,iqtl,j) = get_int(word,is_ok) 504 if ( .not. is_ok ) call stop_on_error (1,file,l,'possible value for model : 0,1 [val:'//trim(word)//'] numQTL '//& 505 '['// trim(str(iqtl))//']') 506 !ALL mode 507 if (all_mode) then 508 do k=i+1,ncar 509 int_qtl_t (k,iqtl,j) = int_qtl_t(i,iqtl,j) 510 end do 511 end if 512 end do 513 else 514 interQtl = .false. 515 end if 516 end do 517 518 if (all_mode) then 519 nb_qtl_def = iqtl 520 else 521 nb_qtl_def(i) = iqtl 522 end if 523 524 call log_mess('Description of ['//trim(carac_t(i))//']... ok',VERBOSE_DEF) 525 526 if (all_mode) then 527 exit ! go out 528 end if 529 end do 530 531 !Lecture des heritabilite et des matrices de correlations genetique et phenotypique 532 ! on lit cette matrice si la ligne "matrice_correlation" existe 533 allocate(h2_t(ncar)) 534 allocate(RhoP_t(ncar,ncar)) 535 allocate(RhoG_t(ncar,ncar)) 536 537 h2_t=0.5d0 538 RhoP_t=0.5d0 539 RhoG_t=0.5d0 540 541 ios=0 542 line_read='' 543 do while( ios == 0 .and. (trim(line_read)=='') ) 544 call GET(unit_mod,line_read,IOSTAT=ios) 545 end do 546 saveLineRead=line_read 547 word = trim(next_word(line_read,is_ok)) 548 549 if ( trim(word) == "correlation_matrix" ) then 550 ios=0 551 ic=0 552 do while( ios == 0 .and. ic < ncar ) 553 call GET(unit_mod,line_read,IOSTAT=ios) 554 if (ios /= 0) cycle 555 if (trim(line_read) =='' ) cycle 556 saveLineRead=line_read 557 ic=ic+1 558 do i=1,ncar 559 !Heritability 560 if ( ic == i ) then 561 word = trim(next_word(line_read,is_ok)) 562 if ( .not. is_ok ) then 563 call stop_application("Bad definition of heritability :"//trim(saveLineRead)) 564 end if 565 h2_t(ic) = get_real(word,is_ok) 566 if ( .not. is_ok ) then 567 call stop_application("Bad definition of heritability :"//trim(saveLineRead)) 568 end if 569 570 !checking heritability 571 if (h2_t(ic)>1) then 572 call stop_application("Heritability for the "//trim(str(ic))//" is greater than 1") 573 end if 574 if (h2_t(ic)<0) then 575 call stop_application("Heritability for the "//trim(str(ic))//" is less than 0") 576 end if 577 end if 578 !Phenotype correlation 579 if ( i < ic ) then 580 word = trim(next_word(line_read,is_ok)) 581 if ( .not. is_ok ) then 582 call stop_application("Bad definition of Phenotypique correlation :"//trim(saveLineRead)) 583 end if 584 RhoP_t(ic,i) = get_real(word,is_ok) 585 if ( .not. is_ok ) then 586 call stop_application("Bad definition of Phenotypique correlation :"//trim(saveLineRead)) 587 end if 588 RhoP_t(i,ic) = RhoP_t(ic,i) 589 end if 590 591 !Genetic correlation 592 if ( i > ic ) then 593 word = trim(next_word(line_read,is_ok)) 594 if ( .not. is_ok ) then 595 call stop_application("Bad definition of Phenotypique correlation :"//trim(saveLineRead)) 596 end if 597 RhoG_t(ic,i) = get_real(word,is_ok) 598 if ( .not. is_ok ) then 599 call stop_application("Bad definition of Phenotypique correlation :"//trim(saveLineRead)) 600 end if 601 RhoG_t(i,ic) = RhoG_t(ic,i) 602 end if 603 end do 604 end do 605 606 if ( ic /= ncar ) then 607 print *,' Format : ' 608 print *,'correlation_matrix' 609 print *,'h2[C1] RhoG[C1,C2] RhoG[C1,C3] ....' 610 print *,'RhoP[C1,C2] h2[C2] RhoG[C2,C3] ....' 611 print *,'RhoP[C1,C3] RhoP[C2,C3] h2[C3] ....' 612 print *,' ** ' 613 call stop_application(" ** Bad definition of Correlation_matrix **") 614 end if 615 616 saveLineRead='' 617 else 618 call log_mess("correlation_matrix keyword are not found. Default value is 0.5 for "//& 619 "heritability and phenotypic,genotypic correlations",WARNING_DEF); 620 621 end if 622 623 ! pour l instant on fait rien si pas de filtre entree 624 if (.not. present(filter_car) ) return 625 626 !AJOUT OFI 627 ! SI L UTILISATEUR MET UNE LISTE, ON CONSIDERE CETTE LISTE COMME ETANT LA LISTE DES CARACTERES A PRENDRE EN COMPTE DANS L ANALYSE 628 ! PAR DEFAUT, donc,ON A ANALYSE TOUS LES CARACTERES DEFINIS 629 allocate (lcar(ncar)) 630 lcar='' 631 ncar2=0 632 line_read=saveLineRead 633 ! Selection of particular trait 634 ios=0 635 ! print *,1,trim(saveLineRead) 636 ! print *,2,trim(line_read) 637 ! stop 638 do while ( trim(line_read)=='' .and. ios == 0 ) 639 call GET(unit_mod,line_read,IOSTAT=ios) 640 end do 641 642 do while( ios == 0) 643 ! liste des caracteres a prendre en consideration lors de l analyse 644 line_orig=trim(line_read) 645 do while ( trim(line_read)/='' ) 646 ncar2 = ncar2 + 1 647 if (ncar2 > ncar ) then 648 call log_mess("Filter trait:"//trim(line_orig),WARNING_DEF) 649 ios=-1 650 exit 651 end if 652 lcar(ncar2)=trim(next_word(line_read,is_ok)) 653 end do 654 line_read='' 655 do while ( trim(line_read)=='' .and. ios == 0 ) 656 call GET(unit_mod,line_read,IOSTAT=ios) 657 end do 658 end do 659 660 if (associated(filter_car)) deallocate(filter_car) 661 allocate (filter_car(ncar2)) 662 663 if (ncar2 == 0 ) return 664 665 if (ncar2 > ncar) then 666 call stop_application ("Too many traits defines in the list trait filter [Model file:"//trim(file)//"]") 667 end if 668 669 close(unit_mod) 670 671 do ii=1,ncar2 672 do i=1,ncar 673 if ( trim(carac_t(i)) == trim(lcar(ii))) exit 674 end do 675 if ( i > ncar ) then 676 !call stop_application("Trait ["//trim(lcar(ii))//"] is not defined in the model file [Model file:"//trim(file)//"]") 677 i=get_int(lcar(ii),is_ok) 678 if (is_ok) then 679 if ( i > ncar .or. i<=0 ) then 680 call stop_application("Bad definition of the index trait ["//trim(lcar(ii))//& 681 "] in the model file [Model file:"//trim(file)//"]") 682 end if 683 684 call log_mess(" ***> get index trait :"//trim(str(i)),INFO_DEF) 685 filter_car(ii)=i 686 else 687 call stop_application("Trait ["//trim(lcar(ii))//"] is not defined in the model file [Model file:"//trim(file)//"]") 688 end if 689 690 else 691 filter_car(ii)=i 692 end if 693 694 end do 695 696 deallocate (lcar) 697 698 call log_mess('END SUBROUTINE read_model',DEBUG_DEF) 699 end subroutine read_model
read_perf_by_column
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
read_perf_by_column
DESCRIPTION
read the phenotypic file with the format :
For each animal, its ID (identical to the ID given in the pedigree file) is followed by information
about nuisance effects (fixed effect levels, covariable value) and then by three information for each trait :
* the performance
* an 0/1 variable IP which indicates if (IP=1) or not (IP=0) the trait was measured for this animal and must be included in the analysis
* and 0/1 variable (IC) which indicates if (IC=0) it was censored or not (IC=1), this IC information being needed for survival analysis (by default IC=1).
Grammar : <animal_id> <fixed_effect1> <fixed_effect2>... <cov1>..<cov2> <IP> <IC> <VALUE>
The following array are filling : bete, niv, cov, valeur, cdt, ndelt
INPUTS
files : list of phenotypic files filter : vector of index to select a subset of trait defined in the model file
NOTES
SOURCE
1296 subroutine read_perf_by_column(files,filter) 1297 character(len=LENGTH_MAX_FILE),dimension(MAX_FILES_INPUT),intent(in) :: files 1298 integer, dimension(:),intent(in) :: filter 1299 integer :: ios,alloc_stat 1300 integer :: nbete 1301 character(len=LEN_DEF) :: vs,temp_vs 1302 character(len=LEN_LINE) ::line_read 1303 character(len=LEN_BUFFER_WORD) :: word_token 1304 integer :: l,i,j,k 1305 logical :: is_ok,filterActive 1306 character(len=LEN_BUFFER_WORD) :: file 1307 integer :: futureNcar 1308 real (kind=dp), dimension (:,:), allocatable :: cdt1,valeur1 1309 !dim : nc,nbete 1310 integer , dimension (:,:),allocatable :: ndelt1 1311 character(len=LEN_DEF), dimension (:), allocatable :: temp 1312 character(len=LEN_DEF), dimension (:,:), allocatable :: temp2 1313 real(kind=dp), dimension (:,:), allocatable :: temp3 1314 1315 call log_mess('SUBROUTINE read_perf_by_column',DEBUG_DEF) 1316 file = files(1) 1317 open(unit_perf,FILE=file,IOSTAT=ios) 1318 call log_mess('reading traits file...',INFO_DEF) 1319 if ( ios /= 0 ) then 1320 call stop_application('Cannot open the traits file ['//trim(file)//']') 1321 endif 1322 1323 futureNcar = ncar 1324 filterActive=.false. 1325 if (size(filter) /=0 ) then 1326 filterActive=.true. 1327 futureNcar= size(filter) 1328 end if 1329 1330 nbete=MAX_ANIMAL 1331 allocate (bete(nbete), stat = alloc_stat) 1332 allocate (niv(nbete,nfix), stat = alloc_stat) 1333 allocate (cov(nbete,ncov), stat = alloc_stat) 1334 allocate (valeur1(ncar,nbete), stat = alloc_stat) 1335 allocate (cdt1(ncar,nbete), stat = alloc_stat) 1336 allocate (ndelt1(ncar,nbete), stat = alloc_stat) 1337 1338 ios = 0 1339 nbete = 0 1340 i=1 1341 rewind(unit_perf) 1342 do while ( ios == 0) 1343 read (unit_perf,*,iostat=ios) bete(i),(niv(i,j),j=1,nfix),(cov(i,j),j=1,ncov),(valeur1(j,i),cdt1(j,i),ndelt1(j,i),j=1,ncar) 1344 if ( trim(bete(i))/='' .and. ( ios == 0) ) then 1345 nbete = nbete+1 1346 ! print *,trim(bete(i)),(niv(i,j),j=1,nfix),(cov(i,j),j=1,ncov),(valeur1(j,i),cdt1(j,i),ndelt1(j,i),j=1,ncar) 1347 i = i + 1 1348 endif 1349 end do 1350 close(unit_perf) 1351 1352 if ( nbete == 0 ) then 1353 call stop_application("Can not read in the trait file:["//trim(file)//"] animal:["//trim(bete(i))//"]") 1354 end if 1355 1356 allocate (valeur(futureNcar,nbete)) 1357 allocate (cdt(futureNcar,nbete)) 1358 allocate (ndelt(futureNcar,nbete)) 1359 1360 do i=1,nbete 1361 1362 ! do j=1,(i-1) 1363 ! if ( bete(i) == bete(j)) then 1364 ! call stop_on_error (1,file,l,'animal ['//trim(bete(i))//'] have two lines definitions of traits !!') 1365 ! end if 1366 ! end do 1367 do j=1,ncar 1368 ! if ( futureNcar /= ncar ) then 1369 if ( filterActive ) then 1370 do k=1,size(filter) 1371 if ( filter(k) == j ) exit 1372 end do 1373 ! the index car is not used 1374 if ( k > size(filter)) cycle 1375 else 1376 k = j 1377 end if 1378 cdt(k,i) = cdt1(j,i) 1379 valeur(k,i) = valeur1(j,i) 1380 ndelt(k,i) = ndelt1(j,i) 1381 end do 1382 end do 1383 deallocate (valeur1) 1384 deallocate (cdt1) 1385 deallocate (ndelt1) 1386 1387 !on retasse les tableau bete,niv,cov 1388 allocate (temp(nbete)) 1389 temp = bete(:nbete) 1390 deallocate (bete) 1391 allocate (bete(nbete)) 1392 bete = temp 1393 1394 allocate (temp2(nbete,nfix)) 1395 temp2 = niv(:nbete,:) 1396 deallocate (niv) 1397 allocate (niv(nbete,nfix)) 1398 niv=temp2 1399 1400 allocate (temp3(nbete,ncov)) 1401 temp3 = cov(:nbete,:) 1402 deallocate (cov) 1403 allocate (cov(nbete,ncov)) 1404 cov=temp3 1405 1406 call log_mess('END SUBROUTINE read_perf_by_column',DEBUG_DEF) 1407 end subroutine read_perf_by_column
read_perf_by_line
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
read_perf_by_line
DESCRIPTION
read the phenotypic file with the expression quantitative trait value format :
The header line is the list of animals phenotyped. The following line are the fixed effects, covariates and finally the phenotype.
The format of the nuisances effects and phenotype line is :
<IDANIMAL> <FIXED_EFFECT1> <FIXED_EFFECT2>...<COV1> <COV2>...<VALUE1> <VALUE2>...
For missing data, insert a character string which is not interpretable as a numeric(e.g. n/a).
The following array are filling : bete, niv, cov, valeur, cdt, ndelt
INPUTS
files : list of phenotypic files filter : vector of index to select a subset of trait defined in the model file
NOTES
SOURCE
1052 subroutine read_perf_by_line(files,filter) 1053 character(len=LENGTH_MAX_FILE), dimension(MAX_FILES_INPUT), intent(in) :: files 1054 integer, dimension(:),intent(in) :: filter 1055 integer :: ios,alloc_stat 1056 integer :: nbete,current_number,l,i,j,futureNcar,k,ii 1057 character(len=LEN_DEF) :: token 1058 character(len=LEN_LINE) :: line_read 1059 logical :: is_ok,finloop 1060 character(len=1) :: buff_read 1061 character(len=LENGTH_MAX_FILE) :: file 1062 character(len=LEN_BUFFER_WORD) :: name 1063 character(len=LEN_BUFFER_WORD),dimension(:),allocatable :: bete_char 1064 character(len=LEN_DEF) ,dimension(:),allocatable :: niv_t 1065 character(len=LEN_DEF) ,dimension(:,:),allocatable ::locvaleur 1066 real(kind=dp) ,dimension(:),allocatable :: cov_t 1067 real(kind=dp) :: v 1068 1069 call log_mess('SUBROUTINE read_perf_by_line',DEBUG_DEF) 1070 call log_mess('reading transcriptom data file...',INFO_DEF) 1071 nbete = 0 1072 futureNcar = ncar 1073 if (size(filter) /=0 ) futureNcar= size(filter) 1074 1075 call log_mess(' ** transcriptomic file : do not write number of animal at the first line ** ',WARNING_DEF) 1076 1077 file = files(1) 1078 if (ncar <= 0 .or. nfix < 0 .or. ncov < 0 ) then 1079 call stop_application('Devel error: call read_model before read_traits') 1080 end if 1081 1082 open(unit_perf,FILE=file,IOSTAT=ios,FORM="formatted",recl=2**20) 1083 1084 if ( ios /= 0 ) then 1085 call stop_application('The application can not open the traits file ['//trim(file)//']') 1086 endif 1087 1088 line_read = '' 1089 l = 1 1090 ios = 0 1091 1092 ! blank line.... 1093 do while ( (ios == 0) .and. trim(line_read)=='') 1094 l = l + 1 1095 call GET(unit_perf,line_read,IOSTAT=ios) 1096 end do 1097 1098 if ( ios /= 0 ) then 1099 call stop_application('The application can not read the first line ['//trim(file)//'].'// & 1100 ' Split the file to resolve this error') 1101 endif 1102 1103 is_ok = .true. 1104 1105 do while ( is_ok ) 1106 token = next_word(line_read,is_ok) 1107 nbete = nbete + 1 1108 end do 1109 1110 nbete = nbete -1 1111 1112 rewind(unit_perf) 1113 1114 call log_mess('Number of animals detected:'//trim(str(nbete)),VERBOSE_DEF) 1115 1116 allocate (bete(nbete), stat = alloc_stat) 1117 call check_allocate(alloc_stat,'bete') 1118 allocate (bete_char(nbete), stat = alloc_stat) 1119 call check_allocate(alloc_stat,'bete_char') 1120 allocate (niv(nbete,nfix), stat = alloc_stat) 1121 call check_allocate(alloc_stat,'niv') 1122 1123 allocate(niv_t(nbete)) 1124 1125 allocate (cov(nbete,ncov), stat = alloc_stat) 1126 call check_allocate(alloc_stat,'cov') 1127 1128 allocate (cov_t(nbete), stat = alloc_stat) 1129 1130 allocate (locvaleur(futureNcar,nbete), stat = alloc_stat) 1131 allocate (valeur(futureNcar,nbete), stat = alloc_stat) 1132 call check_allocate(alloc_stat,'locvaleur') 1133 allocate (cdt(futureNcar,nbete), stat = alloc_stat) 1134 call check_allocate(alloc_stat,'cdt') 1135 allocate (ndelt(futureNcar,nbete), stat = alloc_stat) 1136 call check_allocate(alloc_stat,'ndelt') 1137 1138 l = 1 1139 ! animal name header 1140 read(unit_perf,*,iostat=ios) (bete_char(i),i=1,nbete) 1141 bete = bete_char 1142 if ( ios /= 0 ) then 1143 call log_mess('Problem detected at the line:'//trim(str(l)),ERROR_DEF) 1144 call stop_application('No header (animals id) are finding in traits file:'//file) 1145 end if 1146 1147 i = 1 1148 do while ( i <= nfix ) 1149 l = l + 1 1150 read(unit_perf,*,iostat=ios) name,(niv_t(j),j=1,nbete) 1151 if ( trim(name) == '') cycle 1152 if ( ios /= 0 ) then 1153 call log_mess('Problem detected at the line:'//trim(str(l)),ERROR_DEF) 1154 call stop_application('The application can not initialize the fixed effect ['//& 1155 trim(str(i))//'] file:['//trim(file)//'].') 1156 end if 1157 !on cherche l effet fixe correspondant 1158 do j=1,nfix 1159 if ( trim(namefix(j))==trim(name)) then 1160 exit 1161 end if 1162 end do 1163 1164 if ( j<=nfix) then 1165 niv(:,j)=niv_t(:) 1166 else 1167 call log_mess("QTLMap do not use fixed effect :"//trim(name),WARNING_DEF) 1168 end if 1169 1170 i=i+1 1171 end do 1172 1173 i = 1 1174 do while ( i <= ncov ) 1175 l = l + 1 1176 read(unit_perf,*,iostat=ios) name,(cov_t(j),j=1,nbete) 1177 if ( trim(name) == '') cycle 1178 if ( ios /= 0 ) then 1179 call log_mess('Problem detected at the line:'//trim(str(l)),ERROR_DEF) 1180 call stop_application('The application can not initialize the covariate ['//& 1181 trim(str(i))//'] file:['//trim(file)//'].') 1182 end if 1183 1184 !on cherche l effet fixe correspondant 1185 do j=1,ncov 1186 if ( trim(namecov(j))==trim(name)) then 1187 exit 1188 end if 1189 end do 1190 1191 if ( j<=ncov) then 1192 cov(:,j)=cov_t(:) 1193 else 1194 call log_mess("QTLMap do not use Covariate :"//trim(name),WARNING_DEF) 1195 end if 1196 1197 i=i+1 1198 end do 1199 1200 i = 1 1201 ii=0 1202 ios=0 1203 finloop=.false. 1204 do while ( i <= ncar .and. (.not. finloop) ) 1205 !count line 1206 l = l + 1 1207 call log_mess("setting value with line:"//trim(str(l)),DEBUG_DEF) 1208 1209 if ( futureNcar /= ncar ) then 1210 ii=ii+1 1211 do k=1,size(filter) 1212 if ( filter(k) == ii ) exit 1213 end do 1214 if ( k > size(filter )) then 1215 name="" 1216 do while ( ios == 0 .and. trim(name)=="") 1217 read(unit_perf,*,iostat=ios) name 1218 end do 1219 cycle 1220 else 1221 if ( k == size(filter)) finloop=.true. 1222 end if 1223 else 1224 k = i 1225 end if 1226 name="" 1227 ios=0 1228 do while ( ios == 0 .and. trim(name)=="") 1229 !reading line 1230 read(unit_perf,*,iostat=ios) name,(locvaleur(k,j),j=1,nbete) 1231 end do 1232 if ( ios /= 0 .and. i <= ncar ) then 1233 call stop_application("can not read the trait ["//trim(str(i))//"] [l:"//trim(str(l))//"]-"//trim(file)//".") 1234 end if 1235 !manage blank line 1236 if ( trim(name) == '') cycle 1237 !No enought line defined in the file 1238 if ( ios /= 0 ) then 1239 call log_mess('Problem detected at the line:'//trim(str(l)),ERROR_DEF) 1240 call stop_application('The application can not initialize the trait ['//& 1241 trim(str(i))//'] file:['//trim(file)//'].') 1242 end if 1243 1244 1245 !first colomn : name 1246 carac(k) = trim(name) 1247 1248 do j=1,nbete 1249 valeur(k,j) = get_real(locvaleur(k,j),is_ok) 1250 if ( .not. is_ok ) then 1251 call log_mess('** value unknown for trait ['//trim(carac(i))//& 1252 '] of animal ['//trim(bete(j))//']',WARNING_DEF) 1253 cdt(k,j) = 0 1254 else 1255 cdt(k,j) = 1 1256 end if 1257 end do 1258 1259 ndelt(k,:) = 1 1260 1261 i = i + 1 1262 end do 1263 1264 close(unit_perf) 1265 deallocate (locvaleur) 1266 deallocate(bete_char) 1267 deallocate (niv_t) 1268 deallocate (cov_t) 1269 1270 call log_mess('END SUBROUTINE read_perf_by_line',DEBUG_DEF) 1271 end subroutine
read_traits
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
read_traits
DESCRIPTION
Read the user phenotypic files and the model file
INPUTS
type_data_transcr : true => read the transcriptome format (animal defined by column), otherwise read as classical file (animal defined by line) filter : vector of index to select a subset of trait defined in the model file calculCd : infer the censored data
NOTES
SOURCE
284 subroutine read_traits(type_data_transcr,filter,calculCd) 285 logical,intent(in) :: type_data_transcr 286 integer, dimension(:),intent(in) :: filter 287 logical ,optional ,intent(in) :: calculCd 288 289 character(len=LENGTH_MAX_FILE), dimension(size(in_perfs)) :: traits_files 290 character(len=LENGTH_MAX_FILE) :: model_file 291 integer :: i,j 292 293 character(len=LEN_DEF) :: temp 294 traits_files = in_perfs 295 model_file = in_param_ef 296 call log_mess('SUBROUTINE read_traits',DEBUG_DEF) 297 298 if (type_data_transcr) then ! transcriptome information ... 299 call read_perf_by_line(traits_files,filter) 300 else 301 call read_perf_by_column(traits_files,filter) 302 end if 303 !************* NCAR ACCORDING TO FILTER ********* 304 if ( size(filter)/=0 ) ncar = size(filter) 305 306 if ( present(calculCd) ) then 307 if ( calculCd ) then 308 do i=1,size(bete) 309 call init_perf_animal(datasetUser,bete(i),valeur(:,i),cdt(:,i),ncar) 310 end do 311 end if 312 end if 313 314 call init_model_struct() 315 call initialise_struct_internal() 316 call check_traits_and_fathers() 317 call set_estime() 318 319 deallocate (ndelt) 320 deallocate (cdt) 321 deallocate (valeur) 322 deallocate (niv) 323 deallocate (cov) 324 deallocate (nuis) 325 deallocate (int_qtl) 326 327 call log_mess('END SUBROUTINE read_traits',DEBUG_DEF) 328 end subroutine read_traits
set_count_discrete
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
set_count_discrete
FUNCTION
Initialize structure ydiscretord , nmod and indicemod for discrete data
SYNOPSIS
Manage discrete data to initialized associated structure
INPUTS
RESULT
none.
EXAMPLE
NOTES
SEE ALSO
ydiscretord nmod indicemod
SOURCE
2456 subroutine set_count_discrete() 2457 2458 logical :: found 2459 integer :: i,ic,ideb,m,m1,j,temp 2460 integer , dimension(ncar,nd) :: indicemod_t 2461 2462 if (.not. allocated(ydiscretord)) allocate (ydiscretord(ncar,nd)) 2463 if (.not. allocated(nmod)) allocate (nmod(ncar)) 2464 ydiscretord = 0 2465 2466 !****************************************************************************** 2467 !****************************************************************************** 2468 ! 2469 ! transformation des donn�es discret � des donn�e utilisable par QTLMAP 2470 ! 2471 ! 2472 !***************************************************************************** 2473 !****************************************************************************** 2474 ! 2475 2476 2477 ! Parametres de maximisation 2478 2479 ! nmod est le nombre des modalit�s du caract�re �tudi� 2480 ! on reserve donc la partie du vecteur param�tre qui contiendra le lambda � estimer 2481 nmod=0 2482 2483 do ic=1,ncar 2484 if (natureY(ic) /= 'i') cycle 2485 nmod(ic)=1 2486 i=1 2487 found=.false. 2488 do while(i.le.nd.and..not.found) 2489 if(presentc(ic,i)) then 2490 found=.true. 2491 ideb=i+1 2492 indicemod_t(ic,1) = int(y(ic,i)) 2493 end if 2494 end do 2495 2496 do i=ideb,nd 2497 if(presentc(ic,i))then 2498 m=1 2499 found=.false. 2500 do while (m<=nmod(ic) .and. .not. found) 2501 if (int(y(ic,i))==indicemod_t(ic,m)) then 2502 found=.true. 2503 else 2504 m=m+1 2505 end if 2506 enddo 2507 2508 if (.not. found) then 2509 nmod(ic)=nmod(ic)+1 2510 indicemod_t(ic,nmod(ic))=y(ic,i) 2511 end if 2512 end if 2513 enddo 2514 2515 2516 ! maintenant on tri le tableau des modalit�s du plus petit au plus grand 2517 ! 2518 !************************************************************************* 2519 !**************************tri a bulle************************************ 2520 !************************************************************************* 2521 do j=1, nmod(ic) 2522 do m1=1, j-1 2523 if (indicemod_t(ic,j)<indicemod_t(ic,m1)) then 2524 temp=indicemod_t(ic,m1) 2525 indicemod_t(ic,m1)=indicemod_t(ic,j) 2526 indicemod_t(ic,j)=temp 2527 end if 2528 end do 2529 end do 2530 !************************************************************************** 2531 ! 2532 ! nmod est le nombre de modalit� possible pour le caract�re �tudi� 2533 ! 2534 ! maintenent il faut donc cr�er le nouveau tableau de performances qui sera 2535 !utilis� dans QTLMap, ydiscretord 2536 ! 2537 do i=1,nd 2538 do j=1, nmod(ic) 2539 if (presentc(ic,i)) then 2540 if ( int(y(ic,i))==indicemod_t(ic,j) )ydiscretord(ic,i)=j 2541 end if 2542 end do 2543 end do 2544 2545 end do ! ic 2546 2547 2548 2549 if (.not.allocated(indicemod)) allocate (indicemod(ncar,maxval(nmod))) 2550 indicemod = 0 2551 do ic=1,ncar 2552 do j=1,nmod(ic) 2553 indicemod(ic,j) = indicemod_t(ic,j) 2554 end do 2555 end do 2556 2557 end subroutine set_count_discrete
set_estime
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
set_estime
DESCRIPTION
Initialize the following data :
* estime (dim : ncar,nm) : indicate the estimability of a full sib family (condition ndmin involve full sib family are used in the qtl analysis)
* nmumest (dim ncar) : number of full sib family with enough progenies (ndmin>=)
* namest (dim ncar) : number of unique female (several full sib family are building with the same dam)
* iam (dim : ncar,nfem) : get the index of female (estimable)
NOTES
anciennement dans optinit, on initialisse estime,nmumest,namest,iam dans cette routine
SOURCE
2280 subroutine set_estime() 2281 logical :: estfem(nfem) 2282 integer :: ic,ip,jm,kd,eff,ifem 2283 2284 call log_mess('SUBROUTINE set_estime',DEBUG_DEF) 2285 2286 if ( allocated (estime) ) then 2287 deallocate( estime ) 2288 deallocate( namest ) 2289 deallocate( nmumest ) 2290 end if 2291 2292 if ( allocated (iam) ) deallocate( iam ) 2293 2294 allocate (estime(size(carac),size(mere)) ) 2295 allocate (namest(size(carac)) ) 2296 allocate (nmumest(size(carac))) 2297 allocate (iam(size(carac),size(femelle)) ) 2298 2299 estime=.false. 2300 nmumest=0 2301 namest=0 2302 iam=0.d0 2303 2304 do ic=1,ncar 2305 estfem=.false. 2306 do ip=1,np 2307 do jm=nmp(ip)+1,nmp(ip+1) 2308 eff=0 2309 do kd=ndm(jm)+1,ndm(jm+1) 2310 if(presentc(ic,kd)) eff=eff+1 2311 end do 2312 ifem=repfem(jm) 2313 if(eff.ge.ndmin) then 2314 estime(ic,jm)=.true. 2315 estfem(ifem)=.true. 2316 nmumest(ic)=nmumest(ic)+1 2317 end if 2318 end do 2319 end do 2320 do ifem=1,nfem 2321 iam(ic,ifem)=0 2322 if(estfem(ifem)) then 2323 namest(ic)=namest(ic)+1 2324 iam(ic,ifem)=namest(ic) 2325 end if 2326 end do 2327 end do 2328 2329 call log_mess('END SUBROUTINE set_estime',DEBUG_DEF) 2330 2331 end subroutine set_estime
set_proportion_discrete
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
set_proportion_discrete
FUNCTION
Initialize structure Prop (proportion) and Seuil (Threshold) for discrete data
SYNOPSIS
Mange process according to the type of the data (continue, discrete or categorial)
INPUTS
RESULT
none.
EXAMPLE
NOTES
SEE ALSO
prop seuil MATH_QTLMAP_G01FAF
SOURCE
2583 subroutine set_proportion_discrete 2584 2585 integer :: ic,i,ifail,ii 2586 real(kind=dp) :: cumul,eff 2587 2588 ifail = 0 2589 2590 ! 2591 ! comptage pour la cr�ation du point de d�part 2592 ! 2593 if(.not. allocated(seuil)) allocate (seuil(ncar,maxval(nmod))) 2594 if (.not. allocated(prop)) allocate (prop(ncar,maxval(nmod))) 2595 2596 prop=0.d0 2597 seuil = 0.d0 2598 2599 do ic=1,ncar 2600 if (natureY(ic) /= 'i') cycle 2601 eff=0.d0 2602 do ii=1,nd 2603 if(presentc(ic,ii)) then 2604 prop(ic,ydiscretord(ic,ii))=prop(ic,ydiscretord(ic,ii))+1.d0 2605 eff=eff+1.d0 2606 end if 2607 end do 2608 do ii=1,nmod(ic) 2609 prop(ic,ii) = prop(ic,ii)/eff 2610 end do 2611 2612 cumul=0.d0 2613 do i=1,nmod(ic)-1 2614 cumul=cumul+prop(ic,i) 2615 seuil(ic,i)=MATH_QTLMAP_G01FAF('L',cumul,ifail) 2616 end do 2617 ! 2618 ! fin du comptage 2619 end do ! ic 2620 end subroutine set_proportion_discrete
sim_perf_shuffling
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
sim_perf_shuffling
DESCRIPTION
NOTES
realise un suffling intrafamille:
- de pere et mere lorsque le couple a plus de ndmin descendants
!-de pere, pour tous les descendants des meres ayant moins de ndmin
!descendants avec ce pere (l'ensemble etant > a ndmin)
SOURCE
2013 subroutine sim_perf_shuffling(analyse_is_multi_traits) 2014 2015 logical :: analyse_is_multi_traits 2016 real (kind=dp) ,dimension(:,:),allocatable :: zy,zcd, zcovar 2017 integer ,dimension(:,:),allocatable :: zndelta,znivx 2018 integer ,dimension(:),allocatable :: effpr,iv 2019 logical ,dimension(:),allocatable :: permut_pere 2020 logical ,dimension(:,:),allocatable :: permut_pere_mere, zpresentc,zpresentg 2021 integer ,dimension(:,:),allocatable :: effpm 2022 integer :: ifail,ip,jm,kd,nm1,nm2,nd1,nd2,i,j,id,nc, npresp 2023 integer :: alloc_stat 2024 2025 ALLOCATE (zy(size(carac),nd), stat = alloc_stat) 2026 CALL check_allocate(alloc_stat,'zy') 2027 ALLOCATE (znivx(nd,nfix), stat = alloc_stat) 2028 CALL check_allocate(alloc_stat,'znivx') 2029 ALLOCATE (zcovar(nd,ncov), stat = alloc_stat) 2030 CALL check_allocate(alloc_stat,'zcovar') 2031 ALLOCATE (zndelta(size(carac),nd), stat = alloc_stat) 2032 CALL check_allocate(alloc_stat,'zndelta') 2033 ALLOCATE (zcd(size(carac),nd), stat = alloc_stat) 2034 CALL check_allocate(alloc_stat,'zcd') 2035 ALLOCATE (zpresentc(ncar,nd), stat = alloc_stat) 2036 CALL check_allocate(alloc_stat,'zpresentc') 2037 ALLOCATE (permut_pere(np), stat = alloc_stat) 2038 CALL check_allocate(alloc_stat,'permut_pere') 2039 ALLOCATE (zpresentg(nchr,nd), stat = alloc_stat) 2040 CALL check_allocate(alloc_stat,'zpresentg') 2041 ALLOCATE (permut_pere_mere(np,nm), stat = alloc_stat) 2042 CALL check_allocate(alloc_stat,'permut_pere_mere') 2043 ALLOCATE (effpr(np), stat = alloc_stat) 2044 CALL check_allocate(alloc_stat,'effpr') 2045 ALLOCATE (iv(nd), stat = alloc_stat) 2046 CALL check_allocate(alloc_stat,'iv') 2047 2048 ALLOCATE (effpm(np,nm), stat = alloc_stat) 2049 CALL check_allocate(alloc_stat,'effpm') 2050 2051 call log_mess('simulation of traits by permutation (multitrait mode)...',VERBOSE_DEF) 2052 !***************************************************************************** 2053 !1�boucle sur les peres mere descendants 2054 ! --> comptage du nombre de descendants par couple pere-mere: effpm(ip,jm) 2055 ! --> comptage du nombre minimum de descendants d'un pere parmi 2056 ! les couple ayant plus de nb_des_min descendants: si ce nbre est inf�rieur � ndim 2057 ! alors ont donne un message d'erreur � l'execution 2058 2059 2060 zy=0.d0 2061 zndelta=0 2062 do ip=1,np 2063 permut_pere(ip)=.false. 2064 effpr(ip)=0 2065 nm1=nmp(ip)+1 2066 nm2=nmp(ip+1) 2067 do jm=nm1,nm2 2068 effpm(ip,jm)=0 2069 nd1=ndm(jm)+1 2070 nd2=ndm(jm+1) 2071 permut_pere_mere(ip,jm)=.true. 2072 do kd=nd1,nd2 2073 zy(:,kd)=y(:,kd) 2074 zndelta(:,kd)= ndelta(:,kd) 2075 zcd(:,kd)=cd(:,kd) 2076 zpresentc(:,kd)=presentc(:,kd) 2077 zpresentg(:,kd)=presentg(:,kd) 2078 znivx(kd,:)=nivx(kd,:) 2079 zcovar(kd,:)=covar(kd,:) 2080 if ((analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd))==size(carac))) & 2081 .or.(.not.analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd)).ne.0))) then 2082 effpm(ip,jm)=effpm(ip,jm)+1 2083 effpr(ip)=effpr(ip)+1 2084 endif 2085 end do 2086 ! si une famille de pere-mere est trop petite les permutations se font en utilisant les performances de toute la famille de ce pere 2087 if (effpm(ip,jm)<ndmin.or.effpm(ip,jm)<nb_des_min) then 2088 permut_pere(ip)=.true. 2089 permut_pere_mere(ip,jm)=.false. 2090 ! impression d'un warning lorsque l'utilisateur a permut� dans la famille de p�re au lieu de la famille pere mere 2091 call log_mess ('WARNING: permutation will be performed within sire family for sire ['//trim(pere(ip))//& 2092 '], and dam ['//trim(mere(jm))//']',VERBOSE_DEF) 2093 endif 2094 end do 2095 2096 ! impression d'un warning lorsque des famille de pere sont trop petites pour �tre permut�es 2097 if (effpr(ip).LT.nb_des_min) then 2098 call log_mess ('family size of sire ['//trim(pere(ip))//'] is too small (['//& 2099 str(effpr(ip))//'] animals) to use CONFIDENTLY permutation method, try the simulation method ',WARNING_DEF) 2100 endif 2101 !print *, 'EFFECTIF PERMUTE',effpr(ip) 2102 end do 2103 ! 2104 2105 !***************************************************************************** 2106 !PERMUTATION DES PERFORMANCES 2107 !***************************************************************************** 2108 iv=0 2109 do ip=1, np 2110 npresp=0 2111 nm1=nmp(ip)+1 2112 nm2=nmp(ip+1) 2113 !***************************************************************************** 2114 ! 1--> r�attribution des performances par permutation au hasard intrafamille de pere( issus 2115 ! de MERE AVEC - de nb_des_min DESCENDANTS CHACUNE ) 2116 IF (permut_pere(ip)) then 2117 do jm=nm1,nm2 2118 nd1=ndm(jm)+1 2119 nd2=ndm(jm+1) 2120 do kd=nd1,nd2 2121 if ((analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd))==size(carac))) & 2122 .or.(.not.analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd)).ne.0))) then 2123 npresp=npresp+1 2124 iv(npresp)=kd 2125 endif 2126 end do 2127 end do 2128 ifail=1 2129 !permutation au hasard dans iv 2130 call MATH_QTLMAP_G05EHF(iv,npresp,ifail) 2131 !permutation du vecteur de performances complet (sans donn�es manquantes) lorsqu'on fait une analyse multicaract�re 2132 !permutation du vecteur de performances incluant les donn�es manquantes lorsqu'on fait une analyse unicaract�re 2133 !(seuls les animaux avec aucune performance sont �cart�s 2134 id=0 2135 do jm=nm1,nm2 2136 nd1=ndm(jm)+1 2137 nd2=ndm(jm+1) 2138 do kd=nd1,nd2 2139 if (((analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd))==size(carac))) & 2140 .or.(.not.analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd)).ne.0))) & 2141 .and..not.permut_pere_mere(ip,jm)) then 2142 id=id+1 2143 ! print *, 'kd', kd, 'Y', y(:,kd) 2144 y(:,kd)=zy(:,iv(id)) 2145 ! print *, 'kd', kd, 'Y', y(:,kd) 2146 ndelta(:,kd)=zndelta(:,iv(id)) 2147 cd(:,kd)=zcd(:,iv(id)) 2148 nivx(kd,:)=znivx(iv(id),:) 2149 covar(kd,:)=zcovar(iv(id),:) 2150 presentc(:,kd)=zpresentc(:,iv(id)) 2151 presentg(:,kd)=zpresentg(:,iv(id)) 2152 endif 2153 end do 2154 end do 2155 ! print*, 'EFFECTIF FINAL', id 2156 ENDIF 2157 !***************************************************************************** 2158 ! 2--> permutation par famille de pere-mere 2159 do jm=nm1,nm2 2160 npresp=0 2161 id=0 2162 nd1=ndm(jm)+1 2163 nd2=ndm(jm+1) 2164 IF (permut_pere_mere(ip,jm)) THEN 2165 do kd=nd1,nd2 2166 if ((analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd))==size(carac))) & 2167 .or.(.not.analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd)).ne.0))) then 2168 npresp=npresp+1 2169 iv(npresp)=kd 2170 endif 2171 end do 2172 ifail=1 2173 !permutation au hasard dans iv 2174 call MATH_QTLMAP_G05EHF(iv,npresp,ifail) 2175 !permutation du vecteur de performances complet (sans donn�es manquantes) lorsqu'on fait une analyse multicaract�re 2176 !permutation du vecteur de performances incluant les donn�es manquantes lorsqu'on fait une analyse unicaract�re 2177 !(seuls les animaux avec aucune performance sont �cart�s 2178 do kd=nd1,nd2 2179 if ((analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd))==size(carac))) & 2180 .or.(.not.analyse_is_multi_traits.and.(count(zpresentg(:,kd))>0).and.(count(zpresentc(:,kd)).ne.0))) then 2181 id=id+1 2182 y(:,kd)=zy(:,iv(id)) 2183 ndelta(:,kd)=zndelta(:,iv(id)) 2184 cd(:,kd)=zcd(:,iv(id)) 2185 nivx(kd,:)=znivx(iv(id),:) 2186 covar(kd,:)=zcovar(iv(id),:) 2187 presentc(:,kd)=zpresentc(:,iv(id)) 2188 presentg(:,kd)=zpresentg(:,iv(id)) 2189 endif 2190 end do 2191 ENDIF 2192 end do 2193 !***************************************************************************** 2194 end do 2195 !***************************************************************************** 2196 ! IMPRESSION DES DONNEES PERMUTEES 2197 !***************************************************************************** 2198 ! do ic=1,ncar 2199 ! call log_mess('Trait ['//CHAR(carac(ic))//'] qtlmap normalize data with means:'//& 2200 ! str(xmut(ic))//'+-'//str(sigt(ic)),INFO_DEF) 2201 ! enddo 2202 ! do ip=1,np 2203 ! nm1=nmp(ip)+1 2204 !nm2=nmp(ip+1) 2205 !do jm=nm1,nm2 2206 ! nd1=ndm(jm)+1 2207 ! nd2=ndm(jm+1) 2208 ! print *, 'permute pere',permut_pere(ip) 2209 ! print *, 'permute pere-mere',permut_pere_mere(ip,jm), effpm(ip,jm), ndmin,nb_des_min 2210 ! do kd=nd1,nd2 2211 ! print *, trim(animal(kd)), ' ', trim(pere(ip)),' ', trim(mere(jm)),& 2212 ! ' data ',(trim(carac(ic)),' ',zy(ic,kd), zpresentc(ic,kd),ic=1,ncar),(znivx(kd,i), i=1,nfix), (zcovar(kd, j), j=1, ncov), & 2213 ! ' sim1 ', (trim(carac(ic)),' ',y(ic,kd), presentc(ic,kd),ic=1,ncar), (nivx(kd,i), i=1,nfix), (covar(kd, j), j=1, ncov) 2214 ! enddo 2215 !enddo 2216 ! enddo 2217 DEALLOCATE (zy) 2218 DEALLOCATE (znivx) 2219 DEALLOCATE (zcovar) 2220 DEALLOCATE (zndelta) 2221 DEALLOCATE (zpresentc) 2222 DEALLOCATE (effpr) 2223 DEALLOCATE (iv) 2224 DEALLOCATE (permut_pere) 2225 DEALLOCATE (permut_pere_mere) 2226 DEALLOCATE (zcd) 2227 DEALLOCATE (effpm) 2228 return 2229 2230 end subroutine sim_perf_shuffling
sim_perf_tirage
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
sim_perf_tirage
DESCRIPTION
NOTES
Tirage des performances pour deux caracteres correles , avec ou sans qtl
SOURCE
1753 subroutine sim_perf_tirage(ncar,rho,h2,ys) 1754 integer , intent(in) :: ncar 1755 !heritabilities of traits , dim : ncar 1756 real (kind=dp),dimension (ncar), intent(in) :: h2 1757 ! correlation matrix , dim : ncar,ncar 1758 real (kind=dp),dimension (ncar,ncar) ,intent(in):: rho 1759 !dim : ncar 1760 real ,dimension(ncar) :: u0 1761 real (kind=dp) ,dimension(:),allocatable :: upv,umv,udv,ugpv,ugmv 1762 1763 !dim ncar,ncar 1764 real (kind=dp) ,dimension(:,:),allocatable :: varcovg,varcovgd 1765 real (kind=dp) ,dimension(:,:),allocatable :: varcove 1766 !dim ncar*ncar 1767 real ,dimension(ncar,ncar) :: ccovvar 1768 !dim : ncar,nr 1769 real (kind=dp) ,dimension(:,:),allocatable :: urv 1770 !dim : nd 1771 !integer ,dimension(:),allocatable :: itest 1772 1773 !dim : ncar,nd 1774 real (kind=dp) ,dimension(ncar,size(animal)),intent(out) :: ys 1775 ! 1776 ! dim : nrsx=((ncar+1)*(ncar+2))/2) 1777 real ,dimension(:),allocatable :: zr,zd,zh,zg 1778 1779 real ,dimension(ncar) :: WORK 1780 1781 ! dim : ncar*(ncar+3)/2 + 1 1782 real ,dimension(:),allocatable :: refgg,refg,refe 1783 1784 integer :: i,j,ngm1,ngm2,igm,nr1,nr2,ir,nm2 1785 integer :: jm,nd1,nd2,kd,ifail,igp,ic,ip,nm1 1786 real (kind=dp) :: vare,vargd,varg,varej 1787 character(len=LEN_DEF) :: tempvs 1788 1789 call log_mess('simulation of traits...',VERBOSE_DEF) 1790 1791 if (size(rho) <= 0 ) then 1792 call stop_application("DEVEL ERROR: sim_perf_tirage RHO is not initialized") 1793 end if 1794 1795 if (size(h2) <= 0 ) then 1796 call stop_application("DEVEL ERROR: sim_perf_tirage H2 is not initialized") 1797 end if 1798 1799 allocate (varcove(ncar,ncar)) 1800 allocate (varcovg(ncar,ncar)) 1801 allocate (varcovgd(ncar,ncar)) 1802 allocate (upv(ncar)) 1803 allocate (umv(ncar)) 1804 allocate (udv(ncar)) 1805 allocate (ugpv(ncar)) 1806 allocate (ugmv(ncar)) 1807 allocate (urv(ncar,nr)) 1808 1809 1810 allocate (zr( ((ncar+1)*(ncar+2))/2) ) 1811 allocate (zd( ((ncar+1)*(ncar+2))/2) ) 1812 allocate (zh( ((ncar+1)*(ncar+2))/2) ) 1813 allocate (zg( ((ncar+1)*(ncar+2))/2) ) 1814 1815 allocate (refgg(ncar*(ncar+3)/2 + 1)) 1816 allocate (refg(ncar*(ncar+3)/2 + 1)) 1817 allocate (refe(ncar*(ncar+3)/2 + 1)) 1818 1819 1820 1821 ! call log_mess("Simulation trait only for real value",WARNING_DEF); 1822 ! natureY='r' 1823 ! 1824 !***************************************************************************** 1825 ! 1826 ! Initialisation des matrices de variance covariance residuelles et polygen 1827 ! 1828 do i=1,ncar 1829 u0(i)=0.d0 1830 varg = sqrt(h2(i)) 1831 vargd =sqrt(0.5d0*h2(i)) 1832 vare=sqrt(1.d0-h2(i)) 1833 varcovg(i,i)=varg*varg 1834 varcovgd(i,i)=vargd*vargd 1835 varcove(i,i)=vare*vare 1836 do j=i+1,ncar 1837 varej = sqrt(1.d0-h2(j)) 1838 ! varcovg(i,j)=varg(i)*varg(j) 1839 varcovg(i,j)=0.0d0 1840 varcovg(j,i)=varcovg(i,j) 1841 ! varcovgd(i,j)=vargd(i)*vargd(j) 1842 varcovgd(i,j)=0.0d0 1843 varcovgd(j,i)=varcovgd(i,j) 1844 varcove(i,j)=rho(i,j)*vare*varej 1845 varcove(j,i)=varcove(i,j) 1846 end do 1847 end do 1848 ! 12 format(1x,f10.5,3x,10(f10.5,1x)) 1849 ! 1850 ! 1851 !***************************************************************************** 1852 ! Initialisation des vecteurs de tirages des performances dans une binormale 1853 ! 1854 do i=1,ncar 1855 do j=1,ncar 1856 ccovvar(i,j) = varcovg(i,j) 1857 end do 1858 end do 1859 1860 call setgmn(u0,ccovvar,ncar,refgg) 1861 ! 1862 do i=1,ncar 1863 do j=1,ncar 1864 ccovvar(i,j) = varcovgd(i,j) 1865 end do 1866 end do 1867 1868 call setgmn(u0,ccovvar,ncar,refg) 1869 ! 1870 do i=1,ncar 1871 do j=1,ncar 1872 ccovvar(i,j) = varcove(i,j) 1873 end do 1874 end do 1875 1876 1877 call setgmn(u0,ccovvar,ncar,refe) 1878 ! 1879 ! 1880 !***************************************************************************** 1881 ! Tirage des valeurs g�n�tiques des gparents et des reproducteurs 1882 ! 1883 ifail=0 1884 ! Grd peres 1885 do igp=1,ngp 1886 call genmn(refgg,zg,WORK) 1887 do ic=1,ncar 1888 ugpv(ic)=0.5d0*dble(zg(ic)) 1889 end do 1890 ! 1891 ! Grd meres 1892 ngm1=ngmgp(igp)+1 1893 ngm2=ngmgp(igp+1) 1894 do igm=ngm1,ngm2 1895 call genmn(refgg,zg,WORK) 1896 do ic=1,ncar 1897 ugmv(ic)=0.5d0*dble(zg(ic)) 1898 end do 1899 ! 1900 ! Reproducteurs 1901 nr1=nrgm(igm)+1 1902 nr2=nrgm(igm+1) 1903 1904 do ir=nr1,nr2 1905 call genmn(refg,zr,WORK) 1906 do ic=1,ncar 1907 urv(ic,ir)=dble(zr(ic))+ugpv(ic)+ugmv(ic) 1908 end do 1909 end do 1910 end do 1911 end do 1912 !***************************************************************************** 1913 ! Tirage des performances des descendants 1914 ! 1915 do ip=1,np 1916 if (reppere(ip).eq.9999) call genmn(refgg,zg,WORK) 1917 do ic=1,ncar 1918 if (reppere(ip).eq.9999) then 1919 upv(ic)=0.5d0*dble(zg(ic)) 1920 else 1921 upv(ic)=0.5d0*urv(ic,reppere(ip)) 1922 end if 1923 end do 1924 nm1=nmp(ip)+1 1925 nm2=nmp(ip+1) 1926 do jm=nm1,nm2 1927 if (repmere(jm).eq.9999) call genmn(refgg,zg,WORK) 1928 do ic=1,ncar 1929 if (repmere(jm).eq.9999) then 1930 umv(ic)=0.5d0*dble(zg(ic)) 1931 else 1932 umv(ic)=0.5d0*urv(ic,repmere(jm)) 1933 end if 1934 end do 1935 nd1=ndm(jm)+1 1936 nd2=ndm(jm+1) 1937 do kd=nd1,nd2 1938 call genmn(refg,zd,WORK) 1939 do ic=1,ncar 1940 udv(ic)=upv(ic)+umv(ic)+dble(zd(ic)) 1941 end do 1942 ! 1943 ! Residuelle 1944 call genmn(refe,zh,WORK) 1945 1946 do ic=1,ncar 1947 ys(ic,kd)=udv(ic)+dble(zh(ic)) 1948 end do 1949 end do 1950 end do 1951 end do 1952 1953 deallocate (urv) 1954 deallocate (varcove) 1955 deallocate (varcovg) 1956 deallocate (varcovgd) 1957 deallocate (upv) 1958 deallocate (umv) 1959 deallocate (udv) 1960 deallocate (ugpv) 1961 deallocate (ugmv) 1962 1963 deallocate (zr) 1964 deallocate (zd) 1965 deallocate (zh) 1966 deallocate (zg) 1967 1968 deallocate (refgg) 1969 deallocate (refg) 1970 deallocate (refe) 1971 1972 end subroutine sim_perf_tirage
sim_QTL
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
sim_QTL
DESCRIPTION
NOTES
Simulation des typages sur trois generations , reperage des alleles au qtl des descendants.
SOURCE
1461 subroutine sim_QTL(icar,nqtl,pas,qtl,croisement) 1462 1463 integer ,intent(in) :: icar,nqtl,pas 1464 ! 1465 ! Tableaux dimensionnes selon n, le nombre d'individus dans le pedigree 1466 integer ,intent(out) :: qtl(np+nm+nd,nqtl,2) 1467 character(len=LEN_BUFFER_WORD) , intent(in) :: croisement 1468 ! 1469 !C Divers 1470 real :: x_rand 1471 integer :: geno,iq,ip,ifem,ipos,nm1,nm2,jm,ngeno1,ngeno2,nd1,nd2,chr 1472 integer :: kd,kkd,iph,n 1473 real (kind=dp) :: pp,pm 1474 real,external :: ranf 1475 logical :: ok 1476 1477 ! 1478 ! 1479 !*********************************************************************** 1480 ! Simulation des alleles au QTL des parents 1481 ! 1482 !*********************************************************************** 1483 do iq=1,nqtl 1484 do ip=1,np 1485 ! xlim(iq)=frequence de allele Q2 chez les GP et de Q1 chez les GM 1486 ! si xlim(iq)=1, tous GP Q1Q1 et toutes GM Q2Q2 1487 qtl(ip,iq,1)=2 1488 qtl(ip,iq,2)=1 1489 x_rand=ranf() 1490 if (dble(x_rand).lt.xlim(iq)) qtl(ip,iq,1)=1 1491 x_rand=ranf() 1492 if (dble(x_rand).lt.xlim(iq)) qtl(ip,iq,2)=2 1493 end do 1494 ! 1495 do ifem=1,nfem 1496 ! xlim(iq)=frequence de allele Q2 chez les GP et de Q1 chez les GM 1497 ! si xlim(iq)=1, tous GP Q1Q1 et toutes GM Q2Q2 1498 qtl(np+ifem,iq,1)=2 1499 qtl(np+ifem,iq,2)=1 1500 if (croisement /= BC_KEYWORD) then 1501 x_rand=ranf() 1502 if (dble(x_rand).gt.xlim(iq)) qtl(np+ifem,iq,1)=1 1503 x_rand=ranf() 1504 if (dble(x_rand).gt.xlim(iq)) qtl(np+ifem,iq,2)=2 1505 else 1506 qtl(np+ifem,iq,2)=2 1507 end if 1508 1509 end do 1510 end do 1511 1512 1513 ! 1514 ! 1515 !********************************************************************** 1516 ! Simulation de la transmission des alleles QTL des parents 1517 ! aux descendants : les estimations de phase des meres sont correctes 1518 ! la phase la plus probable est systematiquement transmise 1519 ! ==> simulation de la deviation a la moyenne nulle 1520 !********************************************************************** 1521 ! 1522 do iq=1,nqtl 1523 chr=0 1524 ok=.false. 1525 do while ( .not. ok ) 1526 chr = chr + 1 1527 if ( chr > nchr ) call stop_application("none chromosome "//trim(chrqtl(iq))//" are defined !"); 1528 ok = chrqtl(iq) == chromo(chr) 1529 end do 1530 1531 n=get_pos((posiqtl(iq)-posi(chr,1))) 1532 if ( n > get_ilong(chr)) n = get_ilong(chr) 1533 if ( n <= 0 ) n = 1 1534 do ip=1,np 1535 nm1=nmp(ip)+1 1536 nm2=nmp(ip+1) 1537 do jm=nm1,nm2 1538 ifem=repfem(jm) 1539 ngeno1=ngenom(chr,jm)+1 1540 ngeno2=ngenom(chr,jm+1) 1541 geno=ngeno1 1542 nd1=ngend(chr,geno)+1 1543 nd2=ngend(chr,geno+1) 1544 do kd=nd1,nd2 1545 kkd=ndesc(chr,kd) 1546 ! Transmission de l'allele paternel: a priori re�oit allele 1 (grand p�re) du p�re 1547 ! si pt� transmission pp non d�pass�e, re�oit all�le 2 1548 qtl(np+nfem+kkd,iq,1)=qtl(ip,iq,1) 1549 ! pp=(-pdd(chr,kd,1,n)-pdd(chr,kd,2,n)+pdd(chr,kd,3,n)+pdd(chr,kd,4,n))/2.d0+0.5d0 1550 ! MODIF OFI : proba qu'on est recu le 2eme allele du pere 1551 pp=pdd(chr,kd,3,n)+pdd(chr,kd,4,n) 1552 x_rand=ranf() 1553 if(dble(x_rand).lt.pp)qtl(np+nfem+kkd,iq,1)=qtl(ip,iq,2) 1554 ! 1555 ! Transmission de l'allele maternel: a priori re�oit allele 1 (grand p�re) de la m�re 1556 qtl(np+nfem+kkd,iq,2)=qtl(np+ifem,iq,1) 1557 ! pm=(-pdd(chr,kd,1,n)+pdd(chr,kd,2,n)-pdd(chr,kd,3,n)+pdd(chr,kd,4,n))/2.d0+0.5d0 1558 ! MODIF OFI : proba qu'on est recu le 2eme allele de la mere 1559 pm=pdd(chr,kd,2,n)+pdd(chr,kd,4,n) 1560 x_rand=ranf() 1561 if(dble(x_rand).lt.pm)qtl(np+nfem+kkd,iq,2)=qtl(np+ifem,iq,2) 1562 ! 1563 ! Mise � jour des perf simul�es sous H0 1564 1565 if (presentc(icar,kkd)) then 1566 if (qtl(np+nfem+kkd,iq,1).eq.qtl(np+nfem+kkd,iq,2)) then 1567 if (qtl(np+nfem+kkd,iq,1).eq.1) then 1568 y(icar,kkd)=y(icar,kkd)+2*ue(icar,iq) 1569 else 1570 y(icar,kkd)=y(icar,kkd)-2*ue(icar,iq) 1571 end if 1572 end if 1573 end if 1574 end do 1575 end do 1576 end do 1577 end do 1578 ! call log_mess('DEPLACER LA GENERATION DU FICHIER PERFOSIM1 SIMULER...',WARNING_DEF) 1579 1580 return 1581 end subroutine sim_QTL
sim_transform_discret
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
sim_transform_discret
DESCRIPTION
NOTES
dans cette subroutine, on genere les donnees discretes (a mettre dans le tableau y) a partir des valeurs de la sous jacente (trouvees dans le tableau y)
SOURCE
1984 subroutine sim_transform_discret(ic) 1985 integer , intent(in) :: ic 1986 integer :: kd,m 1987 1988 do kd=1,size(y,2) 1989 if(presentc(ic,kd))then 1990 m=1 1991 do while (m < nbmodsimultrait(ic) .and. y(ic,kd) > soglia(ic,m)) 1992 m=m+1 1993 enddo 1994 y(ic,kd)=m 1995 end if 1996 end do 1997 1998 end subroutine sim_transform_discret
write_perf
[ Top ] [ m_qtlmap_traits ] [ Subroutines ]
NAME
write_perf
DESCRIPTION
NOTES
SOURCE
2242 subroutine write_perf(file_name) 2243 character(len=*),intent(in) :: file_name 2244 integer :: kd,ic 2245 integer :: myunit=99999 2246 integer ,dimension(ncar) :: cd 2247 open (myunit,file=file_name) 2248 do kd=1,nd 2249 cd = 1 2250 do ic=1,ncar 2251 if(.not. presentc(ic,kd)) then 2252 y(ic,kd)=-99.d0 2253 cd(ic) = 0 2254 end if 2255 end do 2256 2257 write(myunit,FMT='(a12,'// trim(str(ncar)) //'(1x,f9.5,1x,i1,1x,"1"))')& 2258 trim(animal(kd)),(y(ic,kd),cd(ic),ic=1,ncar) 2259 end do 2260 close(myunit) 2261 2262 end subroutine write_perf