m_qtlmap_cd

[ Top ] [ TOOLS ] [ Modules ]

NAME

    m_qtlmap_cd

DESCRIPTION

NOTES

SEE ALSO


add_animal_genea

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    add_animal_genea

DESCRIPTION

   Add a animal-triplet (sire-dam-progeny) inside a list of type GENEALOGY_TYPE

INPUTS

   sire   : name of the sire
   dam    : name of the dam
 progeny  : name of the progeny
  nd      : the index corresponding to the animal array

INPUTS/OUTPUTS

  dataset : anima-triplet list

SOURCE

 73      subroutine add_animal_genea(dataset,sire,dam,progeny,nd)
 74           type (GENEALOGY_TYPE)  , intent(inout)  :: dataset
 75           character(LEN=LEN_DEF) , intent(in)     :: sire,dam,progeny
 76           integer  ,optional     , intent(in)     :: nd             ! correspond a l index du tableau animal =>progeny
 77 
 78           integer  ,  parameter :: BLOCK_ALLOC = 1000
 79           logical :: findSire,findDam,findProg
 80           integer :: i,idSire,idDam,associated
 81 
 82           if ( .not. associated(dataset%listNA) ) then
 83              call extend(dataset%listNA,BLOCK_ALLOC)
 84              call extend(dataset%trioKd,BLOCK_ALLOC)
 85              call extend(dataset%corranim,BLOCK_ALLOC)
 86           end if
 87 
 88           ! Le pere et la mere existe ?
 89           findSire=.false.
 90           findDam=.false.
 91 
 92          ! print *,sire,dam,progeny
 93 
 94           do i=1,dataset%na
 95             if ( trim(dataset%listNA(i)) == trim(sire)) then
 96                findSire=.true.
 97                idSire = i
 98             else if ( trim(dataset%listNA(i)) == trim(dam)) then
 99                findDam=.true.
100                idDam = i
101             end if
102             if ( findSire .and. findDam ) exit
103           end do
104 
105           if ( .not. findSire ) then
106               !il faut creer une entree pour le pere
107               if ( size(dataset%listNA) <= dataset%na ) then
108                   call extend(dataset%listNA,BLOCK_ALLOC)
109                   call extend(dataset%trioKd,BLOCK_ALLOC)
110               end if
111               dataset%na = dataset%na + 1
112               dataset%listNA(dataset%na) = trim(sire)
113               idSire = dataset%na
114               !dataset%trioKd(dataset%na) => null()
115           end if
116 
117           if ( .not. findDam ) then
118               !il faut creer une entree pour la mere
119               if ( size(dataset%listNA) <= dataset%na ) then
120                   call extend(dataset%listNA,BLOCK_ALLOC)
121                   call extend(dataset%trioKd,BLOCK_ALLOC)
122               end if
123               dataset%na = dataset%na + 1
124               dataset%listNA(dataset%na) = trim(dam)
125               idDam = dataset%na
126               !dataset%trioKd(dataset%na) => null()
127           end if
128 
129           !check pour savoir si la progeniture n existe pas deja !
130           findProg = .false.
131 
132 !          do i=1,dataset%na
133 !             if ( trim(dataset%listNA(i)) == trim(progeny)) then
134 !                print *,"double entry for kd:",progeny
135 !                stop
136 !             end if
137 !          end do
138 
139           if ( size(dataset%listNA) <= dataset%na ) then
140             call extend(dataset%listNA,BLOCK_ALLOC)
141             call extend(dataset%trioKd,BLOCK_ALLOC)
142           end if
143 
144           dataset%na = dataset%na + 1
145 
146           dataset%listNA(dataset%na) = trim(progeny)
147 
148           dataset%trioKd(dataset%na)%idSire = idSire
149           dataset%trioKd(dataset%na)%idDam  = idDam
150           dataset%trioKd(dataset%na)%idKd   = dataset%na
151 
152           if (present(nd)) then
153              if ( nd > size(dataset%corrAnim) ) then
154                call extend(dataset%corranim,BLOCK_ALLOC)
155              end if
156              dataset%corranim(nd) = dataset%trioKd(dataset%na)%idKd
157           end if
158 
159      end subroutine add_animal_genea

calcul_corcd

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    calcul_corcd

