Question: LIMMA and Agilent single color
1
gravatar for Gordon Smyth
8.8 years ago by
Gordon Smyth38k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth38k wrote:

Dear Patrick,

Before I answer your questions, can I make a comment on the practice of flagging bad spots? Your code is flagging spots as bad just because they are give low intensities. I don't agree with this practice, because it may be a perfectly good observation of a gene that is not expressed in a particular condition. While it is good practice to filter genes that are not expressed in any condition in your experiment, this means filtering whole rows of your data. I don't think that the practice of filtering individual spots on the basis of intensity is a useful one.

Have said that, I've interpolated answers to your questions below.

> Date: Fri, 17 Dec 2010 17:32:14 +0100
> From: De Boever Patrick <patrick.deboever at vito.be>
> To: "bioconductor at r-project.org" <bioconductor at r-project.org>
> Subject: [BioC] LIMMA and Agilent single color
>
> Dear list members,
>
> My question relates to processing of Agilent single-color arrays. With
> data generated using Agilent's feature extraction software.

> My script so-for is mentioned below.

> The flag information contained in raw data files is transformed to
> weights. I have used 1 for 'good' and 0 for 'bad' in myFlagFun.
>
> If this weight information is loaded:

> -where does LIMMA use this information?

Pretty much everywhere.

> Does the lmFit,contrastsFit need additional statements?

No.  The weights are found automatically, provided that you pass the data objects to the functions. Of course, if you extract the expression matrix from the data object, then you have to pass the weights explicitly to lmFit().

> -what happens if a gene is flagged 'bad' on one array, but not on the
> other?

Remember that it is spots that are being flagged, not genes.  For gene (probe), the good spots are used and the bad spots are not.

> -is there a way to verify that only 'good' genes have been used?

Not directly.  But summary(fit$df.residual) will show the residual degrees of freedoms unequal if variable numbers of values have been used.

> Second question,

> My E.avg (MAlist) contains E.avg$genes with gene information, But after
> applying fit->this gene information is lost? No gene annotation my
> topTable.

This is because you didn't pass the information about the genes to lmFit().  If you use

fit <- lmFit(E.avg, design)

then it will automatically find both the weights and the gene annotation. But if you use,

fit <- lmFit(E.avg$M, design)

then you're passing only the expression log-ratios to lmFit, and it has no way to find the gene annotation.

> Am I missing an argument?

You're not making full use of the data objects.  limma now has better features to handle single-channel Agilent arrays.  The following code would be more direct:

targets <- readTargets("targets.txt")
rawObj <- read.maimages(targets,source="agilent",
                        green.only=TRUE,wt.fun=myFlagFun)
Obj.corrected <- backgroundCorrect(rawObj, method="normexp", offset=1)
E <- normalizeBetweenArrays(Obj.corrected, method="quantile")
E.avg <- avereps(E, ID=E$genes$ProbeName)
fit <- lmFit(E.avg,design)
cont.matrix <- makeContrasts(group1vsgroup2=Group1-Group2,levels=design)
fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)
topTable(fit2)

Note that avereps() doesn't use the weights when making averages.  If you average two spots, one with weight 0 and one with weight 1, you'll get the average of the two with weight 0.5.

Best wishes
Gordon

