/*---read in the data; assumes files are stored as text files
     with names sudarsanam1.txt - sudarsanam12.txt; change the
     datafile= pathname below to match your directory---*/
%macro readdata;

proc datasets;
   delete s2;
run;

%do a = 1 %to 12;

   proc import out=s
      datafile="C:\My Documents\microarray\Sudarsanam\sudarsanam&a..txt"
      dbms=tab replace;
      getnames=yes;
   run;

   data s;
      length gene $ 13 name $ 13;
      set s;
      name = upcase(name);
      type = upcase(type);
      gene = upcase(gene);
      array = &a;
      if (name=" ") then name = type;
      if (gene=" ") then gene = name;
      spot = _n_;
      if (flag=0) then do;
         if (array <= 3) then strain = "snf2rich";
         else if (array <= 6) then strain = "snf2mini";
         else if (array <= 9) then strain = "swi1rich";
         else if (array <= 12) then strain = "swi1mini";
         diff = ch1i-ch1b;
         if (diff > 0) then logi = log2(diff);
         else logi = .;
         output;
         strain = "wildtype";
         diff = ch2i-ch2b;
         if (diff > 0) then logi = log2(diff);
         else logi = .;
         output;
      end;
      keep array gene name spot strain logi;
   run;
 
   proc append base=s2 data=s;
   run;

%end;

%mend;

%readdata
run;

libname d 'D:';
run;

/*---overall normalization---*/
proc mixed data=s2 covtest cl lognote;
   class array strain;
   model logi = strain / 
      outp=d.sudarp(keep=array
      gene name spot strain resid);
   random array array*strain;
   lsmeans strain / diff cl;
run;
/*---removes some genes that slow the analysis---*/

data d.sudarp;
   set d.sudarp;
   where gene not in ('EMPTY', 'NORF');
run;

proc sort data=d.sudarp;
   by gene array spot;
run;

/*---turn off printing---*/
ods listing close;
run;

/*---individual incomplete block models for each clone using
     residuals---*/
proc mixed data=d.sudarp;
   by gene;
   class array spot strain;
   model resid = strain / outp=d.sudarr;
   random spot(array);
   lsmeans strain / diff;
   ods output covparms=d.sudarc tests3=d.sudart 
      diffs=d.sudard;
run;

ods listing;
run;

/*---Bonferroni adjustment---*/
data d.sudard1;
   set d.sudard;
   if (probt = .) then delete;
   logp = log10(probt);
   neglogp = -logp - log10(6917*10);
   if (neglogp <= 0) then delete;
run;