# HW 4 ***********************  Cheok et al data **********************************************

setwd("C:/cwru/classes/epbi473_05/lec6(biocAffyLimma)")
library(limma)
library(Biobase)
library(cheokEset)
cheok
p=pData(cheok);p
e=exprs(cheok)
meds=apply(e,2,median,na.rm=T);meds
zmeds=rep(1000,12625)%o%(1/meds)
e=e*zmeds
meds=apply(e,2,median,na.rm=T);meds



MP=factor(p$TX %in% c("MP","LDMTX.MP","HDMTX.MP"),labels=c("absent","present"))
MP

LD=(p$TX =="LDMTX.MP")
LD

HD=(p$TX %in% c("HDMTX","HDMTX.MP"))

MTX=rep(0,120)
MTX[LD==T]=1
MTX[HD==T]=2
mtx=factor(MTX,labels=c("absent","low","high"))
mtx


np=data.frame(p,MP=MP,MTX=mtx)
np
ii=order(p$timePoint,MP,mtx);ii
pdat=np[ii,];pdat

ne=e[,rownames(pdat)]
vars=c(cheok@phenoData@varLabels,MP="6-mercaptopurine",MTX="methotrexate")
vars
pD <- new("phenoData", pData = pdat, varLabels = vars)
neset <- new("exprSet", exprs = as.matrix(ne), phenoData = pD)
neset
dnp=pData(neset);dnp
timeM=as.numeric(dnp$timePoint=="post")
hi=as.numeric(dnp$MTX=="high")
low=as.numeric(dnp$MTX=="low")
mp=as.numeric(dnp$MP=="present")
intHiMp=hi*mp

nMatrix=cbind(low,hi,mp,intHiMp)
nMatrix
nMatrix[1:60,]=0
nMatrix

M060=matrix(0, ncol=60, nrow=60)
diag(M060)=1
	
M120=rbind(M060,M060)
	
nM120=cbind(nMatrix,M120)#part1

fit <- lmFit(neset, nM120) 
ebfit <- eBayes(fit )

tops=NULL #part1
 for (i in 1:4) tops=rbind(tops,topTable(ebfit,coef=i,n=100))

affyID=tops[,"ID"]
Pvals=tops[,"P.Value"]
tstats=tops[,"t"]
length(affyID)

library(hgu95av2)
llnames <- mget(affyID, env = hgu95av2GENENAME)
gnames <- as.character(llnames)
llnames <- mget(affyID, env = hgu95av2SYMBOL)
sym <- as.character(llnames)
write.table(cbind(affyID,sym,gnames,tstats,Pvals),file="hw4.csv",sep=",",row.names=FALSE)


# double check B versus T cells. 
BT=(pD$ALL.Type=="T.Cell")
design=model.matrix(~BT)
design
TBfit <- lmFit(neset, design) 
TBebfit <- eBayes(TBfit ) 
TBebfit
tops=NULL
tops=topTable(TBebfit,coef=2,n=100)
affyID=tops[,"ID"]
Pvals=tops[,"P.Value"]
tstats=tops[,"t"]
length(affyID)
library(hgu95av2) 
llnames <- mget(affyID, env = hgu95av2GENENAME)
gnames <- as.character(llnames)
llnames <- mget(affyID, env = hgu95av2SYMBOL)
sym <- as.character(llnames)
write.table(cbind(affyID,sym,gnames,tstats,Pvals),file="cellType.csv",sep=",",row.names=FALSE)
data.frame(affyID,sym,gnames,tstats,Pvals)