Code test - MA analysis

illumina_test.R

##

library(limma)
library(illuminaHumanv4.db)

args=commandArgs(TRUE)
sam=read.ilmn(args[1],ctrlfiles=args[2],other.columns="Detection")

# Normalization - neqc

sam.norm=neqc(sam)

#Remove low qual probes

ids=as.character(rownames(sam.norm))
qual=unlist(mget(ids,illuminaHumanv4PROBEQUALITY, ifnotfound=NA))
rem=qual=="No match" | qual=="Bad"
sam.filt=sam.norm[!rem,]

#DE analysis

design=model.matrix(~0+factor(c(1,1,1,1,1,1,2,2,2,2,2,2)))
colnames(design)=c("HighE","LowE")
aw=arrayWeights(sam.filt,design)
fit=lmFit(sam.filt,design,weights=aw)
contrasts=makeContrasts(HighE-LowE, levels=design)
contr.fit=eBayes(contrasts.fit(fit,contrasts))

pdf("VolcanoPLot.pdf",width=7,height=5)
volcanoplot(contr.fit,main="HighE - LowE")
dev.off()

write.table(topTable(contr.fit,coef=1,number=Inf,file="DEresult.txt"))