**0**wrote:

**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

You say that you need to analyze your own microarray data and that this current dataset is just for practice. Do you have your own data already and is it from an Affymetrix BeadChip?

36k