sageseqr-report.Rmd

---
date: "`r date()`"
output: 
  html_document:
    toc: true
    toc_float: true
---

```{css, echo = FALSE}
/* Whole document: */
body{
  font-family: Lato;
  font-size: 12pt;
}
/* Headers */
h1,h2,h3,h4,h5,h6{
  font-size: 24pt;
}
```

```{r title, echo = FALSE}
project <- config::get("report")
```
---
title: `r project`
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, 
                      fig.align = 'left',
                      fig.width = 10, 
                      fig.height = 10)
library(sageseqr)
```
# Data


# Explore Metadata

```{r sample_summary, message = FALSE}
targets::tar_load(clean_md)
var <- config::get("x_var")
summary <- dplyr::group_by_at(clean_md, var)
summary <- dplyr::summarise(summary, count = dplyr::n())
knitr::kable(summary)
```

Visualize the distribution of data.
```{r boxplot}
targets::tar_read(boxplots)
```

```{r boxplot_sex_chromosomes, results = "asis", fig.dim = c(5,5)}
if (!is.null(config::get("sex check"))) {
  cat("Visualize gene expression of sex chromosomes using XIST 
      as a X-chromosome marker and UTY as a Y-chromosome marker.")
  targets::tar_read(sex_plot)$plot
}
```

```{r sex_chromosome_pca, results = "asis",  fig.dim = c(10,5)}
if (!is.null(config::get("sex check"))) {
   cat("Visualize gene expression across X and Y chromosomes by principal component analysis (PCA).")
   targets::tar_read(sex_plot)$plot
}
```

Sample swaps identified by discordant sex label, are `r glue::glue_collapse(targets::tar_read(sex_plot_pca)$discordant, ", ", last = " and ")`.

```{r outlier_list, results = "asis"}
if (!is.null(targets::tar_read(sex_plot_pca)$outliers)) {
  list <- glue::glue_collapse(targets::tar_read(sex_plot_pca)$outliers, ", ", last = " and ")
  glue::glue("Possible outliers, identified as 2 standard deviations from the mean of samples n split by {config::get('sex check')} are {list}.")
} else {
  cat("No outliers identified.")
}
```

Visualize the relationships between covariates.
```{r heatmap_covariates}
correlation_input <- targets::tar_read(
  correlation_plot
  )$plot
col2 <- grDevices::colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
                                 "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                                 "#4393C3", "#2166AC", "#053061")))
corrplot::corrplot(correlation_input, col = col2(200), tl.col = "#000000")
```

# Filter Genes

Remove genes that have less than `r config::get("cpm threshold")` counts per million (CPM) in at least `r config::get("percent threshold")` of samples per specified condition.

`r dim(targets::tar_read(filtered_counts))[1]` genes are used in this analysis.
```{r filter_genes}
knitr::kable(targets::tar_read(biotypes))
```

Check distribution of correlation between genes.
```{r histogram}
targets::tar_read(gene_coexpression)
```

# Identify Outliers

```{r plot_outliers}
targets::tar_read(outliers)$plot
```
Outliers, based on logCPM expression, are `r glue::glue_collapse(targets::tar_read(outliers)$outliers, ", ", last = " and ")`.

# Identify Dropped Gene Features

Gene features dropped, because they are not annotated in BioMart, are `r glue::glue_collapse(targets::tar_read(dropped), ", ", last = " and ")`.

# Significant Covariates

Significant covariates are identified by the pearson correlation (p-value of 1%) between principal component analysis (PCA) of normalized transcripts and variables that meet a 0.1 false discovery rate (FDR) threshold. Significant covariates to adjust for are `r glue::glue_collapse(targets::tar_read(significant_covariates_plot)$significant_covariates, ", ", last = " and ")`.

```{r pca_and_significant_covariates}
targets::tar_read(
  significant_covariates_plot
  )$pc_results
```

# Model Identification
Covariates are added as fixed and random effects iteratively if model improvement by  Bayesian Information Criteria (BIC) was observed.
```{r model, results = "asis"}
if (!isTRUE(config::get("skip model"))) {
  summary <- targets::tar_read(model)
  glue::glue('The best model is: {glue::glue_collapse(as.character(summary$formula), " ")}')
} else {
  glue::glue(
    
    '\n\nStepwise regression model not quantified.\n\n'
    
    )
}
```
\n
```{r model_description}
if (!isTRUE(config::get("skip model"))) {
  knitr::kable(summary$to_visualize)
}
```
# Differential Expression

```{r differential_expression}
plots <- targets::tar_read(plot_de_volcano)
for(p in plots) {
  print(p)
}
```