Entering edit mode
@santy-marques-ladeira-6408
Last seen 10.5 years ago
Dear Limma users,
I'm currently working on Affymetrix microarray (Hu-gene 2.1) with only
one
channel. I got 8 experimental stimuli conducted on 7 donors. During
the
experimental process each stimulus was conducted at the same time for
all
the donors so I assume there are no technical replicates, but I need
to
take in account the "donor effect" (or biological replicates).
In the "Limma" documentation there are 2 ways of including this
"effect"
but it is still not clear to me what is the best choice, and the 2
methods
gave me 86% of correlation for a pairwise comparison test.
*1/The first method I used was including the donors in my design model
:*
# I build a matrix with 2 columns (Stimulus + Donor)
# Each row is an array
sample = matrix(0, nrow=55, ncol=2)
rownames(sample) <- colnames("normalised data I use to perform the
analysis")
colnames(sample) <- c("Stimulus", "Donor")
# In the column Stimulus there is directly the name of the stimulus
tested
sample[1:7,1] = "Stim1"
sample[8:14,1] = "Stim2"
sample[15:20,1] = "Stim3"
sample[21:27,1] = "Stim4"
sample[28:34,1] = "Stim5"
sample[35:41,1] = "Stim6"
sample[42:48,1] = "Stim7"
sample[49:55,1] = "Stim8"
# In the column Donor there is the number of the donor (from 2 to 8)
sample[,2] =
c(rep(c(2,3,4,5,6,7,8),2),3,4,5,6,7,8,rep(c(2,3,4,5,6,7,8),5))
targets = as.data.frame(sample)
Stimulus = factor(targets$Stimulus, levels = c("Stim1", "Stim2",
"Stim3",
"Stim4", "Stim5", "Stim6", "Stim7", "Stim8"))
Donor = factor(targets$Donor, levels = c("2", "3", "4", "5", "6", "7",
"8"))
# I use the model.matrix() function to infer the design based on my
matrix
design <- model.matrix(~-1+Stimulus+Donor)
rownames(design) <- colnames("normalised data I use to perform the
analysis")
colnames(design)[1:8] <- c("Stim1", "Stim2", "Stim3", "Stim4",
"Stim5",
"Stim6", "Stim7", "Stim8")
#This gives me :
# Stim1 Stim2 Stim3 ... ... ...
Donor3
Donor4 ... ... Donor8
#experience1 1 0 0
0 0 0
#experience2 1 0 0
1 0 0
#experience3 1 0 0
0 1 0
#...
#...
#...
First thing is that because I remove for a stimulus the donor2, this
donor
is removed from the design model (maybe because of dimension for
regression
?)
Then I do the usual commands :
fit <- lmFit("normalised data I use to perform the analysis", design)
cont.matrix<-makeContrasts(Stim2vsStim4=Stim2-Stim4, levels=design)
fit2<-contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
anadiff.limma.Stim2<-topTable(fit2, number=1000, adjust="BH")
*2/The second method I used is by treating donors as replicates
separately:*
I used the exact same "targets" data.frame, except there is no more
the
Donor column but only the Stimulus column
[...]
design <- model.matrix(~-1+Stimulus)
[...]
# List of donors saved as a variable for blocks
biolrep =
c(rep(c(2,3,4,5,6,7,8),2),3,4,5,6,7,8,rep(c(2,3,4,5,6,7,8),5))
# I used the duplicateCorrelation() function to take in account the
link
between donors
corfit<-duplicateCorrelation(c031, design, ndups=1, block=biolrep)
Then I do the usual commands (incorporating the donor information) :
fit<-lmFit(c031, design, block=biolrep, cor=corfit$consensus)
cont.matrix<-makeContrasts(Stim2vsStim4=Stim2-Stim4, levels=design)
fit2<-contrasts.fit(fit, cont.matrix)
fit2<-eBayes(fit2)
anadiff.limma.Stim2<-topTable(fit2, number=1000, adjust="BH")
When I check the results and compare the list of the 1000 genes given
by
each method, there are 857 genes in common. The fact that I found the
use
of these two methods on this site doesn't helps me much about choosing
one.
Does someone got an idea about what method I should use for biological
replicates ? (And by the way if someone got an answer more precise
about
why Donor2 is excluded, that should be nice)
Best regards,
Santy
[[alternative HTML version deleted]]