prior_precision: Infer an informative prior for the precision

In SimoneTiberi/BANDITS: BANDITS: Bayesian ANalysis of DIfferenTial Splicing

Infer an informative prior for the precision

Description

prior_precision uses DRIMSeq's pipeline to infer an informative prior for the precision parameter of the Dirichlet-Multinomial distribution. The function computes the genewise estimates for the precision via DRIMSeq::dmPrecision, and calculates the mean and standard deviation of the log-precision estimates.

Usage

prior_precision(
  gene_to_transcript,
  transcript_counts,
  n_cores = 1,
  transcripts_to_keep = NULL,
  max_n_genes_used = 100
)

Arguments

gene_to_transcript	a matrix or data.frame with a list of gene-to-transcript correspondances. The first column represents the gene id, while the second one contains the transcript id.
transcript_counts	a matrix or data.frame, with 1 column per sample and 1 row per transcript, containing the estimated abundances for each transcript in each sample.
n_cores	the number of cores to parallelize the tasks on.
transcripts_to_keep	a vector containing the list of transcripts to keep. Ideally, created via filter_transcripts.
max_n_genes_used	the maximum number of genes to compute the prior on. First, genes with at least 2 transcripts are selected. Then, if more than 'max_n_genes_used' such genes are available, 'max_n_genes_used' of these genes are sampled at random and used to calculate the prior of the precision parameter. A smaller 'max_n_genes_used' (minimum 100) will lead to faster but more approximate prior estimates.

Value

A list with 2 objects containing:

• prior: a vector containing the mean and standard deviation of the log-precision, used to formulate an informative prior in test_DTU;

• genewise_log_precision: a numeric vector with the individual genewise estimates for the log-precision.

Author(s)

Simone Tiberi [e-mail]

See Also

test_DTU, plot_precision

Examples

# specify the directory of the internal data:
data_dir = system.file("extdata", package = "BANDITS")

# load gene_to_transcript matching:
data("gene_tr_id", package = "BANDITS")

# Load the transcript level estimated counts via tximport:
library(tximport)
quant_files = file.path(data_dir, "STAR-salmon", paste0("sample", seq_len(4)), "quant.sf")
txi = tximport(files = quant_files, type = "salmon", txOut = TRUE)
counts = txi$counts

# Optional (recommended): transcript pre-filtering
transcripts_to_keep = filter_transcripts(gene_to_transcript = gene_tr_id,
                                         transcript_counts = counts,
                                         min_transcript_proportion = 0.01,
                                         min_transcript_counts = 10,
                                         min_gene_counts = 20)

# Infer an informative prior for the precision parameter
# Use the same filtering criteria as in 'create_data', by choosing the same argument for 'transcripts_to_keep'.
# If transcript pre-filtering is not performed, leave 'transcripts_to_keep' unspecified.
set.seed(61217)
precision = prior_precision(gene_to_transcript = gene_tr_id, transcript_counts = counts,
                            n_cores = 2, transcripts_to_keep = transcripts_to_keep)
precision$prior
head(precision$genewise_log_precision)

SimoneTiberi/BANDITS documentation built on Nov. 15, 2023, 2:35 p.m.