[BioC] R: R: design matrix with technical and biologial replicates
Manuela Di Russo
manuela.dirusso at for.unipi.it
Fri Apr 20 16:06:34 CEST 2012
Dear James,
if I have correctly interpreted, your suggestion is to fit the model to all
probes on the array and then filter. But the approach usually adopted is
based on filtering before fitting the linear model in order to reduce the
number of hypotheses to test and increase the power. May I filter after
fitting the model using the information about higher-variance/lower
intensity probesets obtained applying genefilter() (for example) to the
ExpressionSet object created by gcrma/rma function?
Sorry but you are the first person who tells me to filter the data after
fitting the model in order not to bias the estimates of the prior variance.
Thank you very much for your help and your time.
Manuela
----------------------------------------------------------------------------
----------
Manuela Di Russo, Ph.D. Student
Department of Experimental Pathology, MBIE
University of Pisa
Pisa, Italy
e-mail: manuela.dirusso at for.unipi.it
mobile: +393208778864
phone: +39050993538
-----Messaggio originale-----
Da: James W. MacDonald [mailto:jmacdon at uw.edu]
Inviato: mercoledì 18 aprile 2012 16:19
A: Manuela Di Russo
Cc: bioconductor at r-project.org
Oggetto: Re: R: [BioC] design matrix with technical and biologial replicates
Hi Manuela,
On 4/18/2012 10:11 AM, Manuela Di Russo wrote:
> Thank you James!
> I had already applied both methods you suggested but I wanted to see
> if there was a better way to handle this kind of experimental design.
> I have another question: I computed means of the duplicates arrays
> after preprocessing and filtering but before fitting the linear model,
> is this correct?
I would argue that you should compute means, then fit the model (to the
eBayes step), and then filter. This is because the eBayes() step estimates a
prior variance estimate over all probes. Here is an example:
Say you filter your probes, selecting only those with an inter-quartile
range greater than some threshold. In that case you are removing all probes
with low variance. Now when you compute the prior in the eBayes() step, you
are basing it on the subset of higher-variance probes. One could argue that
you have biased that estimate upward, which will in turn decrease your power
to detect differences.
Best,
Jim
> Thank you!
> Manuela
> ----------------------------------------------------------------------
> ------
> ----------
> Manuela Di Russo, Ph.D. Student
> Department of Experimental Pathology, MBIE University of Pisa Pisa,
> Italy
> e-mail: manuela.dirusso at for.unipi.it
> mobile: +393208778864
> phone: +39050993538
>
> -----Messaggio originale-----
> Da: James W. MacDonald [mailto:jmacdon at uw.edu]
> Inviato: mercoledì 18 aprile 2012 15:43
> A: Manuela Di Russo
> Cc: bioconductor at r-project.org
> Oggetto: Re: [BioC] design matrix with technical and biologial
> replicates
>
> Hi Manuela,
>
> On 4/18/2012 7:52 AM, Manuela Di Russo wrote:
>> Dear list,
>>
>> I'm working with microarray expression data and I am using limma to
>> detect differentially expressed genes. I have some questions about
>> the design matrix and the handling of biological and technical
replicates.
>>
>> The target file is:
>>
>> Sample_name sample_type
> sample_replicate
>> disease_status
>>
>> MPM_07 1
>> 1 1
>>
>> MPM_08 1
>> 2 1
>>
>> MPM_09 1
>> 3 1
>>
>> MPM_10_a 1
>> 4 1
>>
>> MPM_10_b 1
>> 4 1
>>
>> MPM_11 1
>> 5 1
>>
>> MPM_12 1
>> 6 1
>>
>> PP_01_a 2
>> 7 0
>>
>> PP_01_b 2
>> 7 0
>>
>> PP_02 2
>> 8 0
>>
>> PP_03 2
>> 9 0
>>
>> PP_04 2
>> 10 0
>>
>> PP_05 2
>> 11 0
>>
>> PP_06 2
>> 12 0
>>
>> PV_02 3
>> 13 0
>>
>> PV_03 3
>> 14 0
>>
>> PV_04 3
>> 15 0
>>
>> PV_05 3
>> 16 0
>>
>> Each sample is hybridized on an Affymetrix HG-U133-Plus2 array.
>>
>> So I have 7 mesothelioma samples (sample_type=1) where 2 were from
>> the same patient (MPM_10 a e b)), 7 parietal pleural samples
>> (sample_type=
>> 2) where 2 were from the same patient (PP_01 a e b) and 4 visceral
>> pleural samples (sample_type= 3). In reality 4 parietal pleural
>> samples (PP_02,PP_03,PP_04 and PP_05) and 4 visceral pleural samples
>> (PV_02,PV_03,PV_04 and PV_05) come from the same patients.
>>
>> pd<- data.frame(sample_type= c(rep(1,7),rep(2,7),rep(3,4)),
>> sample_replicate = c(1:4,4,5,6,7,7,8:12,13:16),
>> disease_status=c(rep(1,7),rep(0,11)))
>>
>> biolrep<-pd$sample_replicate
>>
>> f<- factor(pd$sample_type)
>>
>> design<- model.matrix(~0+f)
>>
>> colnames(design)<- c("MPM", "PP", "PV")
>>
>> I tried to handle technical replicates using the block argument of
>> function
>> duplicatecorrelation() as follows:
> I don't think you can use duplicateCorrelation() here, as you don't
> have duplicates for all samples. I believe lmFit() with a cor argument
> will fit a block diagonal correlation matrix, which is clearly not
> applicable here. I may be in error however, in which case Gordon Smyth
> will surely post a correction around 5-6 pm EDT or so.
>
> With a mixture of duplicated and not duplicated samples, you will
> likely have to do one of two less than ideal things. First, you could
> simply ignore the duplication, and analyze as if the duplicates were
independent samples.
> This is less than ideal because there will be a correlation between
> these samples, which will tend to lower your estimate of intra-sample
variation.
>
> Second, you could compute means of the duplicates and then use those
> in lieu of the original data. Again, this is not ideal, as the means
> will have an intrinsically lower variance than individual samples. All
> things equal, this is probably the better way to go.
>
> Best,
>
> Jim
>
>
>> corfit<- duplicateCorrelation(eset_norm_genes_ff_filtered, design,
>> ndups=1, block= biolrep) # eset_norm_genes_ff_filtered is an
>> ExpressionSet object containing pre-processed and filtered data
>>
>> I am interested in identifying differentially expressed genes between
>> MPM and PP and between PV and PP.
>>
>> contrast.matrix_all.contrasts<-
>> makeContrasts(MPMvsPP=MPM-PP,PVvsPP=PV-PP,levels=design)
>>
>> fit_ff<-lmFit(eset_norm_genes_ff_filtered, design,block=biolrep,
>> ndups=1,cor=corfit$consensus)
>>
>> fit2_ff<- contrasts.fit(fit_ff, contrast.matrix_all.contrasts)
>>
>> fit2e_ff<-eBayes(fit2_ff)
>>
>> I think that my approach is correct for the first contrast (MPM vs
>> PP) but not for the second one because biolrep doesn't consider the
>> fact that some samples between PP and PV are paired.
>>
>> Am I correct?
>>
>> What about defining biolrep<-c(1:4,4,5,6,7,7,8:12,8:11)?
>>
>> Is there a method to handle such an experimental design?
>>
>> Sorry for my long post!
>>
>> Any suggestion/comment is welcome.
>>
>> Cheers,
>>
>> Manuela
>>
>>
>>
>> ---------------------------------------------------------------------
>> -
>> ------
>> ----------
>>
>> Manuela Di Russo, Ph.D. Student
>> Department of Experimental Pathology, MBIE University of Pisa Pisa,
>> Italy e-mail:<mailto:manuela.dirusso at for.unipi.it>
>> manuela.dirusso at for.unipi.it
>> mobile: +393208778864
>>
>> phone: +39050993538
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
> --
> James W. MacDonald, M.S.
> Biostatistician
> University of Washington
> Environmental and Occupational Health Sciences
> 4225 Roosevelt Way NE, # 100
> Seattle WA 98105-6099
>
--
James W. MacDonald, M.S.
Biostatistician
University of Washington
Environmental and Occupational Health Sciences
4225 Roosevelt Way NE, # 100
Seattle WA 98105-6099
More information about the Bioconductor
mailing list