[BioC] Problem running eBayes/lmFit

[elliott harrison]

Hi BioC,
I'm trying to do some differential expression analysis with one colour
agilent data.
What I have is 4 CY5 arrays in duplicate. They were both hybridised
slightly differently and I want to check that the data generated was
unaffected by the change. The working premise is that the differentially
expressed genes between 2 arrays done under the first condition will be
replicated under the second.

These are the 8 slides. A9802 = SAMPLE1, A9811 = SAMPLE2 and so on.
A9802 A9811 A9813 A9842 SAMPLE1 SAMPLE2 SAMPLE3 SAMPLE4

The data was loaded like this

> RTestR <- read.maimages(targets2, columns = list(G = "rMeanSignal", Gb
= "rBGUsed", R = "rMeanSignal", Rb = "rBGUsed"),annotation= c("Row",
"Col", "FeatureNum", "ProbeUID","ControlType","ProbeName", "GeneName",
"SystematicName"), wt.fun=myFlagFun)

myFlagFun looks like this

function(x) {
   #Weight only strongly positive spots 1, everything else 0
   present <- x$rIsPosAndSignif == 1
   probe <- x$ControlType == 0
   manual <- x$IsManualFlag == 0
   strong <- x$rIsWellAboveBG == 1
   y <- as.numeric(present & probe & manual & strong)

   #Weight weak spots 0.5

   weak <- strong == FALSE
   weak <- (present & probe & manual & weak)
   weak <- grep(TRUE,weak)
   y[weak] <- 0.5

   #Weight flagged spots 0.5

   sat <- x$rIsSaturated == 0
   xdr <- x$rIsLowPMTScaledUp == 0
   featureOL1 <- x$rIsFeatNonUnifOL == 0
   featureOL2 <- x$rIsFeatPopnOL == 0
   flagged <- (sat & xdr & featureOL1 & featureOL2)
   flagged <- grep(FALSE, flagged)
   good <- grep(TRUE, y==1)
   flagged <- intersect(flagged, good)
   y[flagged] <- 0.5
   y
   }

Combined the arrays

> RTotal =
cbind(RTestR[,1],RTestR[,2],RTestR[,3],RTestR[,4],RTestOriR[,1],RTestOri
R[,2],RTestOriR[,3],RTestOriR[,4])

I've checked to see I'm not just weighting everything out as 0
> summary(RTotal$weight>0)
 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\SingleChanneltest\\SAMPLE
1\\test01_251485017912_S01_GE2-v5_91_0806_1_1
 Mode :logical

 FALSE:8796

 TRUE :36219

 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\SingleChanneltest\\SAMPLE
2\\test01_251485017912_S01_GE2-v5_91_0806_1_2
 Mode :logical

 FALSE:10350

 TRUE :34665

 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\SingleChanneltest\\SAMPLE
3\\test01_251485017912_S01_GE2-v5_91_0806_1_3
 Mode :logical

 FALSE:8785

 TRUE :36230

 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\SingleChanneltest\\SAMPLE
4\\test01_251485017912_S01_GE2-v5_91_0806_1_4
 Mode :logical

 FALSE:7541

 TRUE :37474

 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\980\\test01_251485020980_S01_GE2-v5_91_0806_2_1_2
 Mode :logical

 FALSE:10373

 TRUE :34642

 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\981\\981\\test01_251485020981_S02_GE2-v5_91_0806_2_1_1
 Mode :logical

 FALSE:12630

 TRUE :32385

 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\981\\981\\test01_251485020981_S02_GE2-v5_91_0806_2_1_3
 Mode :logical

 FALSE:12492

 TRUE :32523

 C:\\Program Files\\Agilent\\GeneSpring
GX\\data\\CRX\\984\\984\\test01_251485020984_S02_GE2-v5_91_0806_2_1_2
 Mode :logical

 FALSE:9496

 TRUE :35519     


Bg correct and normalise

> RTotalbg = backgroundCorrect(RTotal, method="none")
> RTRN <- normalizeBetweenArrays(RTotalbg$R, method="quantile")

Comparing the slides to one another so get slide numbers from targets

> f <- paste(RTotal$targets$SlideNumber,sep="")
> f <- factor(f)
> design <- model.matrix(~0+f)
> colnames(design) <- levels(f)
> design
  A9802 A9811 A9813 A9842 SAMPLE1 SAMPLE2 SAMPLE3 SAMPLE4
1     0     0     0     0       1       0       0       0
2     0     0     0     0       0       1       0       0
3     0     0     0     0       0       0       1       0
4     0     0     0     0       0       0       0       1
5     1     0     0     0       0       0       0       0
6     0     1     0     0       0       0       0       0
7     0     0     1     0       0       0       0       0
8     0     0     0     1       0       0       0       0
attr(,"assign")
[1] 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$f
[1] "contr.treatment"

>fit <- lmFit(RTRN, design)
>cont.matrix <-
makeContrasts(OneVOne="A9802-SAMPLE1",OneVTwo="A9802-A9811",TwoVOne="SAM
PLE1-SAMPLE2",levels=design)

The matrix seems to be doing what I want

> cont.matrix
         Contrasts
Levels    OneVOne OneVTwo TwoVOne
  A9802         1       1       0
  A9811         0      -1       0
  A9813         0       0       0
  A9842         0       0       0
  SAMPLE1      -1       0       1
  SAMPLE2       0       0      -1
  SAMPLE3       0       0       0
  SAMPLE4       0       0       0

> fit2  <- contrasts.fit(fit, cont.matrix)
> fit2  <- eBayes(fit2)
Error in ebayes(fit = fit, proportion = proportion, stdev.coef.lim =
stdev.coef.lim) : 
        No residual degrees of freedom in linear model fits

I've found a post that says this error message occurs because all data
is weighted out. I've checked the data after it is loaded, after
backgroundCorrect and it does not appear to be. Beyond that I doesn't
look like the normalizeBetweenArrays of RTotalbg$R RTRN has any weights.
So I must not be setting up the design matrix correctly?

Any and all clues as to where I'm going wrong greatly appreciated.


Elliott Harrison




This message has been scanned for viruses by BlackSpider...{{dropped:3}}