###1###
#
#Read a subset of the barley data and examine first few rows.
#
#######

d=read.csv("barley.csv",header=T)
d[1:4,]

boxplot(d[,-1])

logd=log(d[,-1])
boxplot(logd)

###2###
#
#Median center the data.
#
#######

medians=apply(logd,2,median)
norm.d=sweep(logd,2,medians)

boxplot(norm.d)

###5###
#
#Export the median-centered data.
#
#######

write.csv(norm.d,"nd.csv",row.names=F)

###6###
#
#Now use SAS to change the format of the data and analyze it (see barley.sas).
#Read the results file created by SAS.
#
#######

results=read.delim("results.txt",header=T)
results[1:6,]


###7###
#
#Separate the genotype, time, and genotype-by-time p-values and analyze.
#
#######

p=matrix(results[,6],byrow=T,ncol=3)
p.geno=p[,1]
p.time=p[,2]
p.gtime=p[,3]

#Examine the histograms of p-values.

hist(p.geno)
hist(p.time)
hist(p.gtime)

#Estimate the number of genes whose expression is
#constant accross genotypes.

source("http://www.public.iastate.edu/~dnett/microarray/rfunctions.txt")

estimate.m0(p.geno)

#Estimate the number of genes whose expression is
#constant accross the time points profiled.

estimate.m0(p.time)

#Estimate the number of genes with no genotype-by-time interaction

estimate.m0(p.gtime)

#Compute q-values for genotype, time, and genotype-by-time interaction

q.geno=q.value(p.geno)
q.time=q.value(p.time)
q.gtime=q.value(p.gtime)

#Find number of genes with q-values <= 0.05.

sum(q.geno<=0.05)

sum(q.time<=0.05)

sum(q.gtime<=0.05)

#Obtain IDs for the most significant genes for genotype-by-time interaction.

(d[order(p.gtime),1])[1:10]

#Plot lsmeans for most significant genes.

lsmeans=read.delim("lsmeans.txt",header=T)
source("http://www.public.iastate.edu/~dnett/microarray/plotfunction.txt")

plot.lsmeans(lsmeans,456)
plot.lsmeans(lsmeans,1856)
plot.lsmeans(lsmeans,1111)
plot.lsmeans(lsmeans,332)

########################################################################

SAS CODE: barley.sas


options nodate nonumber;
/*import the data file*/
proc import datafile="d:/barley/nd.csv" out=d replace dbms=csv;
run;
quit;

/*look at the first 36 rows of the data*/
proc print data=d (obs=4);
run;
quit;

/* Prepare data for mixed model analysis */
data ndat;
set d;
run;

data ndat;
set ndat;
gene= _n_;
rep=1;
geno="ML6";
time=0;
intensity=V1;
output;
gene= _n_;
rep=1;
geno="ML6";
time=8;
intensity=V2;
output;
gene= _n_;
rep=1;
geno="ML6";
time=16;
intensity=V3;
output;
gene= _n_;
rep=1;
geno="ML6";
time=20;
intensity=V4;
output;
gene= _n_;
rep=1;
geno="ML6";
time=24;
intensity=V5;
output;
gene= _n_;
rep=1;
geno="ML6";
time=32;
intensity=V6;
output;
gene= _n_;
rep=1;
geno="ML13";
time=0;
intensity=V7;
output;
gene= _n_;
rep=1;
geno="ML13";
time=8;
intensity=V8;
output;
gene= _n_;
rep=1;
geno="ML13";
time=16;
intensity=V9;
output;
gene= _n_;
rep=1;
geno="ML13";
time=20;
intensity=V10;
output;
gene= _n_;
rep=1;
geno="ML13";
time=24;
intensity=V11;
output;
gene= _n_;
rep=1;
geno="ML13";
time=32;
intensity=V12;
output;

