testForDiffUsage {JunctionSeq}R Documentation

Test Junctions for Differential Junction Usage

Description

This function runs the hypothesis tests for differential junction usage.

This function is called internally by the runJunctionSeqAnalyses function, and thus for most purposes users should not need to call this function directly. It may be useful to advanced users performing non-standard analyses.

Usage

testForDiffUsage( jscs, 
    test.formula0 = formula(~ sample + countbin), 
    test.formula1 = formula(~ sample + countbin + condition : countbin),
    method.GLM = c(c("advanced","DESeq2-style"), 
                   c("simpleML","DEXSeq-v1.8.0-style")),
    dispColumn="dispersion", nCores=1, 
    keep.hypothesisTest.fit = FALSE,
    meanCountTestableThreshold = "auto",
    optimizeFilteringForAlpha = 0.01,
    method.cooksFilter = TRUE,
    cooksCutoff,
    pAdjustMethod = "BH",
    verbose = TRUE)

Arguments

jscs

A JunctionSeqCountSet. Usually initially created by readJunctionSeqCounts. Dispersions and size factors must be set, usually using functions estimateJunctionSeqSizeFactors and estimateJunctionSeqDispersions.

test.formula0

The formula for the null hypothesis. Note that the condition to be tested must be named "condition".

test.formula1

The formula for the alternative hypothesis. Note that the condition to be tested must be named "condition".

method.GLM

Character string. Can be used to apply alternative methodologies or implementations. Intended for advanced users who have strong opinions about the underlying statistical methodologies.

The default is "advanced" or, equivalently, "DESeq2-style". This uses the model test methodology used by DESeq2 and DEXSeq v1.12.0 or higher. The alternative method is "simpleML" or, equivalently, "DEXSeq-v1.8.0-style". This uses a simpler maximum-likelihood-based method used by the original DESeq and by some earlier versions of DEXSeq (v1.8.0 or less).

dispColumn

Character value. The name of the fData(jscs) column in which the model dispersion is stored.

nCores

Either an integer or a BiocParallelParam object. Either way, this determines The number of cores to use. Note that multicore functionality may not be available on all platforms. If parallel execution is not available then JunctionSeq will automatically fallback to single-core execution. See the BiocParallel package for more information.

keep.hypothesisTest.fit

Logical value. If TRUE, save both complete hypothesis test model fits for every gene. This will require a lot of memory, but may be useful for statistical diagnostics. Default: FALSE.

meanCountTestableThreshold

"auto" or Numeric value. Features with a total mean normalized count of less than this value will be excluded from the analyses. If left as the default ("auto"), then the cutoff threshold will be determined automatically using the DESeq2 independent filtering method.

optimizeFilteringForAlpha

Numeric value between 0 and 1. If meanCountTestableThreshold is set to "auto" then this sets the adjusted-p-value threshold to optimize against.

method.cooksFilter

Logical value. if TRUE, use the cook's filter to detect and remove outliers.

cooksCutoff

The cook's cutoff threshold to use.

pAdjustMethod

The p-adjustment method to use with the p.adjust function.

verbose

if TRUE, send debugging and progress messages to the console / stdout.

Value

A JunctionSeqCountSet, with hypothesis test results included.

Examples

data(exampleDataSet,package="JctSeqData");
jscs <- testForDiffUsage(jscs);

## > Using single-core execution.

## -------> testJunctionsForDiffUsage: Starting hypothesis test iteration. (Thu Mar 30 17:21:25 2017)

## using supplied model matrix

## found results columns, replacing these

## -------> testJunctionsForDiffUsage: Finished hypothesis test iteration. (Thu Mar 30 17:21:26 2017)

## -------> testJunctionsForDiffUsage: Finished compiling hypothesis test results. (Thu Mar 30 17:21:26 2017)

## ---> tJfDU(): No non-NA maxCooks values. Ignoring cooks.

## > Performing final p.adjust filtering.

## >      No cook's cutoffs found.

## >      Automatically selecting a filtering threshold of 0.113841548108117 to optimize results at the alpha < 0.01 significance level.

## >         (Filtering 0 out of 482 "testable" features, using baseMean < 0.113841548108117)

## >         (Rejected H0 for 13 out of 482 features at alpha < 0.01)

## > Final p.adjust filtering complete.

## Not run: 
##D ########################################
##D #Set up example data:
##D decoder.file <- system.file(
##D                   "extdata/annoFiles/decoder.bySample.txt",
##D                   package="JctSeqData");
##D decoder <- read.table(decoder.file,
##D                   header=TRUE,
##D                   stringsAsFactors=FALSE);
##D gff.file <- system.file(
##D             "extdata/cts/withNovel.forJunctionSeq.gff.gz",
##D             package="JctSeqData");
##D countFiles <- system.file(paste0("extdata/cts/",
##D      decoder$sample.ID,
##D      "/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz"),
##D      package="JctSeqData");
##D ########################################
##D #Advanced Analysis:
##D 
##D #Make a "design" dataframe:
##D design <- data.frame(condition = factor(decoder$group.ID));
##D #Read the QoRTs counts.
##D jscs = readJunctionSeqCounts(countfiles = countFiles,
##D            samplenames = decoder$sample.ID,
##D            design = design,
##D            flat.gff.file = gff.file
##D );
##D #Generate the size factors and load them into the JunctionSeqCountSet:
##D jscs <- estimateJunctionSeqSizeFactors(jscs);
##D #Estimate feature-specific dispersions:
##D jscs <- estimateJunctionSeqDispersions(jscs);
##D #Fit dispersion function and estimate MAP dispersion:
##D jscs <- fitJunctionSeqDispersionFunction(jscs);
##D #Test for differential usage:
##D jscs <- testForDiffUsage(jscs);
##D #Estimate effect sizes and expression estimates:
##D jscs <- estimateEffectSizes( jscs);
##D 
## End(Not run)
[Package JunctionSeq version 1.5.4 Index]