DESCRIPTION

   Compute the cd correlation (used by the multi trait analysis)

   ic,jc 2 traits
   kd the sire/dam progenies
   NCK (ic,kd)      : number of progenies from sire/dam kd which theses progenies have a value for the trait ic
   NCiCjK(kd,ic,jc) : number of progenies from sire/dam kd which theses progenies have a value for the trait ic and jc jointly
   RhoG(ic,jc)      : genetic correlation between ic and jc
   RhoP(ic,jc)      : phenotypic correlation between ic and jc


                          RhoG(ic,jc) * SQRT( h2(ic)*h2(jc) )
   v(ic,jc)          =    -----------------------------------
                                     RhoP(ic,jc)

                       NCiCjK(kd,ic,jc) + 0.25d0*v(ic,jc)* [ NCK(kd,ic)*NCK(kd,jc) - NCiCjK(kd,ic,jc) ]
   CORR_CD(ic,jc,kd) = ---------------------------------------------------------------------------------
                                                     NCK(kd,ic) * NCK(kd,jc)

INPUTS/OUTPUTS

 dataset     : anima-triplet list

SOURCE

354      subroutine calcul_corcd(dataset)
355        type (GENEALOGY_TYPE)         ,intent(inout)     :: dataset
356 
357        integer :: ip,jm,kd,s,i,j,ic,k,jc
358        real(kind=dp) :: r,v(ncar,ncar)
359 
360        allocate (dataset%corcd(ncar,ncar,nd))
361 
362        do ic=1,ncar
363         do jc=ic+1,ncar
364              v(ic,jc) = (RhoG(ic,jc)/RhoP(ic,jc))*sqrt(h2(ic)*h2(jc))
365         end do
366        end do
367 
368        dataset%corcd = 1.d0
369        do ip=1,np
370          do jm=nmp(ip)+1,nmp(ip+1)
371            do kd=ndm(jm)+1,ndm(jm+1)
372              do ic=1,ncar
373                if ( dataset%NCK(kd,ic) > 0 ) then
374                  do jc=ic+1,ncar
375                    if ( dataset%NCK(kd,jc) > 0 ) then
376                      r = real( dataset%NCK(kd,ic)*dataset%NCK(kd,jc) - dataset%NCiCjK(kd,ic,jc))
377                      dataset%corcd(ic,jc,kd) = real(dataset%NCiCjK(kd,ic,jc)) + 0.25d0*v(ic,jc)*r
378                      dataset%corcd(ic,jc,kd) = dataset%corcd(ic,jc,kd) / ( real( dataset%NCK(kd,ic)) * real( dataset%NCK(kd,jc)))
379                      dataset%corcd(jc,ic,kd) = dataset%corcd(ic,jc,kd)
380                    end if
381                  end do
382                end if
383              end do
384            end do
385          end do
386        end do
387 
388      end subroutine calcul_corcd

calcul_y_cd

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    calcul_y_cd

DESCRIPTION

   Compute the new value phenotypic and a cd value of progenies who have a progeniture list. The old phenotypic value is removed if exist.
   depends the heritability of trait h2

   NCK (ic,kd)      : number of progenies from sire/dam kd which theses progenies have a value for the trait ic
   NCiCjK(kd,ic,jc) : number of progenies from sire/dam kd which theses progenies have a value for the trait ic and jc jointly

   SUM Yprog(kd)    : sum of trait values of progenies of sire/dam kd
   h2(ic)           : heritability of the trait ic

   Y(ic,kd)  = SYM Yprog(kd) / NCK (ic,kd)
   CD(ic,kd) = 1.d0 + [h2(ic) * (NCK(kd,ic) - 1.d0) / 4] /  NCK (ic,kd)

INPUTS/OUTPUTS

 dataset     : anima-triplet list

SOURCE

