## ----setup, echo=FALSE, results="hide"---------------------------------------- knitr::opts_chunk$set(tidy=FALSE, cache=TRUE, dev="png", message=TRUE, error=FALSE, warning=TRUE) ## ----Bioconductor_installation, eval=FALSE------------------------------------ # if (!requireNamespace("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # BiocManager::install("DifferentialRegulation") ## ----vignettes, eval=FALSE---------------------------------------------------- # browseVignettes("DifferentialRegulation") ## ----citation----------------------------------------------------------------- citation("DifferentialRegulation") ## ----load, message=FALSE------------------------------------------------------ library(DifferentialRegulation) ## ----specify_data-dir--------------------------------------------------------- data_dir = system.file("extdata", package = "DifferentialRegulation") ## ----specify_directories------------------------------------------------------ # specify samples ids: sample_ids = paste0("organoid", c(1:3, 16:18)) # set directories of each sample input data (obtained via alevin-fry): base_dir = file.path(data_dir, "alevin-fry", sample_ids) # Note that alevin-fry needs to be run with `--use-mtx` option to store counts in a `quants_mat.mtx` file. path_to_counts = file.path(base_dir,"/alevin/quants_mat.mtx") path_to_cell_id = file.path(base_dir,"/alevin/quants_mat_rows.txt") path_to_gene_id = file.path(base_dir,"/alevin/quants_mat_cols.txt") ## ----specify_directories_EC--------------------------------------------------- path_to_EC_counts = file.path(base_dir,"/alevin/geqc_counts.mtx") path_to_EC = file.path(base_dir,"/alevin/gene_eqclass.txt.gz") ## ----load_USA_counts---------------------------------------------------------- sce = load_USA(path_to_counts, path_to_cell_id, path_to_gene_id, sample_ids) sce ## ----cell-type---------------------------------------------------------------- path_to_DF = file.path(data_dir,"DF_cell_types.txt") DF_cell_types = read.csv(path_to_DF, sep = "\t", header = TRUE) matches = match(colnames(sce), DF_cell_types$cell_id) sce$cell_type = DF_cell_types$cell_type[matches] table(sce$cell_type) ## ----load_EC_counts----------------------------------------------------------- EC_list = load_EC(path_to_EC_counts, path_to_EC, path_to_cell_id, path_to_gene_id, sample_ids) ## ----samples_design----------------------------------------------------------- design = data.frame(sample = sample_ids, group = c( rep("3 mon", 3), rep("6 mon", 3) )) design ## ----compute_PB_counts-------------------------------------------------------- PB_counts = compute_PB_counts(sce = sce, EC_list = EC_list, design = design, sample_col_name = "sample", group_col_name = "group", sce_cluster_name = "cell_type") ## ----rm_EC_list--------------------------------------------------------------- rm(EC_list) ## ----EC-test------------------------------------------------------------------ # EC-based test: set.seed(1609612) results = DifferentialRegulation(PB_counts, n_cores = 2, traceplot = TRUE) ## ----names-results------------------------------------------------------------ names(results) ## ----visualize_gene_results--------------------------------------------------- head(results$Differential_results, 3) ## ----sort_gene_results-UP----------------------------------------------------- ordering_UP = order(results$Differential_results[,6], decreasing = TRUE) head(results$Differential_results[ordering_UP,], 3) ## ----sort_gene_results-DOWN--------------------------------------------------- ordering_DOWN = order(results$Differential_results[,6], decreasing = FALSE) head(results$Differential_results[ordering_DOWN,], 3) ## ----visualize_US_results----------------------------------------------------- head(results$US_results, 3) ## ----visualize_USA_results---------------------------------------------------- head(results$USA_results, 3) ## ----visualize_convergence_results-------------------------------------------- results$Convergence_results ## ----plot_pi------------------------------------------------------------------ plot_pi(results, type = "US", gene_id = results$Differential_results$Gene_id[1], cluster_id = results$Differential_results$Cluster_id[1]) plot_pi(results, type = "US", gene_id = results$Differential_results$Gene_id[2], cluster_id = results$Differential_results$Cluster_id[2]) ## ----plot_pi_USA-------------------------------------------------------------- plot_pi(results, type = "USA", gene_id = results$Differential_results$Gene_id[1], cluster_id = results$Differential_results$Cluster_id[1]) plot_pi(results, type = "USA", gene_id = results$Differential_results$Gene_id[2], cluster_id = results$Differential_results$Cluster_id[2]) ## ----plot_pi-traceplot-------------------------------------------------------- plot_traceplot(results, gene_id = results$Differential_results$Gene_id[1], cluster_id = results$Differential_results$Cluster_id[1]) plot_traceplot(results, gene_id = results$Differential_results$Gene_id[2], cluster_id = results$Differential_results$Cluster_id[2]) ## ----load-bulk, message=FALSE------------------------------------------------- library(DifferentialRegulation) ## ----specify_data-dir-bulk---------------------------------------------------- data_dir = system.file("extdata", package = "DifferentialRegulation") ## ----specify_directories-bulk------------------------------------------------- # specify samples ids: sample_ids = paste0("sample", seq_len(6)) # US estimates: quant_files = file.path(data_dir, "salmon", sample_ids, "quant.sf") file.exists(quant_files) # Equivalence classes: equiv_classes_files = file.path(data_dir, "salmon", sample_ids, "aux_info/eq_classes.txt.gz") file.exists(equiv_classes_files) ## ----load_US_counts-bulk------------------------------------------------------ SE = load_bulk_US(quant_files, sample_ids) ## ----load_EC_counts-bulk------------------------------------------------------ EC_list = load_bulk_EC(path_to_eq_classes = equiv_classes_files, n_cores = 2) ## ----samples_design-bulk------------------------------------------------------ group_names = rep(c("A", "B"), each = 3) design = data.frame(sample = sample_ids, group = group_names) design ## ----EC-test-bulk------------------------------------------------------------- # EC-based test: set.seed(1609612) results = DifferentialRegulation_bulk(SE = SE, EC_list = EC_list, design = design, n_cores = 2, traceplot = TRUE) ## ----names-results-bulk------------------------------------------------------- names(results) ## ----visualize_gene_results-bulk---------------------------------------------- head(results$Differential_results, 3) ## ----sort_gene_results-UP-bulk------------------------------------------------ ordering_UP = order(results$Differential_results[,4], decreasing = TRUE) head(results$Differential_results[ordering_UP,], 3) ## ----sort_gene_results-DOWN-bulk---------------------------------------------- ordering_DOWN = order(results$Differential_results[,4], decreasing = FALSE) head(results$Differential_results[ordering_DOWN,], 3) ## ----visualize_convergence_results-bulk--------------------------------------- results$Convergence_results ## ----plot_pi-bulk------------------------------------------------------------- plot_bulk_pi(results, transcript_id = results$Differential_results$Transcript_id[1]) plot_bulk_pi(results, transcript_id = results$Differential_results$Transcript_id[2]) ## ----plot_pi-bulk-traceplot--------------------------------------------------- plot_bulk_traceplot(results, transcript_id = results$Differential_results$Transcript_id[1]) plot_bulk_traceplot(results, transcript_id = results$Differential_results$Transcript_id[2]) ## ----sessionInfo-------------------------------------------------------------- sessionInfo()