StoreyLab / subSeq

Subsampling of high-throughput sequencing count data
MIT License
19 stars 2 forks source link

Unable to use it without treatment #11

Open llrs opened 5 years ago

llrs commented 5 years ago

Nice package and nice idea!

I tried running subsample but I found an error:

library(subSeq)
data(hammer)

hammer.counts = Biobase::exprs(hammer)[, 1:4]
hammer.design = Biobase::pData(hammer)[1:4, ]
hammer.counts = hammer.counts[rowSums(hammer.counts) >= 5, ]

ss = subsample(hammer.counts, c(.01, .1, 1),
               method=c("edgeR", "DESeq2", "voomLimma"))
#> Error in factor(treatment): argument "treatment" is missing, with no default

Created on 2019-05-14 by the reprex package (v0.2.1)

I cannot run the subsample function without treatment, and I cannot provide a treatment that it is not a single vector (I tried with a data.frame). How can we test if we have undersequenced if there are several variables relevant to the condition being tested?

By the way, it was nice to see that it is able to plot directly the summary of ss plot(summary(ss)).

ajbass commented 5 years ago

Thanks for the feedback. In the vignette, we discuss options for users to input their own function. For example, here is one way to input general experimental design matrices in the software:

# create function to take in any design matrix (limma used for example)
my_voom <-
  function(count.matrix, design_matrix, ...) {
    v = limma::voom(counts = count.matrix, design = design_matrix, ...)
    f = limma::lmFit(v, design_matrix)
    e = limma::eBayes(f)
    data.frame(coefficient=f$coefficients[, 2], pvalue=e$p.value[, 2], t.test = e$t[,2])
  }

# Load dataset/filtering
library(subSeq)
data(hammer)
hammer.counts = Biobase::exprs(hammer)[, 1:4]
hammer.design = Biobase::pData(hammer)[1:4, ]
hammer.counts = hammer.counts[rowSums(hammer.counts) >= 5, ]

# subsampling proportions (small for illustration)
proportions = 10^seq(-2, 0, .5)
# Create random covariate and add it to design
random_covariate <- rnorm(4)
design <- model.matrix(~ hammer.design$protocol + random_covariate)
# run subsample using my_voom method
subsamples = subsample(hammer.counts,
                       proportions,
                       method=c("my_voom"),
                       design_matrix=design)

plot(summary(subsamples))

Thus you can even do more complicated testing procedures within the user-defined function (e.g., an F-test between an alternative and null model). Take a look at the vignette to see more details on creating functions within subSeq. Let me know if that works.

llrs commented 5 years ago

I can't read the vignette right now (a problem on my side) but I cannot find that this section on the subSeq.Rnw file on github. Nor in the documentation of the subsample function. Is this explained somewhere what must return that custom function ?

I suppose this solves my problem but I cannot test it now (Close the issue if I am not back to you in a couple of weeks)

Thanks for the tool!

ajbass commented 5 years ago

We briefly mention it in the "Writing Your Own Analysis Method" section of the vignette. We will provide additional implementation details in a future update!