limma: duplicate spots, technical replicates and biological replicates
0
0
Entering edit mode
@ana-staninska-3957
Last seen 10.2 years ago
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 array 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+MU4vsW T4+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]]
limma limma • 775 views
ADD COMMENT

Login before adding your answer.

Traffic: 701 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6