Negative s-values

[rfarouni]

What causes the s-values to be negative? In the plot below I show s-values that are close to zero vs logFold change. Is that an identifiability issue?

Thanks in advance!

[stephens999]

this looks like numerical error due to finite precision. They are all basically 0. We should probably do something to have these returned as 0 rather than negative, but you can treat them as 0.

[rfarouni]

Thanks for the getting back to me Matthew. By the way, any downstream filtering using transformed log values will miss these hits -which happened to me initially before I examined the results in more details.

[stephens999]

hmmm - that could be true for 0s too, right?

[stephens999]

this line was intended to ensure all lfsr values are non-negative https://github.com/stephens999/ashr/blob/ce75b471379b9d7f84c2346880311c08c6e518aa/R/ash.R#L524

And this should in principle ensure the s values are non-negative (they are computed using cumsum of sorted lfsr)

@rfarouni Are your lfsr values also negative? Can you share a reproducible code example so we can see why this isn't working?

[rfarouni]

I am using the DESeq2 function lfcShrink(dds, contrast = contrast, type = "ashr", svalue = T ) which does not return the lfsr values. As for the 0s, I never got exact zeros in the s-values so they always end up log transformed unless they are negative.

[stephens999]

oh, then can you check in DESeq2 how they compute the s values? or do they return the svalues direct from ash?

[rfarouni]

This is the relevant code in lfcShrink

        if (!requireNamespace("ashr", quietly = TRUE)) {
            stop("type='ashr' requires installing the CRAN package 'ashr'")
        }
        if (!quiet) 
            message("using 'ashr' for LFC shrinkage. If used in published research, please cite:\n    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.\n    https://doi.org/10.1093/biostatistics/kxw041")
        if (lfcThreshold > 0) 
            message("lfcThreshold is not used by type='ashr'")
        betahat <- res$log2FoldChange
        sebetahat <- res$lfcSE
        fit <- ashr::ash(betahat, sebetahat, mixcompdist = "normal", 
            method = "shrink", ...)
        res$log2FoldChange <- fit$result$PosteriorMean
        res$lfcSE <- fit$result$PosteriorSD
        mcols(res)$description[2] <- sub("MLE", "MMSE", mcols(res)$description[2])
        mcols(res)$description[3] <- sub("standard error", "posterior SD", 
            mcols(res)$description[3])
        if (svalue) {
            coefAlphaSpaces <- sub(".*p-value: ", "", mcols(res)$description[5])
            res <- res[, 1:3]
            res$svalue <- fit$result$svalue
            mcols(res)[4, ] <- DataFrame(type = "results", description = paste("s-value:", 
                coefAlphaSpaces))
        }
        else {
            res <- res[, c(1:3, 5:6)]
        }
        metadata(res)[["lfcThreshold"]] <- lfcThreshold
        priorInfo(res) <- list(type = "ashr", package = "ashr", 
            version = packageVersion("ashr"), fitted_g = fit$fitted_g)
        res <- resultsFormatSwitch(object = dds, res = res, format = format)
        if (returnList) {
            return(list(res = res, fit = fit))
        }
        else {
            return(res)
        }
    }

Maybe @mikelove can comment on the code https://github.com/mikelove/DESeq2/blob/master/R/lfcShrink.R#L453

[mikelove]

Right, so we are passing MLE LFC and SE from DESeq2 and then extracting the posterior mean, SD, and s-value columns.

[pcarbo]

@rfarouni Do you have the latest version of ashr (available from GitHub)?

[rfarouni]

I have version 2.2-47

[stephens999]

ok, good, that version does not have the code line that ensures they are non-negative. Update your version from github and you should be good