diff --git a/scripts/DRIMSeq_dtu.Rmd b/scripts/DRIMSeq_dtu.Rmd index 480d90ac9..b665b7ab6 100755 --- a/scripts/DRIMSeq_dtu.Rmd +++ b/scripts/DRIMSeq_dtu.Rmd @@ -12,7 +12,7 @@ output: keep_md: true --- -```{r setup, include=FALSE} +```{r DRIMSeq-setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` @@ -23,7 +23,7 @@ on abundance estimates from Salmon. It supports testing one or more contrasts. # Load packages -```{r} +```{r DRIMSeq-load-pkg} suppressPackageStartupMessages({ library(dplyr) library(tximport) @@ -39,14 +39,14 @@ suppressPackageStartupMessages({ We will use the transcript-level quantifications. -```{r} +```{r DRIMSeq-print-se} sg <- se$sg st <- se$st ``` # Create dmDSdata object -```{r dmDSdata} +```{r DRIMSeq-dmDSdata} counts <- data.frame(feature_id = rownames(st), gene_id = unlist(rowData(st)$gene_id), assays(st)[["counts"]], @@ -62,7 +62,7 @@ plotData(d) # Filter ************** MODIFY ************** -```{r data-filter} +```{r DRIMSeq-data-filter} d <- dmFilter(d, min_samps_gene_expr = 3, min_samps_feature_expr = 3, min_gene_expr = 10, min_feature_expr = 5) plotData(d) @@ -70,14 +70,14 @@ plotData(d) # Define design. ************** MODIFY ************** -```{r define-design} +```{r DRIMSeq-define-design} # (des <- model.matrix(~ XXXX, data = samples(d))) (des <- model.matrix(~ 0 + group, data = samples(d))) ``` # Calculate precision -```{r calculate-precision} +```{r DRIMSeq-calculate-precision} set.seed(123) d <- dmPrecision(d, design = des) plotPrecision(d) @@ -85,20 +85,20 @@ plotPrecision(d) # Fit model -```{r fit-model} +```{r DRIMSeq-fit-model} d <- dmFit(d, design = des, verbose = 1) ``` # Define contrasts. ************** MODIFY ************** -```{r define-contrasts} +```{r DRIMSeq-define-contrasts} #(contrasts <- as.data.frame(makeContrasts(XXXX, levels = des))) (contrasts <- as.data.frame(makeContrasts("groupN61311-groupN052611", levels = des))) ``` # Perform tests ************** MODIFY ************** -```{r result-genes} +```{r DRIMSeq-result-genes} level <- "gene" ## set to "feature" if transcript-level results are desired signif3 <- function(x) signif(x, digits = 3) DRIMSeq_fits <- lapply(contrasts, function(cm) { @@ -113,7 +113,7 @@ DRIMSeq_res <- lapply(DRIMSeq_fits, function(dr) { ``` -```{r result-transcripts} +```{r DRIMSeq-result-transcripts} level <- "feature" ## set to "feature" if transcript-level results are desired DRIMSeq_feature_fits <- lapply(contrasts, function(cm) { dr <- dmTest(d, contrast = cm, verbose = 1) @@ -128,7 +128,7 @@ DRIMSeq_feature_res <- lapply(DRIMSeq_feature_fits, function(dr) { # Write results to text files -```{r save-result} +```{r DRIMSeq-save-result} if (class(DRIMSeq_res) == "data.frame") { write.table(DRIMSeq_res %>% dplyr::arrange(pvalue), file = "DRIMSeq_dtu_results.txt", @@ -147,7 +147,7 @@ if (class(DRIMSeq_res) == "data.frame") { The `SummarizedExperiment` object on the transcript level -```{r se-gene} +```{r DRIMSeq-se-gene} ## add rows (NA) for genes that are filtered out (if any) DRIMSeq_resA <- lapply(seq_along(DRIMSeq_res), FUN = function(x) { @@ -217,7 +217,7 @@ for (i in seq_along(DRIMSeq_resA)) { The `SummarizedExperiment` object on the transcript level -```{r se-tx} +```{r DRIMSeq-se-tx} ## add rows (NA) for genes that are filtered out (if any) DRIMSeq_resB <- lapply(seq_along(DRIMSeq_feature_res), FUN = function(x) { @@ -289,14 +289,14 @@ for (i in seq_along(DRIMSeq_resB)) { The output is saved as a list. -```{r save-se} +```{r DRIMSeq-save-se} analysis_se <- list(sg = sg, st = st) saveRDS(analysis_se, file = "DRIMSeq_dtu.rds") ``` # Session info -```{r} +```{r DRIMSeq-session-info} date() sessionInfo() ``` diff --git a/scripts/edgeR_dge.Rmd b/scripts/edgeR_dge.Rmd index dad9dd8da..bb1658e59 100755 --- a/scripts/edgeR_dge.Rmd +++ b/scripts/edgeR_dge.Rmd @@ -12,7 +12,7 @@ output: keep_md: true --- -```{r setup, include=FALSE} +```{r edgeR-setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, dev = c("png", "pdf")) ``` @@ -23,7 +23,7 @@ abundance estimates from Salmon. # Load packages -```{r load-pkg} +```{r edgeR-load-pkg} suppressPackageStartupMessages({ library(dplyr) library(tximport) @@ -40,7 +40,7 @@ We load the `SummarizedExperiment` objects prepared using `tximeta`, containing gene- and transcript-level counts and feature lengths. In this report, we will use the gene-level quantifications. -```{r print-se} +```{r edgeR-print-se} ## List of SummarizedExperiment objects (gene/transcript level) se @@ -53,7 +53,7 @@ sg # Create DGEList and include average transcript length offsets -```{r dge-generate} +```{r edgeR-dge-generate} ## Extract count matrix and average transcript lengths cts <- assay(sg, "counts") normMat <- assay(sg, "length") @@ -74,7 +74,7 @@ dge0$genes <- as.data.frame(rowRanges(sg)) # Calculate logCPMs and add as an assay -```{r add-logcpm} +```{r edgeR-add-logcpm} logcpms <- edgeR::cpm(dge0, log = TRUE, prior.count = 2) stopifnot(all(rownames(logcpms) == rownames(sg)), all(colnames(logcpms) == colnames(sg))) @@ -83,7 +83,7 @@ assay(sg, "logcpm") <- logcpms # Define design. ************** MODIFY ************** -```{r define-design} +```{r edgeR-define-design} stopifnot(all(colnames(dge0) == metadata$names)) #(des <- model.matrix(~ XXXX, data = metadata)) #e.g. @@ -92,7 +92,7 @@ stopifnot(all(colnames(dge0) == metadata$names)) # Filter out lowly expressed genes -```{r filter-genes} +```{r edgeR-filter-genes} dim(dge0) keep <- edgeR::filterByExpr(dge0, design = des) dge <- dge0[keep, ] @@ -102,7 +102,7 @@ dim(dge) # Estimate dispersion -```{r estimate-disp} +```{r edgeR-estimate-disp} ## Estimate dispersion and fit model dge <- estimateDisp(dge, design = des) qlfit <- glmQLFit(dge, design = des) @@ -113,7 +113,7 @@ plotBCV(dge) # Define contrasts ************** MODIFY ************** -```{r define-contrasts} +```{r edgeR-define-contrasts} #(contrasts <- as.data.frame(makeContrasts(XXXX, levels = des))) #e.g. (contrasts <- as.data.frame(makeContrasts( @@ -123,7 +123,7 @@ plotBCV(dge) # Perform tests -```{r perform-tests} +```{r edgeR-perform-tests} signif3 <- function(x) signif(x, digits = 3) edgeR_res <- lapply(contrasts, function(cm) { qlf <- glmQLFTest(qlfit, contrast = cm) @@ -134,7 +134,7 @@ edgeR_res <- lapply(contrasts, function(cm) { # Make MA plots -```{r} +```{r edgeR-ma-plots} if (is(edgeR_res, "data.frame")) { print(ggplot(edgeR_res, aes(x = logCPM, y = logFC, color = FDR <= 0.05)) + geom_point() + theme_bw() + @@ -151,7 +151,7 @@ if (is(edgeR_res, "data.frame")) { # Write results to text files -```{r save-results} +```{r edgeR-save-results} ## Write results to text files and make MA plots if (is(edgeR_res, "data.frame")) { write.table(edgeR_res %>% dplyr::arrange(PValue), @@ -168,7 +168,7 @@ if (is(edgeR_res, "data.frame")) { # Output results as list of `SummarizedExperiment` objects -```{r se} +```{r edgeR-se} ## add rows (NA) for genes that are filtered out (if any) edgeR_resA <- lapply(seq_along(edgeR_res), FUN = function(x) { @@ -227,14 +227,14 @@ for (i in seq_along(edgeR_resA)) { The output is saved as a list. -```{r save-se} +```{r edgeR-save-se} analysis_se <- list(sg = sg, st = se$st) saveRDS(analysis_se, file = "edgeR_dge.rds") ``` # Session info -```{r session-info} +```{r edgeR-session-info} date() sessionInfo() ``` diff --git a/scripts/prepare_shiny.Rmd b/scripts/prepare_shiny.Rmd index 34bae8557..4931797a9 100755 --- a/scripts/prepare_shiny.Rmd +++ b/scripts/prepare_shiny.Rmd @@ -14,7 +14,7 @@ editor_options: chunk_output_type: console --- -```{r setup, include=FALSE} +```{r shiny-setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, dev = c("png", "pdf"), root.dir = "../") #knitr::opts_knit$set(root.dir = "../") ``` @@ -25,7 +25,7 @@ This script prepares the data to be used for shiny app. # Load packages -```{r load-pkg} +```{r shiny-load-pkg} suppressPackageStartupMessages({ library(dplyr) library(tidyr) @@ -40,7 +40,7 @@ suppressPackageStartupMessages({ # Load data -```{r load-data} +```{r shiny-load-data} options(ucscChromosomeNames = FALSE) sg <- se$sg st <- se$st @@ -48,7 +48,7 @@ st <- se$st # Gene models -```{r gene-model} +```{r shiny-gene-model} create_genemodels <- function(genemodels) { idx <- match(c("transcript_id", "gene_id", "exon_id"), colnames(mcols(genemodels))) @@ -66,7 +66,7 @@ if (!is.null(gtffile)) { # Vector with bigWig file names and condition information -```{r bigwig} +```{r shiny-bigwig} metadata <- colData(sg) if (!is.null(bigwigdir)) { bwfiles <- normalizePath(list.files(bigwigdir, pattern = "\\.bw$", @@ -85,7 +85,7 @@ if (!is.null(bigwigdir)) { # edgeR - gene-level MDS -```{r MDS} +```{r shiny-MDS} logcpms <- assay(sg, "logcpm") mds <- limma::plotMDS(logcpms, top = 500, labels = NULL, pch = NULL, cex = 1, dim.plot = c(1, 2), ndim = min(7, ncol(logcpms) - 1), @@ -108,7 +108,7 @@ and condition information The multidimensional scale data is stored in `reducedDims`. -```{r sce_gene} +```{r shiny-sce-gene} ## column data nam <- colData(sg)$names colData(sg)$bwFiles <- bwfiles[nam] @@ -137,7 +137,7 @@ and condition information. The multidimensional scale data is stored in `reducedDims`. -```{r sce_tx} +```{r shiny-sce-tx} ## column data nam <- colData(st)$names colData(st)$bwFiles <- bwfiles[nam] @@ -151,7 +151,7 @@ sce_tx <- SingleCellExperiment(assays = assays(st), # Output results -```{r output} +```{r shiny-save-sce} saveRDS(list(sce_tx = sce_tx, sce_gene = sce_gene), file = "shiny_sce.rds") @@ -159,7 +159,7 @@ saveRDS(list(sce_tx = sce_tx, # Session info -```{r session-info} +```{r shiny-session-info} date() sessionInfo() ```