c MCMC code for Birth-Death-Mutation model c Code written in Fortran 95 c c Code for fixed population size NN (NN=10000). c c To run code:- c gfortran -o run BDM.f c ./run c c replace gfortran by appropriate Fortran compiler. c parameter(NN=10000,Nrun=1100000,ND=100,NS=100) integer cia, cib, cic, cid, cie, TOT, G, gpig integer dia, dib double precision UP(200000), UPZ(200000), UQ(326), UQZ(326) double precision UR(200000), URZ(200000) integer UZZ(200000), UZN(200000) double precision y, z, b, d, m, JJ(3) double precision ppig, pb, pd, pm double precision dig, UG(326), prob, proz double precision ca, cb, cc, cd, ce, cf, cg, NB double precision pa, pc, pe, pf, cow, ch, pi double precision ba, bb, bc, sigb, sigd integer X(NN), XO(NN), U(326), MX, GX(30), CONN(4) integer count, runz, tock, nv, cox, CONX(4) integer XS(NN), XT(NN), cop, coa, cob, caz open(14,file='BDMoutput.txt') c Output from the BDM algorithm. call init_random_seed() pi=3.14159265358979323846 caz=0 sigb=0.025 sigd=0.025 c Proposal standard deviations for RWM U(1)=30 U(2)=23 U(3)=15 U(4)=10 U(5)=8 U(6)=5 U(7)=5 do 318 i=8,11 U(i)=4 318 continue do 328 i=12,24 U(i)=3 328 continue do 338 i=25,44 U(i)=2 338 continue do 348 i=45,326 U(i)=1 348 continue MX=326 c Data input (326 total number of genotypes. Code could easily run on different data by changing the above. NB=NN cow=0 do 876 i=1,473 ca=i cb=NN+1 cb=cb-ca cow=cow+log(cb) 876 continue write(*,*) cow c cow is a term used in the log-likelihood count=0 cop=0 coa=0 cob=0 c The above are counts for checking acceptance rates runz=0 tock=0 507 ppig=-5000 do 510 i=1,200000 call random_number(UP(i)) UPZ(i)=UP(i) call random_number(UR(i)) URZ(i)=UR(i) 510 continue do 520 i=1,MX call random_number(UQ(i)) UQZ(i)=UQ(i) 520 continue c The above generates U, W and V for the paper (UP, UR, UQ). UPZ, etc are used when updating. b=0.65 d=0.2 m=1-b-d prob=log(b)-log(b+d)-cow c Initital parameter values and initial conditioning birth occurs before death X(1)=2 g=1 tot=2 c Initial population configuration after first birth. cox=0 c Event count do 301 i=1,4 CONN(i)=0 301 continue CONN(4)=1 37 cox=cox+1 ca=UP(cox) cib=floor(ca*tot)+1 do 47 i=1,g cib=cib-X(i) if (cib.le.0) then cic=i goto 57 endif 47 continue 57 cb=UR(cox) if (cb.le.b) then X(cic)=X(cic)+1 tot=tot+1 CONN(1)=CONN(1)+1 UZZ(cox)=1 elseif ((cb.gt.b).and.(cb.le.(b+d))) then X(cic)=X(cic)-1 tot=tot-1 CONN(2)=CONN(2)+1 if (X(cic).eq.0) then X(cic)=X(g) g=g-1 endif if (tot.eq.1) then prob=prob+log(b)-log(b+d) X(1)=2 tot=2 CONN(4)=CONN(4)+1 endif UZZ(cox)=2 else CONN(3)=CONN(3)+1 if (X(cic).gt.1) then X(cic)=X(cic)-1 g=g+1 X(g)=1 endif UZZ(cox)=3 endif UZN(cox)=UZZ(cox) if((tot.gt.0).and.(tot.lt.NN)) goto 37 c Simulation of process if(g.lt.MX) goto 507 c if number of genotypes is less than 326 simulation is rejected. do 67 i=1,g cia=1 do 77 j=2,g if (X(j).gt.X(cia)) then cia=j endif 77 continue XO(i)=X(cia) XT(i)=XO(i) X(cia)=0 67 continue do 187 i=1,MX if (XO(i).lt.U(i)) then goto 507 endif 187 continue GX(1)=g do 87 i=2,30 do 97 j=1,g if (XO(j).lt.i) then cib=j-1 goto 107 endif 97 continue 107 GX(i)=cib 87 continue do 117 i=1,MX if (U(i).eq.1) then UG(i)=NN else UG(i)=0 do 127 j=1,GX(U(i)) UG(i)=UG(i)+XO(j) 127 continue endif 117 continue do 137 ia=1,MX prob=prob+log(UG(ia)) ca=UQ(ia) cia=floor(ca*UG(ia))+1 do 147 ja=1,GX(U(ia)) cia=cia-XO(ja) if (cia.le.0) then cib=ja goto 157 endif 147 continue 157 joe=1 if (U(ia).gt.1) then do 167 i=1,(U(ia)-1) ca=XO(cib)-i prob=prob+log(ca) NB=NB-1 167 continue endif do 177 i=1,MX UG(i)=UG(i)-XO(cib) 177 continue XO(cib)=0 137 continue c Does the sampling (473 individuals from x) and computes prob the log-likelihood. do 302 i=1,4 CONX(i)=CONN(i) 302 continue write(*,*) b, d, m, prob, CONN do 777 ix=1,Nrun caz=caz+1 count=count+1 517 call random_number(ba) call random_number(bb) ba=2*pi*ba ba=sin(ba) bc=0.5 bb=log(bb) bb=-2*bb bb=bb**bc pb=b+sigb*bb*ba call random_number(ba) call random_number(bb) ba=2*pi*ba ba=sin(ba) bc=0.5 bb=log(bb) bb=-2*bb bb=bb**bc pd=d+sigd*bb*ba pd=abs(pd) pm=1-pb-pd ca=0.5 if(pb.lt.ca) goto 517 if(pm.lt.0) goto 517 do 303 i=1,3 CONX(i)=0 303 continue CONX(4)=1 c Proposes values for b, d, m (birth, death, mutation) probabilities using RWM. proz=log(pb)-log(pb+pd)-cow c proz is the updated log-likelihood which will be compared with prob X(1)=2 g=1 tot=2 cox=0 33 cox=cox+1 ca=UP(cox) cib=floor(ca*tot)+1 do 43 i=1,g cib=cib-X(i) if (cib.le.0) then cic=i goto 53 endif 43 continue 53 cb=UR(cox) if (cb.le.pb) then X(cic)=X(cic)+1 tot=tot+1 CONX(1)=CONX(1)+1 UZN(cox)=1 elseif ((cb.gt.pb).and.(cb.le.(pb+pd))) then X(cic)=X(cic)-1 tot=tot-1 CONX(2)=CONX(2)+1 if (X(cic).eq.0) then X(cic)=X(g) g=g-1 endif if (tot.eq.1) then prob=prob+log(pb)-log(pb+pd) X(1)=2 tot=2 CONX(4)=CONX(4)+1 endif UZN(cox)=2 else CONX(3)=CONX(3)+1 if (X(cic).gt.1) then X(cic)=X(cic)-1 g=g+1 X(g)=1 endif UZN(cox)=3 endif if((tot.gt.0).and.(tot.lt.NN)) goto 33 if(g.lt.MX) goto 193 do 63 i=1,g cia=1 do 73 j=2,g if (X(j).gt.X(cia)) then cia=j endif 73 continue XO(i)=X(cia) XT(i)=XO(i) X(cia)=0 63 continue do 183 i=1,MX if (XO(i).lt.U(i)) then goto 193 endif 183 continue GX(1)=g do 83 i=2,30 do 93 j=1,g if (XO(j).lt.i) then cib=j-1 goto 103 endif 93 continue 103 GX(i)=cib 83 continue do 113 i=1,MX if (U(i).eq.1) then UG(i)=NN else UG(i)=0 do 123 j=1,GX(U(i)) UG(i)=UG(i)+XO(j) 123 continue endif 113 continue do 133 ia=1,MX proz=proz+log(UG(ia)) ca=UQ(ia) cia=floor(ca*UG(ia))+1 do 143 ja=1,GX(U(ia)) cia=cia-XO(ja) if (cia.le.0) then cib=ja goto 153 endif 143 continue 153 joe=1 if (U(ia).gt.1) then do 163 i=1,(U(ia)-1) ca=XO(cib)-i proz=proz+log(ca) NB=NB-1 163 continue endif do 173 i=1,MX UG(i)=UG(i)-XO(cib) 173 continue XO(cib)=0 133 continue goto 203 193 proz=-5000 c The above simulates the BDM process and the sampling. 203 pa=exp(proz-prob) call random_number(pc) c Calculate acceptamnce probability if (pc.lt.pa) then b=pb d=pd m=pm prob=proz cop=cop+1 do 305 i=1,4 CONN(i)=CONX(i) 305 continue do 310 i=1,cox UZZ(i)=UZN(i) 310 continue else do 315 i=1,4 CONX(i)=CONN(i) 315 continue do 320 i=1,cox UZN(i)=UZZ(i) 320 continue endif c Update parameters as appropriate. do 530 i=1,4000 call random_number(ca) cib=floor(50*ca)+1 cic=i-1 cic=cic*50 call random_number(UPZ(cic+cib)) call random_number(ca) cib=floor(50*ca)+1 call random_number(URZ(cic+cib)) 530 continue c Updating of U and W, 2% of values. proz=log(b)-log(b+d)-cow X(1)=2 g=1 tot=2 cox=0 do 321 i=1,3 CONX(i)=0 321 continue CONX(4)=1 34 cox=cox+1 ca=UPZ(cox) cib=floor(ca*tot)+1 do 44 i=1,g cib=cib-X(i) if (cib.le.0) then cic=i goto 54 endif 44 continue 54 cb=URZ(cox) if (cb.le.b) then X(cic)=X(cic)+1 tot=tot+1 CONX(1)=CONX(1)+1 UZN(cox)=1 elseif ((cb.gt.b).and.(cb.le.(b+d))) then X(cic)=X(cic)-1 tot=tot-1 CONX(2)=CONX(2)+1 if (X(cic).eq.0) then X(cic)=X(g) g=g-1 endif if (tot.eq.1) then proz=proz+log(b)-log(b+d) X(1)=2 tot=2 CONX(4)=CONX(4)+1 endif UZN(cox)=2 else CONX(3)=CONX(3)+1 if (X(cic).gt.1) then X(cic)=X(cic)-1 g=g+1 X(g)=1 endif UZN(cox)=3 endif if((tot.gt.0).and.(tot.lt.NN)) goto 34 if(tot.eq.0) goto 194 if(g.lt.MX) goto 194 do 64 i=1,g cia=1 do 74 j=2,g if (X(j).gt.X(cia)) then cia=j endif 74 continue XO(i)=X(cia) XT(i)=XO(i) X(cia)=0 64 continue do 184 i=1,MX if (XO(i).lt.U(i)) then goto 194 endif 184 continue GX(1)=g do 84 i=2,30 do 94 j=1,g if (XO(j).lt.i) then cib=j-1 goto 104 endif 94 continue 104 GX(i)=cib 84 continue do 114 i=1,MX if (U(i).eq.1) then UG(i)=NN else UG(i)=0 do 124 j=1,GX(U(i)) UG(i)=UG(i)+XO(j) 124 continue endif 114 continue do 134 ia=1,MX proz=proz+log(UG(ia)) ca=UQ(ia) cia=floor(ca*UG(ia))+1 do 144 ja=1,GX(U(ia)) cia=cia-XO(ja) if (cia.le.0) then cib=ja goto 154 endif 144 continue 154 joe=1 if (U(ia).gt.1) then do 164 i=1,(U(ia)-1) ca=XO(cib)-i proz=proz+log(ca) NB=NB-1 164 continue endif do 174 i=1,MX UG(i)=UG(i)-XO(cib) 174 continue XO(cib)=0 134 continue goto 204 194 proz=-5000 204 pa=exp(proz-prob) call random_number(pc) if (pc.lt.pa) then do 214 i=1,200000 UP(i)=UPZ(i) UR(i)=URZ(i) 214 continue prob=proz coa=coa+1 do 325 i=1,4 CONN(i)=CONX(i) 325 continue do 330 i=1,cox UZZ(i)=UZN(i) 330 continue else do 224 i=1,200000 UPZ(i)=UP(i) URZ(i)=UR(i) 224 continue do 335 i=1,4 CONX(i)=CONN(i) 335 continue do 340 i=1,cox UZN(i)=UZZ(i) 340 continue endif c The following updates V (UQ) cia=0 cb=0 do 351 i=1,3 JJ(i)=0 do 352 ia=1,(CONN(i)+1) call random_number(ca) JJ(i)=JJ(i)-log(ca) 352 continue cb=cb+JJ(i) 351 continue pb=JJ(1)/cb pd=JJ(2)/cb pm=1-pb-pd cc=CONN(4) cd=cc*(log(pb)-log(pb+pd)) cd=cd-cc*(log(b)-log(b+d)) call random_number(ce) cf=exp(cd) cic=0 if (ce.le.cf) then b=pb d=pd m=pm prob=prob+cd do 353 i=1,cox call random_number(cg) if (UZZ(i).eq.1) then UR(i)=cg*b elseif (UZZ(i).eq.2) then UR(i)=b+d*cg elseif (UZZ(i).eq.3) then UR(i)=b+d+cg*m endif URZ(i)=UR(i) 353 continue endif do 540 i=1,10 call random_number(ca) cic=floor(MX*ca)+1 call random_number(UQZ(cic)) 540 continue proz=log(b)-log(b+d)-cow X(1)=2 g=1 tot=2 cox=0 35 cox=cox+1 ca=UP(cox) cib=floor(ca*tot)+1 do 45 i=1,g cib=cib-X(i) if (cib.le.0) then cic=i goto 55 endif 45 continue 55 cb=UR(cox) if (cb.le.b) then X(cic)=X(cic)+1 tot=tot+1 elseif ((cb.gt.b).and.(cb.le.(b+d))) then X(cic)=X(cic)-1 tot=tot-1 if (X(cic).eq.0) then X(cic)=X(g) g=g-1 endif if (tot.eq.1) then proz=proz+log(b)-log(b+d) X(1)=2 tot=2 endif else if (X(cic).gt.1) then X(cic)=X(cic)-1 g=g+1 X(g)=1 endif endif if((tot.gt.0).and.(tot.lt.NN)) goto 35 if(g.lt.MX) goto 195 do 65 i=1,g cia=1 do 75 j=2,g if (X(j).gt.X(cia)) then cia=j endif 75 continue XO(i)=X(cia) XT(i)=XO(i) X(cia)=0 65 continue do 185 i=1,MX if (XO(i).lt.U(i)) then goto 195 endif 185 continue GX(1)=g do 85 i=2,30 do 95 j=1,g if (XO(j).lt.i) then cib=j-1 goto 105 endif 95 continue 105 GX(i)=cib 85 continue do 115 i=1,MX if (U(i).eq.1) then UG(i)=NN else UG(i)=0 do 125 j=1,GX(U(i)) UG(i)=UG(i)+XO(j) 125 continue endif 115 continue do 135 ia=1,MX proz=proz+log(UG(ia)) ca=UQZ(ia) cia=floor(ca*UG(ia))+1 do 145 ja=1,GX(U(ia)) cia=cia-XO(ja) if (cia.le.0) then cib=ja goto 155 endif 145 continue 155 joe=1 if (U(ia).gt.1) then do 165 i=1,(U(ia)-1) ca=XO(cib)-i proz=proz+log(ca) NB=NB-1 165 continue endif do 175 i=1,MX UG(i)=UG(i)-XO(cib) 175 continue XO(cib)=0 135 continue goto 205 195 proz=-5000 205 pa=exp(proz-prob) call random_number(pc) if (pc.lt.pa) then do 215 i=1,MX UQ(i)=UQZ(i) 215 continue prob=proz cob=cob+1 else do 225 i=1,MX UQZ(i)=UQ(i) 225 continue endif C Stores every NS value (thinning) and prints to screen every ND value. (In this case these are the same.) c if (caz.eq.NS) then write(14,*) b, d, m c caz=0 c endif if (count.eq.ND) then write(*,*) ix, b, d, m, prob write(*,*) cop, coa, cob count=0 endif 777 continue end c Code for setting initial random seed c Downloaded from: www4.ncsu.edu/~cbkenned/rndm.doc SUBROUTINE init_random_seed() INTEGER :: i, n, clock INTEGER, DIMENSION(:), ALLOCATABLE :: seed CALL RANDOM_SEED(size = n) ALLOCATE(seed(n)) CALL SYSTEM_CLOCK(COUNT=clock) seed = clock + 37 * (/ (i - 1, i = 1, n) /) CALL RANDOM_SEED(PUT = seed) DEALLOCATE(seed) END SUBROUTINE