Comparative Time Series Analysis of RNA-Seq data?

[ashuchawla]

I need help with comparative time-series analysis of RNA-Seq Data. I have an untreated sample, then 3 treatments at day 3, the same 3 treatments at day 6 and then at day 9. I have tag-counts(generated from HT-Seq counts) for all of these samples across time. Could you please suggest some tools to get Differentially expressed genes. I do not wish to perform pair-wise comparisons.

Thank you so much

Ashu

[dpryan]

DESeq(2), edgeR, limma/voom, etc. Anything that can handle GLMs will work.

[ashuchawla]

Great, than you so much.

Could you or anybody explain this to me with an example? I am new to this and its tough to understand what functions/methods to use.

DESeq(2), edgeR, limma/voom, etc. Anything that can handle GLMs will work.

[dpryan]

I expect there's an example of this in one of the vignettes/user guides for a couple of the packages I listed. In general, you just want to fit your data with a GLM of the form "counts ~ time*treatment", assuming that there can be a meaningful time:treatment interaction in your experiment (otherwise, just exchange a plus sign for the asterisk).

[ashuchawla]

The pdf(http://bioconductor.org/packages/rel.../doc/DESeq.pdf) from bioconductor does this ->
> fit1 = fitNbinomGLMs( cdsFull, count ~ libType + condition )
> fit0 = fitNbinomGLMs( cdsFull, count ~ libType )
> pvalsGLM = nbinomGLMTest( fit1, fit0 )
> padjGLM = p.adjust( pvalsGLM, method="BH" )

I dont understand what are they trying to do here. Below is what I tried with my data

>glm_design

samples libType time treatment

Samp_1 single-end d0 Untreated
Samp_2 single-end d3 treatment1
Samp_3 single-end d6 treatment1
Samp_4 single-end d9 treatment1
Samp_5 single-end d3 treatment2
Samp_6 single-end d6 treatment2
Samp_7 single-end d9 treatment2
Samp_8 single-end d3 treatment3
Samp_9 single-end d6 treatment3
Samp_10 single-end d9 treatment3

glm_cds = newCountDataSet(combined_file , glm_design)

glm_cds = estimateSizeFactors(glm_cds)

# I dont have replicates
glm_cds = estimateDispersions(glm_cds, method="blind", sharingMode="fit-only" )

fit1 = fitNbinomGLMs( glm_cds, count ~ time*treatment)

fit0 = fitNbinomGLMs( glm_cds, count ~ time)

pvalsGLM = nbinomGLMTest( fit1, fit0 )

padjGLM = p.adjust( pvalsGLM, method="BH" )


IS THIS WHAT I AM SUPPOSED TO DO?

unique(padjGLM)
[1] 1 NA

padj values do not seem right??
*************************************

I expect there's an example of this in one of the vignettes/user guides for a couple of the packages I listed. In general, you just want to fit your data with a GLM of the form "counts ~ time*treatment", assuming that there can be a meaningful time:treatment interaction in your experiment (otherwise, just exchange a plus sign for the asterisk).

[dpryan]

# I dont have replicates
glm_cds = estimateDispersions(glm_cds, method="blind", sharingMode="fit-only" )

fit1 = fitNbinomGLMs( glm_cds, count ~ time*treatment)
fit0 = fitNbinomGLMs( glm_cds, count ~ time)

Well, you sort of have replicates, just not full replicates. You might try the default dispersion method and see how well that fits. For fit0, try "fit- <- fitNbinomGLMs(glm_cds, count~time+time:treatment)". Otherwise, you getting output for things that can vary by both treatment and a time:treatment interaction, which you indicated not wanting.

BTW, an adjusted p-value of 1 is normal and to be expected. That all of them are 1 is too bad, though doing as I suggested above my change that. You can also try doing some independent filtering (see the genefilter package and the accompanying paper in PNAS), which might improve things further.

[ashuchawla]

> glm_cds = estimateDispersions(glm_cds)
Error in .local(object, ...) :
None of your conditions is replicated. Use method='blind' to estimate across conditions, or 'pooled-CR', if you have crossed factors.

Well, you sort of have replicates, just not full replicates. You might try the default dispersion method and see how well that fits. For fit0, try "fit- <- fitNbinomGLMs(glm_cds, count~time+time:treatment)". Otherwise, you getting output for things that can vary by both treatment and a time:treatment interaction, which you indicated not wanting.

BTW, an adjusted p-value of 1 is normal and to be expected. That all of them are 1 is too bad, though doing as I suggested above my change that. You can also try doing some independent filtering (see the genefilter package and the accompanying paper in PNAS), which might improve things further.

[dpryan]

Ah, right, you might try pooled-CR.

[ashuchawla]

Thank you so much for your help.

This is what it tells me when I try pooled-CR

glm_cds = estimateDispersions(glm_cds, method="pooled-CR", sharingMode="fit-only" , modelFormula=count~time+time:treatment)
Error in FUN(newX[, i], ...) :
No residual degrees of freedom. Most likely the design is lacking sufficient replication.

Ah, right, you might try pooled-CR.