I'm trying to study gene expression across the life cycle of a nematode using DESeq2. We've sequenced eggs, L1, L2, L3, L4, L5, and adults. I've figured out how to do stage v. stage comparisons (e.g., egg v. L1, L1 v. L2, L2 v. L3, etc). Now, I'd like to generate a list of genes whose expression fluctuates at some point in the life cycle -- between any of the stages. I'd appreciate some advice on how to tease this out.
Seqanswers Leaderboard Ad
Collapse
Announcement
Collapse
No announcement yet.
X
-
hi,
I think what you want to use is the likelihood ratio test. You can either do this with DESeq():
dds <- DESeq(dds, test="LRT", full=~stage, reduced=~1)
or you can also use the nbinomLRT() function:
dds <- nbinomLRT(dds, full=~stage, reduced=~1)
The likelihood ratio test compares the difference in log likelihood using the stage information vs removing the stage information altogether. This then gives you a significance test for whether the gene expression changes with stage, at any point in the life cycle like you say.
-
Thanks so much for your reply!
Just to be 100% sure...
I run the DESeq function just as you demonstrated, then I grab the results:
res <- results(dds)
Any gene with an adjusted p-value less than 1e-5 (or whatever confidence cutoff we choose) will be differentially expressed at some stage, whereas the others show stable expression levels over the life cycle. Is this correct?
Comment
-
yes. And with the LRT, note that the log2FoldChange and lfcSE columns in the results object give the log2 fold change and standard error of the last level over the first, which is not related to the LRT statistic, p-value or adjusted p-value in the other columns of the results object.
Also, that is a very low adjusted p-value. We usually use a false discovery rate of 10% for examples. The optimal FDR to aim for depends on the cost of follow-up experiments and the cost of missing a true discovery, of course, but 1/10 or 1/20 is often reasonable. For exploratory experiments 1/5 might make sense.
If you meant 1e-5 for a p-value threshold, we recommend to focus on adjusted p-values for building lists, as these are much more interpretable.
Comment
-
Michael, I was hoping you might be able to answer one more related question...
I've decided that I need to control for batch variation since some of the parasites I'm studying were derived from different host animals at different times of year, etc. I still want to generate a list of genes whose expression varies from stage to stage -- I just want to disregard variation due to other factors as best I can.
Would changing this:
dds <- DESeq(dds, test="LRT", full=~stage, reduced=~1)
to this:
dds <- DESeq(dds, test="LRT", full=~Stage+Date, reduced=~Date)
accomplish my goal? I wanted to double check that I'm not making the reverse comparison and accidentally measuring variation due to collection date.
Also, would you expect to see more DE genes controlling for batch variation or less? I seem to be finding more, and that surprised me (which is what inspired this post -- I'm nearly certain that I'm doing something wrong). I figured that some of the variation I was seeing before was likely due to differences in season or host animal and that pulling those factors out would reduce my DE gene count.
Comment
-
hi,
yes i would just change this to:
dds <- DESeq(dds, test="LRT", full=~Date + Stage, reduced=~Date)
in general always put the variable of interest at the end. while this won't make a difference for the LRT statistic, or p-values, the results() also prints out a log2 fold change based on the last variable in the design formula, and comparing the last level of this variable to the base level.
this is admittedly a bit confusing when you have more than 2 levels for Stage, as all of the levels are taken into account when calculating the LRT statistic and p-values, but we wanted to keep the results() object standardized, and adding every possible LFC contrast would get messy.
you can extract other LFCs, if you are interested, like so:
lfcBvsA <- results(dds, contrast=c("Stage","B","A")$log2FoldChange
lfcCvsA <- results(dds, contrast=c("Stage","C","A")$log2FoldChange
lfcCvsB <- results(dds, contrast=c("Stage","C","B")$log2FoldChange
etc
regarding your second question, it can go either way, but often accounting for batch will end up with more DE genes. If the batches are balanced, the unaccounted-for batch effect contributes more to the denominator of the Wald statistic than the top. And similarly in the case of the LRT, unaccounted-for batch effect is likely to reduce the LRT and raise p-values. Not accounting for batch effect increases the estimates of dispersion at the gene-wise level, and additionally increases the fitted values for dispersion over all genes.
for example, with normal linear models:
> condition <- factor(rep(1:2,each=10))
> batch <- factor(rep(rep(1:2,each=5),2))
> condition
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
Levels: 1 2
> batch
[1] 1 1 1 1 1 2 2 2 2 2 1 1 1 1 1 2 2 2 2 2
Levels: 1 2
> y <- rnorm(20,mean=as.numeric(condition) + 3*as.numeric(batch))
> y
[1] 3.563869 3.065336 4.289239 1.650792 4.086886 6.564326 4.452499
[8] 7.612532 7.720049 6.091950 6.644862 5.132749 5.328683 5.825075
[15] 4.834425 6.843204 6.743685 9.943101 8.852742 10.021774
> coef(summary(lm(y ~ condition)))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.909748 0.6246445 7.860068 3.148028e-07
condition2 2.107282 0.8833807 2.385475 2.825605e-02
> coef(summary(lm(y ~ batch + condition)))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.388551 0.4597157 7.370970 1.092649e-06
batch2 3.042395 0.5308340 5.731349 2.446249e-05
condition2 2.107282 0.5308340 3.969757 9.900114e-04
Note the jump in t-statistic for condition from 2.38 to 3.96. This is because the unexplained variance was reduced by adding the batch variable. The estimate stayed at 2.10, while the SE of the estimate went down from 0.88 to 0.53.
Comment
Latest Articles
Collapse
-
by seqadmin
Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...-
Channel: Articles
04-04-2024, 04:25 PM -
-
by seqadmin
Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...-
Channel: Articles
03-22-2024, 06:39 AM -
ad_right_rmr
Collapse
News
Collapse
Topics | Statistics | Last Post | ||
---|---|---|---|---|
Started by seqadmin, 04-11-2024, 12:08 PM
|
0 responses
31 views
0 likes
|
Last Post
by seqadmin
04-11-2024, 12:08 PM
|
||
Started by seqadmin, 04-10-2024, 10:19 PM
|
0 responses
33 views
0 likes
|
Last Post
by seqadmin
04-10-2024, 10:19 PM
|
||
Started by seqadmin, 04-10-2024, 09:21 AM
|
0 responses
28 views
0 likes
|
Last Post
by seqadmin
04-10-2024, 09:21 AM
|
||
Started by seqadmin, 04-04-2024, 09:00 AM
|
0 responses
53 views
0 likes
|
Last Post
by seqadmin
04-04-2024, 09:00 AM
|
Comment