m_qtlmap_phase_offspring

[ Top ] [ Modules ]

NAME

    m_qtlmap_phase_offspring

DESCRIPTION

NOTES

SEE ALSO


haplotype_offspring_v1

[ Top ] [ m_qtlmap_phase_offspring ] [ Subroutines ]

NAME

    haplotype_offspring_v1

DESCRIPTION

  Prediction of haplotypes in offspring:
  assuming a known sire phase 
  predicting the phase of offspring using flanking markers for marker with a unknown phase
  (inspired from the linkphase software)

NOTES

   creation : 14/12/2010 : carole moreno

SOURCE

 42 subroutine haplotype_offspring_v1
 43  implicit none
 44 !A paramétricer par OLIVIER --> message de CMO
 45 
 46 !A paramétricer par OLIVIER --> message de CMO
 47 integer             :: mktot1, mktot2
 48 real   (kind=dp)     :: ph1,ph2, ph
 49  integer :: i, j, kkd, lk, mk1, mk2, ori1, ori2, k1, k2,c, ip,jm, k,kd,nd1, nd2, nm1, nm2
 50   logical  , dimension (:,:),   allocatable        :: haplotyped, genotyped
 51  integer  :: stat
 52 
 53  allocate (haplotyped(size(numero),maxval(nmk)),STAT=stat)
 54  call check_allocate(stat,'haplotyped')
 55  allocate (genotyped(size(numero),maxval(nmk)),STAT=stat)
 56  call check_allocate(stat,'genotyped')
 57  
 58  allocate (reconstructed(nchr,size(numero),maxval(nmk)),STAT=stat)
 59  call check_allocate(stat,'reconstructed')
 60 
 61  DO c=1,nchr
 62 !A paramétricer par OLIVIER --> message de CMO
 63   mktot1=1
 64   mktot2=nmk(c)
 65 !A paramétricer par OLIVIER --> message de CMO
 66   haplotyped=.false.
 67   genotyped=.false.
 68   reconstructed=.false.
 69   do ip=1,np
 70          !print *, ip, ' pere', ' ', correp(ip)
 71     DO k=mktot1,mktot2 !boucle sur les marqueurs
 72 ! HYPOTHESE: on suppose que genotyp est phasé pour les pères
 73 ! Definition des variables: * genotyped pour les pères 
 74 !               * haplotyped pour les pères
 75       if (pheno(c,k,correp(ip),1)/=nmanque.and.pheno(c,k,correp(ip),2)/=nmanque            )  genotyped(correp(ip),k) =.true.
 76       if (genotyp(c,k,correp(ip),1)/=nmanque.and.genotyp(c,k,correp(ip),2)/=nmanque        )  haplotyped(correp(ip),k)=.true.
 77     ENDDO
 78     nm1=nmp(ip)+1
 79     nm2=nmp(ip+1)
 80     do jm=nm1,nm2
 81       nd1=ndm(jm)+1
 82       nd2=ndm(jm+1)
 83       do kd=nd1,nd2
 84         ph=0.d0;ph1=0.d0;ph2=0.d0
 85         if (corred(kd) == 9999 ) cycle
 86         genotyp(c,:,corred(kd),:)=nmanque
 87         IF(count(presentg(:,kd))>0) then !condition pour sélectionner les individus génotypés
 88 
 89 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 90 !!1) PHASAGE DES PERES ET DES DESCENDANTS LORSQUE LE DESCENDANT ET LE PERE SONT GENOTYPES
 91 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 92 
 93   DO k=mktot1,mktot2 !boucle sur les marqueurs
 94 ! Definition de la variable genotyped pour les descendants
 95    if (pheno(c,k,corred(kd),1)/=nmanque.and.pheno(c,k,corred(kd),2)/=nmanque) genotyped(corred(kd),k)=.true. 
 96    !print *, k, animal(kd), pheno(c,k,i,1),pheno(c,k,i,2), pere(ip), pheno(c,k,correp(ip),1), pheno(c,k,correp(ip),2)
 97 
 98    IF (genotyped(corred(kd),k).and.genotyped(correp(ip),k)) then ! condition sur les couple père-des génotypés
 99  ! Le père est homozygote
