if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("affy",update=F)
BiocManager::install("CLL",update=F)
BiocManager::install("affyQCReport",update=F)
BiocManager::install("limma",update=F)
BiocManager::install("ALL",update=F)
BiocManager::install("annotate",update=F)
BiocManager::install("annaffy",update=F)
BiocManager::install("KEGG.db",update=F)
BiocManager::install("hgu95av2.db",update=F)

library(affy)
library(CLL)
data("CLLbatch")
sampleNames(CLLbatch)
data(disease)
rownames(disease)=disease$SampleID
sampleNames(CLLbatch)=sub("\\.CEL$","",sampleNames(CLLbatch))
mt=match(rownames(disease),sampleNames(CLLbatch))
vmd=data.frame(labelDescription=c("Sample ID", "Disease status: progressive or stable"))
phenoData(CLLbatch)=new("AnnotatedDataFrame",data=disease[mt,], varMetadata=vmd)
CLLbatch=CLLbatch[,!is.na(CLLbatch$Disease)]

# QA/QC
library("affyQCReport")
# when installing using 3.3.1 on a windows machine there were problems with writing the package to the appropriate location:
# probably windows protecting writing files to certain locations-find out download locations from R console, copy and paste that
# into top of window in windows explorer and directly move the uncompressed directory to the location in documents that holds R packages
affyQAReport(CLLbatch)

badArray=match("CLL1", sampleNames(CLLbatch))
CLLB=CLLbatch[,-badArray]

#RMA
CLLrma=rma(CLLB)

# filtering
CLLf=nsFilter(CLLrma, remove.dupEntrez=F,var.cutof=0.5)$eset

# testing for DE
CLLtt=rowttests(CLLf, "Disease")

library(limma)
design=model.matrix(~CLLf$Disease)
CLLlim=lmFit(CLLf, design)
CLLeb=eBayes(CLLlim)

# smaller example-resetup ALL example
library(ALL)
data(ALL)
bcell=grep("^B",as.character(ALL$BT))
types=c("NEG", "BCR/ABL")
moltyp=which(as.character(ALL$mol.biol) %in% types)
ALL_bcrneg=ALL[,intersect(bcell,moltyp)]
ALL_bcrneg$mol.biol=factor(ALL_bcrneg$mol.biol)
ALL_bcrneg$BT=factor(ALL_bcrneg$BT)

subs=c(35,65,75,1,69,71)
ALLset2=ALL_bcrneg[,subs]

tt2=rowttests(ALLset2, "mol.biol")
cl=as.numeric(ALL_bcrneg$mol.biol=="BCR/ABL")
design=cbind(mean=1,diff=cl)
fit2=eBayes(lmFit(exprs(ALLset2),design=design[subs,]))

pdf("limmaCompare.pdf")
plot(-log10(tt2$p.value), -log10(fit2$p.value[,"diff"]),
  xlab="t-test",ylab="limma",pch=".")
abline(0,1,col="red")
dev.off()

pdf("volcanoPlt.pdf")
plot(CLLtt$dm, -log10(CLLtt$p.value), pch=".",xlab="log-ratio",ylab=expression(-log[10]~p))
dev.off()

# annotation
library(annotate)
annotation(CLLf)
library(hgu95av2.db)
genenames=featureNames(CLLf)[p.adjust(CLLeb$p.value[,2],method="BH")<.05]

ll=getEG(genenames,"hgu95av2")
sym=getSYMBOL(genenames,"hgu95av2")

tab=topTable(CLLeb, coef=2, adjust.method="BH", n=9)
tab=data.frame(sym,signif(tab[,-1],3))
htmlpage(list(ll),othernames=tab,, filename="GeneList1.html",
  title="HTML Report", table.center=T,
  table.head=c("Entrez ID",colnames(tab)))

library("annaffy")
library("KEGG.db")

atab=aafTableAnn(genenames, "hgu95av2.db", aaf.handler())
saveHTML(atab, file="GeneList2.html")