Skip to content

Commit

Permalink
Consistent naming of Rmd code chunks
Browse files Browse the repository at this point in the history
  • Loading branch information
csoneson committed Jan 7, 2019
1 parent f210dd2 commit fe3430b
Show file tree
Hide file tree
Showing 3 changed files with 41 additions and 41 deletions.
32 changes: 16 additions & 16 deletions scripts/DRIMSeq_dtu.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ output:
keep_md: true
---

```{r setup, include=FALSE}
```{r DRIMSeq-setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Expand All @@ -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)
Expand All @@ -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"]],
Expand All @@ -62,43 +62,43 @@ 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)
```

# 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)
```

# 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) {
Expand All @@ -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)
Expand All @@ -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",
Expand All @@ -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) {
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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()
```
Expand Down
30 changes: 15 additions & 15 deletions scripts/edgeR_dge.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
```

Expand All @@ -23,7 +23,7 @@ abundance estimates from Salmon.

# Load packages

```{r load-pkg}
```{r edgeR-load-pkg}
suppressPackageStartupMessages({
library(dplyr)
library(tximport)
Expand All @@ -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
Expand All @@ -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")
Expand All @@ -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)))
Expand All @@ -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.
Expand All @@ -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, ]
Expand All @@ -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)
Expand All @@ -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(
Expand All @@ -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)
Expand All @@ -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() +
Expand All @@ -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),
Expand All @@ -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) {
Expand Down Expand Up @@ -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()
```
Expand Down
20 changes: 10 additions & 10 deletions scripts/prepare_shiny.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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 = "../")
```
Expand All @@ -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)
Expand All @@ -40,15 +40,15 @@ suppressPackageStartupMessages({

# Load data

```{r load-data}
```{r shiny-load-data}
options(ucscChromosomeNames = FALSE)
sg <- se$sg
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)))
Expand All @@ -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$",
Expand All @@ -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),
Expand All @@ -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]
Expand Down Expand Up @@ -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]
Expand All @@ -151,15 +151,15 @@ 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")
```

# Session info

```{r session-info}
```{r shiny-session-info}
date()
sessionInfo()
```
Expand Down

0 comments on commit fe3430b

Please sign in to comment.