Limma with replicates at different time points
1
0
Entering edit mode
alessandro • 0
@0010fcc2
Last seen 6 weeks ago
Italy

Hello,

I am performing an RNAseq analysis using the limma package.
From the same cell culture 3 samples(lane) were extracted, these samples were analysed at 4 different time points(Time).
No treatment was applied to any of the samples.
The goal of the project is to compare the 4 time points to each other, initially with a DEG.
The result of my analysis produces far too many differentially expressed genes with very low adj_p_values (even more than 100 genes with an adj_p_value < 0.0001).
I am also having difficulty understanding how I should consider the replicates, whether as technical or biological replicates.
Am I doing something wrong in setting up my model?

This is an example of my data:

Samples Lane Time_Point
S1 A T1
S2 B T1
S3 C T1
S4 A T2
S5 B T2
S6 C T2
S7 A T3
S8 B T3
S9 C T3
S10 A T4
S11 B T4
S12 C T4

This is my code:

dge = DGEList(counts = exp.f, genes = rownames(exp.f))

### Preprocessing
dge$counts <- edgeR::cpm(dge, log = F)
dge = calcNormFactors(dge, method = "TMM") 

### Define model
f = as.factor(dataset$Time)
l = as.factor(dataset$Lane)
design = model.matrix(~0+f)
colnames(design) <- c(levels(f))

v = voomWithQualityWeights(dge, design = design, plot = TRUE)
corfit <- duplicateCorrelation(v, design, block = l, ndups = 1)
v = voomWithQualityWeights(dge, design = design, plot = TRUE, correlation = corfit$consensus.correlation)
corfit <- duplicateCorrelation(v, design, block = l, ndups = 1)

### Fit model
fit = lmFit(v, design, block = l, correlation = corfit$consensus.correlation)
contrast.matrix = makeContrasts(T4-T3,T3-T2,T2-T1,T4-T1, levels = design)
fit2 = contrasts.fit(fit, contrast.matrix)
fit2 = eBayes(fit2)

Thanks,

Alessandro

RNAs limma Differen • 297 views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 7 hours ago
WEHI, Melbourne, Australia

I don't entirely understand your experiment, but the analysis you have done looks reasonable.

If the difference between the lanes is large, you could alternatively consider using design <- model.matrix(~0+f+l) instead of the duplicateCorrelation approach.

An equivalent but more convenient way to do the analysis you've done would be by:

fit <- voomLmFit(dge, design, block=l, sample.weights=TRUE)

This would replace 5 lines of code with one.

ADD COMMENT
0
Entering edit mode

Hi gordon,

Sorry for only replying now.
Thank you very much for your reply, I have tried to do as you suggested.
If I can give you more information about what you are not clear about the experiment I can rewrite it to better understand if my method is correct.

Thank you very much.

ADD REPLY

Login before adding your answer.

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