Fit dispersion function to share dispersion information between
features across the genome.
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.
A JunctionSeqCountSet. Usually initially created by
readJunctionSeqCounts. Size factors must be
set, usually using functions estimateSizeFactors and
estimateJunctionSeqDispersions.
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
dispersion estimation methodology used by DESeq2 and DEXSeq v1.12.0 or higher to generate the initial (feature-specific) dispersion estimates. 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 DEXSeq v1.8.0 or less.
method.dispFit
Character string. Can be used to apply alternative methodologies or implementations. Intended for advanced users who have strong opinions about the underlying statistical methodologies.
Determines the method used to generated "fitted" dispersion estimates. One of "parametric" (the default), "local", or "mean". See the DESeq2 documentation for more information.
method.dispFinal
Character string. Can be used to apply alternative methodologies or implementations. Intended for advanced users who have strong opinions about the underlying statistical methodologies.
Determines the method used to arrive at a "final" dispersion estimate. The default, "shrink" uses the maximum a posteriori estimate, combining information from both
the fitted and feature-specific dispersion estimates. This is the method used by DESeq2 and DEXSeq v1.12.0 and above.
fitDispersionsForExonsAndJunctionsSeparately
When running a "junctionsAndExons" type analysis in which both exons and splice junctions are being tested simultaniously, this parameter
determines whether a single fitted dispersion model should be fitted for both exons and splice junctions, or if separate
fitted dispersions should be calculated for each. By default the dispersions are run separately.
verbose
if TRUE, send debugging and progress messages to the console / stdout.
...
If using the depreciated fitDispersionFunction command, use the same syntax as above.
Value
A JunctionSeqCountSet, with dispersion results included.
Examples
data(exampleDataSet,package="JctSeqData");
jscs <- fitJunctionSeqDispersionFunction(jscs);
## Not run:
#Full example (from scratch):
########################################
#Set up example data:
decoder.file <- system.file(
"extdata/annoFiles/decoder.bySample.txt",
package="JctSeqData");
decoder <- read.table(decoder.file,
header=TRUE,
stringsAsFactors=FALSE);
gff.file <- system.file(
"extdata/cts/withNovel.forJunctionSeq.gff.gz",
package="JctSeqData");
countFiles <- system.file(paste0("extdata/cts/",
decoder$sample.ID,
"/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz"),
package="JctSeqData");
########################################
#Advanced Analysis:
#Make a "design" dataframe:
design <- data.frame(condition = factor(decoder$group.ID));
#Read the QoRTs counts.
jscs = readJunctionSeqCounts(countfiles = countFiles,
samplenames = decoder$sample.ID,
design = design,
flat.gff.file = gff.file
);
#Generate the size factors and load them into the JunctionSeqCountSet:
jscs <- estimateJunctionSeqSizeFactors(jscs);
#Estimate feature-specific dispersions:
jscs <- estimateJunctionSeqDispersions(jscs);
#Fit dispersion function and estimate MAP dispersion:
jscs <- fitJunctionSeqDispersionFunction(jscs);
#Test for differential usage:
jscs <- testForDiffUsage(jscs);
#Estimate effect sizes and expression estimates:
jscs <- estimateEffectSizes( jscs);
## End(Not run)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(JunctionSeq)
Loading required package: SummarizedExperiment
Loading required package: GenomicRanges
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:parallel':
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from 'package:stats':
IQR, mad, xtabs
The following objects are masked from 'package:base':
Filter, Find, Map, Position, Reduce, anyDuplicated, append,
as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
rbind, rownames, sapply, setdiff, sort, table, tapply, union,
unique, unsplit
Loading required package: S4Vectors
Loading required package: stats4
Attaching package: 'S4Vectors'
The following objects are masked from 'package:base':
colMeans, colSums, expand.grid, rowMeans, rowSums
Loading required package: IRanges
Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/JunctionSeq/fitJunctionSeqDispersionFunction.Rd_%03d_medium.png", width=480, height=480)
> ### Name: fitJunctionSeqDispersionFunction
> ### Title: Fit Shared Dispersion Function
> ### Aliases: fitJunctionSeqDispersionFunction
>
> ### ** Examples
>
> data(exampleDataSet,package="JctSeqData");
> jscs <- fitJunctionSeqDispersionFunction(jscs);
> fitDispersionFunction() Starting (Tue Jul 5 23:18:00 2016)
> (fitType = parametric)
> (finalDispersionMethod = shrink)
> (fitDispersionsForExonsAndJunctionsSeparately = TRUE)
min(means[useForFit], na.rm=T)=0.146009972023612
> fdf: Fitting dispersions:
> (Iteration 1) Parametric Dispersion Coefs: [0.00472009572814427,0.218460947334769]
> (Iteration 2) Parametric Dispersion Coefs: [0.00256247314248661,0.234164421253589]
> (FINAL) Parametric Dispersion Coefs: [0.00256236453981216,0.234167798819113]
> fdf: Fitting dispersions of exons and junctions to separate fitted trends.
> fdf: Fitting exon dispersions:
> (Iteration 1) Parametric Dispersion Coefs: [0.00471170193076281,0.279913188264851]
> (Iteration 2) Parametric Dispersion Coefs: [0.00195226100860601,0.298767665275019]
> (Iteration 3) Parametric Dispersion Coefs: [0.00187840115925201,0.287865989524349]
> (FINAL) Parametric Dispersion Coefs: [0.00187811534730775,0.287879543718636]
> fdf: Fitting splice-junction dispersions:
> (Iteration 1) Parametric Dispersion Coefs: [0.00389781116922728,0.181506534119097]
> (Iteration 2) Parametric Dispersion Coefs: [0.00376412469359818,0.182903659298507]
> (FINAL) Parametric Dispersion Coefs: [0.0037642302236494,0.182901650968558]
> fdf(): 'Shrinking' fitted and feature-specific dispersion estimates.
> fdf() Dispersion estimate failed for 0 out of 482 features.
> fitDispersionFunction() Done. (Tue Jul 5 23:18:01 2016)
>
> ## Not run:
> ##D #Full example (from scratch):
> ##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)
>
>
>
>
>
> dev.off()
null device
1
>