Dear all,

I am just a beginner of R program. I need to analyze my microarray data. So, I have just started (just 2 days only) to practice with one data set available in GEO called "GSE58231"

I have compiled the script following the wiki page "Analyze your own microarray data in R/Bioconductor". However, a lot of functions don't work properly, and unfortunately, as I am a beginner, I am unable to figure out the problem. So, I am requesting experts' opinion on this issue. Here is the script. I have made bold where I am facing the problem

# Install bioconductor
source("http://bioconductor.org/biocLite.R")

#Install requried packages of bioconductor for Microarray
source("https://bioconductor.org/biocLite.R")
biocLite("affy")
source("https://bioconductor.org/biocLite.R")
biocLite("affydata")
source("https://bioconductor.org/biocLite.R")
biocLite("ArrayExpress")
source("https://bioconductor.org/biocLite.R")
biocLite("limma")
source("https://bioconductor.org/biocLite.R")
biocLite("Biobase")
source("https://bioconductor.org/biocLite.R")
biocLite("Biostrings")
source("https://bioconductor.org/biocLite.R")
biocLite("genefilter")

# Load required libraries
library(affy)
library(affydata)
library(ArrayExpress)
library(limma)
library(Biobase)
library(Biostrings)
library(genefilter)

#Open and import raw CEL File data in Affay
celpath = "C:/Users/Ornia/Documents/R/Microarray Practice"

data = ReadAffy (celfile.path=celpath)

# Retrieve sample annotation
ph = data@phenoData
ph
ph@data

## If there is no sample annatation in previous section then rewrite the sample annotaion.
ph@data[ ,1] = c("wt1","wt2","wt3","erf96oe1","erf96oe2","erf96oe3")
ph

#Plot histogram of each sample
for (i in 1:6)
{
  name = paste("histogram",i,".tiff",sep="")
  tiff(name)
  hist(data[,i],lwd=2,which='pm',ylab='Density',xlab='Log2 intensities',main=ph@data$sample[i])
  dev.off()
}

#Plot histogram of all sample combined (this one is not working)
**op = par(mfrow = c(2,3))
for(i in 1:6){hist(data[,i],lwd=2,which='pm',ylab='Density',xlab='Log2 intensities',main=ph@data$sample[i])}**

#Plot histogram of all sample in one graph (this one is not working)
**color=c('green','green','green','red','red','red')
hist(data[,1:6],lwd=2,which='pm',col=color,ylab='Density',xlab='Log2 intensities',main='Histogram of raw data')**

#Quantificatioin of Quality (None of them are working)
**data.qc = qc(data)
avbg(data.qc)
sfs(data.qc)
percent.present(data.qc)
ratios(data.qc)
Showed there is no qc**

# Data normalization
data.rma = rma(data)
data.matrix = exprs(data.rma)

#MA plot for normalized data (this one is not working)
**for (i in 1:6)
{
  name = paste("MAplotnorm",i,".tiff",sep="")
  tiff(name)
  MAplot(data.rma,which=i)
  dev.off()
}**

#Principal Component Analysis (this one is not working)
**color=c('green','green','green','red','red','red')
data.PC = prcomp(t(data.matrix),scale.=TRUE)
plot(data.PC$x[1:2],col=color)**


# Identification of DEGs
ph@data[ ,2] = c("wt","wt","wt","ox","ox","ox")
colnames(ph@data)[2]="source"
ph@data
groups = ph@data$source
f = factor(groups,levels=c("wt","ox"))
design = model.matrix(~ 0 + f)
colnames(design) = c("wt","ox")

data.fit = lmFit(data.matrix,design)
data.fit$coefficients[1:10,]

contrast.matrix = makeContrasts(ox-wt,levels=design)
data.fit.con = contrasts.fit(data.fit,contrast.matrix)

data.fit.eb = eBayes(data.fit.con)
names(data.fit.eb)
data.fit.eb$coefficients[1:10,]
data.fit.eb$t[1:10,]
data.fit.eb$p.value[1:10,]

#List of IDs that are DEG (this one is not working)
**IDs.up = rownames(topups)
IDs.down = rownames(topdowns)
ups = subset(DEresults, DEresults[,1]==1)
downs = subset(DEresults, DEresults[,1]==-1)
IDs.up = rownames(ups)
IDs.down = rownames(downs)**

I believe any expert can figure out where I am making the mistakes. I am eagerly looking forward to hearing from you all.

Sincerely, Mahmudul Korea University