> Thank you for providing insight !
>
> Patrick
>
>
> setwd('D:/VITO/R')
> library('Biobase')
> library('convert')
> library('limma')
> AgilentFiles <- list.files(pattern="US")
> myFlagFun <- function(x) {
> #Weight only strongly positive spots 1, everything else 0
> present <- x$gIsPosAndSignif == 1
> probe <- x$ControlType == 0
> manual <- x$IsManualFlag == 0
> strong <- x$gIsWellAboveBG == 1
> y <- as.numeric(present & probe & manual & strong)
>
> #Weight weak spots 0
>
> weak <- strong == FALSE #with values not well above background
> weak <- (present & probe & manual & weak)
> weak <- grep(TRUE,weak)
> y[weak] <- 0
>
> #Weight flagged spots 0
>
> sat <- x$gIsSaturated == 0
> xdr <- x$gIsLowPMTScaledUp == 0
> featureOL1 <- x$gIsFeatNonUnifOL == 0
> featureOL2 <- x$gIsFeatPopnOL == 0
> flagged <- (sat & xdr & featureOL1 & featureOL2)
> flagged <- grep(FALSE, flagged)
> good <- grep(TRUE, y==1)
> flagged <- intersect(flagged, good)
> y[flagged] <- 0
> y
> }
>
> targets <- readTargets("targets.txt")
> rawObj<-read.maimages(AgilentFiles, columns = list(G ="gMeanSignal", Gb = "gBGUsed",
                                            R ="gProcessedSignal", Rb ="gBGMedianSignal"),
> annotation= c("Row", "Col", "FeatureNum", "ProbeUID", "ControlType","ProbeName",
                      "GeneName", "SystematicName"), wt.fun=myFlagFun)
> Obj.corrected <- backgroundCorrect(rawObj, method="normexp",offset=1)
> Obj<-Obj.corrected
> Obj$R <- normalizeBetweenArrays(Obj.corrected$R, method="quantile")
> Obj$R <- log2(Obj$R)
>
> E <- new("MAList", list(targets=Obj$targets, genes=Obj$genes,
weights=Obj$weights, source=Obj$source, M=Obj$R, A=Obj$G))
> E.avg <- avereps(E, ID=E$genes$ProbeName)
>
> design<- as.matrix(read.table("targets.txt",row.names="FileName",header=T))
> fit<-lmFit(E.avg$M,design, weights=E.avg$weights)
> cont.matrix <- makeContrasts(group1vsgroup2=Group1-Group2,levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
> topTable(fit2, genelist=fit2$genes, adjust="BH")
>
>
> Patrick De Boever, PhD, MSc
> Flemish Institute for Technological Research (VITO)
> Unit Environmental Risk and Health, Toxicology group
> Industriezone Vlasmeer 7, 2400 Mol, Belgium
> Tel.   + 32 14 33 51 45
> Fax.  + 32 14 58 05 23
> patrick.deboever at vito.be
>
> Visit our website: www.vito.be/english

annotation probe limma agilent • 3.1k views
ADD COMMENTlink modified 2.7 years ago • written 8.8 years ago by Gordon Smyth38k
Answer: LIMMA and Agilent single color
1
gravatar for De Boever Patrick
8.8 years ago by
De Boever Patrick60 wrote:

Dear Gordon,

Thank you for the support and good remarks. You are right low intensity genes should be kept in the analysis. Also, very useful to know that a direct load for Agilent single color exists.

Best wishes,

Patrick

ADD COMMENTlink modified 2.7 years ago by Gordon Smyth38k • written 8.8 years ago by De Boever Patrick60
Answer: LIMMA and Agilent single color
0
gravatar for amalbouzid.ing
2.7 years ago by
amalbouzid.ing0 wrote:

Hi,

I'm studying Agilent Human Gene Expression one color microarray using Limma. It results significants probesets with topTable function but then I didn't succed to do the conversion of my probenames with "ids2indices" function? nor the annotation of my genes?

Please, there is someone who can help me to solve my problem??

Thanks,

Best,

Amal

ADD COMMENTlink written 2.7 years ago by amalbouzid.ing0
1

Dear Amal,

You've posted an "answer", but really you are asking a new question rather than answering the above question by Patrick de Boever. Your question apparently relates to gene set testing in limma rather than to any issue raised by Patrick, so please post your query to Bioconductor as a fresh question rather than continuing this thread.

Also, we can't tell you what you've done incorrectly if you don't tell us what you've done. Please show us the code you have used, otherwise no one will be able to help you. Have a read of the posting guide at

  https://www.bioconductor.org/help/support/posting-guide/

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by Gordon Smyth38k

Dear Gordon,

Thanks for your reply and your consideration.

I've already posted my question to Bioconductor! 

Looking forward to your help :)

Thanks,

Best regards,

ADD REPLYlink written 2.7 years ago by amalbouzid.ing0
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: 141 users visited in the last hour