SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
cuffdiff : time series and cummeRbund analysis apratap Bioinformatics 1 04-28-2014 03:31 AM
can DESeq handle time series data now? Jeremy RNA Sequencing 0 12-17-2012 09:15 PM
How to analyse time-series data with edgeR ca85 RNA Sequencing 0 10-06-2011 05:04 AM
time series analysis of count data? greigite Bioinformatics 0 06-14-2011 05:11 PM

Reply
 
Thread Tools
Old 07-22-2013, 10:05 PM   #1
ashuchawla
Member
 
Location: san diego

Join Date: Jan 2012
Posts: 38
Default Comparative Time Series Analysis of RNA-Seq data?

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
ashuchawla is offline   Reply With Quote
Old 07-23-2013, 12:07 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

DESeq(2), edgeR, limma/voom, etc. Anything that can handle GLMs will work.
dpryan is offline   Reply With Quote
Old 07-23-2013, 06:39 AM   #3
ashuchawla
Member
 
Location: san diego

Join Date: Jan 2012
Posts: 38
Default

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.

Quote:
Originally Posted by dpryan View Post
DESeq(2), edgeR, limma/voom, etc. Anything that can handle GLMs will work.
ashuchawla is offline   Reply With Quote
Old 07-23-2013, 06:58 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

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 is offline   Reply With Quote
Old 07-23-2013, 10:08 AM   #5
ashuchawla
Member
 
Location: san diego

Join Date: Jan 2012
Posts: 38
Default

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??
*************************************

Quote:
Originally Posted by dpryan View Post
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).

Last edited by ashuchawla; 07-23-2013 at 10:13 AM.
ashuchawla is offline   Reply With Quote
Old 07-23-2013, 12:55 PM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
# 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.
dpryan is offline   Reply With Quote
Old 07-23-2013, 01:11 PM   #7
ashuchawla
Member
 
Location: san diego

Join Date: Jan 2012
Posts: 38
Default

> 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.



Quote:
Originally Posted by dpryan View Post
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 is offline   Reply With Quote
Old 07-23-2013, 01:39 PM   #8
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Ah, right, you might try pooled-CR.
dpryan is offline   Reply With Quote
Old 07-25-2013, 12:12 PM   #9
ashuchawla
Member
 
Location: san diego

Join Date: Jan 2012
Posts: 38
Default

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.

Quote:
Originally Posted by dpryan View Post
Ah, right, you might try pooled-CR.
ashuchawla is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:31 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO