time course experiment with many time points in limmc for Rna-seq data
2
0
Entering edit mode
@joanna-zyprych-3333
Last seen 8.0 years ago

Dear all,

I have a question concerning the RNA-seq experiment. 

My file targets is as follows:

   Filename   Group Time
1   A2780_1 Control    0
2   A2780_2 Control    0
3    4PTX_1   Treat    1
4    4PTX_2   Treat    1
5    8PTX_1   Treat    2
6    8PTX_2   Treat    2
7   16PTX_1   Treat    3
8   16PTX_2   Treat    3
9   32PTX_1   Treat    4
10  32PTX_2   Treat    4
11  64PTX_1   Treat    5
12  64PTX_2   Treat    5
13 128PTX_1   Treat    6
14 128PTX_2   Treat    6

and I'm trying to find genes that are significant across these time points and groups. Then I'm using code from limma guide for time course experiment with many time points and at the end I'm getting some warnings. Could someone help me what is wrong in my code...
> X <- ns(targets$Time, df=5)
> Group <- factor(targets$Group)
> design <- model.matrix(~Group*X)
> x <- DGEList(counts=counts$counts[,c(9,41,13:17,45:49,65,66)],
+              genes=counts$annotation[,c("GeneID","Length")])
> isexpr <- rowSums(cpm(x) > 5) >= 2
> x <- x[isexpr,]
> y <- voom(x, design)
Coefficients not estimable: GroupTreat:X1 GroupTreat:X2 GroupTreat:X3 GroupTreat:X4 GroupTreat:X5
Partial NA coefficients for 11610 probe(s)
> fit <- lmFit(y, design)

Coefficients not estimable: GroupTreat:X1 GroupTreat:X2 GroupTreat:X3 GroupTreat:X4 GroupTreat:X5
Partial NA coefficients for 11610 probe(s)

 fit <- eBayes(fit)
  Estimation of var.prior failed - set to default value
> topTableF(fit, adjust="BH")
F-statistics not found in fit

I will be grateful for any help.

Best regards, Joanna

 

 

 

limma • 1.1k views
ADD COMMENT
0
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 7 hours ago
The city by the bay

Remove the Group factor in your design matrix, because it's confounded with Time when you set it up as an interaction model. That is to say, the only control sample you have is that at the zero time point, while all other time points are treated. This means that each of your interaction terms is identical to one of the spline coefficients, such that they're both unestimable. If, for example, you had DE between A2780 and 4PTX, you wouldn't be able to tell whether that was due to the time effect or due to the time:treatment interaction effect.

I suspect your "control" sample is really a "treatment at time zero" sample, in which case you don't need the Group factor. Just leave it out and test for the effect of time by dropping all the spline coefficients. However, if that's not the case, then you've designed your experiment poorly; you'd need control samples at several time points to do a proper analysis, e.g., to compare the time response between control and treatment.

I suppose a compromise would be to set it up like this:

design <- model.matrix(~Group + X)

This would allow you to test for the effect of time within the treated group, and it would allow you to test for DE between control and treatment at time zero only. I suspect the latter won't get a lot of DE, though, due to the confounding issue; you wouldn't be able to tell whether a difference between two groups was due to the treatment effect or due to the time effect.

ADD COMMENT
0
Entering edit mode
@joanna-zyprych-3333
Last seen 8.0 years ago

Many thanks for your advise! I will try to do this in this way. Best wishes, Joanna

 

 

 

ADD COMMENT

Login before adding your answer.

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