What Limma model to use when corfit$consensus is negative or nearly zero
1
0
Entering edit mode
ialwkcrfm • 0
@ialwkcrfm-13801
Last seen 7.3 years ago

I am analyzing a microarray experiment. The design of my study is as follows:

Each mouse has 3 measurements (baseline, 14 days, 30 days) There are NO replicates (technical or biological). I have 28 samples in total that corresponds to 10 mice and there are 3 treatments. So per treatment I have 3 mice.

 

Group   Treatment   MiceID  Status
GR1 TrtA    1   baseline
GR1 TrtA    1   Day14
GR1 TrtA    1   Day30
GR3 Combo   2   baseline
GR3 Combo   2   Day14
GR3 Combo   2   Day30
GR3 Combo   3   baseline
GR3 Combo   3   Day30
GR2 TrtB    4   baseline
GR2 TrtB    4   Day14
GR2 TrtB    4   Day30
GR2 TrtB    5   baseline
GR2 TrtB    5   Day14
GR2 TrtB    6   baseline
GR2 TrtB    6   Day14
GR2 TrtB    6   Day30
GR1 TrtA    7   baseline
GR1 TrtA    7   Day14
GR1 TrtA    7   Day30
GR3 Combo   8   baseline
GR3 Combo   8   Day14
GR3 Combo   8   Day30
GR3 Combo   9   baseline
GR3 Combo   9   Day14
GR3 Combo   9   Day30
GR1 TrtA    10  baseline
GR1 TrtA    10  Day14
GR1 TrtA    10  Day30

Goal: The question is what are the differentially expressed genes for day 14 vs baseline and day30 vs day14 per treatment?

To account for the correlation within each mice, I used the following Limma model and my corfit$consensus = -0.0068

1) Given this is a negative and nearly zero corfit correlation, should I be including correlation in my lmFit model? If not, then how do I account for the correlation per mice?

2) Does it make sense to run two Limma, one comparing baseline vs day14 and another comparing day30 vs day 14 samples

LIMMA R Code: annot : sampleInfo shown above, data.log= log2(gene expression signals)

Treat <- factor(paste(annot$Treatment,annot$Status,sep=".")) 
design <- model.matrix(~0+Treat) 
colnames(design) <- c("TrtA.baseline", "TrtA.Day30", "TrtA.Day14", "Combo.baseline", "Combo.Day30", "Combo.Day14", "TrtB.baseline", "TrtB.Day30", "TrtB.Day14") 
corfit <- duplicateCorrelation(data.log, design, block=annot$MiceID)corfit$consensus  
fit1 <- lmFit(data.log,design,block=annot$miceID,correlation=corfit$consensus)
Differences <- makeContrasts(TrtALatevsTrtAwk2= TrtA.Day30-TrtA.Day14, ComboLatevsCombowk2=Combo.Day30-Combo.Day14, TrtBLate=TrtB.Day30-TrtB.Day14, TrtAwk2vsTrtA0= TrtA.Day14-TrtA.baseline, Combowk2vsCombo.0=Combo.Day14-Combo.baseline, TrtB=TrtB.Day14-TrtB.baseline,levels=design) 
fit2 <- contrasts.fit(fit1,Differences) fit <- eBayes(fit2)

 

 

Limma corfit • 1.7k views
ADD COMMENT
1
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 2 hours ago
The city by the bay

For your first question: limma and duplicateCorrelation have no problems with a consensus correlation of zero. Just proceed with the DE analysis like you would do for any other value of the correlation.

For your second question: use all the samples when running lmFit and eBayes, as this provides more information for variance estimation. You can use makeContrasts to specify particular contrasts.

More generally, it seems you want to compare between time points within each treatment. In that case, you can use an alternative approach whereby you just block on MiceID in the design matrix, i.e.,

design <- model.matrix(~factor(annot$MiceID)+Treat)
design <- design[,!grepl("baseline", colnames(design))] # get to full rank

... which avoids the anticonservativeness associated with duplicateCorrelation.

ADD COMMENT
0
Entering edit mode

Thanks Aaron. Let me run this model and update you on the outcome.

ADD REPLY

Login before adding your answer.

Traffic: 796 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