Error when using gamma for delta prior with spiked-in data

[lilythepooh]

Dear BASiCS team,

Thanks a lot for reading this message.

I was trying to run BASiCS "WithSpikes=TRUE, Regression=FALSE, PriorDelta="gamma"" option with spiked-in data. My code is as the following:

sce <- SingleCellExperiment::SingleCellExperiment(
  assays = list(counts = Counts4[!Tech1,])
)
spike_sce<- SingleCellExperiment::SingleCellExperiment(
  assays = list(counts = Counts4[Tech1,])
)
altExp(sce, "spike-ins") <- spike_sce
rowData(altExp(sce, "spike-ins")) <- data.frame(Name = rownames(Counts4[Tech1,]), Molecules = SpikesInfo1$SpikeInput, row.names = rownames(Counts4[Tech1,]))
ZData<-sce
prior.param<-BASiCS_PriorParam(ZData,
                               PriorDelta = "gamma")

MCMC_Output<-BASiCS_MCMC(Data=ZData,N=17000,Thin=5,Burn=12000,
                         Regression=FALSE,WithSpikes = TRUE,
                         PriorParam = prior.param)

and I got the following error:

altExp 'spike-ins' is assumed to contain spike-in genes.
see help(altExp) for details. 

-------------------------------------------------------------
NOTE: default choice PriorDelta = 'log-normal'  (recommended value). 
Vallejos et al (2015) used a 'gamma' prior instead.
-------------------------------------------------------------

Error in .BASiCS_MCMC_NoSpikesParam(GPar$BioCounts, PriorParam$StochasticRef,  : 
  PriorDelta = 'gamma' is not supported for the no-spikes case

I am not sure why it is giving me a "no-spikes case" error when "WithSpikes=TRUE". The same data worked when "PriorDelta="log-normal", WithSpikes=TRUE". I also tried to use example data from data<-makeExampleBASiCS_Data(WithBatch = FALSE, WithSpikes = TRUE) , but it gave the same error when "PriorDelta="gamma"".

Could you please tell me how to fix this problem? Thanks very much!

Best regards, Sijia Li

[alanocallaghan]

Hi Sijia!

Is there a reason you want to use a gamma prior for delta in the no-spikes case? The gamma prior for delta is considered a legacy feature for compatibility reasons and so is not implemented for the no-spikes case, as the error message states. Further, it is not compatible at all for the regression case (where we use a joint prior on mu and delta). We consider a log-normal prior on delta to be a more natural framing. As you can see eq 5 in the genome biology paper, a log-normal prior induces a gaussian prior on log fold changes in over-dispersion when comparing between groups. https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0930-3

Cheers Al

[lilythepooh]

Hi Sijia!

Is there a reason you want to use a gamma prior for delta in the no-spikes case? The gamma prior for delta is considered a legacy feature for compatibility reasons and so is not implemented for the no-spikes case, as the error message states. Further, it is not compatible at all for the regression case (where we use a joint prior on mu and delta). We consider a log-normal prior on delta to be a more natural framing. As you can see eq 5 in the genome biology paper, a log-normal prior induces a gaussian prior on log fold changes in over-dispersion when comparing between groups. https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0930-3

Cheers Al

Dear Alan,

Thanks a lot for the quick reply.

Sorry if I was not clear in my last post. My question is that the gamma prior was not working even with Regression=FALSE, WithSpikes=TRUE. I have spiked-in genes in my data, but BASiCS_MCMC is giving me an error about no spikes data. for example:

data<-makeExampleBASiCS_Data(WithBatch=False,WithSpikes=TRUE)
prior.param<-BASiCS_PriorParam(data,PriorDelta="gamma")
chain<-BASiCS_MCMC(data,N=5000,Thin=5,Burn=2500, WithSpikes=TRUE,Regression=FALSE)

Here the data contains spiked-ins, and I am not using regression. But I got error

altExp 'spike-ins' is assumed to contain spike-in genes.
see help(altExp) for details. 

-------------------------------------------------------------
NOTE: default choice PriorDelta = 'log-normal'  (recommended value). 
Vallejos et al (2015) used a 'gamma' prior instead.
-------------------------------------------------------------

Error in .BASiCS_MCMC_NoSpikesParam(GPar$BioCounts, PriorParam$StochasticRef,  : 
  PriorDelta = 'gamma' is not supported for the no-spikes case

The session info is

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)

Thank you a lot!

Best regards, Sijia Li

[alanocallaghan]

My question is that the gamma prior was not working even with Regression=FALSE, WithSpikes=TRUE. I have spiked-in genes in my data, but BASiCS_MCMC is giving me an error about no spikes data.

Yes I understood. My comment on the regression case aside, the rest of the comment still applies. We have not implemented the ability to use a gamma prior for the no-spikes case. Can I ask why you would like to use it here? As stated in the NOTE, the lognormal prior is recommended.

[lilythepooh]

We have not implemented the ability to use a gamma prior for the no-spikes case.

Dear Alan,

Sorry if I was not clear. I am trying to use the gamma prior for a with-spikes case, not no-spikes case. Is the gamma prior feature disabled completely regardless of spikes or no-spikes? (It gives no-spikes error even when there are spikes.)

I was just curious about the different result from using different priors...

Thanks a lot for your time. Have a good evening!

Best regards, Sijia Li

[alanocallaghan]

Hi Sijia! Sorry I misunderstood and failed to parse your code properly. Totally my fault! That is just a bug probably introduced when we refactored the code slightly.

This issue is fixed in branch RELEASE_3_12. You can install it directly from github before we push it to Bioconductor if you like:

devtools::install_github("catavallejos/BASiCS", ref="RELEASE_3_12")

It's slightly embarrassing for me that you stated the problem plainly and politely twice and I insisted that you misunderstood. I hope you can find the funny side - let's blame it on tiredness this time :)

[lilythepooh]

Hi Sijia! Sorry I misunderstood and failed to parse your code properly. Totally my fault! That is just a bug probably introduced when we refactored the code slightly.

This issue is fixed in branch RELEASE_3_12. You can install it directly from github before we push it to Bioconductor if you like:

devtools::install_github("catavallejos/BASiCS", ref="RELEASE_3_12")

It's slightly embarrassing for me that you stated the problem plainly and politely twice and I insisted that you misunderstood. I hope you can find the funny side - let's blame it on tiredness this time :)

Hello Alan,

Thanks a lot for the timely update! I appreciate your patience and your time very much.

Have a wonderful week! :-)

Best wishes, Sijia