Agilent single color array

Question

[Leite]

Hello everyone,

I would like your opinion about this code, to analyze Agilent single color array

# Reference: http://matticklab.com/index.php?title=Single_channel_analysis_of_Agilent_microarray_data_with_Limma


#Configure o diretório de trabalho
setwd("C:/User/___")

rm(list=ls(all=TRUE)) 

library(Biobase)
library(GEOquery)
library(limma)

#Read in tab-delimited targets file
targets <- readTargets("targets.txt")

#load data into an RGList object
x <- read.maimages( files=targets, path=".", source= "agilent", green.only= T, columns= list( G= "gProcessedSignal") )

#Pre-processing
ex <- x$E
rownames(ex) <- x$genes$ProbeName
ex <- log(ex,2)
y <- x
y$E <- ex
y <- avereps.EList(y, ID=y$genes$ProbeName) #avereps_H.EList
y$E <- 2 ^ y$E 
y <- backgroundCorrect(y, method="normexp", offset=16)
y <- normalizeBetweenArrays(y, method="quantile") #normalizeBetweenArrays normalizeWithinArrays quantile 
#average replicate spots
y.ave <- avereps(y, ID=y$genes$ProbeName)

#Build the design matrix for the linear modelling function
f <- factor(targets$Target, levels = unique(targets$Target))
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)

#Apply the intensity values to lmFit
fit <- lmFit(y.ave, design)
write.table(fit, file="fit.txt", sep="\t", quote=FALSE)

#Create a contrast matrix
contrast.matrix <- makeContrasts("Disease-Control", levels=design)

#Apply this contrast matrix to the modeled data and compute statistics for the data
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

#Output the statistics for the dataset and write them to disk
output <- topTable(fit2, coef=1, genelist=y.ave$genes, number=Inf, adjust.method="BH", lfc=1.5)
write.table(output, file="Disease-Control.txt", sep="\t", quote=FALSE)

Best regards, Leite

[GenoMax]

Not sure what you mean by your opinion. As noted on the link included above it is code that has worked for this person and you will need to try it to see if it works for you.

@Kevin Blighe has his version of how to analyze Agilent data that you can check out: build the expression matrix step by step from GEO raw data

[Leite]

Thank you @genomax,

What I would like to know is if the code looks correct, or if it has something that looks strange on it, I'm new working with microarray. Sometimes lacking confidence in the codes even if it works correctly, there is something that can be wrong that only the most experienced people can see.

The @kevin tutorial is very good, but it's for two-colour.

Best regards

[GenoMax]

Well try it out. Take a look at limma vignette. Best way to learn is by trying things out.

[Leite]

Thanks again @genomax I will study again Limma users guide!

[Kevin Blighe]

Hi Leite, that pipeline (if you have copied it straight from Mattick's web-site) is almost definitely fine to use. John Mattick is a highly reputable figure who is now based in London - on his web-page (to which you have linked), he additionally mentions assistance from Gordon Smyth, i.e., the person who co-developed limma. Cannot argue with that.

Remember that there are different ways to process microarray data, just as here are different ways to process RNA-seq.

[Leite]

Hey @Kevin

Thank you very much for the comment, I am more "comfortable" in knowing the reputation of this laboratory.

I will test it with other data for further knowledge, as suggested by @genomax.

Leite

Answer 1

[andrew.j.skelton73]

I suggest you follow @Kevin Blighe and @Genomax's helpful suggestions, one thing I did notice is that when you call topTable, your lfc parameter is 1.5, which in log2 space is 2^1.5 = 2.82, so you're actually only looking for genes that are nearly 3 fold different. Remember to wrap your linear space value in log2, so it would be lfc = log2(1.5), and secondly, always visualise your results where you can! - A volcano plot might be a very helpful visualisation for you to understand what's going on in differential expression - here's a good tutorial to get you started

[Leite]

hey @andrew.j.skelton73 Thanks for the reply, I use too p-value = 0.05 as a cutoff, but I do it manually! I I appreciated the website about volcano plot it is very helpful.

Leite