m_qtlmap_analyse_gen
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_gen
DESCRIPTION
NOTES
SEE ALSO
cary
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
cary
DESCRIPTION
sum of square of the trait value in the full sib family jm
DIMENSIONS
nm
NOTES
carydf
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
carydf
DESCRIPTION
sum of square of the trait value in the half sib family ip (only with dam estimable)
DIMENSIONS
np
NOTES
eff
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
eff
DESCRIPTION
number of progenies with a trait value in the full sib family jm
DIMENSIONS
nm
NOTES
effdf
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
effdf
DESCRIPTION
number of progenies (only with dam estimable) with a trait value in the half sib family ip
DIMENSIONS
np
NOTES
effp
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
effp
DESCRIPTION
number of progenies with a trait value in the half sib family ip
DIMENSIONS
np
NOTES
estmum
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
estmum
DESCRIPTION
number of female estimable in a half sib family ip
DIMENSIONS
np
NOTES
sig0
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
sig0
DESCRIPTION
variance of traits values for the half sib family ip
DIMENSIONS
np
NOTES
somcd
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
somcd
DESCRIPTION
sum of the censured data in the full sib family jm
DIMENSIONS
nm
NOTES
somcddf
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
somcddf
DESCRIPTION
sum of the censured data in the half sib family ip (only with dam estimable)
DIMENSIONS
np
NOTES
somy
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
somy
DESCRIPTION
sum of the trait value in the full sib family jm
DIMENSIONS
nm
NOTES
somydf
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
somydf
DESCRIPTION
sum of the trait value in the half sib family ip (only with dam estimable)
DIMENSIONS
np
NOTES
xmu0m
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
xmu0p
DESCRIPTION
mean of traits values for the full sib family jm
DIMENSIONS
nm
NOTES
xmu0p
[ Top ] [ m_qtlmap_analyse_gen ] [ Variables ]
NAME
xmu0p
DESCRIPTION
mean of traits values for the half sib family ip
DIMENSIONS
np
NOTES
courbe_lin_ldla
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
courbe_lin_ldla
DESCRIPTION
Print the information about result (curve, maximum finded, solution of the estimation)
INPUTS
ic : index of the trait
chr : index of the chromosome
est_moy : True if the mean is given
est_var : True if the variance is given
xlrmax : LRT maximum
nbfem : number of female estimable
nbniv : number of level fixed effect
nbef : number of fixed effect
nbco : number of covariate
ntlevp : number of level for qtl sire interaction with fixed effect
ntlevm : number of level for qtl dam interaction with fixed effect
par0 : solution (bestim) under H0
par1 : solution (bestim) under H1
ntniv : number of level in the incidence matrix
vecsol0 : logical array : estimability of level in the incidence matrix under H0
vecsol1 : logical array : estimability of level in the incidence matrix under H1
precis0 : precision information about each level under H0
precis1 : precision information about each level under H1
ordo : LRT curve under H1
npo : number of position tested
nposx : the maximum position finded
printRiskFactor : print cox case with Risk factor
qtl : qtl effect are estimated
hsire : haplotype effect are estimated
hdam : haplotype effect with consideration of female estimable are estimated
NOTES
En fonction du profil de stat de test trace la courbe correspondante Sous programme appele par analyse Version 1 Date 1er Mars 2010. JM Elsen
SOURCE
643 subroutine courbe_lin_ldla(chr,est_moy,est_var,xlrmax,ic, & 644 nbfem,nbniv,nbef,nbco,ntlevp,ntlevm,par0,par1,& 645 ntniv,vecsol0,precis0,vecsol1,precis1, & 646 ordo,npo,nposx,printRiskFactor,qtl,hsire,hdam) 647 logical , intent(in) :: est_moy,est_var 648 integer ,intent(in) :: ic 649 real (kind=dp) ,intent(in) :: xlrmax 650 integer ,intent(in) :: chr,nbfem 651 integer ,intent(in) :: nbniv 652 integer ,intent(in) :: nbef 653 integer ,intent(in) :: nbco 654 integer ,intent(in) :: ntlevp 655 integer ,intent(in) :: ntlevm 656 real (kind=dp),dimension(:),intent(in) :: par0 657 real (kind=dp),dimension(:),intent(in) :: par1 658 integer ,intent(in) :: ntniv 659 real (kind=dp),dimension(:),intent(in) :: precis0 660 logical ,dimension(:),intent(in) :: vecsol0 661 real (kind=dp),dimension(:),intent(in) :: precis1 662 logical ,dimension(:),intent(in) :: vecsol1 663 integer ,intent(in) :: npo,nposx 664 real (kind=dp) ,intent(in) ,dimension(npo) :: ordo 665 logical , intent(in) , optional :: printRiskFactor 666 real (kind=dp) :: par_refm, par_refp 667 real (kind=dp), dimension(nbef) :: par_refef 668 logical,intent(in) :: qtl,hsire,hdam 669 670 ! ************* deb chgt ************* 671 integer :: n_est_par,i_haplo,lll 672 ! ************* fin chgt ************* 673 674 675 logical :: estRiskFactor=.false. 676 character*40 chrom 677 678 ! 679 ! Pour la courbe 680 real :: x_p(npo), y_p(npo) 681 integer :: effp(np) 682 integer :: efft,ipos,ip,nm1,nm2,jm,ipar,indest,ntot,km 683 integer :: ilev,ief,iniveau,ico,lp,ilp,lm,ilm,i 684 logical :: impfem 685 real (kind=dp),dimension(size(par0)) :: par0_t 686 real (kind=dp),dimension(size(par1)) :: par1_t 687 real(kind=dp) :: dxmax 688 689 if ( present(printRiskFactor) ) then 690 estRiskFactor=printRiskFactor 691 end if 692 dxmax=absi(chr,nposx) 693 694 ! comptages 695 ! 696 efft=0 697 impfem=.false. 698 do ip=1,np 699 effp(ip)=0 700 nm1=nmp(ip)+1 701 nm2=nmp(ip+1) 702 do jm=nm1,nm2 703 if(estime(ic,jm))impfem=.true. 704 efft=efft+eff(jm) 705 effp(ip)=effp(ip)+eff(jm) 706 end do 707 end do 708 709 ! write(nficout,3003) 710 ! 3003 format(//1x,'Maximum likelihood ratio test :'/) 711 ! write(nficout,3000)dxmax*1.d2,xlrmax 712 ! 3000 format(1x,'The maximum is reached at position ',f7.2,' cM, with value ',f8.3/) 713 ! 714 ! impression des r�sultats sous H0 715 ! 716 write(nficout,600) 717 600 format(//80('*'),/,'Estimation of parameters under H0'//) 718 719 par0_t=par0 720 ! ************* deb chgt ************ 721 if (est_var) then 722 ! remise � l'�chelle des variables 723 ! 724 ! ipar-> 1 to np : std 725 ! np+1: General mean (les autres analyses on un mean / pere) 726 ! 727 par0_t=par0*sigt(ic) 728 end if 729 ! ************* fin chgt ************ 730 731 indest=0 732 if (est_var) then 733 indest=np 734 if (est_moy.and.vecsol0(1)) par0_t(np+1)=par0(np+1)+xmut(ic) 735 736 write(nficout,601) 737 601 format('Within sire standard deviation',/) 738 do ip = 1,np 739 write(nficout,602)trim(pere(ip)),par0_t(ip) 740 602 format(' sire ',associated, ' s.d. :',f10.3) 741 end do 742 endif 743 744 if (.not. estRiskFactor) then 745 write(nficout,603) 746 603 format(//,' parameter ',' estimable ? value ','precision'/) 747 if (est_moy.and.vecsol0(1)) then 748 indest=indest+1 749 write(nficout,604) par0_t(indest),precis0(1) 750 604 format('General mean yes ',2f10.3,1x) 751 else 752 write(nficout,6041) 753 6041 format('General mean no ') 754 end if 755 else 756 write(nficout,6030) 757 6030 format(//,' parameter ',' estimable ? risk factor '/) 758 759 endif 760 761 write(nficout,606) 762 606 format(/,'Sire polygenic effects') 763 par_refp=par0_t(indest+1) 764 do ip =1,np 765 i=ip 766 if (est_moy) i=ip+1 767 if (vecsol0(i)) then 768 indest=indest+1 769 if (.not. estRiskFactor) then 770 write(nficout,607)trim(pere(ip)), par0_t(indest),precis0(i) 771 607 format(' sire ',associated, 15x,'yes ',2f10.3,1x) 772 else 773 write(nficout,6070)trim(pere(ip)), dexp(par0_t(indest)-par_refp) 774 6070 format(' sire ',associated, 15x,'yes ',f10.3) 775 endif 776 else 777 write(nficout,608)trim(pere(ip)) 778 608 format(' sire ',associated, 15x,'no ') 779 end if 780 781 end do 782 783 ntot=np 784 par_refm=par0_t(indest+1) 785 if (est_moy) ntot=np+1 786 if(impfem) write(nficout,609) 787 609 format(/,'Dam polygenic effects') 788 km=0 789 do jm=1,nm 790 if (estime(ic,jm))then 791 km=km+1 792 if(vecsol0(ntot+km)) then 793 indest=indest+1 794 if (.not. estRiskFactor) then 795 write(nficout,610)trim(mere(jm)), par0_t(indest),precis0(ntot+km) 796 610 format(' dam ',associated, 15x,'yes ',2f10.3,1x) 797 else 798 write(nficout,6100)trim(mere(jm)), dexp(par0_t(indest)-par_refm) 799 endif 800 6100 format(' dam ',associated, 15x,'yes ',f10.3) 801 else 802 write(nficout,611)trim(mere(jm)) 803 end if 804 611 format(' dam ',associated, 15x,'no ') 805 end if 806 end do 807 808 ntot=ntot+nbfem 809 if(nbniv.ne.0)then 810 ilev=0 811 write(nficout,612) 812 612 format(/,'fixed effects') 813 do ief=1,nbef 814 par_refef(ief)=par0_t(indest+1) 815 do iniveau=1,nlev(modele(ic,3+ief)) 816 ilev=ilev+1 817 if (vecsol0(ntot+ilev)) then 818 indest=indest+1 819 if (.not. estRiskFactor) then 820 write(nficout,613)trim(namefix(modele(ic,3+ief))),iniveau,par0_t(indest),precis0(ntot+ilev) 821 613 format(a15,' level',i3,' yes ',2f10.3,1x) 822 else 823 write(nficout,6130)trim(namefix(modele(ic,3+ief))),iniveau,dexp(par0_t(indest)-par_refef(ief)) 824 6130 format(a15,' level',i3,' yes ',f10.3) 825 endif 826 else 827 write(nficout,614)trim(namefix(modele(ic,3+ief))),iniveau 828 614 format(a15,' level',i3,' no ') 829 end if 830 end do 831 end do 832 end if 833 834 ntot=ntot+nbniv 835 if(nbco.ne.0)then 836 write(nficout,615) 837 615 format(/,'covariables') 838 do ico=1,nbco 839 if (vecsol0(ntot+ico)) then 840 indest=indest+1 841 if (.not. estRiskFactor) then 842 write(nficout,616)trim(namecov(modele(ic,3+nbef+ico))),par0_t(indest),precis0(ntot+ico) 843 616 format(a15,12x,'yes ',2f10.3,1x) 844 else 845 write(nficout,6160)trim(namecov(modele(ic,3+nbef+ico))),dexp(par0_t(indest)) 846 6160 format(a15,12x,'yes ',f10.3) 847 endif 848 else 849 write(nficout,617)trim(namecov(modele(ic,3+nbef+ico))) 850 617 format(a15,12x,'no ') 851 end if 852 end do 853 end if 854 ! 855 ! impression des resultats sous H1 856 write(nficout,618) 857 618 format(//80('*'),/,'Estimation of parameters under H1'//) 858 859 par1_t=par1 860 ! 861 ! remise � l'�chelle des variables 862 ! 863 ! ************* deb chgt ************* 864 n_est_par=0 865 if (est_moy) n_est_par=1 866 ntot=n_est_par 867 868 if(est_var)then 869 ! ************* fin chgt ************* 870 par1_t=par1*sigt(ic) 871 end if 872 873 indest=0 874 if (est_var) then 875 indest=np 876 if(est_moy.and.vecsol1(1))par1_t(np+1)=par1(np+1)+xmut(ic) 877 878 write(nficout,619) 879 619 format('Within sire standard deviation ',/) 880 do ip = 1,np 881 write(nficout,620)trim(pere(ip)),par1_t(ip) 882 620 format(' sire ',associated, ' s.d. :',f10.3) 883 end do 884 endif 885 886 887 888 if (.not. estRiskFactor) then 889 write(nficout,621) 890 621 format(//,' parameter ',' estimable ? value ','precision'/) 891 if (est_moy.and.vecsol1(1)) then 892 indest=indest+1 893 write(nficout,604) par1_t(indest),precis1(1) 894 else 895 write(nficout,6041) 896 end if 897 else 898 write(nficout,6210) 899 6210 format(//,' parameter ',' estimable ? risk factor '/) 900 end if 901 902 if( hsire .or. hdam ) then 903 904 write(nficout,640) 905 640 format(/,'Haplotypes effects') 906 907 ! do i_haplo=1,nb_max_haplo(nposx) 908 ! do i_haplo=1,nb_max_haplo-1 909 do i_haplo=1,nb_haplo_reduit 910 ! 911 ! TRIM non disponible au CSIRO 912 ! 913 if(vecsol1(n_est_par+i_haplo)) then 914 indest=indest+1 915 if (.not. estRiskFactor)then 916 ! write(nficout,641)name_haplo_reduit(i_haplo+1)(:longhap),count_haplo(i_haplo+1),& 917 ! write(nficout,641)name_haplo_reduit(i_haplo),pb_haplo_reduit(i_haplo),& 918 ! par1_t(indest),precis1(n_est_par+i_haplo) 919 ! print*,name_haplo_reduit(i_haplo),pb_haplo_reduit(i_haplo),& 920 ! par1_t(indest),precis1(n_est_par+i_haplo) 921 write(nficout,641)name_haplo(i_haplo)(:longhap),count_haplo(i_haplo),& 922 par1_t(indest),precis1(n_est_par+i_haplo) 923 ! 641 format(' haplotype ',a, 14x,' yes ',2f10.3,9x) 924 641 format(' haplotype ',associated, ' freq = ',f4.2,2x,' yes ',2f10.3,9x) 925 else 926 ! write(nficout,642)name_haplo(nposx,i_haplo)(:longhap), count_haplo(nposx,i_haplo),& 927 ! write(nficout,642)name_haplo_reduit(i_haplo+1)(:longhap), count_haplo(i_haplo+1),& 928 write(nficout,642)name_haplo_reduit(i_haplo),pb_haplo_reduit(i_haplo),& 929 dexp(par1_t(indest)) 930 642 format(' haplotype ',associated, ' freq = ',f4.2,2x,' yes ',f10.3,9x) 931 endif 932 else 933 ! write(nficout,643)name_haplo(nposx,i_haplo)(:longhap),count_haplo(nposx,i_haplo) 934 ! write(nficout,643)name_haplo_reduit(i_haplo+1)(:longhap),count_haplo(i_haplo+1) 935 write(nficout,643)name_haplo_reduit(i_haplo),pb_haplo_reduit(i_haplo) 936 643 format(' haplotype ',associated, ' freq = ',f4.2,2x,' no ') 937 end if 938 end do 939 940 ! n_est_par=n_est_par+nb_max_haplo(nposx) 941 ! n_est_par=n_est_par+nb_max_haplo-1 942 n_est_par=n_est_par+nb_haplo_reduit 943 944 end if ! option 945 946 if(qtl) then 947 write(nficout,624) 948 624 format(/,'Sire QTL effects') 949 950 if(ntlevp.eq.1) then 951 write(nficout,6241) 952 6241 format(59x,'allelic origin'/) 953 do ip = 1,np 954 955 i=n_est_par+ip 956 957 ! ****** POUR JM ******** 958 call liste_haplo(chr,dxmax,nposx,hsire,hdam) 959 call sort_haplo(chr,nposx,hsire,hdam) 960 ! print *,'num_haplo_pere:',num_haplo_pere(ip,1),num_haplo_pere(ip,2) 961 chrom=name_haplo(num_haplo_pere(ip,1,1))(:longhap)//' '//& 962 name_haplo(num_haplo_pere(ip,2,1))(:longhap) 963 964 if (vecsol1(i)) then 965 indest=indest+1 966 if(phasp(chr,ip)) then 967 if (.not. estRiskFactor)then 968 write(nficout,6251)trim(pere(ip)), par1_t(indest),precis1(i),chrom 969 6251 format(' sire ',associated, 14x,' yes ',2f10.3,9x,'known ',associated) 970 else 971 write(nficout,62510)trim(pere(ip)), dexp(par1_t(indest)),chrom 972 62510 format(' sire ',associated, 14x,' yes ',f10.3,9x,'known ',associated) 973 endif 974 else 975 if (.not. estRiskFactor)then 976 write(nficout,6252)trim(pere(ip)), par1_t(indest),precis1(i),chrom 977 6252 format(' sire ',associated, 14x,' yes ',2f10.3,9x,'unknown ',associated) 978 else 979 write(nficout,62520)trim(pere(ip)), dexp(par1_t(indest)),chrom 980 62520 format(' sire ',associated, 14x,' yes ',f10.3,9x,'unknown ',associated) 981 endif 982 end if 983 else 984 write(nficout,626)trim(pere(ip)),chrom 985 626 format(' sire ',associated, 14x,' no ',associated) 986 end if 987 end do 988 end if 989 990 if(ntlevp.ne.1) then 991 write(nficout,6242) 992 6242 format(59x,'allelic origin'/) 993 do ip = 1,np 994 do lp = 1,ntlevp 995 ! ************* deb chgt ************* 996 ilp=ntlevp*(ip-1)+lp-1+ n_est_par 997 ! ilp=ntlevp*(ip-1)+lp-1 998 ! if (est_moy) ilp=ntlevp*(ip-1)+lp 999 ! ************* fin chgt ************* 1000 if (vecsol1(1+ilp)) then 1001 indest=indest+1 1002 if (.not. estRiskFactor)then 1003 if(phasp(chr,ip)) then 1004 write(nficout,6271)trim(pere(ip)), lp,par1_t(indest),precis1(1+ilp) 1005 6271 format(' sire ',associated, 'level ',i3,5x,' yes ',2f10.3,9x,'known') 1006 else 1007 write(nficout,6272)trim(pere(ip)), lp,par1_t(indest),precis1(1+ilp) 1008 6272 format(' sire ',associated, 'level ',i3,5x,' yes ',2f10.3,9x,'unknown') 1009 end if 1010 else 1011 if(phasp(chr,ip)) then 1012 write(nficout,62710)trim(pere(ip)), lp,dexp(par1_t(indest)) 1013 62710 format(' sire ',associated, 'level ',i3,5x,' yes ',f10.3,9x,'known') 1014 else 1015 write(nficout,62720)trim(pere(ip)), lp,dexp(par1_t(indest)) 1016 62720 format(' sire ',associated, 'level ',i3,5x,' yes ',f10.3,9x,'unknown') 1017 end if 1018 endif 1019 else 1020 write(nficout,628)trim(pere(ip)),lp 1021 628 format(' sire ',associated, 'level ',i3,5x,' no ') 1022 end if 1023 end do 1024 end do 1025 end if 1026 1027 ! ************* deb chgt ************* 1028 ntot=ntlevp*np+ n_est_par 1029 ! ntot=ntlevp*np 1030 ! if (est_moy) ntot=1+ntlevp*np 1031 ! ************* fin chgt ************* 1032 if(impfem) write(nficout,629) 1033 629 format(/,'Dam QTL effects') 1034 1035 if(ntlevm.eq.1) then 1036 1037 if(impfem) write(nficout,6241) 1038 do jm=1,nm 1039 if (estime(ic,jm))then 1040 km=iam(ic,repfem(jm)) 1041 if(vecsol1(ntot+km)) then 1042 indest=indest+1 1043 1044 if (.not. estRiskFactor)then 1045 if(phasm(chr,jm)) then 1046 write(nficout,6301)trim(mere(jm)), par1_t(indest),precis1(ntot+km) 1047 6301 format(' dam ',associated, 14x,' yes ',2f10.3,9x,'known') 1048 else 1049 write(nficout,6302)trim(mere(jm)), par1_t(indest),precis1(ntot+km) 1050 6302 format(' dam ',associated, 14x,' yes ',2f10.3,9x,'unknown') 1051 end if 1052 else 1053 if(phasm(chr,jm)) then 1054 write(nficout,63010)trim(mere(jm)), dexp(par1_t(indest)) 1055 63010 format(' dam ',associated, 14x,' yes ',f10.3,9x,'known') 1056 else 1057 write(nficout,63020)trim(mere(jm)), dexp(par1_t(indest)) 1058 63020 format(' dam ',associated, 14x,' yes ',f10.3,9x,'unknown') 1059 end if 1060 endif 1061 1062 end if 1063 end if 1064 end do 1065 end if 1066 1067 if(ntlevm.ne.1) then 1068 1069 if(impfem) write(nficout,6242) 1070 do jm=1,nbfem 1071 if(estime(ic,jm))then 1072 km=iam(ic,repfem(jm)) 1073 do lm=1,ntlevm 1074 ilm=ntlevm*(km-1)+lm 1075 if(vecsol1(ntot+ilm))then 1076 indest=indest+1 1077 1078 if (.not. estRiskFactor)then 1079 if(phasm(chr,jm)) then 1080 write(nficout,6321)trim(mere(jm)),lm,par1_t(indest),precis1(ntot+ilm) 1081 6321 format(' dam ',associated, 'level ',i3,6x,' yes ',2f10.3,9x,'known') 1082 else 1083 write(nficout,6322)trim(mere(jm)),lm,par1_t(indest),precis1(ntot+ilm) 1084 6322 format(' dam ',associated, 'level ',i3,6x,' yes ',2f10.3,9x,'unknown') 1085 end if 1086 else 1087 if(phasm(chr,jm)) then 1088 write(nficout,63210)trim(mere(jm)),lm,dexp(par1_t(indest)) 1089 63210 format(' dam ',associated, 'level ',i3,6x,' yes ',f10.3,9x,'known') 1090 else 1091 write(nficout,63220)trim(mere(jm)),lm,dexp(par1_t(indest)) 1092 63220 format(' dam ',associated, 'level ',i3,6x,' yes ',f10.3,9x,'unknown') 1093 end if 1094 endif 1095 1096 end if 1097 end do 1098 end if 1099 end do 1100 1101 end if 1102 1103 write(nficout,*) 1104 write(nficout,*)' NOTE: known allelic origin means QTL effect = maternal - paternal allele effects' 1105 1106 ntot=ntot+namest(ic)*ntlevm 1107 end if ! option 1108 1109 1110 write(nficout,606) 1111 par_refp=par1_t(indest+1) 1112 do ip = 1,np 1113 if (vecsol1(ntot+ip)) then 1114 indest=indest+1 1115 if (.not. estRiskFactor)then 1116 write(nficout,607)trim(pere(ip)),par1_t(indest),precis1(ntot+ip) 1117 else 1118 write(nficout,6070)trim(pere(ip)),dexp(par1_t(indest)-par_refp) 1119 endif 1120 else 1121 write(nficout,608)trim(pere(ip)) 1122 end if 1123 end do 1124 1125 ntot=ntot+np 1126 write(nficout,609) 1127 km=0 1128 par_refm=par1_t(indest+1) 1129 do jm=1,nm 1130 if (estime(ic,jm)) then 1131 km=km+1 1132 if(vecsol1(ntot+km)) then 1133 indest=indest+1 1134 if (.not. estRiskFactor)then 1135 write(nficout,610)trim(mere(jm)),par1_t(indest),precis1(ntot+km) 1136 else 1137 write(nficout,6100)trim(mere(jm)),dexp(par1_t(indest)-par_refm) 1138 endif 1139 else 1140 write(nficout,611)trim(mere(jm)) 1141 end if 1142 end if 1143 end do 1144 1145 ntot=ntot+nbfem 1146 if(nbniv.ne.0)then 1147 ilev=0 1148 write(nficout,612) 1149 1150 do ief=1,nbef 1151 par_refef(ief)=par0_t(indest+1) 1152 do iniveau=1,nlev(modele(ic,3+ief)) 1153 ilev=ilev+1 1154 if (vecsol1(ntot+ilev)) then 1155 indest=indest+1 1156 if (.not. estRiskFactor)then 1157 write(nficout,613)trim(namefix(modele(ic,3+ief))),iniveau,par1_t(indest),precis1(ntot+ilev) 1158 else 1159 write(nficout,6130)trim(namefix(modele(ic,3+ief))),iniveau,dexp(par1_t(indest)-par_refef(ief)) 1160 endif 1161 1162 else 1163 write(nficout,614)trim(namefix(modele(ic,3+ief))),iniveau 1164 end if 1165 end do 1166 end do 1167 end if 1168 1169 ntot=ntot+nbniv 1170 if(nbco.ne.0)then 1171 write(nficout,615) 1172 do ico=1,nbco 1173 if (vecsol1(ntot+ico)) then 1174 indest=indest+1 1175 if (.not. estRiskFactor)then 1176 write(nficout,616)trim(namecov(modele(ic,3+nbef+ico))),par1_t(indest),precis1(ntot+ico) 1177 else 1178 write(nficout,6160)trim(namecov(modele(ic,3+nbef+ico))),dexp(par1_t(indest)) 1179 endif 1180 else 1181 write(nficout,617)trim(namecov(modele(ic,3+nbef+ico))) 1182 end if 1183 end do 1184 end if 1185 1186 1187 return 1188 end subroutine courbe_lin_ldla
end_analyse_gen
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
end_analyse_gen
DESCRIPTION
deallocation of arrays provide by this module.
NOTES
SOURCE
247 subroutine end_analyse_gen 248 249 deallocate (estmum) 250 deallocate (eff) 251 deallocate (somcd) 252 deallocate (somy) 253 deallocate (cary) 254 deallocate (effp) 255 deallocate (effdf) 256 deallocate (somcddf) 257 deallocate (somydf) 258 deallocate (carydf) 259 deallocate (sig0) 260 deallocate (xmu0p) 261 deallocate (xmu0m) 262 263 end subroutine end_analyse_gen
get_eff_paternal_and_total
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
get_eff_paternal_and_total
DESCRIPTION
Get the number of femal with enough progeny with perf
OUTPUTS
effp : number of progenies in the half sib family ip with a trait value efft : number of progenies in the dataset with a trait value
NOTES
SOURCE
465 subroutine get_eff_paternal_and_total(ic,effp,efft) 466 integer , intent(in) :: ic 467 integer , dimension(np) , intent(out) :: effp 468 integer , intent(out) :: efft 469 470 integer :: ip, nm1,nm2,jm 471 472 efft=0 473 do ip=1,np 474 effp(ip)=0 475 nm1=nmp(ip)+1 476 nm2=nmp(ip+1) 477 do jm=nm1,nm2 478 efft=efft+count(presentc(ic,ndm(jm)+1:ndm(jm+1))) 479 effp(ip)=effp(ip)+eff(jm) 480 print *,'get-eff:',efft 481 end do 482 end do 483 end subroutine get_eff_paternal_and_total
init_analyse_gen
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
init_analyse_gen
DESCRIPTION
allocation of arrays provide by this module.
NOTES
SOURCE
207 subroutine init_analyse_gen 208 integer :: stat 209 210 allocate (estmum(np),STAT=stat) 211 call check_allocate(stat,'estmum [m_qtlmap_analyse_gen]') 212 allocate (eff(nm),STAT=stat) 213 call check_allocate(stat,'eff [m_qtlmap_analyse_gen]') 214 allocate (somcd(nm),STAT=stat) 215 call check_allocate(stat,'somcd [m_qtlmap_analyse_gen]') 216 allocate (somy(nm),STAT=stat) 217 call check_allocate(stat,'somy [m_qtlmap_analyse_gen]') 218 allocate (cary(nm),STAT=stat) 219 call check_allocate(stat,'cary [m_qtlmap_analyse_gen]') 220 allocate (effp(np)) 221 !call check_allocate(stat,'estime [m_qtlmap_analyse_gen]') 222 allocate (effdf(np),STAT=stat) 223 call check_allocate(stat,'effdf [m_qtlmap_analyse_gen]') 224 allocate (somcddf(np),STAT=stat) 225 call check_allocate(stat,'somcddf [m_qtlmap_analyse_gen]') 226 allocate (somydf(np),STAT=stat) 227 call check_allocate(stat,'somydf [m_qtlmap_analyse_gen]') 228 allocate (carydf(np),STAT=stat) 229 call check_allocate(stat,'carydf [m_qtlmap_analyse_gen]') 230 allocate (sig0(np),STAT=stat) 231 call check_allocate(stat,'sig0 [m_qtlmap_analyse_gen]') 232 allocate (xmu0p(np),STAT=stat) 233 call check_allocate(stat,'xmu0p [m_qtlmap_analyse_gen]') 234 allocate (xmu0m(nm),STAT=stat) 235 call check_allocate(stat,'xmu0m [m_qtlmap_analyse_gen]') 236 end subroutine init_analyse_gen
optinit
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
optinit
DESCRIPTION
initialisation for half sib, full sib family :
* number of progenies
* sum of cd
* sum of y
* sum of square y
* sig0
* mean0
INPUTS
ic : index of the trait
NOTES
Initialise les ecart types et moyennes intra-famille Sous programme appele par analyse Version XXX _ HEG030507
SOURCE
285 subroutine optinit(ic) 286 integer , intent(in) :: ic 287 288 ! Divers 289 integer :: ifem,ip,nm1,nm2,jm,nd1,nd2,kd,imumest,namest 290 real (kind=dp) :: somyp 291 logical :: estfem(nfem) 292 ! 293 !****************************************************************************** 294 effp=0 295 imumest=0 296 estfem=.false. 297 do ip=1,np 298 somyp=0.d0 299 sig0(ip)=0.d0 300 effdf(ip)=0 301 somcddf(ip)=0.d0 302 somydf(ip)=0.d0 303 carydf(ip)=0.d0 304 estmum(ip)=0 305 nm1=nmp(ip)+1 306 nm2=nmp(ip+1) 307 do jm=nm1,nm2 308 eff(jm)=0 309 somcd(jm)=0.d0 310 somy(jm)=0.d0 311 cary(jm)=0.d0 312 nd1=ndm(jm)+1 313 nd2=ndm(jm+1) 314 ifem=repfem(jm) 315 do kd=nd1,nd2 316 if(presentc(ic,kd)) then 317 eff(jm)=eff(jm)+1 318 somcd(jm)=somcd(jm)+cd(ic,kd) 319 somy(jm)=somy(jm)+y(ic,kd)*cd(ic,kd) 320 cary(jm)=cary(jm)+(y(ic,kd)*y(ic,kd))*cd(ic,kd) 321 end if 322 end do 323 324 effp(ip)=effp(ip)+eff(jm) 325 somyp=somyp+somy(jm) 326 if( estime(ic,jm) ) then 327 imumest=imumest+1 328 estmum(ip)=estmum(ip)+1 329 xmu0m(imumest)=somy(jm)/dble(eff(jm)) 330 estfem(ifem)=.true. 331 else 332 effdf(ip)=effdf(ip)+eff(jm) 333 somcddf(ip)=somcddf(ip)+somcd(jm) 334 somydf(ip)=somydf(ip)+somy(jm) 335 carydf(ip)=carydf(ip)+cary(jm) 336 end if 337 end do 338 if (effp(ip) == 0.d0) then 339 call stop_application('Father ['//trim(pere(ip))//'] has got no child with trait value') 340 else 341 xmu0p(ip)=somyp/dble(effp(ip)) 342 end if 343 do jm=nm1,nm2 344 nd1=ndm(jm)+1 345 nd2=ndm(jm+1) 346 do kd=nd1,nd2 347 if(presentc(ic,kd)) then 348 sig0(ip)=sig0(ip)+(y(ic,kd)-xmu0p(ip))*(y(ic,kd)-xmu0p(ip)) 349 end if 350 end do 351 end do 352 if ( effp(ip) == 1) then 353 call stop_application('Father ['//trim(pere(ip))//'] has got only one child with trait value') 354 else 355 sig0(ip)=dsqrt(sig0(ip)/dble(effp(ip)-1)) 356 end if 357 358 if(estmum(ip).gt.0) then 359 estmum(ip)=estmum(ip)-1 360 nmumest(ic)=nmumest(ic)-1 361 end if 362 end do 363 end subroutine optinit
optinit_da
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
optinit_da
DESCRIPTION
INPUTS
yda :
NOTES
SOURCE
375 subroutine optinit_da(yda) 376 ! integer , intent(in) :: ic 377 real (kind=dp), dimension (nd) :: yda 378 379 ! Divers 380 integer :: effp,ifem,ip,nm1,nm2,jm,nd1,nd2,kd,ic,imumest 381 real (kind=dp) :: somyp 382 logical :: estfem(nfem) 383 ! 384 !****************************************************************************** 385 ic=1 386 imumest=0 387 388 do ifem=1,nfem 389 estfem(ifem)=.false. 390 end do 391 do ip=1,np 392 effp=0 393 somyp=0.d0 394 sig0(ip)=0.d0 395 effdf(ip)=0 396 somydf(ip)=0.d0 397 carydf(ip)=0.d0 398 estmum(ip)=0 399 nm1=nmp(ip)+1 400 nm2=nmp(ip+1) 401 do jm=nm1,nm2 402 eff(jm)=0 403 somy(jm)=0.d0 404 cary(jm)=0.d0 405 nd1=ndm(jm)+1 406 nd2=ndm(jm+1) 407 ifem=repfem(jm) 408 do kd=nd1,nd2 409 if(presentc(ic,kd)) then 410 eff(jm)=eff(jm)+1 411 somy(jm)=somy(jm)+yda(kd) 412 cary(jm)=cary(jm)+(yda(kd)*yda(kd)) 413 end if 414 end do 415 416 effp=effp+eff(jm) 417 somyp=somyp+somy(jm) 418 if(estime(ic,jm)) then 419 imumest=imumest+1 420 estmum(ip)=estmum(ip)+1 421 xmu0m(imumest)=somy(jm)/dble(eff(jm)) 422 estfem(ifem)=.true. 423 else 424 effdf(ip)=effdf(ip)+eff(jm) 425 somydf(ip)=somydf(ip)+somy(jm) 426 carydf(ip)=carydf(ip)+cary(jm) 427 end if 428 end do 429 if (effp == 0.d0) then 430 call log_mess('optinit : effp <- 0.d0',WARNING_DEF) 431 xmu0p(ip)= 0.d0 432 else 433 xmu0p(ip)=somyp/dble(effp) 434 end if 435 do jm=nm1,nm2 436 nd1=ndm(jm)+1 437 nd2=ndm(jm+1) 438 do kd=nd1,nd2 439 if(presentc(ic,kd)) then 440 sig0(ip)=sig0(ip)+(yda(kd)-xmu0p(ip))*(yda(kd)-xmu0p(ip)) 441 end if 442 end do 443 end do 444 sig0(ip)=dsqrt(sig0(ip)/dble(effp-1)) 445 if(estmum(ip).gt.0) then 446 estmum(ip)=estmum(ip)-1 447 nmumest(:)=imumest-1 448 end if 449 end do 450 451 end subroutine optinit_da
scale_variables_dam
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
scale_variables_dam
DESCRIPTION
INPUTS
chr : index of the chromosome
ic : index of the trait
am : qtleffect
sigt : trait variance
OUTPUTS
xmoym1 :
am1 :
subphasm :
nmsub :
submere :
impfem :
NOTES
SOURCE
506 subroutine scale_variables_dam(chr,ic,am,xmoym,sigt,submere,nmsub,am1,xmoym1,subphasm,impfem) 507 integer ,intent(in) :: chr 508 integer ,dimension(:), intent(in) :: ic 509 real (kind=dp) , dimension(nm),intent(in) :: xmoym 510 real (kind=dp) , dimension(nm),intent(in) :: am 511 real (kind=dp) ,intent(in) :: sigt 512 real (kind=dp) , dimension(nm),intent(out) :: xmoym1 513 real (kind=dp) , dimension(nm),intent(out) :: am1 514 logical , dimension(nm),intent(out) :: subphasm 515 integer , intent(out) :: nmsub 516 character(len=LEN_DEF) , dimension(nm),intent(out) :: submere 517 logical , intent(out) :: impfem 518 519 520 integer :: ip, nm1,nm2,jm,i 521 logical ,dimension(size(ic)) :: estime_ic 522 523 impfem=.false. 524 nmsub = 0 525 526 do ip=1,np 527 nm1=nmp(ip)+1 528 nm2=nmp(ip+1) 529 do jm=nm1,nm2 530 do i=1,size(ic) 531 estime_ic(i) = estime(ic(i),jm) 532 end do 533 if( count(estime_ic(:))== size(ic) )then 534 impfem=.true. 535 nmsub=nmsub+1 536 submere(nmsub)=mere(jm) 537 subphasm(nmsub)=phasm(chr,jm) 538 xmoym1(nmsub)=xmoym(jm)*sigt 539 am1(nmsub)=am(jm)*sigt 540 end if 541 end do 542 end do 543 end subroutine
set_ntnivmax
[ Top ] [ m_qtlmap_analyse_gen ] [ Subroutines ]
NAME
set_ntnivmax
DESCRIPTION
get the maximum number parameter to build the incidence matrix
INPUTS
ic : index of the trait nqtl : the hypothesis
OUTPUTS
ntnivmax : maximum number of parameter nteffmax : maximum number of effect ntlev : number of level interaction fixed effect-qtl nbniv : number of level fixed effect
NOTES
SOURCE
562 subroutine set_ntnivmax(ic,nqtl,ntnivmax,nteffmax,ntlev,nbniv) 563 integer ,intent(in) :: ic 564 integer ,intent(in) :: nqtl 565 integer ,intent(out) :: ntnivmax 566 integer ,intent(out) :: nteffmax 567 integer ,intent(out) :: ntlev 568 integer ,intent(out) :: nbniv 569 570 integer :: nbtef,nbtco,nbtint 571 integer :: nbtp,nbtm,jef 572 573 nbtef=modele(ic,1) 574 nbtco=modele(ic,2) 575 nbtint=modele(ic,3) 576 ntlev=1 577 nbtp=3+nbtef+nbtco 578 nbtm=nbtp+nbtint 579 580 !compute number of level for qtl*interaction 581 if(nbtint > 0) then 582 do jef=1,nbtint 583 ntlev=ntlev*nlev(modele(ic,nbtp+jef)) 584 end do 585 end if 586 587 ! number of level for fixed effect 588 nbniv=0 589 do jef=1,nbtef 590 nbniv=nbniv+nlev(modele(ic,3+jef)) 591 end do 592 593 594 !General mean + polygenic sire + polygenic dam + 2*qtl*interaction+nbniv effet fixed + number of covariate 595 ! namest : nombre de fem estim 596 ! count(estime) : nombre de dam estimable 597 ntnivmax = 1 + np + count(estime(ic,:)) + nqtl*(ntlev*np+ntlev*namest(ic))+nbniv+nbtco + NB_HAPLO_PRIOR 598 nteffmax = 3+2*nqtl+nbtef+nbtco + NB_HAPLO_PRIOR 599 600 call log_mess('Value NTNIV MAX:'//trim(str(ntnivmax)),VERBOSE_DEF) 601 call log_mess('Value NTEFF MAX:'//trim(str(nteffmax)),VERBOSE_DEF) 602 603 end subroutine set_ntnivmax