gosianow / DRIMSeq_before_BioC3.6

Package for differential analyses using the Dirichlet Multinomial distribution
3 stars 0 forks source link

warnings produced by dmPrecision() #4

Open kvittingseerup opened 7 years ago

kvittingseerup commented 7 years ago

I ran in the a problem where the dmPrecision() function produces a lot of warnings that looks like a problem. The issue can be reproduced by the following code example:

rm(list=ls())
library('DRIMSeq')
library('PasillaTranscriptExpr')

### Lets used the pasilla data
data_dir  <- system.file("extdata", package = "PasillaTranscriptExpr")
pasilla_metadata <- read.table(file.path(data_dir, "metadata.txt"), header = TRUE, as.is = TRUE)
pasilla_counts <- read.table(file.path(data_dir, "counts.txt"), header = TRUE, as.is = TRUE)
gene_id_subset <- readLines(file.path(data_dir, "gene_id_subset.txt"))

### pasilla_samples design matrix with 3 groups
pasilla_samples <- data.frame(sample_id = pasilla_metadata$SampleName, group = pasilla_metadata$condition)
# add a bath effect
pasilla_samples$batch <- c(1,1,2,1,2,2,2)

### Create and subset dmData
d <- dmDSdata(counts = pasilla_counts, samples = pasilla_samples)
d <- d[names(d) %in% gene_id_subset, ]
### Filtering
d <- dmFilter(d, min_samps_gene_expr = 7, min_samps_feature_expr = 3, min_gene_expr = 10, min_feature_expr = 10)

## Create the design matrix without batch effect
designWithoutBatch <- model.matrix(~ group        , data = pasilla_samples)
designWithBatch    <- model.matrix(~ group + batch, data = pasilla_samples)

designWithoutBatch
designWithBatch

### Calculate precision
dWithoutBatch <- dmPrecision(d, design = designWithoutBatch)
dWithBatch    <- dmPrecision(d, design = designWithBatch   )

The dmPrecision() works fine when I analyze the one design without the batch effect but when the batch effect is added I get the the following warning:

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

! Using common_precision = 21.2862 as prec_init !

! Using 0 as a shrinkage factor !

There were 50 or more warnings (use warnings() to see the first 50)
> head(warnings())
Warning messages:
1: In rowSums(lgamma(y + prec * prop) - lgamma(prop * prec)) :
  value out of range in 'lgamma'
2: In rowSums(lgamma(y + prec * prop) - lgamma(prop * prec)) :
  value out of range in 'lgamma'
3: In rowSums(lgamma(y + prec * prop) - lgamma(prop * prec)) :
  value out of range in 'lgamma'
4: In rowSums(lgamma(y + prec * prop) - lgamma(prop * prec)) :
  value out of range in 'lgamma'
5: In rowSums(lgamma(y + prec * prop) - lgamma(prop * prec)) :
  value out of range in 'lgamma'
6: In rowSums(lgamma(y + prec * prop) - lgamma(prop * prec)) :
  value out of range in 'lgamma'

And I get the same warning if I just ignore them and continue with the dmFit and dmTest.

Is this a serious problem?

/Kristoffer

gosianow commented 7 years ago

Hello,

no, that is not a serious problem. This usually happens in the Hessian calculation when the regression framework is used. First, I try to calculate the exact Hessian but when this fails then the approximation from the optim package is used. The warnings occur during the exact calculation. For the next version I will try to prevent to display those warnings.

kvittingseerup commented 7 years ago

Thanks :-)