/*rep2*/
gene= _n_;
rep=2;
geno="ML6";
time=0;
intensity=V13;
output;
gene= _n_;
rep=2;
geno="ML6";
time=8;
intensity=V14;
output;
gene= _n_;
rep=2;
geno="ML6";
time=16;
intensity=V15;
output;
gene= _n_;
rep=2;
geno="ML6";
time=20;
intensity=V16;
output;
gene= _n_;
rep=2;
geno="ML6";
time=24;
intensity=V17;
output;
gene= _n_;
rep=2;
geno="ML6";
time=32;
intensity=V18;
output;
gene= _n_;
rep=2;
geno="ML13";
time=0;
intensity=V19;
output;
gene= _n_;
rep=2;
geno="ML13";
time=8;
intensity=V20;
output;
gene= _n_;
rep=2;
geno="ML13";
time=16;
intensity=V21;
output;
gene= _n_;
rep=2;
geno="ML13";
time=20;
intensity=V22;
output;
gene= _n_;
rep=2;
geno="ML13";
time=24;
intensity=V23;
output;
gene= _n_;
rep=2;
geno="ML13";
time=32;
intensity=V24;
output;

/*rep3*/
gene= _n_;
rep=3;
geno="ML6";
time=0;
intensity=V25;
output;
gene= _n_;
rep=3;
geno="ML6";
time=8;
intensity=V26;
output;
gene= _n_;
rep=3;
geno="ML6";
time=16;
intensity=V27;
output;
gene= _n_;
rep=3;
geno="ML6";
time=20;
intensity=V28;
output;
gene= _n_;
rep=3;
geno="ML6";
time=24;
intensity=V29;
output;
gene= _n_;
rep=3;
geno="ML6";
time=32;
intensity=V30;
output;
gene= _n_;
rep=3;
geno="ML13";
time=0;
intensity=V31;
output;
gene= _n_;
rep=3;
geno="ML13";
time=8;
intensity=V32;
output;
gene= _n_;
rep=3;
geno="ML13";
time=16;
intensity=V33;
output;
gene= _n_;
rep=3;
geno="ML13";
time=20;
intensity=V34;
output;
gene= _n_;
rep=3;
geno="ML13";
time=24;
intensity=V35;
output;
gene= _n_;
rep=3;
geno="ML13";
time=32;
intensity=V36;
output;
drop V1-V36;
run;

proc print data=ndat (obs=36);
run;

/*the data for the first gene*/
data temp;  set ndat;
  if gene=1;
run;

ods trace on; /* turns on the writing of the trace record */
run;
/*fitting mixed linear model for the first gene*/
proc mixed data=temp covtest IC;
 class gene rep geno time intensity; /*classification variables*/
 model intensity=geno time geno*time/DDFM=KENWARDROGER;
 random rep rep*geno rep*time;
 lsmeans geno*time;
run;
quit;


ods trace off; /* turns off the writing of the trace record */
run;

ods noresults;
run;

ods listing close; /*do not print in the output window*/
run;

options nomlogic nonotes; /*Using SAS System Options to Suppress Log Output*/

proc mixed data=ndat covtest IC;
class gene rep geno time intensity;
by gene;
ods output InfoCrit=IC;
ods output Tests3=results;
ods output CovParms=varcomps;
ods output lsmeans=means;
model intensity=geno time geno*time/DDFM=KENWARDROGER outp=resids;
random rep rep*geno rep*time;
lsmeans geno*time;
run;
quit;
ods listing;  /*to print in the output window*/
run;

proc print data=IC (obs=10);
run;

proc print data=varcomps (obs=10);
run;

proc print data=results (obs=10);
run;

proc print data=means (obs=10);
run;


/*export estimates of the variance components*/
proc export data=varcomps outfile="D:\barley\vars.txt" dbms=TAB replace;
run;

/*export information criteria*/
proc export data=IC outfile="D:\barley\IC.txt" dbms=TAB replace;
run;

/*export p-values for the test of fixed effects, genotype, time, and genotype-time interaction */
data results; set results;
 pval=1*probf;
 drop probf;
run;
proc export data=results outfile="D:\barley\results.txt" dbms=TAB replace;
run;

/*export lsmeans*/
proc export data=means outfile="D:\barley\lsmeans.txt" dbms=TAB replace;
run;


/*export residuals*/
proc export data=resids outfile="D:\barley\resids.txt" dbms=TAB replace;
run;