Question: paired microarray samples
0
Chris Fjell • 60 wrote:
I've a question on how to properly deal with paired samples in limma.
For example, I've got an experiment with 3 blood donors (D1, D2, D3)
and
4 treatments (T1, T2, T3, T4) of the blood cells. One singlechannel
array is done for each donortreatment (e.g. D1 cells treated with
T2).
In the past I've treated these as simple biological replicates, so the
code flow is like:
treatments=pData(BSData.quantile)$Sample_Group # T1, T2, T3, T4, from
Illumina beadarray data
treatments=as.factor(treatments)
design=model.matrix(~0+treatments)
colnames(design)=levels(treatments)
fit = lmFit(exprs(BSData.quantile), design)
cont.matrix = makeContrasts(T1vT2 = T1  T2, T3vT1 = T3  T1, levels =
design)
fit = contrasts.fit(fit, cont.matrix)
ebFit = eBayes(fit)
for( comp in colnames( cont.matrix ) ) {
tt = topTable(ebFit, coef = comp)
...
[write results to file...]
}
But donors are responding quite differently so I don't want to simple
pool all the treatments, I want to consider response of each donor,
i.e.
"Paired Samples".
I found this snippet on the BioC newsgroup at
https://stat.ethz.ch/pipermail/bioconductor/2006May/012969.html
Calculate the logFC ratios manually per donor, then pass it to lmFit()

seems like the standard paired ttest.
For my example, for the comparison of T1 to T2
exM = exprs(BSData.quantile)
diffD1 = exM[, which( donors == "D1" & samples == "T1" )]  exM[,
which(
donors == "D1" & samples == "T2" )]
diffD2 = exM[, which( donors == "D2" & samples == "T1" )]  exM[,
which(
donors == "D2" & samples == "T2" )]
diffD3 = exM[, which( donors == "D3" & samples == "T1" )]  exM[,
which(
donors == "D3" & samples == "T2" )]
paired_ex = cbind( diffD1, diffD2, diffD3)
fit = lmFit(paired_ex)
The example cited above also uses "block = origin, correlation =
dupcor$cons" parameters to lmFit() but this is for technical
replicates (?)
This works but I get dramatically worse pvalues (due I assume to
using
half the number of "samples": 3 ratios tested against 0 versus
testing
two groups of 3).
But the limma user's guide (Oct 22, 2008 version, Chapter 8, p.40)
gives
another example, using model.matrix() not makeContrasts():
> SibShip < factor(targets$SibShip)
> Treat < factor(targets$Treatment, levels=c("C","T"))
> design < model.matrix(~SibShip+Treat)
> fit < lmFit(eset, design)
> fit < eBayes(fit)
> topTable(fit, coef="TreatT")
I do the similar thing for my data ...
design_paired_treatments = model.matrix( ~donors+treatments )
fit_ps = lmFit( exprs(BSData.quantile),design_paired_treatments)
fit_ps = eBayes( fit_ps )
and look at the columns of my design matrix and it's like
colnames( design_paired_treatments )
[1] "(Intercept)" "donorsD1" "donorsD2" "treatmentsT1"
[5] "treatmentsT2" "treatmentsT3"
Donor D3 and treatment T4 don't appear in the design.
and I get some results...
tt = topTable( fit_ps, coef="donorsD1")
tt
ID logFC AveExpr t P.Value adj.P.Val B
6370315 6.58 8.74 119.4 4.26e17 2.08e12 17.6
7160474 3.59 8.77 56.5 7.43e14 1.37e09 16.6
5260484 5.69 8.71 55.8 8.39e14 1.37e09 16.6
This is a regression model, I think... I can see why the two terms
(donor D3 and treatment T4) disappeared  they have to be restated in
terms of the rest, if I remember right for doing ANOVA with indicator
variables.
Is there a recommended way to do this with limma? Or another package?
Thanks for any advice.
Chris
ADD COMMENT
• link
•
modified 9.2 years ago
by
Naomi Altman • 6.0k
•
written
9.2 years ago by
Chris Fjell • 60