/*
Copyright (C) Jing Hua Zhao, MRC-Epid 4-2-2009

The software is freely available for academic use provided appropriate
acknowledgement is made.
*/

%macro wtl(data, map, trait, snpid=rsn, vlist=, inc=);
data out.long (keep=&snpid id &vlist a1a2 add n);
     set &data;
     fid=open("&data");
     length id $11. add 3. a1a2 $3.;
     format add 1.;
     set &map point=_n_;
     n=0;
     do col=2 to attrn(fid,"nvars");
        iid=col-1;
        set &trait (keep=&vlist &inc) point=iid;
        if &inc=1 then do;
           id=varname(fid,col);
           a1a2=vvaluex(id);
           add=.;
           if a1a2 ne "   " then do;
              a1=substr(a1a2,1,1);
              a2=substr(a1a2,3,1);
              add=(a1=b)+(a2=b);
              n+1;
           end;
           output;
        end;
     end;
     rc=close(fid);
run;
%mend wtl;

libname in ("." "/genetics/data/GWA/EPIC/6-5-7/wide2");
libname out '/tmp/scratch';
options compress=yes mprint ps=max;
proc sql;
     create table ind as select id, age from in.id where cohort="1" order by id;
     create table bmi as select id, bmi from in.trait order by id;
     create table map as select chr, rsn, a, b, pos from in.map order by chr, rsn;
quit;
data trait;
     merge ind (in=id_flag) bmi;
     format bmi 5.2 obesity 1.;
     if bmi ne . then obesity=(bmi>=30);
     by id;
     one=id_flag;
run;
%wtl(in.a1a2, map, trait, snpid=rsn, vlist=age bmi obesity, inc=one);
ods select none;
data out.cr (keep=rsn n);
     set out.long;
     by rsn notsorted;
     if last.rsn;
run;
proc allele data=out.long genocol;
     by rsn notsorted;
     var a1a2;
     ods output markersumm=out.ms allelefreq=out.af;
run;
proc reg data=out.long;
     by rsn notsorted;
     ods output parameterestimates=out.bmipm;
     model bmi = age add / b stb;
quit;
proc logistic data=out.long descending;
     by rsn notsorted;
     ods output parameterestimates=out.obpm CLOddsPL=out.obclpm;
     model obesity = age add / expb clodds=pl;
run;
proc casecontrol data=out.long tall marker=rsn indiv=id genocol delimiter="/" outstat=out.cc;
     by rsn notsorted;
     trait obesity;
     var a1a2;
run;
ods select all;

%macro w2l(data, map, trait, snpid=rsn, vlist=, inc=);
%let fid=%sysfunc(open(&data));
%let nvars=%sysfunc(attrn(&fid,nvars));
%let i1=%sysfunc(varname(&fid,2));
%let in=%sysfunc(varname(&fid,&nvars));
%let rc=%sysfunc(close(&fid));
data out.long (keep=&snpid id &vlist a1a2 add n);
     set &data;
     array ids $ &i1.--&in;
     length add 3.;
     format add 1.;
     set &map point=_n_;
     n=0;
     do col=2 to &nvars;
        iid=col-1;
        set &trait (keep=&vlist &inc) point=iid;
        if &inc=1 then do;
           id=vname(ids(iid));
           a1a2=ids(iid);
           add=.;
           if a1a2 ne "   " then do;
              a1=substr(a1a2,1,1);
              a2=substr(a1a2,3,1);
              add=(a1=b)+(a2=b);
              n+1;
           end;
           output;
        end;
     end;
run;
%mend w2l;
