The support.bioconductor.org editor has been updated to markdown! Please see more info at: Tutorial: Updated Support Site Editor

Question: Need expert suggestion about Bioconductor package for microarray data anlaysis
0
gravatar for mahmudornob
4 weeks ago by
mahmudornob0 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

affy limma • 95 views
ADD COMMENTlink modified 4 weeks ago • written 4 weeks ago by mahmudornob0

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?

ADD REPLYlink written 4 weeks ago by Gordon Smyth36k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 199 users visited in the last hour