Question: DESEq2 warning while using GLM fitting
2
gravatar for bioinforupesh2009.au
4.4 years ago by
Spain
bioinforupesh2009.au100 wrote:

Hello all,

Using time series (two condition with 4 time point) data in DESeq2, i got warning while fitting by GLM ??

bckCDS <- DESeqDataSetFromMatrix(countData = x, colData=ExpDesign, design= ~condition.label + condition.Time + condition.label:condition.Time)

> ddsLRT= DESeq (bckCDS, test ="LRT", reduced = ~condition.label + condition.Time)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
NOTE: fitType='parametric', but the dispersion trend was not well captured by the
  function: y = a/x + b, and a local regression fit was automatically substituted.
  specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing

 

May i have fit it by "local", indeed its work by this however, my question is why its not working by "fitType='parametric'". What is advantage and disadvantage of it.

Thanks

rna-seq deseq2 next-gen R • 4.3k views
ADD COMMENTlink modified 4.4 years ago by Michael Love1.8k • written 4.4 years ago by bioinforupesh2009.au100
2
gravatar for Michael Love
4.4 years ago by
Michael Love1.8k
United States
Michael Love1.8k wrote:

Depending on the version you are using, this used to be a warning but is now just a message/note.

The printed message here is literally all you need to know. The software cannot find parameters a and b to fit the dispersion over mean trend using that formula (read the manuscript for more information about the point of modeling dispersion over mean).

So a local regression was automatically substituted. There is no further action required by the user.

ADD COMMENTlink written 4.4 years ago by Michael Love1.8k

Dear sir,
Thanks for your kind answer. Its done. But now i am trying to use Time series analysis for one condition (one condition with 4 time points).... Indeed i tried and got error.

bckCDS <- DESeqDataSetFromMatrix(countData = wt, colData=ExpDesign, design= ~condition.Time)
>head(ExpDesign)
    condition.label condition.Time
WT.T1.1              WT             T1
WT.T1.2              WT             T1
WT.T1.3              WT             T1
WT.T2.1              WT             T2
....
.......

> ddsLRT= DESeq(bckCDS, test ="LRT", reduced = ~condition.Time, fitType='local')
Error in nbinomLRT(object, full = full, reduced = reduced, betaPrior = betaPrior,  :
  less than one degree of freedom, perhaps full and reduced models are not in the correct order

> ?DESeq
> ddsLRT= DESeq(bckCDS, test ="LRT", fitType='local')
Error in checkLRT(full, reduced) :
  provide a reduced formula for the likelihood ratio test, e.g. nbinomLRT(object, reduced = ~ 1)

ADD REPLYlink written 4.4 years ago by bioinforupesh2009.au100

Now i tried.....

ddsLRT= DESeq(bckCDS, test ="LRT", fitType='local',full=~condition.Time, reduced=~1)

and it seem its work, is it okey what i am doing ??

ADD REPLYlink written 4.4 years ago by bioinforupesh2009.au100

Yes, that's likely what you want. What you're asking there is, "Which model fits a gene's expression better, one where expression can change with each timepoint (~condition.Time), or one where expression is constant or at least doesn't vary at the timepoints we looked at (~1)."

If you know someone locally who's taken more statistics (or even just linear algebra...so any physicicist, computational biologist, engineer, etc.), you might ask him/her to draw out what ~condition.Time and ~1 are doing. I know that these things largely look meaningless, but they're probably pretty easy to understand if someone just goes through a simple example on a white board.

ADD REPLYlink modified 4.4 years ago • written 4.4 years ago by Devon Ryan89k

You used the same model for the full and reduced fitting...that won't work very well. I'm assuming with the LRT you're trying to ask the question, "Is there an effect of time?" In that case, the reduced design should be ~1.

ADD REPLYlink written 4.4 years ago by Devon Ryan89k

Thanks Devon for your kind reply,
that mean i should not use full = ~condition, right ?
now i end up with..

bckCDS <- DESeqDataSetFromMatrix(countData = wt, colData=ExpDesign, design= ~condition.Time)

ddsLRT= DESeq(bckCDS, test ="LRT", fitType='local', reduced=~1)

resLRT=results(ddsLRT)

ADD REPLYlink written 4.4 years ago by bioinforupesh2009.au100

You didn't define a condition column in ExpDesign, so using ~condition would just lead to an error message.

ADD REPLYlink modified 4.4 years ago • written 4.4 years ago by Devon Ryan89k

I didnt get any error bcz i have condition as name condition.label in ExpDesign and its only one condition, therefore i can't use contrast model, however i used only condition.Time in my model matrix.

bckCDS <- DESeqDataSetFromMatrix(countData = wt, colData=ExpDesign, design= ~condition.Time)

ddsLRT= DESeq(bckCDS, test ="LRT", fitType='local', reduced=~1)


>head(ExpDesign)
    condition.label condition.Time
WT.T1.1              WT             T1
WT.T1.2              WT             T1
WT.T1.3              WT             T1
WT.T2.1              WT             T2
....
.......

ADD REPLYlink written 4.4 years ago by bioinforupesh2009.au100

If condition.label is only ever WT then there's no point in including it. It can't change anything.

ADD REPLYlink written 4.4 years ago by Devon Ryan89k

Thanks Dev,
I got my results, just last answer i would like to have from your side. I got comparison between
T1 vs T2, T1 vs T3 and T1 vs T4. however, is there way to compare all 6 combinations. like
T1 vs T2, T1 vs T3, T1 vs T4, T2 vs T3, T2 vs T4 and T3 vs T4 ??? i guess, i have to use nbinomLRT () ??

thanks, i got it.

 

ADD REPLYlink modified 4.4 years ago • written 4.4 years ago by bioinforupesh2009.au100

See help(results) and note the contrast parameter. Also, you can just use a Wald test, which is the default.

ADD REPLYlink written 4.4 years ago by Devon Ryan89k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 936 users visited in the last hour