C C Simulate quantitative trait due to multiple unlinked QTLs for C nuclear families including twins C C Marker information is available for the first QTL which may C be X-linked. The other oligenes are autosomal C A recombination can occur 0 or 1 times within an interval on that C chromosome C C 1|2 1|1 1 2|1 C 5|3 3|3 2 3|3 C D|d d|d 3 ---------------------> D|d C 2|2 1|3 4 2|3 C 3|4 4|1 5 3|4 C integer CSTRM,OSTR parameter(CSTRM=4,OSTR=7) C C linelength, max no of alleles for marker, max words per line, C largest pedigree, missing values, max number of markers, C max no of QTLs C integer LINSIZ, MAXALL, MAXCAT, MAXCOL, MAXSIZ, MISS, MAXLOC, & MAXQTL parameter(LINSIZ=1024, MAXALL=10, MAXCAT=10, MAXCOL=512, & MAXSIZ=20, MISS=-9999, MAXLOC=60, MAXQTL=10) C C nfam Total number of nuclear families, C nmzfam Number of families containing MZ twins C nhsfam Number of families containing halfsibs C nsibs Number of children per family C integer nfam, nhsfam, nmzfam, nsibs C C Variance components C integer ncat double precision mu, prev(MAXCAT), thresh(MAXCAT) double precision va, vc, vq, vqd, ve, vt C Pedigree structure integer pedigree integer num,nfound,id(MAXSIZ),fa(MAXSIZ),mo(MAXSIZ) integer sex(MAXSIZ), twin(MAXSIZ) C C Trait loci -- map position of first QTL (posqtl=rank order, mapqtl=cM) double precision pqtl(2) integer posqtl double precision mapqtl double precision mu1, mu2, mu3, mmu1, mmu2 integer nqtl, qtl(MAXSIZ, MAXQTL, 2) double precision trait(MAXSIZ) integer showqtl C C Markers integer nmark, nall double precision cumfrq(MAXALL) integer set(MAXSIZ,MAXLOC,2) C C Ascertainment integer n_high, n_low double precision hi_thresh, lo_thresh C C Meiosis double precision map(MAXLOC), theta(MAXLOC) C C Commands character*3 keyword, keyw2 character*(LINSIZ) lin character*20 words(MAXCOL) character*80 bigwor(MAXCOL/4) equivalence (words(1),bigwor(1)) character*80 locfil, infil, pedfil character*8 version C integer i,ilevel,nlin,t0, typ logical echo, filexist, prompt, xlinkd character*7 mtyp(2) double precision frq, het, tot C C random number generator seeds integer ix,iy,iz common /rand/ ix,iy,iz C C functions C intrinsic isatty, time integer eow, ival, time logical isatty double precision fval, ppnd external fdate character*24 fdate C C main loop -- reads stdin, parses C echo=.false. ilevel=1 nlin=0 prompt=isatty(5) showqtl=0 t0=time() ix=mod(t0,29282) if (t0.gt.29282) then iy=mod(t0/29282,29282) else iy=mod(419*ix+6173,29282) end if iz=mod(419*iy+6173,29282) mtyp(1)='marker ' mtyp(2)='xmarker' infil=' ' locfil=' ' pedfil=' ' version='20040817' C stamp and initialize C write(*,'(a/2a/a/2a/)') 2 '|||| TWINSIM: A program for simulating genetic data', 3 '|\\/| Version: ',version, 4 '|/\\| Author : David L Duffy (c) 1995-2004', 5 '|||| Job run: ',fdate() 999 continue xlinkd=.FALSE. typ=1 n_high=0 n_low=0 hi_thresh=0.0d0 lo_thresh=0.0d0 nfam=1000 nsibs=2 nmzfam=500 nhsfam=0 nqtl=4 nmark=4 nall=4 call mkfreq(nall, cumfrq, het) map(1)=0.0d0 map(2)=10.0d0 map(3)=20.0d0 map(4)=30.0d0 mapqtl=15.0d0 posqtl=0 call place(mapqtl,posqtl,nmark,nloci,map,theta) mu=0.0d0 pqtl(1)=0.5d0 pqtl(2)=1.0d0 ncat=MISS prev(1)=0.0d0 thresh(1)=MISS va=0.0d0 vc=0.0d0 vq=0.5d0 vqd=0.0d0 ve=0.5d0 mu1=-1.0d0 mu2=0.0d0 mu3=1.0d0 mmu1=-1.0d0 mmu2=1.0d0 C C read-eval loop C 1 continue if (prompt) write(*,'(/a,$)') '>> ' i=1 3 continue if (ilevel.eq.1) then read(*,'(a)',end=2) lin(i:LINSIZ) else read(CSTRM,'(a)',end=2) lin(i:LINSIZ) end if nlin=nlin+1 i=eow(lin) if (i.gt.1 .and. lin(i-1:i).eq.' \\') goto 3 if (echo) then write(*,'(/2a/)') '-> ',lin(1:75) end if narg=MAXCOL call args(lin,narg,words,1) keyword=words(1)(1:3) keyw2=words(2)(1:3) C C Locus description if (keyword.eq.'set' .and. keyw2.eq.'nal') then if (words(3)(1:3).eq.'per') then nall=0 else nall=fval(words(3)) end if call mkfreq(nall, cumfrq, het) elseif (keyword.eq.'set' .and. keyw2.eq.'fre') then nall=narg-2 tot=0.0d0 het=1.0d0 write(*,'(a,i2,a,$)') & 'Marker allele frequencies (Nall=',nall,'):' do 40 i=3,narg frq=fval(words(i)) het=het-frq*frq tot=tot+frq cumfrq(i-2)=tot write(*,'(1x,f5.3,$)') frq 40 continue write(*,*) cumfrq(nall)=1.0d0 elseif (keyword.eq.'set' .and. words(2)(1:2).eq.'va') then va=fval(words(3)) if (narg.ge.4) then vq=fval(words(4)) call qtomu(pqtl, vq, mu, mu1, mu2, mu3) end if if (narg.ge.5) then vc=fval(words(5)) end if if (narg.ge.6) then ve=fval(words(6)) end if call showvc(6,xlinkd,ncat,prev,mu,va,vq,vc,ve,vt, & pqtl,mu1,mu2,mu3,mmu1,mmu2) elseif (keyword.eq.'set' .and. words(2).eq.'vc') then vc=fval(words(3)) call showvc(6,xlinkd,ncat,prev,mu,va,vq,vc,ve,vt, & pqtl,mu1,mu2,mu3,mmu1,mmu2) elseif (keyword.eq.'set' .and. words(2).eq.'ve') then ve=fval(words(3)) call showvc(6,xlinkd,ncat,prev,mu,va,vq,vc,ve,vt, & pqtl,mu1,mu2,mu3,mmu1,mmu2) elseif (keyword.eq.'set' .and. words(2).eq.'vq') then vq=fval(words(3)) call qtomu(pqtl, vq, mu, mu1, mu2, mu3) call showvc(6,xlinkd,ncat,prev,mu,va,vq,vc,ve,vt, & pqtl,mu1,mu2,mu3,mmu1,mmu2) elseif (keyword.eq.'set' .and. keyw2.eq.'qtl') then mapqtl=fval(words(3)) if (narg.ge.4) pqtl(1)=fval(words(4)) if (narg.eq.5) vq=fval(words(5)) call qtomu(pqtl, vq, mu, mu1, mu2, mu3) call place(mapqtl,posqtl,nmark,nloci,map,theta) elseif (keyword.eq.'set' .and. keyw2.eq.'pen') then mu1=fval(words(3)) mu2=fval(words(4)) mu3=fval(words(5)) mmu1=fval(words(6)) mmu2=fval(words(7)) call mutoq(pqtl, mu1, mu2, mu3, mu, vq, vqd) elseif (keyword.eq.'set' .and. keyw2.eq.'xli') then xlinkd=.not.(narg.eq.3 .and. words(3).eq.'off') typ=1 if (xlinkd) then typ=2 write(*,'(a)') 'First QTL and markers now X-linked.' end if elseif (keyword.eq.'set' .and. (keyw2.eq.'pri' .or. & keyw2.eq.'out' .or. keyw2.eq.'ple')) then if (narg.eq.2) then words(3)='on' showqtl=1 elseif (words(3).eq.'on') then showqtl=1 else showqtl=ival(words(3)) end if write(*,'(a,i4)') 'Print QTL genos = ',showqtl elseif (keyword.eq.'set' .and. keyw2.eq.'pre') then ncat=narg-2 if (ncat.lt.2) then ncat=2 prev(2)=fval(words(3)) prev(1)=1.0d0-prev(2) else do 55 i=1, ncat prev(i)=fval(words(i+2)) 55 continue end if mu=0.0d0 do 56 i=1,ncat-1 mu=mu+prev(i) thresh(i)=ppnd(mu) 56 continue hi_thresh=1.5d0 lo_thresh=1.5d0 mu=0.0d0 if (ncat.eq.2) then write(*,'(a,f8.3)') 'Prevalence = ',prev(2) write(*,'(a,f8.3)') 'Threshold = ',thresh(1) else write(*,'(a,10(1x,f8.3):)') & 'Proportions = ',(prev(i), i=1, ncat) write(*,'(a,9(1x,f8.3):)') & 'Thresholds = ',(thresh(i), i=1, ncat-1) end if elseif (keyword.eq.'set' .and. keyw2.eq.'mea') then mu=fval(words(3)) if (mu.eq.MISS) mu=0.0d0 write(*,'(a,f8.3)') 'Mean = ',mu elseif (keyword.eq.'set' .and. keyw2.eq.'nqt') then nqtl=ival(words(3)) if (nqtl.gt.MAXQTL) nqtl=MAXQTL elseif (keyword.eq.'set' .and. keyw2.eq.'map') then nmark=narg-2 do 70 i=1,nmark map(i)=fval(words(i+2)) 70 continue posqtl=0 call place(mapqtl,posqtl,nmark,nloci,map,theta) elseif ((keyword.eq.'sho' .and. keyw2.eq.'map') & .or. keyword.eq.'ls' .or. keyword.eq.'lis') then call showvc(6,xlinkd,ncat,prev,mu,va,vq,vc,ve,vt, & pqtl,mu1,mu2,mu3,mmu1,mmu2) call showmap(nloci, posqtl, pqtl, nall, het, map, theta) elseif (keyword.eq.'set' .and. keyw2.eq.'nfa') then nfam=ival(words(3)) write(*,'(a,i6)') 'No. of families = ',nfam elseif (keyword.eq.'set' .and. keyw2.eq.'nsi') then nsibs=ival(words(3)) write(*,'(a,i6)') 'No. of kids/fam = ',nsibs elseif (keyword.eq.'set' .and. keyw2(1:2).eq.'mz') then nmzfam=ival(words(3)) write(*,'(a,i6)') 'No. of MZ fams = ',nmzfam elseif (keyword.eq.'set' .and. keyw2.eq.'hal') then nhsfam=ival(words(3)) write(*,'(a,i6)') 'No. of HS fams = ',nhsfam elseif (keyword.eq.'set' .and. keyw2.eq.'hig') then n_high=ival(words(3)) if (narg.eq.4) hi_thresh=fval(words(4)) write(*,'(a,i6)') 'No. of probands = ',n_high write(*,'(a,f10.3)') 'Threshold = ',hi_thresh elseif (keyword.eq.'set' .and. keyw2.eq.'low') then n_low=ival(words(3)) if (narg.eq.4) lo_thresh=fval(words(4)) write(*,'(a,i6)') 'No. of probands = ',n_low write(*,'(a,f10.3)') 'Threshold = ',lo_thresh elseif (keyword.eq.'wri' .and. keyw2.eq.'loc') then locfil=words(3) if (narg.eq.4) pedfil=words(4) if (locfil.eq.' ') locfil='twinsim.in' if (pedfil.eq.' ') pedfil='twinsim.ped' write(*,'(3a)') 'Writing Sib-pair control file "', & locfil(1:eow(locfil)),'"' open(OSTR,file=locfil,status='unknown') call showvc(OSTR,xlinkd,ncat,prev,mu,va,vq,vc,ve,vt, & pqtl,mu1,mu2,mu3,mmu1,mmu2) write(OSTR,'(a)') 'set locus mztwin qua' if (ncat.eq.2) then write(OSTR,'(a)') 'set locus trait aff' else write(OSTR,'(a)') 'set locus trait qua' end if j=0 do 95 i=1,nloci if (i.ne.posqtl) then j=j+1 write(OSTR,'(a,i3.3,1x,a,1x,f7.2)') & 'set loc marker',j,mtyp(typ),map(i) else if (showqtl.eq.0) then write(OSTR,'(a,f7.2)') '# qtl at ',map(i) else write(OSTR,'(a,1x,a,1x,f7.2)') & 'set locus qtl',mtyp(typ),map(i) end if end if 95 continue if (showqtl.gt.1) then do 100 i=1,nqtl write(OSTR,'(a,i3.3,a,f7.2)') & 'set locus qtl',i,' marker ', 1000.0*i 100 continue end if write(OSTR,'(2a)') 'read pedigree ',pedfil write(OSTR,'(a)') 'set ple -1','run' close(OSTR,status='keep') elseif (keyword.eq.'sim' .or. keyword.eq.'wri') then pedfil=words(2) if (pedfil.eq.' ') pedfil='twinsim.ped' if (narg.ge.3) nfam=ival(words(3)) if (narg.ge.4) nsibs=ival(words(4)) if (narg.ge.5) nmzfam=ival(words(5)) write(*,'(/a,i6,a/a,i6/a,i6,a,f5.1,a)') 2 'Simulating ', nfam, ' nuclear families:', 3 ' No. children/fam = ',nsibs, 4 ' No. MZ families = ',nmzfam, 5 ' (',100.0d0*dfloat(nmzfam)/dfloat(nfam),'%)' if (nhsfam.gt.0) then write(*,'(/a,i6,a,f5.1,a)') 2 ' No. HS families = ',nhsfam, 3 ' (',100.0d0*dfloat(nhsfam)/dfloat(nfam),'%)' end if if (n_high.gt.0 .or. n_low.gt.0) then write(*,'(/a,i3,a/25x,a,i3,a)') 2 'Families selected to contain at least ',n_high, 3 ' high probands,','and at least ', n_low, ' low probands.' end if write(*,'(a,i6,a,f7.1,a/3a)') 2 ' No. markers = ',nmark, 3 ' (spanning ',map(nloci)-map(1),' cM)', 4 'Writing to "',pedfil(1:eow(pedfil)),'".' open(OSTR,file=pedfil,status='unknown') call simped(OSTR, xlinkd,posqtl,pqtl,nqtl, mu1, mu2, mu3, 2 mmu1, mmu2, mu, va, vc, ve, 3 ncat, thresh, n_high, hi_thresh, n_low, lo_thresh, 4 nfam, nmzfam, nhsfam, nsibs, nloci, 5 pedigree, num,nfound,id,fa,mo,sex,twin, 6 theta, nall, cumfrq, qtl, trait, set, showqtl) close(OSTR,status='keep') elseif (keyword.eq.'set'.and.keyw2.eq.'see') then ix=ival(words(3)) iy=ival(words(4)) iz=ival(words(5)) elseif (keyword.eq.'set'.and.keyw2.eq.'ech') then echo=.true. if (words(3).eq.'off') then echo=.false. end if elseif (keyword.eq.'set'.and.keyw2.eq.'pro') then prompt=.true. if (words(3).eq.'off') then prompt=.false. end if elseif (keyword.eq.'tim') then call stamp(t0) elseif (keyword(1:1).eq.'!' .or. keyword(1:1).eq.'#') then write(*,'(a)') lin(1:79) elseif (keyword(1:1).eq.'%'.or.keyword(1:1).eq.'$') then call shell(lin) elseif (keyword.eq.'inc' .and. ilevel.eq.1) then call args(lin,narg,bigwor,1) infil=bigwor(2) inquire(file=infil,exist=filexist) if (filexist) then write(*,'(/3a)') & 'Reading commands from "', infil(1:eow(infil)),'".' ilevel=ilevel+1 open(CSTRM,file=infil,status='unknown') else write(*,'(/3a/)') & 'ERROR: Unable to open "', infil(1:eow(infil)),'".' end if elseif (keyword.eq.'cle') then goto 999 elseif (keyword.eq.'hel') then write(*,'(a)') '* MODEL SETTINGS *', 2 'set mean ','set prevalence [...]', 3 'set va [vq [vc [ve]]]', 4 'set vc ', 5 'set ve ', 6 'set vq ' write(*,'(a)') 2 'set xlinked [on|off]', 3 'set pen [ ] ' write(*,'(a,a)') 2 'set qtl [ ', 3 '[]]' write(*,'(a)') 2 'set nqtl ', 3 'set map ...', 4 'set nall |perfect', 5 'set fre ...' write(*,'(a)') 2 'set nfam ', 3 'set mzfam ', 4 'set halfsibs ', 5 'set nsibs ' write(*,'(a)') '* SUMMARIZE MODEL *', 2 'list|ls|show map', '* SIMULATION *', 'set print|plevel 0|1|2', 3 'set high ', 4 'set low ', 5 'write loc []', 6 'simulate|write [nfam [nsibs [mzfam]]]' write(*,'(a)') '* MISCELLANEOUS *', 2 'set seeds ', 'set echo on|off ', 3 'set prompt on|off', 'time', '!|# ', 4 '%|$ ', 'include ', 5 'clear all settings', 'help', 'stop|quit' elseif (keyword(1:1).eq.' ') then continue elseif (keyword.eq.'sto' .or. keyword.eq.'qui' .or. & keyword.eq.'exi' .or.keyword(1:2).eq.'by') then goto 2 C else write(*,'(a,i3,a/7x,a/)') & 'ERROR: problematic input at line ',nlin,':',lin(1:72) end if goto 1 2 continue if (ilevel.gt.1) then write(*,'(/3a/)') & 'Closing include file "',infil(1:eow(infil)),'".' ilevel=ilevel-1 close(CSTRM,status='keep') goto 1 end if C C end of main loop C call stamp(t0) end C end-of-main C C Print variance components C subroutine showvc(str,xlinkd,ncat,prev,mu,va,vq,vc,ve,vt, & pqtl,mu1,mu2,mu3,mmu1,mmu2) integer CSTRM,OSTR parameter(CSTRM=4,OSTR=7) C integer MAXCAT, MISS parameter(MAXCAT=10, MISS=-9999) integer ncat, str logical xlinkd double precision mu, pqtl(2), prev(MAXCAT) double precision va, vc, vq, ve, vt, & mu1, mu2, mu3, mmu1, mmu2 character*1 com com=' ' if (str.eq.OSTR) com='#' vt=va+vq+vc+ve if (ncat.eq.MISS) then if (xlinkd) then write(str,'(/a,1x,a,f8.3,a,f8.3,a,f8.3,a)') 2 com,'Phenotypic mean = ',mu,' (Female SD=',sqrt(vt), 3 '; Male SD=',sqrt(vt+vq),')' else write(str,'(/a,1x,a,f8.3,a,f8.3,a)') & com,'Phenotypic mean = ',mu,' (SD=',sqrt(vt),')' end if elseif (ncat.eq.2) then write(str,'(/a,1x,a,f8.3))') com,'Trait prevalence= ',prev(2) else write(str,'(a,1x,a,10(1x,f8.3):)') & com,'Proportions = ',(prev(i), i=1, ncat) end if write(str,'(/a,1x,a,1x,f8.3/a,1x,a,3(1x,f8.3))') 2 com,'QTL allele prop = ', pqtl(1), 3 com,'QTL genot means = ', mu1, mu2, mu3 if (xlinkd) then write(str,'(a,1x,a,2(1x,f8.3))') & com,'QTL male means = ', mmu1, mmu2 write(str,'(2(/a,1x,a,f8.3,a,f6.1,a))') 2 com,'VQ female = ',vq,' (',100.0d0*vq/vt,'%)', 3 com,'VQ male = ',2*vq,' (',200.0d0*vq/(vq+vt),'%)' else write(str,'(/a,1x,a,f8.3,a,f6.1,a)') & com,'VQ = ',vq,' (',100.0d0*vq/vt,'%)' end if write(str,'(3(a,1x,a,f8.3,a,f6.1,a/))') 2 com,'VA = ',va,' (',100.0d0*va/vt,'%)', 3 com,'VC = ',vc,' (',100.0d0*vc/vt,'%)', 4 com,'VE = ',ve,' (',100.0d0*ve/vt,'%)' write(str,*) return end C end-of-showvc C C Show map C subroutine showmap(nloci, posqtl, pqtl, nall, het, map, theta) integer MAXLOC parameter(MAXLOC=60) double precision pqtl(2) integer nall, posqtl double precision het C Meiosis integer nloci double precision map(MAXLOC), theta(MAXLOC) integer i,j write(*,'(a//a/a)') 'User specified genetic map:', 2 'Locus Pos (cM) Theta Nallele Heterozyg ', 3 '------ -------- ----- ------- -----------' j=0 do 90 i=1,nloci if (i.ne.posqtl) then j=j+1 write(*,'(a6,i3,1x,f8.2,f8.3,3x,i5,5x,f6.4)') & 'Marker',j,map(i),theta(i),nall,het else write(*,'(a10,f8.2,f8.3,3x,i5,5x,f6.4)') 2 'QTL ',map(i),theta(i),2, 3 2.0d0*pqtl(1)*(1.0d0-pqtl(1)) end if 90 continue write(*,*) return end C end-of-showmap C C incorporate QTL position into marker map C subroutine place(mapqtl,posqtl,nmark,nloci,map,theta) integer MAXLOC parameter(MAXLOC=60) C QTL position on map double precision mapqtl C map integer nmark, nloci, posqtl double precision map(MAXLOC), theta(MAXLOC) double precision prev_pos C local variables integer i,j nloci=nmark+1 j=0 do 5 i=1, nloci if (i.ne.posqtl) then j=j+1 map(j)=map(i) end if 5 continue do 10 i=1, nmark if (map(i).gt.mapqtl) then do 20 j=nmark,i,-1 map(j+1)=map(j) 20 continue map(i)=mapqtl posqtl=i goto 30 end if 10 continue C else map(nmark+1)=mapqtl posqtl=nmark+1 30 continue prev_pos=0.0d0 do 50 i=1,nloci theta(i)=0.5d0*(1-exp(-0.02d0*(map(i)-prev_pos))) prev_pos=map(i) 50 continue return end C end-of-place C C Write out simulated pedigrees C subroutine simped(str,xlinkd,posqtl,pqtl,nqtl, mu1, mu2, mu3, 2 mmu1, mmu2, mu, va, vc, ve, 2 ncat, thresh, n_high, hi_thresh, n_low, lo_thresh, 3 nfam, nmzfam, nhsfam, nsibs, nloci, 4 pedigree, num,nfound,id,fa,mo,sex,twin, 5 theta, nall, cumfrq, qtl, trait, set, showqtl) C integer MAXALL, MAXCAT, MAXSIZ, MISS, MAXLOC, MAXQTL parameter(MAXALL=10, MAXCAT=10, MAXSIZ=20, MISS=-9999, & MAXLOC=60, MAXQTL=10) C C nfam Total number of nuclear families, C nmzfam Number of families containing MZ twins C nhsfam Number of families containing a half-sib to rest of family C nsibs Number of children per family C posqtl Location of 1st qtl on genetic map C showqtl 1=print 1st 2= print all qtl genotypes to pedigree file C 3=print ibd for 1st qtl C integer nfam, nhsfam, nmzfam, nqtl, nsibs, posqtl, showqtl, str logical xlinkd C C Variance components C integer ncat double precision mu, thresh(MAXCAT), va, vc, ve double precision mu1, mu2, mu3, mmu1, mmu2 C Ascertainment integer n_high, n_low double precision hi_thresh, lo_thresh C Pedigree structure integer pedigree integer num,nfound,id(MAXSIZ),fa(MAXSIZ),mo(MAXSIZ) integer sex(MAXSIZ), twin(MAXSIZ) C C Frequency of increasing allele at first qtl (others set at 0.5) C Genotype at all trait loci C double precision pqtl(2), pq(2) integer qtl(MAXSIZ, MAXQTL, 2) double precision trait(MAXSIZ) C C Markers allele frequencies and genotypes integer nloci, nall double precision cumfrq(MAXALL) integer set(MAXSIZ,MAXLOC,2) C C meiosis double precision theta(MAXLOC) C local variables integer c1,c2,f1,f2,first,gtp,hs,idxtwin,m1,m2,i,nhs,nmz,nadd integer hi_cou, lo_cou logical ordinal, xmale character*1 sx(2), aff(2) double precision aconst, aval, cval, eval, qval C functions integer irandom double precision rannor data sx,aff /'m','f','n','y'/ ordinal=(ncat.gt.0) C C Initialize marker positions and allele frequencies C pq(1)=0.5d0 pq(2)=1.0d0 C C Standardize QTL contributions C aconst=sqrt(2*va/dfloat(nqtl)) C C loop over pedigrees C idxtwin=0 pedigree=0 nhs=0 nmz=0 10 continue C C if MZ twins, they are first two siblings: C only simulate genotypes for second twin and subsequent children C hs=0 first=1 nfound=2 if (nmz.lt.nmzfam .and. nsibs.ge.2) then nmz=nmz+1 first=2 end if if (nhs.lt.nhsfam .and. (first.eq.1 .or. nsibs.gt.2)) then nhs=nhs+1 nfound=3 hs=nfound+1 end if num=nfound+nsibs do 6 i=1,nfound id(i)=i fa(i)=0 mo(i)=0 sex(i)=i if (i.gt.2) sex(i)=irandom(1,2) twin(i)=1 6 continue do 7 i=nfound+first,num id(i)=i fa(i)=1 mo(i)=2 if (i.eq.hs) then if (sex(3).eq.1) then fa(i)=3 else mo(i)=3 end if end if sex(i)=irandom(1,2) twin(i)=1 7 continue if (first.gt.1) then idxtwin=nfound+first twin(idxtwin)=2 do 8 i=nfound+1,idxtwin-1 id(i)=i fa(i)=fa(idxtwin) mo(i)=mo(idxtwin) sex(i)=sex(idxtwin) twin(i)=2 8 continue end if C first initialize the founder ibd-genotypes f1=1 f2=2 do 30 i=1,nfound xmale=(sex(i).eq.1 .and. xlinkd) if (.not.xmale) then do 31 j=1,nloci set(i,j,1)=f1 set(i,j,2)=f2 31 continue else do 32 j=1,nloci set(i,j,1)=f1 set(i,j,2)=f1 32 continue end if f1=f1+2 f2=f2+2 30 continue C C then gene drop the nonfounders ibd-genotypes C do 50 i=nfound+first,num if (sex(i).eq.1 .and. xlinkd) then call mumson(i,mo(i),set,nloci,theta) else call genof(i,fa(i),mo(i),set,nloci,theta) end if 50 continue C C now go through and simulate the founder marker genotypes C if (nall.gt.0) then do 60 i=1,nfound xmale=(sex(i).eq.1 .and. xlinkd) do 61 j=1,posqtl-1 call found(cumfrq,set(i,j,1)) if (.not.xmale) then call found(cumfrq,set(i,j,2)) else set(i,j,2)=set(i,j,1) end if 61 continue do 62 j=posqtl+1,nloci call found(cumfrq,set(i,j,1)) if (.not.xmale) then call found(cumfrq,set(i,j,2)) else set(i,j,2)=set(i,j,1) end if 62 continue 60 continue end if C and founder qtl genotypes do 65 i=1,nfound call found(pqtl,set(i,posqtl,1)) if (sex(i).eq.1 .and. xlinkd) then set(i,posqtl,2)=set(i,posqtl,1) else call found(pqtl,set(i,posqtl,2)) end if C the polygenes are autosomal do 63 j=1,nqtl call found(pq,qtl(i,j,1)) call found(pq,qtl(i,j,2)) 63 continue 65 continue C C and update the nonfounders based on the ibd-genotype they have inherited C do 70 i=nfound+first,num do 71 j=1,nloci c1=set(i,j,1) c2=set(i,j,2) f1=(c1-1)/2+1 f2=2-mod(c1,2) m1=(c2-1)/2+1 m2=2-mod(c2,2) set(i,j,1)=set(f1,j,f2) set(i,j,2)=set(m1,j,m2) 71 continue do 72 j=1,nqtl qtl(i,j,1)=qtl(1,j,irandom(1,2)) qtl(i,j,2)=qtl(2,j,irandom(1,2)) 72 continue 70 continue C C Copy genotypes from index twin to cotwin(s) C if (first.gt.1) then do 75 i=nfound+1,idxtwin-1 do 76 j=1,nloci if (j.eq.posqtl) then set(i,j,1)=set(idxtwin,j,1) set(i,j,2)=set(idxtwin,j,2) else set(i,j,1)=0 set(i,j,2)=0 end if 76 continue do 77 j=1,nqtl qtl(i,j,1)=qtl(idxtwin,j,1) qtl(i,j,2)=qtl(idxtwin,j,2) 77 continue 75 continue end if C C Write out simulated data for this family C hi_cou=0 lo_cou=0 cval=sqrt(vc)*rannor(0.0d0, 1.0d0) do 100 i=1,num eval=sqrt(ve)*rannor(0.0d0, 1.0d0) gtp=set(i,posqtl,1)+set(i,posqtl,2)-1 if (gtp.eq.1) then qval=mu1 else if (gtp.eq.2) then qval=mu2 else if (gtp.eq.3) then qval=mu3 end if nadd=0 do 110 j=1, nqtl nadd=nadd+qtl(i,j,1)+qtl(i,j,2)-3 110 continue aval=aconst*nadd trait(i)=mu+aval+qval+cval+eval if (ordinal) then do 115 k=1, ncat-1 if (trait(i).lt.thresh(k)) then trait(i)=dfloat(k) goto 116 end if 115 continue trait(i)=dfloat(ncat) 116 continue end if if (i.gt.nfound .and. trait(i).ge.hi_thresh) hi_cou=hi_cou+1 if (i.gt.nfound .and. trait(i).le.lo_thresh) lo_cou=lo_cou+1 100 continue C C Rejection sampling C if (hi_cou.lt.n_high .or. lo_cou.lt.n_low) goto 10 pedigree=pedigree+1 do 150 i=1,num if (ncat.eq.2) then write(str,'(i6,3(1x,i6),1x,a1,1x,i1,1x,a1,$)') 2 pedigree,id(i),fa(i),mo(i),sx(sex(i)), 3 twin(i), aff(int(trait(i))) else write(str,'(i6,3(1x,i6),1x,a1,1x,i1,1x,f12.6,$)') 2 pedigree,id(i),fa(i),mo(i),sx(sex(i)), 3 twin(i), trait(i) end if do 120 j=1, nloci if (j .ne. posqtl) then write(str,'(1x,i2,1x,i2,$)') set(i,j,1),set(i,j,2) elseif (showqtl.ne.0) then write(str,'(1x,i2,1x,i2,$)') set(i,j,1),set(i,j,2) end if 120 continue if (showqtl.gt.1) then do 130 j=1, nqtl write(str,'(1x,i2,1x,i2,$)') qtl(i,j,1),qtl(i,j,2) 130 continue end if write(str,*) 150 continue if (pedigree.lt.nfam) goto 10 200 continue return end C end-of-simped C C load cumulative frequencies for a marker with nall equifrequent alleles C subroutine mkfreq(nall, cumfrq, het) integer nall double precision cumfrq(*), het integer j double precision frq, tot if (nall.eq.0) then het=1.0d0 return end if frq=1.0d0/dfloat(nall) tot=0.0d0 do 5 j=1,nall-1 tot=tot+frq cumfrq(j)=tot 5 continue cumfrq(nall)=1.0d0 het=1.0d0-frq return end C end-of-mkfreq C C founder frequency C subroutine found(cumfrq,allele) integer allele double precision cumfrq(*) double precision x C functions double precision random x=random() allele=0 10 continue allele=allele+1 if (x.gt.cumfrq(allele)) goto 10 return end C end-of-found C C transmit haplotypes from each parent to child, including possibility C of recombination within haplotype C subroutine genof(indx,fa,mo,set,nloci,theta) integer MAXSIZ, MAXLOC parameter(MAXSIZ=20, MAXLOC=60) integer indx,fa,mo,nloci,set(MAXSIZ,MAXLOC,2) double precision theta(nloci) C local variables integer fagranp,mogranp,j C functions integer irandom double precision random C fagranp=irandom(1,2) mogranp=irandom(1,2) do 50 j=1,nloci if (theta(j).gt.random()) fagranp=3-fagranp if (theta(j).gt.random()) mogranp=3-mogranp set(indx,j,1)=set(fa,j,fagranp) set(indx,j,2)=set(mo,j,mogranp) 50 continue return end C end-of-genof C C transmit X chromosome haplotypes from mother to son, including possibility C of recombination within haplotype C subroutine mumson(indx,mo,set,nloci,theta) integer MAXSIZ, MAXLOC parameter(MAXSIZ=20, MAXLOC=60) integer indx,mo,nloci,set(MAXSIZ,MAXLOC,2) double precision theta(nloci) C local variables integer mogranp,j C functions integer irandom double precision random C mogranp=irandom(1,2) do 50 j=1,nloci if (theta(j).gt.random()) mogranp=3-mogranp set(indx,j,2)=set(mo,j,mogranp) set(indx,j,1)=set(indx,j,2) 50 continue return end C end-of-genof C C convert VQ to genotypic means C subroutine qtomu(pqtl, vq, mu, mu1, mu2, mu3) double precision pqtl(2) double precision mu, mu1, mu2, mu3, vq double precision qconst qconst=sqrt(0.5d0*vq/pqtl(1)/(1.0d0-pqtl(1))) mu1=mu-qconst mu2=mu mu3=mu+qconst return end C end-of-qtomu C C convert genotypic means to VQ (and VQD) C subroutine mutoq(pqtl, mu1, mu2, mu3, mu, vq, vqd) double precision pqtl(2) double precision mu, mu1, mu2, mu3, vq, vqd mu=pqtl(1)**2*mu1+2*(1.0d0-pqtl(1))*pqtl(1)*mu2+ & (1.0d0-pqtl(1))**2*mu3 mu3=mu3-mu mu2=mu2-mu mu1=mu1-mu vq=2*pqtl(1)*pqtl(2)*(pqtl(1)*(mu1-mu2)+pqtl(2)*(mu2-mu3))**2 vqd=pqtl(1)*pqtl(1)*pqtl(2)*pqtl(2)*(mu1-2*mu2+mu3)**2 return end C end-of-mutoq C C Algorithm AS 111, Appl.Statist., vol.26, 118-121, 1977. C Produces normal deviate corresponding to lower tail area = p. C double precision function ppnd(p) double precision p, q, r double precision a0, a1, a2, a3, b1, b2, b3, b4, & c0, c1, c2, c3, d1, d2, split double precision half, one, zero data split/0.42d0/ data a0,a1,a2,a3/2.50662823884d0,-18.61500062529d0, 1 41.39119773534d0,-25.44106049637d0/, b1,b2,b3,b4/ 2 -8.47351093090d0,23.08336743743d0,-21.06224101826d0, 3 3.13082909833d0/, c0,c1,c2,c3/-2.78718931138d0,-2.29796479134d0, 4 4.85014127135d0,2.32121276858d0/, d1,d2/3.54388924762d0, 5 1.63706781897d0/ data zero/0.d0/, one/1.d0/, half/0.5d0/ C q = p-half if (abs(q).gt.split) go to 10 C C 0.08 < p < 0.92 C r = q*q ppnd = q*(((a3*r + a2)*r + a1)*r + a0)/((((b4*r + b3)*r + b2)*r 1 + b1)*r + one) return C C p < 0.08 or p > 0.92, set r = min(p,1-p) C 10 r = p if (q.gt.zero) r = one-p if (r.le.zero) go to 20 r = sqrt(-log(r)) ppnd = (((c3*r + c2)*r + c1)*r + c0)/((d2*r + d1)*r + one) if (q.lt.zero) ppnd = -ppnd return 20 continue ppnd = zero return end C end-of-ppnd (AS111) C C The function RANDN() returns a normally distributed pseudo-random C number with zero mean and unit variance. Calls are made to a C function subprogram RANDOM() which returns independent random C numbers uniform in the interval (0,1). C C The algorithm uses the ratio of uniforms method of A.J. Kinderman C and J.F. Monahan augmented with quadratic bounding curves. C double precision function rannor(mu,sd) double precision mu,sd C functions double precision a,b,r1,r2,s,t C functions double precision random data s,t,a,b / 0.449871d0, -0.386595d0, 0.19600d0, 0.25472d0/ data r1,r2/ 0.27597d0, 0.27846d0/ C C Generate P = (u,v) uniform in rectangle enclosing acceptance region 50 u = random() v = random() v = 1.7156d0 * (v - 0.5d0) C Evaluate the quadratic form x = u - s y = abs(v) - t q = x**2 + y*(a*y - b*x) C Accept P if inside inner ellipse if (q .lt. r1) go to 100 C Reject P if outside outer ellipse if (q .gt. r2) go to 50 C Reject P if outside acceptance region if (v**2 .gt. -4.0d0*log(u)*u**2) go to 50 C Return ratio of P's coordinates as the normal deviate 100 randn = v/u rannor=randn*sd+mu return end C end-of-rannor C C random integer from U(lo..hi) C integer function irandom(lo,hi) integer lo,hi C functions double precision random irandom=lo+int(dfloat(hi-lo+1)*random()) return end C end-of-irandom C C Algorithm AS 183 Appl Stat 1982; 31:188 C Returns a pseudo-random number from U(0,1) C C ix,iy,iz should be "randomly" initialised to 1-3000 C eg via time C double precision function random() integer ix,iy,iz common /rand/ ix,iy,iz ix=171*mod(ix,177)-2*(ix/177) iy=172*mod(iy,176)-35*(iy/176) iz=170*mod(iz,178)-63*(iz/178) C if (ix.lt.0) ix=ix+30269 if (iy.lt.0) iy=iy+30307 if (iz.lt.0) iz=iz+30323 C random=mod(dfloat(ix)/30269.0d0+dfloat(iy)/30307.0d0 + & dfloat(iz)/30323.0d0,1.0d0) return end C end-of-random (AS183) C C pass line to shell -- requires existence of fairly C standard routine system() C subroutine shell(lin) character*(*) lin integer sta,fin external system sta=0 5 continue sta=sta+1 if (lin(sta:sta).ne.'$' .and. lin(sta:sta).ne.'%') goto 5 7 continue sta=sta+1 if (lin(sta:sta).eq.' ') goto 7 fin=len(lin) 10 continue fin=fin-1 if (lin(fin:fin).eq.' ') goto 10 write(*,'(a/3a/a)') '!','! "',lin(sta:min(fin,75)),'"','!' call system(lin(sta:fin)) return end C end-of-shell C C write elapsed time since last asked C subroutine stamp(t0) integer elapsed, t0 C functions C intrinsic time integer time elapsed=time()-t0 if (elapsed.lt.120) then write(*,'(/a,i5,a)') & 'This job took ',elapsed,' seconds' elseif (elapsed.lt.7200) then write(*,'(/a,f5.1,a)') & 'This job took ',float(elapsed)/60.0,' minutes' else write(*,'(/a,f5.1,a)') & 'This job took ',float(elapsed)/3600.0,' hours' end if return end C end-of-stamp C C extracts narg arguments from input string s C C typ=1 whitespace separated C typ=2 whitespace separated or reserved character (id by opchar()) C subroutine args(s,narg,arg,typ) character*(*) s integer narg, typ character*(*) arg(narg) integer eol,i,iarg,n,sarg,sol C functions integer eow, sow logical opchar do 5 i=1,narg arg(i)=' ' 5 continue sol=sow(s) eol=eow(s) n=1 i=sol C C start of main loop 30 continue if (i.gt.eol) goto 40 if (typ.eq.2 .and. opchar(s(i:i))) then arg(n)=s(i:i) n=n+1 elseif (s(i:i).eq.'"') then iarg=-1 i=i+1 sarg=i 45 continue if (i.gt.eol .or. s(i:i).eq.'"') goto 60 iarg=iarg+1 i=i+1 goto 45 elseif (s(i:i).ne.' ' .and. s(i:i).ne.'\t') then iarg=-1 sarg=i 50 continue if (i.gt.eol .or. s(i:i).eq.' '.or. s(i:i).eq.'\t' .or. & s(i:i).eq.'"' .or. (typ.eq.2 .and. opchar(s(i:i)))) then i=i-1 goto 60 end if iarg=iarg+1 i=i+1 goto 50 60 continue C arg(n)=s(sarg:(sarg+iarg)) n=n+1 end if i=i+1 goto 30 40 continue C C return # arguments actually found narg=n-1 return end C end-of-arg C C is a reserved character for primitives? "()*+-/<=>^" C logical function opchar(ch) character*1 ch integer ich ich=ichar(ch) opchar=((ich.ge.40 .and. ich.le.43) .or. ich.eq.45 .or. & ich.eq.47 .or. (ich.ge.60 .and. ich.le.62) .or. ich.eq.94) return end C end-of-opchar C C character to integer conversion via internal read C integer function ival(string) integer MISS parameter(MISS=-9999) character*20 string integer i if (string.eq.' ') then ival=0 elseif (string.eq.'x' .or. string.eq.'X' .or. string.eq.'.') then ival=MISS else read(string,'(i20)',err=10) i ival=i end if return C error -- word is not an integer 10 write(*,'(2a/)') 'ERROR: Unable to read integer ',string ival=0 return end C end-of-ival C C character to float conversion via internal read C double precision function fval(string) integer BLANK, MISS parameter(BLANK=-9999,MISS=-9999) character*20 string double precision v if (string.eq.' ' .or. string.eq.'-') then fval=BLANK elseif (string.eq.'x' .or. string.eq.'X' .or. string.eq.'.') then fval=MISS elseif (string.eq.'y' .or. string.eq.'Y') then fval=2.0d0 elseif (string.eq.'n' .or. string.eq.'N') then fval=1.0d0 else read(string,'(f20.0)',err=10) v fval=v end if return C error -- word is not a number 10 write(*,'(2a/)') 'ERROR: Unable to read number ',string fval=0.0d0 return end C end-of-fval C C skip leading whitespace C integer function sow(string) character*(*) string do 10 i=1,len(string) 10 if (string(i:i).ne.' ' .and. string(i:i).ne.'\t') goto 20 20 sow=i return end C C find end of string C integer function eow(string) character*(*) string do 10 i=len(string),1,-1 10 if (string(i:i).ne.' ') goto 20 20 eow=i return end C end-of-eow C C C end-of-program