Updating final report

assets/_quarto.yml

@@ -0,0 +1,12 @@
+format:
+  html:
+    code-fold: true
+    code-line-numbers: true
+    code-link: true
+    theme:
+      light:
+      - flatly
+      - styles.scss
+    toc: true
+    toc-depth: 3
+    toc-location: left

assets/analysis_de.qmd

@@ -42,6 +42,10 @@ cohort_id <- gsub("^analysis_de\\/(.*)\\.de\\.full\\.csv$", "\\1", de_table_file
 # Get all comparisons
 all_comparisons <- unique(de_table$comparison)
 
+# Get comparisons CSV
+comparisons_csv <- "samplesheets/comparisons.csv"
+comparisons_df <- read_csv(comparisons_csv)
+
 # Order and reduce to necessary columns
 de_table_col_order <- c(
   "Gene",
@@ -79,24 +83,23 @@ p_dist_plot_file <- paste0("analysis_de/results/", cohort_id, ".p_val_dist.png")
 log_fc_dist_plot_file <- paste0("analysis_de/results/", cohort_id, ".log_fc_dist.png")
 ```
 
-#### Dataset name: `{r} cohort_id`
-
-#### All comparisons made:
+This section shows the differential expression (DE) results for the comparisons made for your dataset. Each DE analysis is a simple case vs. control comparison of two [pseudobulk groups](#pseudobulked-groups-availablle-for-comparison) from your dataset. The following table lists the DE comparisons that were made:
 
 ```{r}
-#| echo: false
-#| results: asis
-for (cmp in all_comparisons) {
-  cat("- ")
-  cat(cmp)
-  cat("\n")
-}
+datatable(comparisons_df, colnames = c("Reference group", "Test group"))
 ```
 
-### Differential expression results
+The results are presented in four sections:
+
+- [p-value distribution results](#p-value-distributions)
+- [Log-fold change distribution results](#log-fold-change-distributions)
+- [Volcano plots](#volcano-plots)
+- [Differential expression results table](#differential-expression-table)
 
 #### p-value distributions
 
+The following plots show the distribution of p-values for each comparison made. These can be useful to determine the quality of the results: a well-powered experiment with true positive results should show an enrichment of low p-values (i.e. towards the left side of each plot). Uniformly distributed p-values indicate no true positive results or low power to detect true positives.
+
 ```{r}
 #| results: asis
 #| echo: false
@@ -105,6 +108,10 @@ cat(paste0('<img src="', p_dist_plot_file, '" alt="p-value distributions: ', coh
 
 #### Log-fold-change distributions
 
+Log-fold-change distributions are also useful in determining **how much** your significant genes are differentially expressed by. In the plots below, the red histograms represent the distribution of non-significant ("ns") genes, while the blue histograms represent significant genes ("sig").
+
+These plots can be used to determine whether an additional fold-change threshold should be applied to the data. Such a threshold can help to remove noise from very lowly-perturbed genes that may not represent a significant biological change. However, **keep in mind** that in many biological contexts, small fold-changes can be biologically significant, so apply a fold-change threshold with care.
+
 ```{r}
 #| results: asis
 #| echo: false
@@ -113,6 +120,8 @@ cat(paste0('<img src="', log_fc_dist_plot_file, '" alt="log-fold-change distribu
 
 #### Volcano plots
 
+Volcano plots combine both the p-value and fold-change information into one very useful plot. At a glance, you can identify the most significant results (higher up on the vertical axis, and highlighted in red) and the most differentially expressed genes (further away from the centre). If a fold-change threshold was specified, vertical lines will be drawn at the relevant fold-change positions on the plot and any genes between them will be greyed out and considered *non-significant*. The plots below additionally highlight the top-most significant genes for your convenience.
+
 ```{r}
 #| results: asis
 #| echo: false
@@ -121,6 +130,8 @@ cat(paste0('<img src="', volcano_plot_file, '" alt="volcano plots: ', cohort_id,
 
 #### Differential expression table
 
+Finally, below is a table summarising all of the differential expression results performed. The "Significant results only" tab shows just those genes that pass the p-value and fold-change thresholds provided. The "All results" tab shows all genes from all comparisons, regardless of significance.
+
 :::panel-tabset
 
 ##### Significant results only

assets/analysis_gsea.qmd

@@ -78,27 +78,19 @@ all_comparisons_str <- unique_comparisons %>%
   .$comparison
 names(all_comparisons_str) <- all_comparisons
 
+# Get comparisons CSV
+comparisons_csv <- "samplesheets/comparisons.csv"
+comparisons_df <- read_csv(comparisons_csv)
+
 # Get GSEA plots
 gsea_plots <- sapply(all_comparisons_str, function(cmp) {
   paste0("analysis_gsea/plots/", cohort_id, ".gsea.", cmp, ".png")
 })
 ```
 
-#### Dataset name: `{r} cohort_id`
-
-#### All comparisons made:
-
-```{r}
-#| echo: false
-#| results: asis
-for (cmp in all_comparisons) {
-  cat("- ")
-  cat(cmp)
-  cat("\n")
-}
-```
+This section presents the results from the GSEA analysis conducted on your dataset. GSEA works by first ranking the entire set of differential expression results by their fold-change values; this includes **both** significant **and** non-significant genes. The method then looks at a database of **biological pathways** and the **gene sets** associated with each. For each gene set, GSEA assesses whether these genes appear clustered together within your ranked dataset, more so than would be expected by chance. The degree of clustering towards the upregulated or downregulated ends of your dataset is used to calculate an **enrichment score** for that gene set.
 
-### GSEA expression results
+One GSEA analysis was conducted for [each differential expression comparison](#differential-expression-analysis) that was made. Below, you will find plots showing the enrichment score for each **significant** pathway from each analysis: more upregulated pathways have bars pointing to the right and are highlighted in orange, while downregulated pathways point to the left and are highlighted in blue.
 
 :::panel-tabset
 
@@ -108,11 +100,11 @@ purrr::walk(
   sort(all_comparisons),
   function(n) {
     p <- gsea_plots[[n]]
-    cat(paste("####", n, "\n\n"))
-    cat(paste("#####", "Gene Set Enrichment Analysis Results:", n))
+    cat(paste("#####", n, "\n\n"))
+    cat(paste("######", "Gene Set Enrichment Analysis Results:", n))
     cat(paste0('<img src="', p, '" alt="gsea results: ', n, '" width="100%">\n\n'))
   }
 )
 ```
 
-:::
+:::

assets/analysis_ora.qmd

@@ -78,27 +78,19 @@ all_comparisons_str <- unique_comparisons %>%
   .$comparison
 names(all_comparisons_str) <- all_comparisons
 
+# Get comparisons CSV
+comparisons_csv <- "samplesheets/comparisons.csv"
+comparisons_df <- read_csv(comparisons_csv)
+
 # Get ORA plots
 ora_plots <- sapply(all_comparisons_str, function(cmp) {
   paste0("analysis_ora/plots/", cohort_id, ".ora.", cmp, ".png")
 })
 ```
 
-#### Dataset name: `{r} cohort_id`
-
-#### All comparisons made:
-
-```{r}
-#| echo: false
-#| results: asis
-for (cmp in all_comparisons) {
-  cat("- ")
-  cat(cmp)
-  cat("\n")
-}
-```
+This section presents the results from the ORA analysis conducted on your dataset. ORA works by taking the set of your **significantly differentially expressed** genes and comparing that against each set of genes mapped to **biological pathways** in a curated database. For each pathway, ORA looks at the proportion of those genes that appear in your list of significant genes. If that proportion is higher than what would be expected by chance, the gene set is considered **over-represented** or **enriched** in your dataset. In contrast, some pathways may have very low representation in your dataset, in which case they can be considered **under-represented** or **depleted**. The ORA method generates a list of **enrichment ratios** for each pathway.
 
-### ORA expression results
+One ORA analysis was conducted for [each differential expression comparison](#differential-expression-analysis) that was made. Below, you will find plots showing the log2-transformed *enrichment ratio* for each **significant** pathway from the analysis: more enriched pathways have bars pointing to the right and are highlighted in orange, while depleted or under-represented pathways point to the left and are highlighted in blue.
 
 :::panel-tabset
 
@@ -108,8 +100,8 @@ purrr::walk(
   sort(all_comparisons),
   function(n) {
     p <- ora_plots[[n]]
-    cat(paste("####", n, "\n\n"))
-    cat(paste("#####", "Over-Representation Analysis Results:", n))
+    cat(paste("#####", n, "\n\n"))
+    cat(paste("######", "Over-Representation Analysis Results:", n))
     cat(paste0('<img src="', p, '" alt="ora results: ', n, '" width="100%">\n\n'))
   }
 )

assets/analysis_pseudo.qmd

@@ -37,8 +37,13 @@ pseudo_groups_file <- "pseudo_groups.txt"
 pseudo_groups <- scan(pseudo_groups_file, character())
 ```
 
-#### Dataset name: `{r} cohort_id`
-#### Fields used for pseudobulking:
+With the dataset filtered, integrated, and annotated, the next step is to perform pseudobulking. This is the process of aggregating the RNA counts across many cells within a sample into a single data point. For example, we may have a sample with 5000 B cells and 4000 T cells. During pseudobulking, we would aggregate (i.e. sum together) the expression profiles of all B cells and all T cells, leaving us with one pseudobulked B cell expression profile and one pseudobulked T cell expression profile for that sample. We would then repeat this procedure for every sample.
+
+The process of pseudobulking effectively performs a bulk RNA-seq experiment for each cell type - hence the name "pseudobulk". But, instead of relying on cell sorting methods to isolate each cell type from one another, we can use our single cell resolution data and cell annotations. The main advantage of pseudobulking is that we can use well-established statistical techniques for assessing differentially expressed genes that were developed for bulk RNA-seq and are more robust than newer, more complex scRNA-seq techniques. See the [single-cell best practices guide](https://www.sc-best-practices.org/conditions/differential_gene_expression.html) for more information on the benefits of pseudobulking.
+
+#### Fields used for pseudobulking
+
+The following [metadata fields](#all-annotation-fields-available-for-downstream-analysis) were used for aggregation of cells:
 
 ```{r}
 #| echo: false
@@ -52,6 +57,8 @@ for (g in pseudo_groups) {
 
 ### Pseudobulked groups availablle for comparison
 
+Below is a table summarising all of your final pseudobulk groups and the number of samples in your dataset that are represented in each group. These pseudobulk group names are now available for comparison during differential expression analysis. **Note** that for differential expression analysis, you will ideally have *at least three* samples for each of your case and control groups, if not more.
+
 ```{r}
 #| echo: false
 #| warning: false

assets/annotation.qmd

@@ -27,6 +27,19 @@ available_annotation_files <- dir(
   full.names = TRUE
 )
 available_annotations <- unlist(sapply(available_annotation_files, scan, character()))
+annotation_priority <- c(
+  "custom_cell_type",
+  "custom_cell_type.max_score",
+  "SingleR.annotation",
+  "SingleR.hpca_main",
+  "SingleR.hpca_fine",
+  "Phase"
+)
+annotation_priority <- annotation_priority[annotation_priority %in% available_annotations]
+annotation_used_for_clusters <- scan("cluster_annotation.txt", character())
+if (length(annotation_used_for_clusters) == 0) {
+  annotation_used_for_clusters <- annotation_priority[1]
+}
 ```
 
 ```{r}
@@ -138,9 +151,65 @@ if (is.null(cohort_id)) {
 }
 ```
 
-#### Dataset name: `{r} cohort_id`
+There are three main annotation options supported by the scRNAvigator pipeline:
+
+- Cell cycle annotation
+- Automatic cell type annotation
+- Custom marker gene-based annotation
+
+Cell cycle annotation is useful if you are studying rapidly-dividing cells such as tumor samples and wish to capture information about which stage of the cell cycle each cell is in.
+
+Automatic cell type annotation is performed by comparing your data against curated and pre-annotated public datasets and assigning annotations to your cells based on their similarities to the curated datasets. For human samples, the Human Primary Cell Atlas (HPCA) database is used for this purpose by default, while for non-human species, an external database needs to be provided in order to perform automatic cell type annotation.
+
+Custom marker gene-based annotation takes a list of cell types and associated marker genes and runs a scoring algorithm to determine the closest matching cell type for each cell. Custom cell type annotations creates two annotation fields in the sample metadata: `custom_cell_type` and `custom_cell_type.max_score`. The `custom_cell_type.max_score` annotation simply picks the custom cell type that scored the highest. The `custom_cell_type` annotation, however, also takes into account how much the highest score differs from the second-highest score; when these are too close together for a given cell, the cell will be labelled "Ambiguous". Of these two annotation fields, we recommend using `custom_cell_type` for downstream analysis due to it being more conservative.
+
+After annotation with one or more of the above methods, your cell **clusters** will be annotated as well. This takes one of the annotation categories from above and determines whether that cluster predominantly consists of a single cell type or label. If so, that label is applied to the whole cluster. If there is too much heterogeneity within a cluster, it is labelled "Ambiguous".
+
+#### All annotation fields available for cluster annotation
+
+Your cells were annotated with the following annotation fields:
+
+```{r}
+#| echo: false
+#| warning: false
+#| error: false
+#| message: false
+all_annotations_df <- tibble(
+  annotation_field = c(),
+  description = c()
+)
+annotation_descriptions <- list(
+  custom_cell_type = "Custom marker gene-based cell type annotation (highest score significantly higher than second-highest score)",
+  custom_cell_type.max_score = "Custom marker gene-based cell type annotation (highest score only)",
+  SingleR.annotation = "SingleR annotation with custom cell annotation database",
+  SingleR.hpca_main = "SingleR annotation with HPCA database ('main' cell type category)",
+  SingleR.hpca_fine = "SingleR annotation with HPCA database ('fine' cell type category)",
+  Phase = "Cell cycle phase"
+)
+for (ann in annotation_priority) {
+  all_annotations_df <- all_annotations_df %>%
+    bind_rows(tibble(
+      annotation_field = ann,
+      description = annotation_descriptions[[ann]]
+    ))
+}
+```
+
+```{r}
+#| echo: false
+#| eval: !expr (dim(all_annotations_df)[1] > 0)
+
+datatable(
+  all_annotations_df,
+  colnames = c("Annotation field", "Description")
+)
+```
+
+Cell **clusters** were annotated with the labels from the **`r annotation_used_for_clusters`** annotation field.
+
+#### All annotation fields available for downstream analysis
 
-### All annotation fields available for downstream analyses
+The following list contains all of the annotations listed above, as well as the `cluster_annotation` field and additional metadata fields that were provided in the samplesheet. These are all of the fields available for performing downstream pseudobulking and subsequent differential expression and functional enrichment analyses.
 
 ```{r}
 #| echo: false
@@ -158,26 +227,26 @@ for (ann in available_annotations) {
 #| results: asis
 #| echo: false
 #| eval: !expr cell_cycle_annotated
-cat("### Cell cycle annotation\n\n")
+cat("#### Cell cycle annotation\n\n")
 cat(paste0('<img src="', cc_umap_file, '" alt="cell cycle umap: ', cohort_id, '" width="100%">\n\n'))
 ```
 
 ```{r}
 #| results: asis
 #| echo: false
 #| eval: !expr db_annotated
-cat("### Automatic cell type annotation\n\n")
-cat("#### Cell type UMAP\n\n")
+cat("#### Automatic cell type annotation\n\n")
+cat("##### Cell type UMAP\n\n")
 cat(paste0('<img src="', db_umap_file, '" alt="cell type umap: ', cohort_id, '" width="100%">\n\n'))
-cat("#### Cell type proportions\n\n")
+cat("##### Cell type proportions\n\n")
 cat(paste0('<img src="', cell_type_props, '" alt="cell type proportions: ', cohort_id, '" width="100%">\n\n'))
 ```
 
 ```{r}
 #| results: asis
 #| echo: false
 #| eval: !expr db_annotated && file.exists(cell_type_prop_table_hpca_main)
-cat("#### HPCA Main annotations")
+cat("##### HPCA Main annotations")
 ```
 
 ```{r}
@@ -193,7 +262,7 @@ datatable(read_csv(cell_type_prop_table_hpca_main))
 #| results: asis
 #| echo: false
 #| eval: !expr db_annotated && file.exists(cell_type_prop_table_hpca_fine)
-cat("#### HPCA Fine annotations")
+cat("##### HPCA Fine annotations")
 ```
 
 ```{r}
@@ -209,7 +278,7 @@ datatable(read_csv(cell_type_prop_table_hpca_fine))
 #| results: asis
 #| echo: false
 #| eval: !expr db_annotated && file.exists(cell_type_prop_table_singler)
-cat("#### SingleR annotations")
+cat("##### SingleR annotations")
 ```
 
 ```{r}
@@ -225,10 +294,10 @@ datatable(read_csv(cell_type_prop_table_singler))
 #| results: asis
 #| echo: false
 #| eval: !expr custom_annotated
-cat("### Custom cell type annotations\n\n")
-cat("#### Cell type UMAP\n\n")
+cat("#### Custom cell type annotations\n\n")
+cat("##### Cell type UMAP\n\n")
 cat(paste0('<img src="', custom_umap, '" alt="cell type umap: ', cohort_id, '" width="100%">\n\n'))
-cat("#### Cell type proportions\n\n")
+cat("##### Cell type proportions\n\n")
 cat(paste0('<img src="', custom_cell_type_props, '" alt="cell type proportions: ', cohort_id, '" width="100%">\n\n'))
 ```
 
@@ -245,10 +314,10 @@ datatable(custom_cell_type_props_table)
 #| results: asis
 #| echo: false
 #| eval: !expr custom_annotated
-cat("#### Custom cell type feature plots\n\n")
+cat("##### Custom cell type feature plots\n\n")
 for (cell_type in names(custom_cell_type_feature_plots)) {
   p <- custom_cell_type_feature_plots[[cell_type]]
-  cat(paste0("##### ", cell_type, "\n\n"))
+  cat(paste0("###### ", cell_type, "\n\n"))
   cat(paste0('<img src="', p, '" alt="feature plot: ', cell_type, ' - ', cohort_id, '" width="100%">\n\n'))
 }
 ```
@@ -257,7 +326,7 @@ for (cell_type in names(custom_cell_type_feature_plots)) {
 #| results: asis
 #| echo: false
 #| eval: !expr clusters_annotated
-cat("### Cluster cell type annotation\n\n")
+cat("#### Cluster cell type annotation\n\n")
 cat(paste0('<img src="', cluster_umap_file, '" alt="cluster cell type assignment umap: ', cohort_id, '" width="100%">\n\n'))
 ```