diff --git a/.Rbuildignore b/.Rbuildignore index 8bb0571..ca52a6f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -14,3 +14,4 @@ ucsc.knownGenes.db ^appveyor\.yml$ setupEnv.R +^\.github$ diff --git a/.github/workflows/check-bioc.yml b/.github/workflows/check-bioc.yml new file mode 100644 index 0000000..95605be --- /dev/null +++ b/.github/workflows/check-bioc.yml @@ -0,0 +1,333 @@ +## Read more about GitHub actions the features of this GitHub Actions workflow +## at https://lcolladotor.github.io/biocthis/articles/biocthis.html#use_bioc_github_action +## +## For more details, check the biocthis developer notes vignette at +## https://lcolladotor.github.io/biocthis/articles/biocthis_dev_notes.html +## +## You can add this workflow to other packages using: +## > biocthis::use_bioc_github_action() +## +## Using GitHub Actions exposes you to many details about how R packages are +## compiled and installed in several operating system.s +### If you need help, please follow the steps listed at +## https://github.com/r-lib/actions#where-to-find-help +## +## If you found an issue specific to biocthis's GHA workflow, please report it +## with the information that will make it easier for others to help you. +## Thank you! + +## Acronyms: +## * GHA: GitHub Action +## * OS: operating system + +on: + push: + pull_request: + +name: R-CMD-check-bioc + +## These environment variables control whether to run GHA code later on that is +## specific to testthat, covr, and pkgdown. +## +## If you need to clear the cache of packages, update the number inside +## cache-version as discussed at https://github.com/r-lib/actions/issues/86. +## Note that you can always run a GHA test without the cache by using the word +## "/nocache" in the commit message. +env: + has_testthat: 'true' + run_covr: 'true' + run_pkgdown: 'false' + has_RUnit: 'false' + cache-version: 'cache-v1' + run_docker: 'false' + +jobs: + build-check: + runs-on: ${{ matrix.config.os }} + name: ${{ matrix.config.os }} (${{ matrix.config.r }}) + container: ${{ matrix.config.cont }} + ## Environment variables unique to this job. + + strategy: + fail-fast: false + matrix: + config: + - { os: ubuntu-latest, r: '4.2', bioc: '3.16', cont: "bioconductor/bioconductor_docker:RELEASE_3_16", rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest" } + - { os: ubuntu-latest, r: 'next', bioc: '3.17', cont: "bioconductor/bioconductor_docker:RELEASE_3_17", rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest" } + - { os: macOS-latest, r: '4.2', bioc: '3.16'} + - { os: windows-latest, r: '4.2', bioc: '3.16'} + ## Check https://github.com/r-lib/actions/tree/master/examples + ## for examples using the http-user-agent + env: + R_REMOTES_NO_ERRORS_FROM_WARNINGS: true + RSPM: ${{ matrix.config.rspm }} + NOT_CRAN: true + TZ: UTC + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + + steps: + + ## Set the R library to the directory matching the + ## R packages cache step further below when running on Docker (Linux). + - name: Set R Library home on Linux + if: runner.os == 'Linux' + run: | + mkdir /__w/_temp/Library + echo ".libPaths('/__w/_temp/Library')" > ~/.Rprofile + + ## Most of these steps are the same as the ones in + ## https://github.com/r-lib/actions/blob/master/examples/check-standard.yaml + ## If they update their steps, we will also need to update ours. + - name: Checkout Repository + uses: actions/checkout@v3 + + ## R is already included in the Bioconductor docker images + - name: Setup R from r-lib + if: runner.os != 'Linux' + uses: r-lib/actions/setup-r@v2 + with: + r-version: ${{ matrix.config.r }} + http-user-agent: ${{ matrix.config.http-user-agent }} + + ## pandoc is already included in the Bioconductor docker images + - name: Setup pandoc from r-lib + if: runner.os != 'Linux' + uses: r-lib/actions/setup-pandoc@v2 + + - name: Query dependencies + run: | + install.packages('remotes') + saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2) + shell: Rscript {0} + + - name: Restore R package cache + if: "!contains(github.event.head_commit.message, '/nocache') && runner.os != 'Linux'" + uses: actions/cache@v3 + with: + path: ${{ env.R_LIBS_USER }} + key: ${{ env.cache-version }}-${{ runner.os }}-biocversion-RELEASE_3_16-r-4.2-${{ hashFiles('.github/depends.Rds') }} + restore-keys: ${{ env.cache-version }}-${{ runner.os }}-biocversion-RELEASE_3_16-r-4.2- + + - name: Cache R packages on Linux + if: "!contains(github.event.head_commit.message, '/nocache') && runner.os == 'Linux' " + uses: actions/cache@v3 + with: + path: /home/runner/work/_temp/Library + key: ${{ env.cache-version }}-${{ runner.os }}-biocversion-RELEASE_3_16-r-4.2-${{ hashFiles('.github/depends.Rds') }} + restore-keys: ${{ env.cache-version }}-${{ runner.os }}-biocversion-RELEASE_3_16-r-4.2- + + - name: Install Linux system dependencies + if: runner.os == 'Linux' + run: | + sysreqs=$(Rscript -e 'cat("apt-get update -y && apt-get install -y", paste(gsub("apt-get install -y ", "", remotes::system_requirements("ubuntu", "20.04")), collapse = " "))') + echo $sysreqs + sudo -s eval "$sysreqs" + + - name: Install macOS system dependencies + if: matrix.config.os == 'macOS-latest' + run: | + ## Enable installing XML from source if needed + brew install libxml2 + echo "XML_CONFIG=/usr/local/opt/libxml2/bin/xml2-config" >> $GITHUB_ENV + + ## Required to install magick as noted at + ## https://github.com/r-lib/usethis/commit/f1f1e0d10c1ebc75fd4c18fa7e2de4551fd9978f#diff-9bfee71065492f63457918efcd912cf2 + brew install imagemagick@6 + + ## For textshaping, required by ragg, and required by pkgdown + brew install harfbuzz fribidi + + ## For installing usethis's dependency gert + brew install libgit2 + + ## Required for tcltk + brew install xquartz --cask + + - name: Install Windows system dependencies + if: runner.os == 'Windows' + run: | + ## Edit below if you have any Windows system dependencies + shell: Rscript {0} + + - name: Init tinytex + uses: r-lib/actions/setup-tinytex@v2 + + - name: Install tinytex packages + run: | + system("tlmgr --version") + install.packages("tinytex") + tinytex::tlmgr_install(pkgs = c("bera", "caption", "changepage", "enumitem", "everysel", "fancyhdr", "footmisc", "grfext", "index", "marginfix", "mathtools", + "ms", "nowidow", "parnotes", "parskip", "placeins", "preprint", "ragged2e", "side", "soul", "titlesec", "tocbibind", "xstring")) + shell: Rscript {0} + + - name: Install BiocManager + run: | + message(paste('****', Sys.time(), 'installing BiocManager ****')) + remotes::install_cran("BiocManager") + shell: Rscript {0} + + - name: Set BiocVersion + run: | + BiocManager::install(version = "${{ matrix.config.bioc }}", ask = FALSE, force = TRUE) + shell: Rscript {0} + + - name: Install dependencies pass 1 + run: | + ## Try installing the package dependencies in steps. First the local + ## dependencies, then any remaining dependencies to avoid the + ## issues described at + ## https://stat.ethz.ch/pipermail/bioc-devel/2020-April/016675.html + ## https://github.com/r-lib/remotes/issues/296 + ## Ideally, all dependencies should get installed in the first pass. + + ## Set the repos source depending on the OS + ## Alternatively use https://storage.googleapis.com/bioconductor_docker/packages/ + ## though based on https://bit.ly/bioc2021-package-binaries + ## the Azure link will be the main one going forward. + gha_repos <- if( + .Platform$OS.type == "unix" && Sys.info()["sysname"] != "Darwin" + ) c( + "AnVIL" = "https://bioconductordocker.blob.core.windows.net/packages/3.16/bioc", + BiocManager::repositories() + ) else BiocManager::repositories() + + ## For running the checks + message(paste('****', Sys.time(), 'installing rcmdcheck and BiocCheck ****')) + install.packages(c("rcmdcheck", "BiocCheck"), repos = gha_repos) + + ## Pass #1 at installing dependencies + ## This pass uses AnVIL-powered fast binaries + ## details at https://github.com/nturaga/bioc2021-bioconductor-binaries + ## The speed gains only apply to the docker builds. + message(paste('****', Sys.time(), 'pass number 1 at installing dependencies: local dependencies ****')) + remotes::install_local(dependencies = TRUE, repos = gha_repos, build_vignettes = FALSE, upgrade = TRUE) + continue-on-error: true + shell: Rscript {0} + + - name: Install dependencies pass 2 + run: | + ## Pass #2 at installing dependencies + ## This pass does not use AnVIL and will thus update any packages + ## that have seen been updated in Bioconductor + message(paste('****', Sys.time(), 'pass number 2 at installing dependencies: any remaining dependencies ****')) + remotes::install_local(dependencies = TRUE, repos = BiocManager::repositories(), build_vignettes = TRUE, upgrade = TRUE, force = TRUE) + shell: Rscript {0} + + - name: Install BiocGenerics + if: env.has_RUnit == 'true' + run: | + ## Install BiocGenerics + BiocManager::install("BiocGenerics") + shell: Rscript {0} + + - name: Install covr + if: github.ref == 'refs/heads/master' && env.run_covr == 'true' && runner.os == 'Linux' + run: | + remotes::install_cran("covr") + shell: Rscript {0} + + - name: Install pkgdown + if: github.ref == 'refs/heads/master' && env.run_pkgdown == 'true' && runner.os == 'Linux' + run: | + remotes::install_cran("pkgdown") + shell: Rscript {0} + + - name: Session info + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + shell: Rscript {0} + + - name: Run CMD check + env: + _R_CHECK_CRAN_INCOMING_: false + DISPLAY: 99.0 + run: | + options(crayon.enabled = TRUE) + rcmdcheck::rcmdcheck( + args = c("--no-manual", "--no-vignettes", "--timings"), + build_args = c("--no-manual", "--keep-empty-dirs", "--no-resave-data"), + error_on = "error", + check_dir = "check" + ) + shell: Rscript {0} + + ## Might need an to add this to the if: && runner.os == 'Linux' + - name: Reveal testthat details + if: env.has_testthat == 'true' + run: find . -name testthat.Rout -exec cat '{}' ';' + + - name: Run RUnit tests + if: env.has_RUnit == 'true' + run: | + BiocGenerics:::testPackage() + shell: Rscript {0} + + - name: Run BiocCheck + env: + DISPLAY: 99.0 + run: | + BiocCheck::BiocCheck( + dir('check', 'tar.gz$', full.names = TRUE), + `quit-with-status` = TRUE, + `no-check-R-ver` = TRUE, + `no-check-bioc-help` = TRUE + ) + shell: Rscript {0} + + - name: Test coverage + if: github.ref == 'refs/heads/master' && env.run_covr == 'true' && runner.os == 'Linux' + run: | + covr::codecov() + shell: Rscript {0} + + - name: Install package + if: github.ref == 'refs/heads/master' && env.run_pkgdown == 'true' && runner.os == 'Linux' + run: R CMD INSTALL . + + - name: Build pkgdown site + if: github.ref == 'refs/heads/master' && env.run_pkgdown == 'true' && runner.os == 'Linux' + run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE) + shell: Rscript {0} + ## Note that you need to run pkgdown::deploy_to_branch(new_process = FALSE) + ## at least one locally before this will work. This creates the gh-pages + ## branch (erasing anything you haven't version controlled!) and + ## makes the git history recognizable by pkgdown. + + - name: Install deploy dependencies + if: github.ref == 'refs/heads/master' && env.run_pkgdown == 'true' && runner.os == 'Linux' + run: | + apt-get update && apt-get -y install rsync + + - name: Deploy pkgdown site to GitHub pages 🚀 + if: github.ref == 'refs/heads/master' && env.run_pkgdown == 'true' && runner.os == 'Linux' + uses: JamesIves/github-pages-deploy-action@releases/v4 + with: + clean: false + branch: gh-pages + folder: docs + + - name: Upload check results + if: failure() + uses: actions/upload-artifact@master + with: + name: ${{ runner.os }}-biocversion-RELEASE_3_16-r-4.2-results + path: check + + ## Note that DOCKER_PASSWORD is really a token for your dockerhub + ## account, not your actual dockerhub account password. + ## This comes from + ## https://seandavi.github.io/BuildABiocWorkshop/articles/HOWTO_BUILD_WORKSHOP.html#6-add-secrets-to-github-repo + ## Check https://github.com/docker/build-push-action/tree/releases/v1 + ## for more details. + - uses: docker/build-push-action@v1 + if: "!contains(github.event.head_commit.message, '/nodocker') && env.run_docker == 'true' && runner.os == 'Linux' " + with: + username: ${{ secrets.DOCKER_USERNAME }} + password: ${{ secrets.DOCKER_PASSWORD }} + repository: gagneurlab/outrider + tag_with_ref: true + tag_with_sha: true + tags: latest diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml deleted file mode 100644 index 87cb157..0000000 --- a/.gitlab-ci.yml +++ /dev/null @@ -1,47 +0,0 @@ -variables: - PKGNAME: "OUTRIDER" - PKGDIR: "./" - R_LIBS_U: "tmp_lib" - RVER: "3.6.0-Bioc3.9" - TEXVER: "2016" - # CI_DEBUG_TRACE: "true" - -stages: - - testing - -before_script: - - module load "i12g/R/$RVER" - - module load "i12g/texlive/$TEXVER" - - export R_LIBS_U="`realpath ${CI_JOB_ID}_${R_LIBS_U}`" - - mkdir $R_LIBS_U - - R --vanilla CMD build --no-build-vignettes --no-manual $PKGDIR - - PKG_FILE_NAME=$(ls -1t *.tar.gz | head -n 1) - - R_LIBS_USER=$R_LIBS_U Rscript --vanilla -e "install.packages('${PKG_FILE_NAME}', lib='$R_LIBS_U', repo=NULL)" - -R-check: - stage: testing - when: always - tags: [shell] - script: - - R_LIBS_USER=$R_LIBS_U R --vanilla CMD check --no-vignettes --timings "${PKG_FILE_NAME}" - -R-BiocCheck: - stage: testing - when: always - tags: [shell] - script: - - R_LIBS_USER=$R_LIBS_U R --vanilla CMD BiocCheck --no-check-vignettes "${PKG_FILE_NAME}" - -R-coverage: - stage: testing - when: always - tags: [shell] - script: - - R_LIBS_USER=$R_LIBS_U Rscript --vanilla -e "library('$PKGNAME'); covr::package_coverage('$PKGDIR')" - -R-check_vignettes: - stage: testing - when: always - tags: [shell] - script: - - cd "$PKGDIR/vignettes/" && R_LIBS_USER=$R_LIBS_U R --vanilla CMD Sweave --engine=knitr::knitr --pdf "$PKGNAME.Rnw" diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 978852b..0000000 --- a/.travis.yml +++ /dev/null @@ -1,54 +0,0 @@ -language: r -sudo: required -dist: xenial -cache: - - directories: - - $HOME/R/Library - - -addons: - apt: - packages: - - libxml2-dev - - tk-dev - - libgit2-dev - -r: - - bioc-devel - - bioc-release - - 3.6.0 - -stages: - - Setup Cache - - Test - -# setup cache first -jobs: - include: - - stage: Setup Cache - r: bioc-devel - script: true - - stage: Setup Cache - r: bioc-release - script: true - - stage: Setup Cache - r: 3.6.0 - script: true - -before_install: - - tlmgr install index marginfix bera nowidow parnotes - -install: - - Rscript setupEnv.R - - R --version - - R -e "BiocManager::version()" - -script: - - R CMD build . - - R -e 'BiocCheck::BiocCheck(".", `no-check-vignettes`=TRUE)' - - R CMD check --no-vignettes --timings *tar.gz - - R -e 'devtools::run_examples(run=TRUE, document=TRUE)' - - cd vignettes && R CMD Sweave --engine=knitr::knitr --pdf OUTRIDER.Rnw - -after_success: - - test $TRAVIS_BUILD_STAGE_NAME = "Test" && Rscript -e 'covr::codecov()' diff --git a/DESCRIPTION b/DESCRIPTION index 8fa4dab..03f6ffa 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: OUTRIDER Title: OUTRIDER - OUTlier in RNA-Seq fInDER Type: Package -Version: 1.16.1 +Version: 1.16.2 Date: 2020-10-20 URL: https://github.com/gagneurlab/OUTRIDER BugRepots: https://github.com/gagneurlab/OUTRIDER/issues @@ -9,11 +9,13 @@ Authors@R: c( person("Felix", "Brechtmann", role=c("aut"), email="brechtma@in.tum.de"), person("Christian", "Mertes", role=c("aut", "cre"), - email="mertes@in.tum.de"), + email="mertes@in.tum.de", comment=c(ORCID="0000-0002-1091-205X")), person("Agne", "Matuseviciute", role=c("aut")), person(c("Michaela", "Fee"), "Müller", role=c("ctb")), - person("Vicente", "Yepez", role=c("aut"), email="yepez@in.tum.de"), - person("Julien", "Gagneur", role=c("aut"), email="gagneur@in.tum.de")) + person("Vicente", "Yepez", role=c("aut"), email="yepez@in.tum.de", + comment=c(ORCID="0000-0001-7916-3643")), + person("Julien", "Gagneur", role=c("aut"), email="gagneur@in.tum.de", + comment=c(ORCID="0000-0002-8924-8365"))) Description: Identification of aberrant gene expression in RNA-seq data. Read count expectations are modeled by an autoencoder to control for confounders in the data. Given these expectations, the RNA-seq read counts @@ -27,7 +29,7 @@ biocViews: ImmunoOncology, RNASeq, Transcriptomics, Alignment, Sequencing, License: MIT + file LICENSE NeedsCompilation: yes Encoding: UTF-8 -RoxygenNote: 7.1.1 +RoxygenNote: 7.2.3 Depends: R (>= 3.6), BiocParallel, @@ -53,7 +55,6 @@ Imports: pcaMethods, PRROC, RColorBrewer, - Rcpp, reshape2, S4Vectors, scales, diff --git a/NAMESPACE b/NAMESPACE index aa2cc67..3ab9fc0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -73,6 +73,7 @@ importFrom(BiocGenerics,plotDispEsts) importFrom(BiocParallel,SerialParam) importFrom(BiocParallel,bpisup) importFrom(BiocParallel,bplapply) +importFrom(BiocParallel,bpmapply) importFrom(BiocParallel,bpparam) importFrom(BiocParallel,bpstart) importFrom(BiocParallel,bpstop) diff --git a/NEWS b/NEWS index 237c70f..215e368 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,9 @@ +CHANGES IN VERSION 1.16.2 +----------------------------- + o Add option to restrict FDR correction to sets of genes of interest + o Add opttion to retrieve results based on those FDR values on subsets of + genes + CHANGES IN VERSION 1.7.1 ----------------------------- o Move to S3/S4 methods to be compatible with FRASER diff --git a/R/OUTRIDER.R b/R/OUTRIDER.R index 6e39227..0016258 100644 --- a/R/OUTRIDER.R +++ b/R/OUTRIDER.R @@ -22,6 +22,12 @@ #' @inheritParams controlForConfounders #' @param controlData If TRUE, the default, the raw counts are controled #' for confounders by the autoencoder +#' @param subsets A named list of named lists specifying any number of gene +#' subsets (can differ per sample). For each subset, FDR correction +#' will be limited to genes in the subset, and the FDR corrected +#' pvalues stored as an assay in the ods object in addition to the +#' transcriptome-wide FDR corrected pvalues. See the examples for +#' how to use this argument. #' @param ... Further arguments passed on to \code{controlForConfounders} #' @return OutriderDataSet with all the computed values. The values are stored #' as assays and can be accessed by: \code{assay(ods, 'value')}. @@ -44,9 +50,19 @@ #' plotAberrantPerSample(ods) #' plotVolcano(ods, 1) #' +#' # example of restricting FDR correction to subsets of genes of interest +#' genesOfInterest <- list("sample_1"=sample(rownames(ods), 3), +#' "sample_2"=sample(rownames(ods), 8), +#' "sample_6"=sample(rownames(ods), 5)) +#' genesOfInterest +#' ods <- OUTRIDER(ods, subsets=list("exampleSubset"=genesOfInterest)) +#' padj(ods, subsetName="exampleSubset")[1:10,1:10] +#' res <- results(ods, all=TRUE) +#' res +#' #' @export OUTRIDER <- function(ods, q, controlData=TRUE, implementation='autoencoder', - BPPARAM=bpparam(), ...){ + subsets=NULL, BPPARAM=bpparam(), ...){ checkOutriderDataSet(ods) implementation <- tolower(implementation) @@ -65,7 +81,7 @@ OUTRIDER <- function(ods, q, controlData=TRUE, implementation='autoencoder', } message(date(), ": P-value calculation ...") - ods <- computePvalues(ods, BPPARAM=BPPARAM) + ods <- computePvalues(ods, subsets=subsets, BPPARAM=BPPARAM) message(date(), ": Zscore calculation ...") ods <- computeZscores(ods, diff --git a/R/getNSetterFuns.R b/R/getNSetterFuns.R index 6d14e1b..3cc938a 100644 --- a/R/getNSetterFuns.R +++ b/R/getNSetterFuns.R @@ -4,6 +4,8 @@ #' the values within the OUTRIDER model. #' #' @param ods,object An OutriderDataSet object. +#' @param subsetName Name of a gene subset for which to store or retrieve +#' FDR corrected p values #' @param ... Further arguments passed on to the underlying assay function. #' @return A matrix or vector dependent on the type of data retrieved. #' @@ -58,18 +60,26 @@ pValue <- function(ods){ #' @rdname getter_setter_functions #' @export padj -padj <- function(ods){ - if(!'padjust' %in% assayNames(ods)){ +padj <- function(ods, subsetName=NULL){ + assayName <- 'padjust' + if(!is.null(subsetName)){ + assayName <- paste0(assayName, "_", subsetName) + } + if(!assayName %in% assayNames(ods)){ stop('Please compute first the P-values before retrieving', ' the adjusted ones.') } - assay(ods, 'padjust') + assay(ods, assayName) } -`padj<-` <- function(ods, ..., value){ +`padj<-` <- function(ods, subsetName=NULL, ..., value){ stopifnot(is.matrix(value)) stopifnot(dim(ods) == dim(value)) - assay(ods, 'padjust', ...) <- value + assayName <- 'padjust' + if(!is.null(subsetName)){ + assayName <- paste0(assayName, "_", subsetName) + } + assay(ods, assayName, ...) <- value return(ods) } diff --git a/R/method-evaluation.R b/R/method-evaluation.R index 43a780a..0cdd431 100644 --- a/R/method-evaluation.R +++ b/R/method-evaluation.R @@ -1,8 +1,9 @@ aberrant.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, - by=c("none", "sample", "gene")){ + by=c("none", "sample", "gene"), subsetName=NULL, ...){ checkFullAnalysis(object) - aberrantEvents <- padj(object) <= padjCutoff + aberrantEvents <- padj(object, subsetName=subsetName) <= padjCutoff + if(isScalarNumeric(zScoreCutoff, na.ok=FALSE)){ aberrantEvents <- aberrantEvents & abs(zScore(object)) >= zScoreCutoff } @@ -26,6 +27,10 @@ aberrant.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, #' if NA or NULL no Z-score cutoff is used #' @param by if the results should be summarized by 'sample', #' 'gene' or not at all (default). +#' @param subsetName The name of a subset of genes of interest for which FDR +#' corected pvalues have been previously computed. Those FDR values +#' on the subset will then be used to determine aberrant status. +#' Default is NULL (using transcriptome-wide FDR corrected pvalues). #' @param ... Currently not in use. #' #' @return The number of aberrent events by gene or sample or a TRUE/FALSE diff --git a/R/method-gridSearch.R b/R/method-gridSearch.R index 4d980a0..ee114ab 100644 --- a/R/method-gridSearch.R +++ b/R/method-gridSearch.R @@ -200,7 +200,14 @@ injectOutliers <- function(ods, freq, zScore, inj, lnorm, sdlog){ # only insert outliers if they are different from before # and not too large if(art_out < max_out & counts[idxRow, idxCol] != art_out){ - counts[idxRow, idxCol] <- art_out + + # do not introduce outlier such that all samples have 0 counts for that gene. + # if all other samples (excluding idxCol) are 0 and art_out is 0 do not inject that outlier. + if(all(counts[idxRow,-idxCol] == 0) & art_out == 0){ + index[idxRow, idxCol] <- 0 + }else{ + counts[idxRow, idxCol] <- art_out + } }else{ index[idxRow, idxCol] <- 0 zScore diff --git a/R/method-results.R b/R/method-results.R index 1e35b0c..d1e68c3 100644 --- a/R/method-results.R +++ b/R/method-results.R @@ -1,5 +1,5 @@ compileResults.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, - round=2, all=FALSE, ...){ + round=2, all=FALSE, subsetName=NULL, ...){ # # input check and parsing @@ -19,14 +19,17 @@ compileResults.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, round <- 2 } + meanCorrectedCounts <- rowMeans(counts(object, normalized=TRUE)) if(isFALSE(all)){ abByGene <- aberrant(object, padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, by="gene") abBySample <- aberrant(object, padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, by="sample") object <- object[abByGene > 0, abBySample > 0] + meanCorrectedCounts <- meanCorrectedCounts[abByGene > 0] } + # warning if no rows to return if(nrow(object) == 0){ if(isFALSE(all)){ warning('No significant events: use all=TRUE to print all events.') @@ -34,43 +37,50 @@ compileResults.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, warning('Please provide an object with at least one feature.') } return(data.table( - geneID=NA_character_, - sampleID=NA_character_, - pValue=NA_real_, - padjust=NA_real_, - zScore=NA_real_, - l2fc=NA_real_, - rawcounts=NA_integer_, - normcounts=NA_real_, - meanCorrected=NA_real_, - theta=NA_real_, - aberrant=NA, - AberrantBySample=NA_integer_, - AberrantByGene=NA_integer_, - padj_rank=NA_real_)[0]) + geneID=NA_character_, + sampleID=NA_character_, + pValue=NA_real_, + padjust=NA_real_, + zScore=NA_real_, + l2fc=NA_real_, + rawcounts=NA_integer_, + normcounts=NA_real_, + meanCorrected=NA_real_, + theta=NA_real_, + aberrant=NA, + AberrantBySample=NA_integer_, + AberrantByGene=NA_integer_, + padj_rank=NA_real_, + FDR_set=NA_character_)[0]) } # # extract data - # + # ans <- data.table( geneID = rownames(object), sampleID = rep(colnames(object), each=nrow(object)), pValue = c(pValue(object)), - padjust = c(padj(object)), + padjust = c(padj(object, subsetName=subsetName)), zScore = c(zScore(object)), l2fc = c(assay(object, "l2fc")), rawcounts = c(counts(object)), normcounts = c(counts(object, normalized=TRUE)), - meanCorrected = rowMeans(counts(object, normalized=TRUE)), + meanCorrected = meanCorrectedCounts, theta = theta(object), aberrant = c(aberrant(object, - padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff)), + padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, + subsetName=subsetName)), AberrantBySample = rep(each=nrow(object), aberrant(object, - padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, by="sample")), + padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, by="sample", + subsetName=subsetName)), AberrantByGene = aberrant(object, - padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, by="gene"), - padj_rank = c(apply(padj(object), 2, rank))) + padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, by="gene", + subsetName=subsetName), + padj_rank = c(apply(padj(object), 2, rank)), + FDR_set = ifelse(is.null(subsetName), "transcriptome-wide", + subsetName) + ) # round columns if requested if(is.numeric(round)){ @@ -80,16 +90,65 @@ compileResults.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, } # - # keep only aberrent events and sort by padj value + # keep only aberrant events and sort by padj value # if(isFALSE(all)){ ans <- ans[aberrant == TRUE] + + } else{ + # if return full subset is requested, retrieve those as all non-NA padj + ans <- ans[!is.na(padjust),] } ans <- ans[order(padjust)] return(ans) } +compileResultsAll.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, + round=2, all=FALSE, + returnTranscriptomewideResults=TRUE, ...){ + # + # input check and parsing + # + checkOutriderDataSet(object) + checkFullAnalysis(object) + + # get all padjust assays (transcriptome-wide and on subsets) + padjAssays <- grep('padjust', assayNames(object), value=TRUE) + + # dont retrieve transcriptome-wide results if requested + if(isFALSE(returnTranscriptomewideResults)){ + if(length(padjAssays) == 1 && padjAssays == 'padjust'){ + warning("Retrieving transcriptome-wide results as no other ", + "padjust assays are available in the ods object.") + } else{ + padjAssays <- padjAssays[padjAssays != 'padjust'] + } + } + + # get results for all available padjust assays + resSubsets <- lapply(padjAssays, function(padjAssay){ + if(padjAssay == 'padjust'){ + subsetName <- NULL + } else{ + subsetName <- gsub("padjust_", "", padjAssay) + } + resSub <- compileResults.OUTRIDER(object=object, + padjCutoff=padjCutoff, zScoreCutoff=zScoreCutoff, + all=all, round=round, + subsetName=subsetName, ...) + return(resSub) + }) + res <- rbindlist(resSubsets) + + # sort it if existing + if(length(res) > 0){ + res <- res[order(padjust)] + } + + return(res) +} + #' #' Accessor function for the 'results' object in an OutriderDataSet object. #' @@ -105,9 +164,15 @@ compileResults.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, #' more user friendly #' @param all By default FALSE, only significant read counts are listed in the #' results. If TRUE all results are assembled resulting in a -#' data.table of length samples x genes +#' data.table of length samples x genes. If FDR corrected pvalues +#' for subsets of genes of interest have been calculated, this +#' indicates whether the whole subset will be listed in the results. +#' @param returnTranscriptomewideResults If FDR corrected pvalues for subsets +#' of genes of interest have been calculated, this parameter +#' indicates whether additionally the transcriptome-wide results +#' should be returned as well (default), or whether only results +#' for those subsets should be retrieved. #' @param ... Additional arguments, currently not used -#' #' @return A data.table where each row is an outlier event and the columns #' contain additional information about this event. Eg padj, l2fc #' @@ -122,10 +187,20 @@ compileResults.OUTRIDER <- function(object, padjCutoff=0.05, zScoreCutoff=0, #' res <- results(ods, all=TRUE) #' res #' +#' # example of retrieving results with FDR correction limited to a +#' # set of genes of interest +#' genesOfInterest <- list("sample_1"=sample(rownames(ods), 3), +#' "sample_2"=sample(rownames(ods), 8), +#' "sample_6"=sample(rownames(ods), 5)) +#' genesOfInterest +#' ods <- computePvalues(ods, subsets=list("exampleSubset"=genesOfInterest)) +#' res <- results(ods, all=TRUE, returnTranscriptomewideResults=FALSE) +#' res +#' #' @name results #' @rdname results #' @aliases results results,OutriderDataSet-method #' #' @export -setMethod("results", "OutriderDataSet", compileResults.OUTRIDER) +setMethod("results", "OutriderDataSet", compileResultsAll.OUTRIDER) diff --git a/R/pValMatrix.R b/R/pValMatrix.R index 3dbb05b..f011e80 100644 --- a/R/pValMatrix.R +++ b/R/pValMatrix.R @@ -11,6 +11,12 @@ #' the alternative hypothesis used to calculate the P-values, #' defaults to "two.sided" #' @param method Method used for multiple testing +#' @param subsets A named list of named lists specifying any number of gene +#' subsets (can differ per sample). For each subset, FDR correction +#' will be limited to genes in the subset, and the FDR corrected +#' pvalues stored as an assay in the ods object in addition to the +#' transcriptome-wide FDR corrected pvalues. See the examples for +#' how to use this argument. #' @param BPPARAM Can be used to parallelize the computation, defaults to #' bpparam() #' @param ... additional params, currently not used. @@ -36,6 +42,16 @@ #' #' assays(ods)[['pValue']][1:10,1:10] #' +#' # example of restricting FDR correction to subsets of genes of interest +#' genesOfInterest <- list("sample_1"=sample(rownames(ods), 3), +#' "sample_2"=sample(rownames(ods), 8), +#' "sample_6"=sample(rownames(ods), 5)) +#' ods <- computePvalues(ods, subsets=list("exampleSubset"=genesOfInterest)) +#' padj(ods, subsetName="exampleSubset")[1:10,1:10] +#' ods <- computePvalues(ods, +#' subsets=list("anotherExampleSubset"=rownames(ods)[1:5])) +#' padj(ods, subsetName="anotherExampleSubset")[1:10,1:10] +#' #' @exportMethod computePvalues setGeneric("computePvalues", function(object, ...) standardGeneric("computePvalues")) @@ -44,13 +60,26 @@ setGeneric("computePvalues", #' @export setMethod("computePvalues", "OutriderDataSet", function(object, alternative=c("two.sided", "greater", "less"), method='BY', - BPPARAM=bpparam()){ + subsets=NULL, BPPARAM=bpparam()){ alternative <- match.arg(alternative) object <- pValMatrix(object, alternative, BPPARAM=BPPARAM) if(method != 'None'){ object <- padjMatrix(object, method) + + # calculate FDR for each provided subset and assign to ods + if(!is.null(subsets)){ + stopifnot(is.list(subsets)) + stopifnot(!is.null(names(subsets))) + for(setName in names(subsets)){ + geneListSubset <- subsets[[setName]] + object <- padjOnSubset(ods=object, method=method, + BPPARAM=BPPARAM, + genesToTest=geneListSubset, + subsetName=setName) + } + } } object }) @@ -145,3 +174,84 @@ padjMatrix <- function(ods, method='BH'){ return(ods) } +#' +#' FDR correction on subset of genes +#' +#' Function to correct the computed pValues using BY FDR correction, limited to +#' a subset of genes. +#' +#' @param ods OutriderDataSet +#' @param method adjustment method, by default 'BY' +#' @param genesToTest A named list with the subset of genes to test per sample. +#' The names must correspond to the sampleIDs in the ods object. +#' @param subsetName The name under which the resulting FDR corrected pvalues +#' will be stored as an assay with name 'padjust_'+subsetName. +#' @return OutriderDataSet containing adjusted pValues +#' +#' @noRd +padjOnSubset <- function(ods, genesToTest, subsetName, method='BY', + BPPARAM=bpparam()){ + + # check input + stopifnot(!is.null(genesToTest)) + stopifnot(is.list(genesToTest) || is.vector(genesToTest)) + + # replicate subset genes for each sample if given as vector + if(!is.list(genesToTest)){ + genesToTest <- rep(list(genesToTest), ncol(ods)) + names(genesToTest) <- colnames(ods) + } + + # check that names are present and correspond to samples in ods + stopifnot(!is.null(names(genesToTest))) + if(!all(names(genesToTest) %in% colnames(ods))){ + stop("names(genesToTest) need to be sampleIDs in the given ods object.") + } + + # compute FDR on the given subsets of genes + fdrSubset <- bpmapply(colnames(ods), FUN=function(sampleId){ + + # get genes to test for this sample + genesToTestSample <- genesToTest[[sampleId]] + padj <- rep(NA, nrow(ods)) + + # if no genes present in the subset for this sample, return NAs + if(is.null(genesToTestSample)){ + return(padj) + } + + # get idx of junctions corresponding to genes to test + if(is.character(genesToTestSample)){ + rowIdx <- sort(which(rownames(ods) %in% genesToTestSample)) + } else{ + stop("Genes in the list to test must be a character vector ", + "of geneIDs.") + } + + # check that rowIdx is not empty vector + if(length(rowIdx) == 0){ + warning("No genes from the given subset found in the ods ", + "object for sample: ", sampleId) + return(padj) + } + + # retrieve pvalues of genes to test + p <- pValue(ods)[rowIdx, sampleId] + + # FDR correction on subset + padjSub <- p.adjust(p, method=method) + + # return FDR on subset, filled with NA for all other genes + padj[rowIdx] <- padjSub + return(padj) + + }, SIMPLIFY=TRUE, BPPARAM=BPPARAM) + rownames(fdrSubset) <- rownames(ods) + colnames(fdrSubset) <- colnames(ods) + + # add FDR subset info to ods object and return + padj(ods, subsetName=subsetName) <- fdrSubset + validObject(ods) + return(ods) +} + diff --git a/R/package-OUTRIDER.R b/R/package-OUTRIDER.R index 0d39f9e..88be64b 100644 --- a/R/package-OUTRIDER.R +++ b/R/package-OUTRIDER.R @@ -24,6 +24,7 @@ #' isScalarValue isScalarNA chunk seq_col seq_row #' #' @importFrom BiocParallel bplapply bpparam SerialParam bpisup bpstart bpstop +#' bpmapply #' #' @importFrom GenomicFeatures makeTxDbFromGFF exonsBy #' diff --git a/README.md b/README.md index 8f1a2c5..647cda3 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,17 @@ # OUTRIDER # +[![Build](https://github.com/gagneurlab/OUTRIDER/actions/workflows/check-bioc.yml/badge.svg)](https://github.com/gagneurlab/OUTRIDER/actions/workflows/check-bioc.yml) +[![Version](https://img.shields.io/github/v/release/gagneurlab/OUTRIDER)](https://github.com/gagneurlab/OUTRIDER/releases) +[![Coverage status](https://codecov.io/gh/gagneurlab/OUTRIDER/branch/master/graph/badge.svg)](https://codecov.io/github/gagneurlab/OUTRIDER?branch=master) +[![License](https://img.shields.io/github/license/mashape/apistatus.svg?maxAge=2592000)](https://github.com/gagneurlab/OUTRIDER/blob/master/LICENSE) + OUTRIDER is a tool to find aberrantly expressed genes in RNA-seq samples. The method is published in the [AJHG](https://doi.org/10.1016/j.ajhg.2018.10.025) and available through [Bioconductor](http://bioconductor.org/packages/release/bioc/html/OUTRIDER.html). +It is also part of the workflow [DROP](https://github.com/gagneurlab/drop), which is described in [Nature Protocols](https://www.nature.com/articles/s41596-020-00462-5). -[![Pipeline status](https://travis-ci.org/gagneurlab/OUTRIDER.svg?branch=master)](https://travis-ci.org/gagneurlab/OUTRIDER) -[![AppVeyor build status](https://ci.appveyor.com/api/projects/status/a2f6io5isq0jhobf/branch/master?svg=true)](https://ci.appveyor.com/project/c-mertes/outrider/branch/master) -[![Version](https://img.shields.io/badge/Version-1.7.1-green.svg)](https://github.com/gagneurlab/OUTRIDER/tree/master) -[![Coverage status](https://codecov.io/gh/gagneurlab/OUTRIDER/branch/master/graph/badge.svg)](https://codecov.io/github/gagneurlab/OUTRIDER?branch=master) -[![License](https://img.shields.io/github/license/mashape/apistatus.svg?maxAge=2592000)](https://github.com/gagneurlab/OUTRIDER/blob/master/LICENSE) +Please cite our method paper if you use it in a publication: + +> Brechtmann F, Mertes C, Matusevičiūtė A, *et al.* OUTRIDER: A Statistical Method for Detecting Aberrantly Expressed Genes in RNA Sequencing Data. *Am J Hum Genet. 2018;103(6):907-917*. https://doi.org/10.1016/j.ajhg.2018.10.025 ## Installation @@ -43,7 +47,7 @@ install.packages('devtools') devtools::install_github('gagneurlab/OUTRIDER', dependencies=TRUE) # installing a specific version/tag of OUTRIDER -devtools::install_github('gagneurlab/OUTRIDER@1.3.5', dependencies=TRUE) +devtools::install_github('gagneurlab/OUTRIDER@1.16.2', dependencies=TRUE) ``` To check which versions/tags are available you can check the GitHub repo @@ -70,7 +74,7 @@ To install those packages, please run as administrator: For Ubuntu or debian based systems: ``` -sudo apt-get install build-essential libcurl4-gnutls-dev libxml2-dev libssl-dev zlib1g-dev libmysqld-dev +sudo apt-get install build-essential libcurl4-gnutls-dev libxml2-dev libssl-dev zlib1g-dev libmariadbd-dev ``` For centOS or RHEL based systems: diff --git a/appveyor.yml b/appveyor.yml deleted file mode 100644 index 13b9136..0000000 --- a/appveyor.yml +++ /dev/null @@ -1,79 +0,0 @@ -# -# Use the Bootstrap code from the rwinlib project for CRAN -# https://github.com/rwinlib -# -init: - ps: | - $ErrorActionPreference = "Stop" - Invoke-WebRequest https://raw.githubusercontent.com/rwinlib/base/master/scripts/appveyor-tool.ps1 -OutFile "..\appveyor-tool.ps1" - Import-Module '..\appveyor-tool.ps1' - -install: - - ps: Bootstrap - # install ImageMagick and set path correctly at first position - # only needed for building the vignette - #- ps: InstallMiktex - #- choco install -y -PackageParameters LegacySupport=true imagemagick - #- set "PATH=C:\Program Files\ImageMagick-7.0.8-Q16;%PATH%" - -#cache: -# - C:\RLibrary - -environment: - # parameters to setup R environment - # https://github.com/krlmlr/r-appveyor/blob/master/README.md - PKGTYPE: win.binary - BIOC_USE_DEVEL: FALSE - USE_RTOOLS: TRUE - KEEP_VIGNETTES: TRUE - GCC_PATH: mingw_64 - R_ARCH: x64 -# APPVEYOR_SAVE_CACHE_ON_ERROR: TRUE - - matrix: -# - R_VERSION: devel -# BIOC_USE_DEVEL: TRUE -# - R_VERSION: devel -# BIOC_USE_DEVEL: TRUE -# GCC_PATH: mingw_32 -# R_ARCH: i386 - - R_VERSION: release -# - R_VERSION: 3.6.0 - -build_script: - - travis-tool.sh install_bioc BiocManager BiocCheck - - R -e "BiocManager::install(Ncpus=5, ask=FALSE, upgrade=TRUE)" - - Rscript setupEnv.R - -test_script: - - R CMD build --no-build-vignettes . - - R CMD check --install-args=--build --no-vignettes --no-manual --timings *tar.gz - - R -e "BiocCheck::BiocCheck(list.files(path='.', pattern='.tar.gz'))" - - R -e "devtools::run_examples(document=FALSE)" - # runs currently into an error saying: "LaTeX Error: Missing \begin{document}." - # - cd vignettes && R CMD Sweave --engine=knitr::knitr --pdf OUTRIDER.Rnw - -on_failure: - - 7z a failure.zip *.Rcheck\* - - appveyor PushArtifact failure.zip - - appveyor PushArtifact \*_*.tar.gz - - appveyor PushArtifact \*_*.zip - -artifacts: - - path: '*.Rcheck\**\*.log' - name: Logs - - - path: '*.Rcheck\**\*.out' - name: Logs - - - path: '*.Rcheck\**\*.fail' - name: Logs - - - path: '*.Rcheck\**\*.Rout' - name: Logs - - - path: '\*_*.tar.gz' - name: Bits - - - path: '\*_*.zip' - name: Bits diff --git a/man/OUTRIDER.Rd b/man/OUTRIDER.Rd index e9e2768..ec684a6 100644 --- a/man/OUTRIDER.Rd +++ b/man/OUTRIDER.Rd @@ -9,6 +9,7 @@ OUTRIDER( q, controlData = TRUE, implementation = "autoencoder", + subsets = NULL, BPPARAM = bpparam(), ... ) @@ -25,6 +26,13 @@ for confounders by the autoencoder} implementation. Also 'pca' and 'peer' can be used to control for confounding effects} +\item{subsets}{A named list of named lists specifying any number of gene +subsets (can differ per sample). For each subset, FDR correction +will be limited to genes in the subset, and the FDR corrected +pvalues stored as an assay in the ods object in addition to the +transcriptome-wide FDR corrected pvalues. See the examples for +how to use this argument.} + \item{BPPARAM}{A \code{\link[BiocParallel:BiocParallelParam-class]{BiocParallelParam}} instance to be used for parallel computing.} @@ -72,4 +80,14 @@ res plotAberrantPerSample(ods) plotVolcano(ods, 1) +# example of restricting FDR correction to subsets of genes of interest +genesOfInterest <- list("sample_1"=sample(rownames(ods), 3), + "sample_2"=sample(rownames(ods), 8), + "sample_6"=sample(rownames(ods), 5)) +genesOfInterest +ods <- OUTRIDER(ods, subsets=list("exampleSubset"=genesOfInterest)) +padj(ods, subsetName="exampleSubset")[1:10,1:10] +res <- results(ods, all=TRUE) +res + } diff --git a/man/aberrant.Rd b/man/aberrant.Rd index 44b3b8c..5ce6152 100644 --- a/man/aberrant.Rd +++ b/man/aberrant.Rd @@ -11,7 +11,9 @@ aberrant(object, ...) object, padjCutoff = 0.05, zScoreCutoff = 0, - by = c("none", "sample", "gene") + by = c("none", "sample", "gene"), + subsetName = NULL, + ... ) } \arguments{ @@ -26,6 +28,11 @@ if NA or NULL no Z-score cutoff is used} \item{by}{if the results should be summarized by 'sample', 'gene' or not at all (default).} + +\item{subsetName}{The name of a subset of genes of interest for which FDR +corected pvalues have been previously computed. Those FDR values +on the subset will then be used to determine aberrant status. +Default is NULL (using transcriptome-wide FDR corrected pvalues).} } \value{ The number of aberrent events by gene or sample or a TRUE/FALSE diff --git a/man/computePvalues.Rd b/man/computePvalues.Rd index f1f9108..6d48e60 100644 --- a/man/computePvalues.Rd +++ b/man/computePvalues.Rd @@ -12,6 +12,7 @@ computePvalues(object, ...) object, alternative = c("two.sided", "greater", "less"), method = "BY", + subsets = NULL, BPPARAM = bpparam() ) } @@ -26,6 +27,13 @@ defaults to "two.sided"} \item{method}{Method used for multiple testing} +\item{subsets}{A named list of named lists specifying any number of gene +subsets (can differ per sample). For each subset, FDR correction +will be limited to genes in the subset, and the FDR corrected +pvalues stored as an assay in the ods object in addition to the +transcriptome-wide FDR corrected pvalues. See the examples for +how to use this argument.} + \item{BPPARAM}{Can be used to parallelize the computation, defaults to bpparam()} } @@ -50,6 +58,16 @@ ods <- computePvalues(ods) assays(ods)[['pValue']][1:10,1:10] +# example of restricting FDR correction to subsets of genes of interest +genesOfInterest <- list("sample_1"=sample(rownames(ods), 3), + "sample_2"=sample(rownames(ods), 8), + "sample_6"=sample(rownames(ods), 5)) +ods <- computePvalues(ods, subsets=list("exampleSubset"=genesOfInterest)) +padj(ods, subsetName="exampleSubset")[1:10,1:10] +ods <- computePvalues(ods, + subsets=list("anotherExampleSubset"=rownames(ods)[1:5])) +padj(ods, subsetName="anotherExampleSubset")[1:10,1:10] + } \seealso{ p.adjust diff --git a/man/getter_setter_functions.Rd b/man/getter_setter_functions.Rd index 2c7e61d..18a6cbc 100644 --- a/man/getter_setter_functions.Rd +++ b/man/getter_setter_functions.Rd @@ -22,7 +22,7 @@ zScore(ods) pValue(ods) -padj(ods) +padj(ods, subsetName = NULL) \S4method{dispersions}{OutriderDataSet}(object, ...) @@ -31,6 +31,9 @@ theta(ods) \arguments{ \item{ods, object}{An OutriderDataSet object.} +\item{subsetName}{Name of a gene subset for which to store or retrieve +FDR corrected p values} + \item{...}{Further arguments passed on to the underlying assay function.} } \value{ diff --git a/man/results.Rd b/man/results.Rd index 37fee56..9702ab2 100644 --- a/man/results.Rd +++ b/man/results.Rd @@ -13,6 +13,7 @@ results(object, ...) zScoreCutoff = 0, round = 2, all = FALSE, + returnTranscriptomewideResults = TRUE, ... ) } @@ -31,7 +32,15 @@ more user friendly} \item{all}{By default FALSE, only significant read counts are listed in the results. If TRUE all results are assembled resulting in a -data.table of length samples x genes} +data.table of length samples x genes. If FDR corrected pvalues +for subsets of genes of interest have been calculated, this +indicates whether the whole subset will be listed in the results.} + +\item{returnTranscriptomewideResults}{If FDR corrected pvalues for subsets +of genes of interest have been calculated, this parameter +indicates whether additionally the transcriptome-wide results +should be returned as well (default), or whether only results +for those subsets should be retrieved.} } \value{ A data.table where each row is an outlier event and the columns @@ -53,4 +62,14 @@ ods <- OUTRIDER(ods) res <- results(ods, all=TRUE) res +# example of retrieving results with FDR correction limited to a +# set of genes of interest +genesOfInterest <- list("sample_1"=sample(rownames(ods), 3), + "sample_2"=sample(rownames(ods), 8), + "sample_6"=sample(rownames(ods), 5)) +genesOfInterest +ods <- computePvalues(ods, subsets=list("exampleSubset"=genesOfInterest)) +res <- results(ods, all=TRUE, returnTranscriptomewideResults=FALSE) +res + } diff --git a/setupEnv.R b/setupEnv.R deleted file mode 100644 index 209e9ed..0000000 --- a/setupEnv.R +++ /dev/null @@ -1,71 +0,0 @@ -BTYPE <- ifelse(.Platform$OS.type == 'unix', "source", "both") -NCPUS <- 6 -START_TIME <- Sys.time() - -print_log <- function(...){ - hash_line <- paste0(rep("#", 10), collapse="") - message(paste0("\n", hash_line, "\n### ", date(), ": ", ..., "\n", hash_line, "\n")) -} - -installIfReq <- function(p, type=BTYPE, Ncpus=NCPUS, ...){ - for(j in p){ - if(!requireNamespace(j, quietly=TRUE)){ - print_log("Install ", j) - INSTALL(j, type=type, Ncpus=Ncpus, ...) - } - } -} - -# install Bioconductor -if(!requireNamespace("BiocManager", quietly=TRUE)){ - print_log("Install BiocManager") - install.packages("BiocManager", Ncpus=NCPUS) -} -INSTALL <- BiocManager::install - -# since the current XML package is not compatible with 3.6 anymore -if(!requireNamespace("XML", quietly=TRUE) & R.version[['major']] == "3"){ - installIfReq(p="devtools", type=BTYPE, Ncpus=NCPUS) - devtools::install_version("XML", version="3.99-0.3") -} - -# because of https://github.com/r-windows/rtools-installer/issues/3 -if("windows" == .Platform$OS.type){ - print_log("Install XML on windows ...") - BTYPE <- "win.binary" - installIfReq(p=c("XML", "xml2", "RSQLite", "progress", "AnnotationDbi", "BiocCheck")) - - print_log("Install source packages only for windows ...") - INSTALL(c("GenomeInfoDbData", "org.Hs.eg.db", "TxDb.Hsapiens.UCSC.hg19.knownGene"), type="both") -} else { - BTYPE <- "source" -} - -# install needed packages -# add testthat to pre installation dependencies due to: https://github.com/r-lib/pkgload/issues/89 -for(p in c("XML", "xml2", "testthat", "devtools", "covr", "roxygen2", "BiocCheck", "R.utils", - "GenomeInfoDbData", "rtracklayer", "hms")){ - installIfReq(p=p, type=BTYPE, Ncpus=NCPUS) -} - -# install OUTRIDER with its dependencies with a timeout due to -# travis (50 min) and appveyor (60 min) set installation warmup to 30 min max -maxTime <- max(30, (60*30 - difftime(Sys.time(), START_TIME, units="sec"))) -R.utils::withTimeout(timeout=maxTime, { - try({ - print_log("Update packages") - INSTALL(ask=FALSE, type=BTYPE, Ncpus=NCPUS) - - print_log("Install OUTRIDER") - devtools::install(".", dependencies=TRUE, upgrade=TRUE, - type=BTYPE, Ncpus=NCPUS) - }) -}) - -# fix knitr for 3.6 for more details see BiocStyle issue 78 -# https://github.com/Bioconductor/BiocStyle/issues/78 -if(R.version[['major']] == "3"){ - options(repos=c(CRAN="http://cran.rstudio.com")) - installIfReq(p="devtools", type=BTYPE, Ncpus=NCPUS) - devtools::install_version("knitr", version="1.29", type="source") -} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 779b0ac..47e4407 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -6,6 +6,11 @@ using namespace Rcpp; +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + // predictMatY arma::mat predictMatY(arma::mat x, arma::mat E, arma::mat D, arma::vec b); RcppExport SEXP _OUTRIDER_predictMatY(SEXP xSEXP, SEXP ESEXP, SEXP DSEXP, SEXP bSEXP) {