Hello,

Thank you very much for your package.
I have questions about model and K value settings in zinbwave function when combining with DESeq2.

I’m working on SMART-seqHT dataset of a specific cell type, human sample, to compare gene expressions under two conditions(positive vs negative).
The samples are from 3 patients.
My interest is in the DEGs between positive and negative cells, and the normalized values for heatmap visualization.
With Seurat package, I failed to find enough DEGs for further functional analysis, so I tried the zinbwave and DESeq2, in which I found few hundred DEGs with the Wald test.

These are the tutorials I followed:
https://github.com/mikelove/zinbwave-deseq2/blob/master/zinbwave-deseq2.knit.md
https://bioconductor.org/packages/release/bioc/vignettes/zinbwave/inst/doc/intro.html
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

My questions are:

1. About the model:

My current code is:

zinb <- zinbwave(zinb, K=0, BPPARAM=SerialParam(), epsilon=1e12, normalizedValues=TRUE, observationalWeights = TRUE)
dds <- DESeqDataSet(zinb, design= ~ condition + patient)

But should I include the model in the zinbwave function too?:

zinb <- zinbwave(zinb, K=0, X=“~ condition + patient”, BPPARAM=SerialParam(), epsilon=1e12, normalizedValues=TRUE, observationalWeights = TRUE)
dds <- DESeqDataSet(zinb, design= ~ condition + patient)

The second codes gives me about 40 less DEGs than the first codes.
According to issue #36( #36) , should I use the second codes ?
And where does this different number of DEGs come from?

2. Is the K value setting (“K=0“) appropriate, for both DEGs and normalized values?

Here is my current code:

sce <- SingleCellExperiment(assays = list(counts = as.matrix(count_matrix)), colData = sample)

filter <- rowSums(assay(zinb)>5)>5
zinb <- sce[filter,]

zinb <- zinbwave(zinb, K=0, BPPARAM=SerialParam(), epsilon=1e12, normalizedValues=TRUE, observationalWeights = TRUE)

write.table(assay(zinb,"normalizedValues"), file="normalizedValues.txt",sep="\t")

dds <- DESeqDataSet(zinb, design= ~ condition + patient)
dds <- estimateSizeFactors(dds, type="poscounts")

scr <- scran::calculateSumFactors(dds)
sizeFactors(dds) <- scr

dds <- DESeq(dds, sfType="poscounts", useT=TRUE, minmu=1e-6, minRep=Inf) 
res <- lfcShrink(dds, contrast=c(“condition”, "positive, "negative”), type = "normal")

I apologize for my poor understanding and duplicated question.
Thank you very much in advance.