[BioC] (no subject)
Naomi Altman
naomi at stat.psu.edu
Fri Mar 12 18:28:06 CET 2010
The estimated error variance used for the test denominator will be an
average of technical and biological replication, and therefore not
really appropriate for your analysis. However, you could average the
2 technical replicates prior to running limma which would give you
the right error structure.
--Naomi
At 12:04 PM 3/12/2010, Ana Staninska wrote:
>Dear Bioconductor,
>I have a simple experiment that I have to analyze in order to find
>differentially expressed genes. I have 10 biological replicates, and
>each biological replicate has two technical replicates which appear
>as dye swapped. So in total I have 20 arrays. Each of the probes are
>spotted twice on the array (on the left and on the right hand side).
>I use limma to do my analysis. I know at the moment it is not
>possible to treat duplicate spots, technical replicates and
>biological replicates, but I though if I use the
>duplicateCorrelation function on my duplicate spots, and then to use
>a contrast matrix to average of all of the Treated vs Non-treated
>biological samples, I could address all 3 replications. Am I correct?
>
>
>I am sending a copy of my code, if someone could look at it at tell
>me whether I made somewhere a mistake.
>Thank you very much in advance,
>Sincerely Ana Staninska
>
>
> library(limma)> library(statmod)> library(marray)>
> library(convert)> library(hexbin)> library(gridBase)>
> library(RColorBrewer)> > targets <-
> readTargets("Lysi_270705.txt")> > ### Only manually removed ot
> absent spots are given 0 weight ###> RGa <- read.maimages(targets,
> source="genepix", wt.fun=wtflags(weight=0,
> cutoff=-75), other.columns=c("F635 SD","B635 SD","F532 SD","B532
> SD","B532 Mean","B635 Mean","F Pixels","B Pixels"))Read
> LYSI270705_1_200905.gpr Read LYSI270705_1dw_200905.gpr Read
> LYSI270705_2_200905.gpr Read LYSI270705_2dw_200905.gpr Read
> LYSI270705_3_121005.gpr Read LYSI270705_3dw_121005.gpr Read
> LYSI270705_4_121005.gpr Read LYSI270705_4dw_121005.gpr Read
> LYSI270705_5_121005.gpr Read LYSI270705_5dw__121005.gpr Read
> LYSI270705_6_121005.gpr Read LYSI270705_6dw__121005.gpr Read
> LYSI270705_7_151001.gpr Read LYSI270705_7dw_151005.gpr Read
> LYSI270705_8_151005.gpr Read LYSI270705_8dw_151005.gpr Read
> LYSI270705_9_151005.gpr Read LYSI270705_9dw_151005.gpr Read LYSI270705!
> _10_151005.gpr Read LYSI270705_10dw_151005.gpr > for(i in
> 1:nrow(RGa)){+ for(j in
> 1:ncol(RGa)){+ if(RGa$Rb[i,j]+RGa$R[i,j]+ RGa$G[i,j]+
> RGa$Gb[i,j] ==0)+ RGa$weights[i,j]<-0+ }+ }> >
> ####################################################> ###
> Background Correction = Normexp + offset 25 ####>
> ####################################################> > RG
> <-backgroundCorrect(RGa, method="normexp", , normexp.method="mle",
> offset=25)Green channelCorrected array 1 Corrected array 2
> Corrected array 3 Corrected array 4 Corrected array 5 Corrected
> array 6 Corrected array 7 Corrected array 8 Corrected array 9
> Corrected array 10 Corrected array 11 Corrected array 12 Corrected
> array 13 Corrected array 14 Corrected array 15 Corrected array 16
> Corrected array 17 Corrected array 18 Corrected array 19 Corrected
> array 20 Red channelCorrected array 1 Corrected array 2 Corrected
> array 3 Corrected array 4 Corrected array 5 Corrected array 6
> Corrected array 7 Corrected array 8 Corrected array !
> 9 Corrected array 10 Corrected array 11 Corrected array 12 Corrected a
>rray 13 Corrected array 14 Corrected array 15 Corrected array 16
>Corrected array 17 Corrected array 18 Corrected array 19 Corrected
>array 20 > ####################################################>
>##### normalize Within arrays #########>
>####################################################> > MA
><-normalizeWithinArrays(RG, method="loess")> >
>####################################################> ######
>Contrast Matrix ############>
>####################################################> >
>design<-cbind( + MU1vsWT1=c(
>1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),+ MU2vsWT2=c(0,0,
>1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),+ MU3vsWT3=c(0,0,0,0,
>1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0),+ MU4vsWT4=c(0,0,0,0,0,0,
>1,-1,0,0,0,0,0,0,0,0,0,0,0,0),+ MU5vsWT5=c(0,0,0,0,0,0,0,0,
>1,-1,0,0,0,0,0,0,0,0,0,0),+ MU6vsWT6=c(0,0,0,0,0,0,0,0,0,0,
>1,-1,0,0,0,0,0,0,0,0),
>+ MU7vsWT7=c(0,0,0,0,0,0,0,0,0,0!
> ,0,0,
> 1,-1,0,0,0,0,0,0),+ MU8vsWT8=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,
> 1,-1,0,0,0,0),+ MU9vsWT9=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
> 1,-1,0,0),+ MU10vsWT10=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
> 1,-1))> > cont.matrix <-
> makeContrasts(MUvsWT=(MU1vsWT1+MU2vsWT2+MU3vsWT3+MU4vsWT4+MU5vsWT5+MU6vsWT6+MU7vsWT7+MU8vsWT8+MU9vsWT9+MU10vsWT10)/10,
> levels=design)> >
> ####################################################>
> ### Duplicate Correlations on duplicate spots ####>
> ####################################################> >
> corfit<-duplicateCorrelation(MA, ndups=2, spacing=192)> >
> ####################################################> ##### Linear
> Fit Model and Contrasts fit #######>
> ####################################################> >
> fit<-lmFit(MA, design, ndups=2, spacing=192,
> cor=corfit$consensus)> > fit<-contrasts.fit(fit, cont.matrix)> >
> ####################################################>
> ######### eBayes Statistics ###############> #################!
> ###################################> > fit<-eBayes(fit)> > ###########
>###################################################> ### Writing
>the Results ######>
>##############################################################>
>TTnew<-topTable(fit,coef=1, number=100, adjust="BH")
>
>
>Ana StaninskaHelmholtz-Zentrum MuenchenDepartment of Scientific
>ComputingNeuherberg, Deutschland+49 (0) 89 3187 2656
>
> [[alternative HTML version deleted]]
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at stat.math.ethz.ch
>https://stat.ethz.ch/mailman/listinfo/bioconductor
>Search the archives:
>http://news.gmane.org/gmane.science.biology.informatics.conductor
Naomi S. Altman 814-865-3791 (voice)
Associate Professor
Dept. of Statistics 814-863-7114 (fax)
Penn State University 814-865-1348 (Statistics)
University Park, PA 16802-2111
More information about the Bioconductor
mailing list