Entering edit mode
De Boever Patrick
▴
60
@de-boever-patrick-3981
Last seen 10.2 years ago
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? Does the lmFit,contrastsFit
need additional statements?
-what happens if a gene is flagged 'bad' on one array, but not on the
other?
-is there a way to verify that only 'good' genes 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.
Am I missing an argument?
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@vito.be<mailto:patrick.deboever@vito.be>
Visit our website: www.vito.be/english<http: www.vito.be="" english="">
---
This e-mail, any attachments and the information it contains are
confidential and meant only for the use of the addressee(s) only.
Access to this e-mail by anyone other than the addressee(s) is
unauthorized. If you are not the intended addressee (or responsible
for delivery of the message to such person), you may not use, copy,
distribute or deliver to anyone this message (or any part of its
contents) or take any action in reliance on it. In such case, you
should destroy this message and notify the sender immediately. If you
have received this e-mail in error, please notify us immediately by
e-mail or telephone and delete the e-mail from any computer.
All reasonable precautions have been taken to ensure no viruses are
present in this e-mail and its attachments. As our company cannot
accept responsibility for any loss or damage arising from the use of
this e-mail or attachments we recommend that you subject these to your
virus checking procedures prior to use.
[[alternative HTML version deleted]]