dmDStest-class {DRIMSeq} | R Documentation |
dmDStest extends the dmDSfit
class by adding the null
model Dirichlet-multinomial (DM) and beta-binomial (BB) likelihoods and the
gene-level and feature-level results of testing for differential
exon/transcript usage. Result of calling the dmTest
function.
## S4 method for signature 'dmDStest' design(object, type = "null_model") results(x, ...) ## S4 method for signature 'dmDStest' results(x, level = "gene")
type |
Character indicating which design matrix should be returned.
Possible values |
x, object |
dmDStest object. |
... |
Other parameters that can be defined by methods using this generic. |
level |
Character specifying which type of results to return. Possible
values |
results(x)
: get a data frame with gene-level or
feature-level results.
design_fit_null
Numeric matrix of the design used to fit the null model.
lik_null
Numeric vector of the per gene DM null model likelihoods.
lik_null_bb
Numeric vector of the per gene BB null model likelihoods.
results_gene
Data frame with the gene-level results including:
gene_id
- gene IDs, lr
- likelihood ratio statistics based on
the DM model, df
- degrees of freedom, pvalue
- p-values and
adj_pvalue
- Benjamini & Hochberg adjusted p-values.
results_feature
Data frame with the feature-level results including:
gene_id
- gene IDs, feature_id
- feature IDs, lr
-
likelihood ratio statistics based on the BB model, df
- degrees of
freedom, pvalue
- p-values and adj_pvalue
- Benjamini &
Hochberg adjusted p-values.
Malgorzata Nowicka
dmDSdata
, dmDSprecision
,
dmDSfit
# -------------------------------------------------------------------------- # Create dmDSdata object # -------------------------------------------------------------------------- ## Get kallisto transcript counts from the 'PasillaTranscriptExpr' package library(PasillaTranscriptExpr) data_dir <- system.file("extdata", package = "PasillaTranscriptExpr") ## Load metadata pasilla_metadata <- read.table(file.path(data_dir, "metadata.txt"), header = TRUE, as.is = TRUE) ## Load counts pasilla_counts <- read.table(file.path(data_dir, "counts.txt"), header = TRUE, as.is = TRUE) ## Create a pasilla_samples data frame pasilla_samples <- data.frame(sample_id = pasilla_metadata$SampleName, group = pasilla_metadata$condition) levels(pasilla_samples$group) ## Create a dmDSdata object d <- dmDSdata(counts = pasilla_counts, samples = pasilla_samples) ## Use a subset of genes, which is defined in the following file gene_id_subset <- readLines(file.path(data_dir, "gene_id_subset.txt")) d <- d[names(d) %in% gene_id_subset, ] # -------------------------------------------------------------------------- # Differential transcript usage analysis - simple two group comparison # -------------------------------------------------------------------------- ## Filtering ## Check what is the minimal number of replicates per condition table(samples(d)$group) d <- dmFilter(d, min_samps_gene_expr = 7, min_samps_feature_expr = 3, min_gene_expr = 10, min_feature_expr = 10) plotData(d) ## Create the design matrix design_full <- model.matrix(~ group, data = samples(d)) ## To make the analysis reproducible set.seed(123) ## Calculate precision d <- dmPrecision(d, design = design_full) plotPrecision(d) head(mean_expression(d)) common_precision(d) head(genewise_precision(d)) ## Fit full model proportions d <- dmFit(d, design = design_full) ## Get fitted proportions head(proportions(d)) ## Get the DM regression coefficients (gene-level) head(coefficients(d)) ## Get the BB regression coefficients (feature-level) head(coefficients(d), level = "feature") ## Fit null model proportions and perform the LR test to detect DTU d <- dmTest(d, coef = "groupKD") ## Plot the gene-level p-values plotPValues(d) ## Get the gene-level results head(results(d)) ## Plot feature proportions for a top DTU gene res <- results(d) res <- res[order(res$pvalue, decreasing = FALSE), ] top_gene_id <- res$gene_id[1] plotProportions(d, gene_id = top_gene_id, group_variable = "group") plotProportions(d, gene_id = top_gene_id, group_variable = "group", plot_type = "lineplot") plotProportions(d, gene_id = top_gene_id, group_variable = "group", plot_type = "ribbonplot")