[BioC] LIMMA on RT-PCR data

[Sandhya Pemmasani Kiran]

Hi Gordon and Ryan,

Thank you so much for all your help.
Relieved after knowing that I can use limma on this gene set, because for microarray I used limma. This is validation experiment with RT PCR. So I wanted to maintain uniformity (if that makes sense)

My experimental design is spread across 85 samples like this-

4 batches- A, B, C and D
In each Batch, 3 treatments - T1, T2, T3
Same patients were taken for treatments of a batch.
I don't have comparisons across the batches, but I have comparisons within batches

SampleName      Batch   Trt     Patient
S1      A       1       1
S2      A       2       1
S3      A       3       1
S4      A       1       2
S5      A       2       2
S6      A       3       2
....................
.....................
S7      D       1       30
S8      D       2       30
S9      D       3       30

Observation are missing on some samples.

TD = factor(paste(Batch, Trt, sep="_"))
Patient <- factor(Patient)
design <- model.matrix(~0+TD+Patient)
fit <- lmFit(dCT, design)

cont.matrix <- makeContrasts(TDA_2 - TDA_1, levels = design)
fit1 = contrasts.fit(fit, cont.matrix)
fit2 = eBayes(fit1)


Is this correct?  I hope I don't need to use duplicateCorrelation as I don't want to compare Batch A with Batch B.

Please advise.


Thanks,
Sandhya






-----Original Message-----
From: Gordon K Smyth [mailto:smyth at wehi.EDU.AU]
Sent: Wednesday, November 13, 2013 5:16 AM
To: Ryan
Cc: Bioconductor mailing list; Sandhya Pemmasani Kiran
Subject: Re: [BioC] LIMMA on RT-PCR data

Hi Ryan,

Here's a simulated example, just to illustrate:

Simulate 50 patients with before and after treatment:

> library(limma)
> y <- matrix(rnorm(4*100),4,100)
> patient <- gl(50,2)
> treatment <- gl(2,1,100)

Drop 15 observations as missing:

> keep <- sample(1:100,85)
> y <- y[,keep]
> patient <- factor(patient[keep])
> treatment <- treatment[keep]

> design <- model.matrix(~treatment)
> dupcor <- duplicateCorrelation(y2,design,block=patient)
> dupcor$consensus
[1] 0.0166576

Theoretically the correlation should be zero.

> fit <- lmFit(y,design,block=patient,correlation=dupcor$consensus)
> fit <- eBayes(fit)
> fit$df.total
[1] 115.5603 115.5603 115.5603 115.5603
> fit$df.prior
[1] 32.56032
> fit$df.residual
[1] 83 83 83 83

The perfect solution here is to pool the variances which would give df.total = 4*83 = 332.  The estimated value is slighty more conservative.
limma takes care never to use df.total greater than the pooled value.

> topTable(fit,coef=2)
         logFC     AveExpr          t    P.Value adj.P.Val         B
4  0.47806544  0.07233048  1.8941941 0.06069901 0.2427960 -3.954093
1 -0.12812909  0.04446857 -0.5713735 0.56885611 0.6946831 -5.112549
3 -0.09334881  0.08623914 -0.4705573 0.63884385 0.6946831 -5.150553
2  0.09463560 -0.00908128  0.3934894 0.69468313 0.6946831 -5.174667

No DE, which is correct.

Cheers
Gordon


On Tue, 12 Nov 2013, Ryan wrote:

> Wow, I wasn't aware that limma could be beneficial for such a small
> number of genes. That's good to know.
>
> On Tue Nov 12 14:50:21 2013, Gordon K Smyth wrote:
>> Hi Sandhya,
>>
>> Yes, limma should work fine on this data, although you are on the
>> lower boundary in terms of number of genes.  Theoretically, the
>> minimum number of genes for the empirical Bayes procedure to be
>> beneficial is three. Four genes is probably the minimum from a
>> practical point of view.
>>
>> You may already know how to use duplicateCorrelation.  If you have a
>> simple before vs after paired comparison with some unmatched samples,
>> then you could proceed like this:
>>
>>  design <- model.matrix(~treatment)
>>  dupcor <- duplicateCorrelation(y, design, block=patient)  fit <-
>> lmFit(y, design, block=patient, correlation=dupcor$consensus)  fit <-
>> eBayes(fit)
>>  topTable(fit,coef=2)
>>
>> Be sure to check that dupcor$consensus is greater than zero.
>>
>> We used this strategy to compare tumour vs normal tissue in the
>> presence of unmatched samples in
>>
>>  http://www.ncbi.nlm.nih.gov/pubmed/17236199
>>
>> although that was microarray data rather than PCR.
>>
>> Best wishes
>> Gordon
>>
>>
>> ---------- original message ---------- BioC] LIMMA on RT-PCR data
>> Sandhya Pemmasani Kiran sandhya.p at ocimumbio.com Tue Nov 12
>> 12:16:25 CET 2013
>>
>> Dear list,
>>
>> I have RT-PCR data on 4 genes and 85 samples.
>> Can I use 'limma' on this small set of genes.
>>
>> I want to use limma rather than usual paired t-test because I have
>> missing values and I don't want to miss the information available on
>> paired samples..
>>
>> Please advise.
>>
>>
>> Thanks,
>> Sandhya

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:25}}