100       if ((pheno(c,k,correp(ip),1)==pheno(c,k,correp(ip),2)) &
101      .and.(pheno(c,k,correp(ip),1)==pheno(c,k,corred(kd),1)))  then
102           genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),1)
103           genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),2)
104           haplotyped(corred(kd),k)=.true.
105           genotyp(c,k,correp(ip),1)=pheno(c,k,correp(ip),1)
106           genotyp(c,k,correp(ip),2)=pheno(c,k,correp(ip),2)
107           haplotyped(correp(ip),k)=.true.
108 ! print *, animal(kd), k,'OK'
109       endif
110       if ((pheno(c,k,correp(ip),1)==pheno(c,k,correp(ip),2)) &
111      .and.(pheno(c,k,correp(ip),1)==pheno(c,k,corred(kd),2)))  then
112           genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),2)
113           genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),1)
114           haplotyped(corred(kd),k)=.true.
115           genotyp(c,k,correp(ip),1)=pheno(c,k,correp(ip),2)
116           genotyp(c,k,correp(ip),2)=pheno(c,k,correp(ip),1)
117           haplotyped(correp(ip),k)=.true.
118       endif
119 !LE père est héterozygote l'allèle de l'haplotype 1 ou 2 est transmis
120       if     ( (genotyp(c,k,correp(ip),1)/=genotyp(c,k,correp(ip),2)        )  &
121         .and.  (pheno(c,k,corred(kd),1)==pheno(c,k,corred(kd),2)) &
122         .and.  (genotyp(c,k,correp(ip),1)==pheno(c,k,corred(kd),1    ))) then
123           genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),1)
124           genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),2)
125           haplotyped(corred(kd),k)=.true.
126       endif
127       if    (  (genotyp(c,k,correp(ip),1)/=genotyp(c,k,correp(ip),2) )   &
128         .and.  (pheno(c,k,corred(kd),1)==pheno(c,k,corred(kd),2)) &
129         .and.  (genotyp(c,k,correp(ip),2)==pheno(c,k,corred(kd),2  ) )) then  
130           genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),2)
131           genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),1)
132           haplotyped(corred(kd),k)=.true.
133       endif
134 
135    ENDIF ! FIN de la condition sur les couple père-des génotypés
136   ENDDO ! FIn de la boucle sur les marqueurs
137 
138  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139 !! 2) PHASAGE POUR LES MARQUEURS sans phase en fonction des phases des marqueurs voisins
140 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141 
142  DO k=mktot1,mktot2 ! boucle marker mktot1 mktot2
143     IF(.not.haplotyped(corred(kd),k))then
144      mk1=0;mk2=0;ori1=0;ori2=0
145      
146 ! search for informative marker on the left side   
147      if (k==mktot1) then
148          mk1=0
149      else
150      do k1=k-1,mktot1,-1
151 !       if(corred(kd)==5)  print *, corred(kd), k, k1,haplotyped(corred(kd),k1), haplotyped(correp(ip),k1)
152         if(haplotyped(corred(kd),k1) .and. haplotyped(correp(ip),k1)) then
153          if(genotyp(c,k1,correp(ip),1)/=genotyp(c,k1,correp(ip),2))then ! search for informative marker
154            if(genotyp(c,k1,corred(kd),1)==genotyp(c,k1,correp(ip),1))ori1=1
155            if(genotyp(c,k1,corred(kd),1)==genotyp(c,k1,correp(ip),2))ori1=2
156            mk1=k1;exit
157          endif
158         endif
159        enddo
160      endif
161 
162 ! search for informative marker on the right side   
163      if (k==mktot2) then
164          mk2=0
165      else
166        do k2=k+1,mktot2
167         if(haplotyped(corred(kd),k2) .and. haplotyped(correp(ip),k2)) then
168          if(genotyp(c,k2,correp(ip),1)/=genotyp(c,k2,correp(ip),2))then ! search for informative marker
169            if(genotyp(c,k2,corred(kd),1)==genotyp(c,k2,correp(ip),1))ori2=1
170            if(genotyp(c,k2,corred(kd),1)==genotyp(c,k2,correp(ip),2))ori2=2
171            mk2=k2;exit
172          endif
173         endif
174        enddo
175      endif
176 
177 !! reconstruction of phases for markers with unknown phase using flanking phased markers
178      if(mk1/=0 .and. mk2/=0)then !mk1, mk2
179         if(ori1==ori2)then ! no recombination
180          ph=(1.00-(xaldane(posim(c,k)-posim(c,k1))))*(1.00-xaldane(posim(c,k2)-posim(c,k)))/(1.00-xaldane(posim(c,k2)-posim(c,k1)))
181          if(ph>PROB_PHASE_DESC)then
182            reconstructed(c,corred(kd),k)=.true.
183            if(genotyp(c,k,correp(ip),ori1)==pheno(c,k,corred(kd),1)) then
184              genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),1);genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),2)
185            else
186              genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),1);genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),2)
187          endif
188         endif
189        endif
190 
191       else if(mk1/=0 .and. mk2==0)then
192         ph1=1.00-xaldane(posim(c,k)-posim(c,k1))
193         if(ph1>PROB_PHASE_DESC)then
194           reconstructed(c,corred(kd),k)=.true.
195          if(genotyp(c,k,correp(ip),ori1)==pheno(c,k,corred(kd),1))then
196            genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),1);genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),2)
197          else
198            genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),1);genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),2)
199          endif
200         endif
201 
202        else if(mk1==0 .and. mk2/=0)then
203          ph2=1.00-xaldane(posim(c,k2)-posim(c,k))
204          if(ph2>PROB_PHASE_DESC)then
205           reconstructed(c,corred(kd),k)=.true.
206           if(genotyp(c,k,correp(ip),ori2)==pheno(c,k,corred(kd),1))then
207             genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),1);genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),2)
208           else
209             genotyp(c,k,corred(kd),2)=pheno(c,k,corred(kd),1);genotyp(c,k,corred(kd),1)=pheno(c,k,corred(kd),2)
210           endif
211         endif
212 
213       endif !mk1, mk2
214 
215      ENDIF ! not.haplotyped
216 
217   END DO
218 
219 ENDIF! fin de la condition pour sélectionner les individus génotypés
220 ENDDO ! end animal
221 ENDDO ! end mere
222 ENDDO ! end pere
223 ENDDO ! end chromosome
224 
225  deallocate(haplotyped)
226  deallocate(genotyped)
227 
228 end subroutine haplotype_offspring_v1