DRIMSeq problems with dmPrecision

Question

[poecile.pal]

Hello all,

I would like to use DRIMSeq (tutorial - https://bioconductor.org/packages/release/bioc/vignettes/DRIMSeq/inst/doc/DRIMSeq.pdf) In one moment I performed step dmPrecision:

> d <- dmPrecision(d, design = design_full)
! Using a subset of 0.1 genes to estimate common precision !

! Using common_precision = 196.046 as prec_init !

! Using loess fit as a shrinkage factor !

There were 50 or more warnings (use warnings() to see the first 50)

Warnings were like this:

50: In rowSums(lgamma(y + prec * prop) - lgamma(prop * prec)) :
  value out of range in 'lgamma'

I found that it is not the serious problem (https://groups.google.com/forum/#!topic/isoformswitchanalyzer/ihb4aK2fjG0), but I had some alarming messages:

> head(mean_expression(d), 3)
             gene_id mean_expression
1 ENSG00000000003.15       4494.7260
2  ENSG00000000005.6        247.1312
3 ENSG00000000419.12       1273.0802
Warning messages:
1: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
2: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
3: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
4: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
5: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
6: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
7: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
8: In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete
> common_precision(d)
[1] 196.046
Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :
  display list redraw incomplete

Could you please help me to deal this problem? By the way, I performed this operation with an unfiltered object. Also I am confused by "! Using loess fit as a shrinkage factor !". In tutorial I saw "! Using loess fit as a shrinkage factor !" . I can't find in the tutorial, when I should use this shrinkage factor and how. I knew about shrinkage before performing of DEG analysis, where I used

> res_tableOEg_ex38 <- DESeq2::results(ddsg_ex38, contrast=contrast_oeg_ex38, alpha = 0.05) 
> res_tableOEg_ex38 <- lfcShrink(ddsg_ex38, contrast=contrast_oeg_ex38, res=res_tableOEg_ex38)

with the table obut in this case I am afraid to do smth without instructions. Thank you!

Best regards, Poecile

Answer 1

[Kristoffer Vitting-Seerup]

This is a problem which occurs when there are to much many lowly expressed features. You can try filtering more strictly.

But I would not recommend using DRIMSeq in the first place. There are many other tools which performs better - try DEXSeq instead.

Also do you know about IsoformSwitchAnalyzeR - my R package which will help you make sense of your differential transcript usage analysis?

[poecile.pal]

Thank you a lot, I tried to repeat with the d3 with 34% of the initial number of genes (it was obtained after the strong filtering), and with d2 (25%), but had the same problem( Possibly, should I filter more and more strictly? I think, leaving 34% and 25% - it's very strict?

> d <- dmDSdata(counts=counts, samples=samps)
> d
An object of class dmDSdata 
with 38577 genes and 7 samples
* data accessors: counts(), samples()
  design()
  mean_expression(), common_precision(), genewise_precision()
  proportions(), coefficients()

> d2 <- dmFilter(d,
+ min_samps_feature_expr=n.small, min_feature_expr=10,
+ min_samps_feature_prop=n.small, min_feature_prop=0.1,
+ min_samps_gene_expr=n, min_gene_expr=10)
> d2
An object of class dmDSdata 
with 9587 genes and 7 samples
* data accessors: counts(), samples()

> d3 <- dmFilter(d,
+ min_samps_feature_expr=n.small, min_feature_expr=10,
+ min_samps_gene_expr=n, min_gene_expr=10)
> d3
An object of class dmDSdata 
with 13222 genes and 7 samples
* data accessors: counts(), samples()

Yes, I saw DEXSeq. Thank you very much, I am the beginner and I haven't know other programs, in addition to those discussed in details (with commands) in this article. I will try to use IsoformSwitchAnalyzeR too.