269      subroutine calcul_y_cd(dataset)
270        type (GENEALOGY_TYPE)         ,intent(inout)     :: dataset
271 
272        type(TRIO_TYPE) ,dimension(:) ,pointer  :: list
273        integer :: ip,jm,kd,s,i,j,ic,k,jc
274        real(kind=dp) :: sumy(ncar)
275 
276        call log_mess("Calcul CD and Y with generation 3",INFO_DEF)
277 
278        allocate (dataset%NCK(nd,ncar))
279        dataset%NCK=0
280        allocate (dataset%NCiCjK(nd,ncar,ncar))
281        dataset%NCiCjK=0
282 
283        do ip=1,np
284          do jm=nmp(ip)+1,nmp(ip+1)
285            do kd=ndm(jm)+1,ndm(jm+1)
286              call get_listProgenies(dataset,kd,animal(kd),s,list)
287              call log_mess(trim(animal(kd))//":Number of progeny :"//str(s),DEBUG_DEF)
288              sumy=0.d0
289              do i=1,s
290               do j=1,ncar
291                if ( list(i)%presentc(j) ) then
292                  dataset%NCK(kd,j) = dataset%NCK(kd,j) + 1
293                  sumy(j) = sumy(j) + list(i)%perfkd(j)
294                  do k=j+1,ncar
295                     if ( list(i)%presentc(k) ) then
296                      dataset%NCiCjK(kd,j,k) = dataset%NCiCjK(kd,j,k) + 1
297                      dataset%NCiCjK(kd,k,j) = dataset%NCiCjK(kd,j,k)
298                     end if
299                  end do
300                end if
301               end do
302              end do
303 
304 
305              do ic=1,ncar
306                presentc(ic,kd) = .false. !on invalide la performave de l individu kd
307                if ( dataset%NCK(kd,ic) > 0 ) then
308                 presentc(ic,kd) = .true.
309                 y(ic,kd)  = sumy(ic) / real(dataset%NCK(kd,ic))
310                 cd(ic,kd) = 1.d0 + 0.25d0*h2(ic) * (real(dataset%NCK(kd,ic)) - 1.d0)
311                 cd(ic,kd) = cd(ic,kd) / real(dataset%NCK(kd,ic))
312                 !cd = 0.25d0*h2(ic)
313        !         print *,ic,kd,cd(ic,kd)
314                end if
315              end do
316 
317 
318              if (s>0) deallocate(list)
319            end do
320          end do
321        end do
322 
323        !cd=0.0d0
324 
325      end subroutine calcul_y_cd

extend

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    extend

DESCRIPTION

    Generic interface to increase an array

SOURCE

47     interface extend
48         module procedure extendListNa
49         module procedure extendListTrio
50         module procedure extendListCorrAnim
51     end interface extend

extendListCorrAnim

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    extendListCorrAnim

DESCRIPTION

   increase the list of type integer with BLOCK_ALLOC values

INPUTS

  BLOCK_ALLOC : the size to extends

INPUTS/OUTPUTS

  list        : the integer list

SOURCE

503      subroutine extendListCorrAnim(list,BLOCK_ALLOC)
504         integer                ,dimension(:) ,pointer  ,intent(inout) :: list
505         integer                                        ,intent(in)    :: BLOCK_ALLOC
506 
507         integer ,dimension(:) ,pointer  :: listTemp => null()
508         integer :: i
509 
510         !start
511         if ( .not. associated(list) ) then
512            allocate (list(BLOCK_ALLOC))
513            list=0
514            return
515         end if
516 
517         allocate ( listTemp(size(list)) )
518         listTemp = list
519         deallocate ( list )
520         allocate ( list(size(listTemp)+BLOCK_ALLOC) )
521         list=0
522         list(:size(listTemp)) = listTemp
523         deallocate (listTemp)
524 
525         return
526      end subroutine extendListCorrAnim

extendListNa

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    extendListNa

DESCRIPTION

   increase the list of type character(LEN=LEN_DEF) with BLOCK_ALLOC values

INPUTS

  BLOCK_ALLOC : the size to extends

INPUTS/OUTPUTS

  list        : the character list

SOURCE

430       subroutine extendListNa(list,BLOCK_ALLOC)
431         character(LEN=LEN_DEF) ,dimension(:) ,pointer  ,intent(inout) :: list
432         integer                                        ,intent(in)    :: BLOCK_ALLOC
433 
434         character(LEN=LEN_DEF) ,dimension(:) ,pointer  :: listTemp => null()
435 
436         !start
437         if ( .not. associated(list) ) then
438            allocate (list(BLOCK_ALLOC))
439            return
440         end if
441 
442         allocate ( listTemp(size(list)) )
443         listTemp = list
444         deallocate ( list )
445         allocate ( list(size(listTemp)+BLOCK_ALLOC) )
446         list(:size(listTemp)) = listTemp
447         deallocate (listTemp)
448 
449         return
450      end subroutine extendListNa

extendListTrio

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    extendListTrio

DESCRIPTION

   increase the list of type type(TRIO_TYPE) with BLOCK_ALLOC values

INPUTS

  BLOCK_ALLOC : the size to extends

INPUTS/OUTPUTS

  list        : the type(TRIO_TYPE) list

SOURCE

466       subroutine extendListTrio(list,BLOCK_ALLOC)
467         type(TRIO_TYPE) ,dimension(:) ,pointer  ,intent(inout) :: list
468         integer                                        ,intent(in)    :: BLOCK_ALLOC
469 
470         type(TRIO_TYPE) ,dimension(:) ,pointer  :: listTemp => null()
471         integer :: i
472 
473         !start
474         if ( .not. associated(list) ) then
475            allocate (list(BLOCK_ALLOC))
476            return
477         end if
478 
479         allocate ( listTemp(size(list)) )
480         listTemp = list
481         deallocate ( list )
482         allocate ( list(size(listTemp)+BLOCK_ALLOC) )
483         list(:size(listTemp)) = listTemp
484         deallocate (listTemp)
485 
486         return
487      end subroutine extendListTrio

get_listProgenies

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    get_listProgenies

DESCRIPTION

   get the list of progeniture which the sire or the dam is the animal indexed by ijkd

INPUTS

 dataset     : anima-triplet list
 ijkd        : index of the parental progeny
 nameAnimal  : name of the parental progeny

OUTPUTS

 sizeList    : size of the result list
 list        : list of progenies

SOURCE

215      subroutine get_listProgenies(dataset,ijkd,nameAnimal,sizeList,list)
216         type (GENEALOGY_TYPE)         ,intent(in)     :: dataset
217         integer                       , intent(in)    :: ijkd
218         character(len=LEN_DEF)        , intent(in)    :: nameAnimal
219         integer                      , intent(out)    :: sizeList
220         type(TRIO_TYPE) ,dimension(:) , intent(out) , pointer  :: list
221 
222         integer  ,  parameter :: BLOCK_ALLOC = 100
223         integer :: i,id
224 
225         list =>null()
226 
227         id = dataset%corranim(ijkd)
228 
229         if ( id <=0 .or. id > size(dataset%corranim)) then
230             print *,'(get_listProgenies) can not get corresponding index of animal :',ijkd,' name:',nameAnimal
231             stop
232         end if
233 
234         call extendListTrio(list,BLOCK_ALLOC)
235         sizeList = 0
236         do i=1,dataset%na
237             if ( dataset%trioKd(i)%idSire == id .or.  dataset%trioKd(i)%idDam == id ) then
238               if ( sizeList >= size(list) ) then
239                 call extendListTrio(list,BLOCK_ALLOC)
240               end if
241               sizeList = sizeList + 1
242               list(sizeList)=dataset%trioKd(i)
243             end if
244         end do
245 
246      end subroutine get_listProgenies

init_perf_animal

[ Top ] [ m_qtlmap_cd ] [ Subroutines ]

NAME

    init_perf_animal

DESCRIPTION

   Add performance and cd information about a progeniture contents in the animal-triplet list

INPUTS

 progeny   : name of the progeny
 ncar      : number of trait (dimension of listPerf,listCdt)
 listPerf  : phenotypic values
 listCdt   : cd values

INPUTS/OUTPUTS

  dataset : anima-triplet list

SOURCE

178      subroutine init_perf_animal(dataset,progeny,listPerf,listCdt,ncar)
179        type (GENEALOGY_TYPE)      ,intent(inout)     :: dataset
180        character(len=LEN_DEF)        , intent(in)    :: progeny
181        integer                    ,intent(in)        :: ncar
182        real(kind=dp)  ,dimension(ncar) , intent(in)  :: listPerf,listCdt
183 
184        integer :: i,j
185 
186 
187        do i=1,dataset%na
188          if (trim(dataset%listNA(i))  == trim(progeny)) then
189              dataset%trioKd(i)%presentc=.false.
190              dataset%trioKd(i)%perfKd(:ncar) = listPerf(:ncar)
191              do j=1,ncar
192                 if ( listCdt(j) /= 0 ) dataset%trioKd(i)%presentc(j)=.true.
193              end do
194              return
195          end if
196        end do
197      end subroutine init_perf_animal