*! Date    : 1 Jun 2005
*! Version : 1.04
*! Author  : Adrian Mander
*! Email   : junk.ade@ntlworld.com

/* DATA is all in allsnp
qui insheet using original.csv,clear comma
qui _genp
qui compress
qui save allsnp,replace
*/

prog def hapblock
version 8.0
syntax [varlist] [using/] [, MV MVDEL HLEN(numlist) Start(integer 1) SModel(integer 1) Block(string) REPLACE]
preserve

global hapmod1 ""
global hapdf1 ""
global hapllhd1 ""
global hapmod0 ""
global hapdf0 ""
global hapllhd0 ""

if "`mvdel'"=="" & "`mv'"=="" local mvdel "mvdel"
if "`mvdel'"~="" & "`mv'"~="" local mvdel ""
if "`hlen'"=="" local hlen "2"
if "`using'"=="" local using "results"
if "`block'"=="" local block "block"

di as text in smcl "{hline}"
di in text "Missing data will be specified in hapipf using the " in res "`mv'`mvdel'" in text " option"
di in text "The window length will be " in res "`hlen'"
di in text "All the results will be stored to file " in res "`using'"
di in text "All the test statistics will be stored to file " in res "`block'"
di in text in smcl "{hline}"

qui drop _all
qui set obs 0
qui gen model=""
qui gen df=.
qui gen double llhd=.
qui gen str80 vlist=""
qui gen str5 miss=""
qui gen sloc = .
qui gen plusi = .
qui gen haplen = .
qui g pctmiss =.
cap qui save "`using'",`replace'
if _rc~=0 {
  if _rc==602 di as error "`using'.dta already exists use -replace- option"
  exit(602)
}
restore

local wc:word count `varlist'
local mxnloc = `wc'/2


foreach haplen of numlist `hlen' { /*loop for the length of haplotype of interest */
  local step 1
  local lastloc = `mxnloc'-`haplen'+1
  forvalues loci=`start'(`step')`lastloc' {
    di as res _continue "`loci',"

    /*Create the varlist*/
    local vlist ""
    local sloc = 2*`loci'-1
    local lloc = 2*(`loci'+`haplen'-1)
    forv var = `sloc'/`lloc' {
       local temp: word `var' of `varlist'
       local vlist "`vlist' `temp'"
    }

    /* Loop through the models more at start and end*/

    local midloc=int(`haplen'/2+1)
    local endloc=`midloc'
    local startloc 1
    if `loci'==1 local modstep 1
    if `loci'~=1 local modstep=`midloc'-1
    if `loci'==`lastloc' {
       local modstep 1
       local endloc = `haplen'+1
       local startloc=`midloc'
    }
    /* Start at a particular model */
    if `smodel'~=1 {
      local startloc `smodel'
      local smodel 1
    }

    foreach i of numlist `startloc'(`modstep')`endloc' {
      local ipf "l1"
      local midloc = `haplen'/2+1
      forvalues loc = 2/`haplen' {
        if `i'==1 local ipf "`ipf'*l`loc'"
        if `i'~=1 & `loc'==`i'  local ipf "`ipf'+l`loc'"
        if `i'~=1 & `loc'~=`i'  local ipf "`ipf'*l`loc'"
      }

      preserve
      qui hapipf `vlist', ipf(`ipf') nolog quiet `mv' `mvdel'  model(1) acc(0.005) ipfacc(0.001)

      qui use "`using'",replace
      local nobs=_N+1
      qui set obs `nobs'
      qui replace sloc = `loci' in `nobs'
      qui replace plusi = `i' in `nobs'
      qui replace haplen = `haplen' in `nobs'
      qui replace model= "$hapmod1" in `nobs'
      qui replace df= $hapdf1 in `nobs'
      qui replace llhd = $hapllhd1 in `nobs'
      qui replace vlist = "`vlist'" in `nobs'
      qui replace miss = "`mv'`mvdel'" in `nobs'
      qui replace pctmiss = `r(nmiss)'/`r(N)' in `nobs'
      qui save "`using'",replace
      restore

    }

  }
}


prores using results, b("`block'") `replace'

end

/* Process the output to find out where the blocks are */

prog def prores
syntax using/ [,Log Block(string) REPLACE]
preserve
qui use "`using'",replace

if _N==0 {
  di as error "No models have been fitted..!"
  di as error "Perhaps adjust the smodel() or start() choice?"
  exit(198)
}

/* Generate the test statistics */
qui compress
sort miss haplen sloc df
qui by miss haplen sloc: gen double chi = 2*(llhd[1]-llhd)
qui by miss haplen sloc: gen ddf = df-df[1]
qui by miss haplen sloc: gen double pv= chi2(chi,ddf)
qui gen position = sloc+plusi-1.5
qui gen mlen = cond(miss=="mv", "Window width "+string(haplen)+", MAR", "Window width "+string(haplen)+", MCAR")
qui gen double lpv = log(1-pv)
qui replace lpv = -40 if lpv==.

/* Calculate the blocks */

qui drop if df==0
qui gen sig=pv<0.05
sort miss haplen sloc position

by miss haplen: gen blocks=cond(sig==1,"1", "") if _n==1

qui by miss haplen: gen blocke= string(position+0.5) if _n==_N
qui by miss haplen: replace blocks=cond(sig==1,                        /*
*/     cond( blocks[_n-1]=="", string(position-0.5), blocks[_n-1])  ,  /*
*/     ""   /*
*/ )
qui by miss haplen: replace blocke=cond(sig[_n+1]==0,     /*
*/     string(position+0.5)   ,                           /*
*/     ""   /*
*/ ) if _n~=_N
qui gen block = blocks+"-"+blocke if blocks~="" & blocke~=""
qui gen blocksize = real(blocke)-real(blocks)+1
qui compress

/*
tab blocksize mlen
tab block mlen
*/

lab var pv "p-value"
lab var lpv "log(1-pvalue)"
lab var position "Model Position of the + Sign"

cap qui save "`block'",`replace'
if _rc~=0 di as error "NOTE: `block'.dta already exists use the -replace- option"

if "`log'"=="" {
  local y "pv"
  local ylab "ylab(0 0.05 0.5 1.0)"
}
else {
  local y "lpv"
  local ylab "yline(-.05129329)"
}

scatter `y' position `if', by(mlen) `ylab' yline(0.05)

restore
end