From d97a2538835d1927a75f1be0d4bfc6549496eef7 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Mon, 30 Mar 2026 16:34:40 +1100 Subject: [PATCH 01/18] WIP: add csv files to final report inputs, start adding text to report --- assets/index.qmd | 19 +++++++++++++++++++ main.nf | 9 +++++++-- modules/report.nf | 3 +++ 3 files changed, 29 insertions(+), 2 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 66ead11..816f794 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -2,3 +2,22 @@ title: "scRNAvigator Report" format: html --- + +Thank you for using scRNAvigator! + +This website contains a report of your analysis. Depending on what sections of the pipeline were ran, it will include one or more of the following: + +- Quality control metrics and plots +- Doublet detection and filtering results +- Dataset integration results +- Cell annotation results +- Pseudobulking and differential expression results +- Functional enrichment analysis results + +## Cite us + +If you use these results in published work, please cite us with: + +``` +# TODO: add citation information +``` \ No newline at end of file diff --git a/main.nf b/main.nf index e224ae6..6975ba7 100644 --- a/main.nf +++ b/main.nf @@ -77,7 +77,7 @@ workflow { gsea_db = !!params.gsea_db ? channel.fromPath(params.gsea_db, checkIfExists: true).first() : channel.value([]) // Read in samplesheet - samplesheet = channel.fromPath(params.input) + samplesheet = channel.fromPath(params.input, checkIfExists: true) .splitCsv( header: true ) .map { row -> { def sample = row.sample @@ -408,6 +408,11 @@ workflow { .map { tmp -> [ tmp ] } report_style = channel.fromPath("${projectDir}/assets/styles.scss") + // Gather CSV inputs + samplesheet_csv = channel.fromPath(params.input, checkIfExists: true) + custom_marker_genes_csv = !!params.custom_marker_genes ? channel.fromPath(params.custom_marker_genes, checkIfExists: true) : channel.value([]) + comparisons_csv = !!params.comparisons ? channel.fromPath(params.comparisons, checkIfExists: true) : channel.value([]) + // Run the report module report_input = cohort_id .merge(report_metadata) @@ -427,5 +432,5 @@ workflow { .merge(analysis_ora_results) .merge(available_annotation_files) .merge(report_style) - REPORT(report_input) + REPORT(report_input, samplesheet_csv, custom_marker_genes_csv, comparisons_csv) } \ No newline at end of file diff --git a/modules/report.nf b/modules/report.nf index 29a4521..7c1b9ae 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -24,6 +24,9 @@ process REPORT { path(analysis_ora_results, stageAs: "analysis_ora/*"), path(available_annotation_files, stageAs: "available_annotations/??.txt"), path(report_style) + path samplesheet + path custom_marker_genes + path comparisons output: tuple val(cohort_name), path("report"), emit: report From 89b5cdd6361fe9faace8dc9256e10d9c33793001 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Tue, 31 Mar 2026 12:18:47 +1100 Subject: [PATCH 02/18] add samplesheet and modules run to index page of report --- assets/index.qmd | 140 ++++++++++++++++++++++++++++++++++++++++++++++ main.nf | 4 +- modules/report.nf | 7 ++- 3 files changed, 146 insertions(+), 5 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 816f794..4306590 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -14,6 +14,146 @@ This website contains a report of your analysis. Depending on what sections of t - Pseudobulking and differential expression results - Functional enrichment analysis results +```{r} +#| echo: false +#| results: hide +#| warning: false +#| error: false +#| message: false +library(tidyverse) +library(DT) + +# Get cohort name +cohort_name <- scan("cohort_name.txt", character())[1] + +# Read sample sheet +samplesheet <- read_csv("samplesheets/samplesheet.csv") + +# Get samplesheet column names +samplesheet_cols <- list( + sample = "Sample name", + rds = "File path", + min_ncount = "Min. UMI count", + max_ncount = "Max. UMI count", + min_nfeature = "Min. gene count", + max_nfeature = "Max. gene count", + min_mt_pct = "Min. MT gene proportion (%)", + max_mt_pct = "Max. MT gene proportion (%)", + cells_to_remove = "Cells to remove file path", + res = "Sample clustering resolution", + multiplet_rate = "Expected multiplet rate" +) + +samplesheet_cols_show <- list( + sample = FALSE, + rds = FALSE, + min_ncount = FALSE, + max_ncount = FALSE, + min_nfeature = FALSE, + max_nfeature = FALSE, + min_mt_pct = FALSE, + max_mt_pct = FALSE, + cells_to_remove = FALSE, + res = FALSE, + multiplet_rate = FALSE +) + +for (col in colnames(samplesheet)) { + samplesheet_cols_show[[col]] <- TRUE + + if (!(col %in% names(samplesheet_cols))) { + samplesheet_cols[[col]] <- col + } +} + +samplesheet_cols_show$rds <- FALSE +samplesheet_cols_show$cells_to_remove <- FALSE +samplesheet_cols_final <- samplesheet_cols[unlist(samplesheet_cols_show)] + +samplesheet_show <- samplesheet %>% + select(all_of(names(samplesheet_cols_final))) + +# Read custom marker genes if present +custom_marker_genes <- NA +custom_marker_genes_file <- "samplesheets/custom_marker_genes.csv" +if (file.exists(custom_marker_genes_file)) { + custom_marker_genes <- read_csv(custom_marker_genes_file) +} + +# Read comparisons if present +comparisons <- NA +comparisons_file <- "samplesheets/comparisons.csv" +if (file.exists(comparisons_file)) { + comparisons <- read_csv(comparisons_file) +} + +# Get list of pipeline sections that were run +modules_run <- list( + qc_filter = TRUE, + doublets = dir.exists("qc_doublets"), + integration = dir.exists("integration_qc"), + annotation_cc = dir.exists("annotation_cc"), + annotation_db = dir.exists("annotation_db"), + annotation_custom = dir.exists("annotation_custom"), + annotation_cluster = dir.exists("annotation_clusters"), + pseudobulk = dir.exists("analysis_pseudo"), + de = dir.exists("analysis_de"), + gsea = dir.exists("analysis_gsea"), + ora = dir.exists("analysis_ora") +) + +module_names <- list( + qc_filter = "Quality control and filtering", + doublets = "Doublet detection and filtering", + integration = "Dataset integration", + annotation_cc = "Cell cycle annotation", + annotation_db = "Automatic cell type annotation", + annotation_custom = "Custom cell type annotation", + annotation_cluster = "Cluster annotation", + pseudobulk = "Pseudobulking", + de = "Differential expression", + gsea = "Gene set enrichment analysis", + ora = "Over-representation analysis" +) + +modules_run_names <- unlist(module_names[unlist(modules_run)]) +modules_run_names_str <- paste("- ", modules_run_names, collapse = "\n") +modules_run_names_str <- paste0(modules_run_names_str, "\n") +``` + +## Dataset summary + +### Cohort name: `{r} cohort_name` + +### Sample sheet + +The following sample sheet table summarises the metadata and potential filtering thresholds for each sample in the dataset: + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false + +datatable( + samplesheet_show, + colnames = unname(unlist(samplesheet_cols_final)) +) +``` + +The following modules were run for this analysis: + +```{r} +#| echo: false +#| results: asis +#| warning: false +#| error: false +#| message: false + +cat(modules_run_names_str) +``` + + ## Cite us If you use these results in published work, please cite us with: diff --git a/main.nf b/main.nf index 6975ba7..84da383 100644 --- a/main.nf +++ b/main.nf @@ -77,7 +77,8 @@ workflow { gsea_db = !!params.gsea_db ? channel.fromPath(params.gsea_db, checkIfExists: true).first() : channel.value([]) // Read in samplesheet - samplesheet = channel.fromPath(params.input, checkIfExists: true) + samplesheet_csv = channel.fromPath(params.input, checkIfExists: true) + samplesheet = samplesheet_csv .splitCsv( header: true ) .map { row -> { def sample = row.sample @@ -409,7 +410,6 @@ workflow { report_style = channel.fromPath("${projectDir}/assets/styles.scss") // Gather CSV inputs - samplesheet_csv = channel.fromPath(params.input, checkIfExists: true) custom_marker_genes_csv = !!params.custom_marker_genes ? channel.fromPath(params.custom_marker_genes, checkIfExists: true) : channel.value([]) comparisons_csv = !!params.comparisons ? channel.fromPath(params.comparisons, checkIfExists: true) : channel.value([]) diff --git a/modules/report.nf b/modules/report.nf index 7c1b9ae..1b81482 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -24,9 +24,9 @@ process REPORT { path(analysis_ora_results, stageAs: "analysis_ora/*"), path(available_annotation_files, stageAs: "available_annotations/??.txt"), path(report_style) - path samplesheet - path custom_marker_genes - path comparisons + path samplesheet_csv, stageAs: "samplesheets/samplesheet.csv" + path custom_marker_genes_csv, stageAs: "samplesheets/custom_marker_genes.csv" + path comparisons_csv, stageAs: "samplesheets/comparisons.csv" output: tuple val(cohort_name), path("report"), emit: report @@ -38,6 +38,7 @@ process REPORT { def meta_fields = (metadata.meta_fields == null) ? '' : metadata.meta_fields.join('\n') """ # Save important metadata to file for reporting + echo "${cohort_name}" > cohort_name.txt mkdir -p available_annotations echo "${int_res}" > integration_resolution.txt echo -e "${pseudo_groups}" > pseudo_groups.txt From c14a662fd4c90f12e2ddcf17595d53127b8ecfee Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Wed, 1 Apr 2026 16:50:38 +1100 Subject: [PATCH 03/18] update qc filter report text --- assets/index.qmd | 9 ------- assets/qc_filter.qmd | 61 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 9 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 4306590..3ca7670 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -152,12 +152,3 @@ The following modules were run for this analysis: cat(modules_run_names_str) ``` - - -## Cite us - -If you use these results in published work, please cite us with: - -``` -# TODO: add citation information -``` \ No newline at end of file diff --git a/assets/qc_filter.qmd b/assets/qc_filter.qmd index 7beba75..674eabe 100644 --- a/assets/qc_filter.qmd +++ b/assets/qc_filter.qmd @@ -61,6 +61,67 @@ post_filter_plots <- sapply(sample_map, function(i) { s <- samples[[i]] paste("qc_filter", i, paste(s, "feature_count_plot.post_filtered.png", sep = "."), sep = "/") }) + +# Read sample sheet +samplesheet <- read_csv("samplesheets/samplesheet.csv") + +filtering_params <- samplesheet %>% + select(any_of(c( + "min_ncount", + "max_ncount", + "min_nfeature", + "max_nfeature", + "min_mt_pct", + "max_mt_pct", + "cells_to_remove" + ))) + +is_filtered <- !all(is.na(filtering_params) | filtering_params == "") +``` + +The following plots show the distributions of UMI counts (called `nCount_RNA`) and unique gene counts (called `nFeature_RNA`) per cell, as well as the mitochondiral gene percentages (`percent.mt`) per cell. These three values are plotted in the "Features vs counts" scatter plots towards the bottom of the page. Ideally, you want to see a roughly linear relationship between these values, up to a plateau point where greater sequencing depth no longer captures additional genes. Poor quality or dying cells may appear as additional clusters of cells towards the bottom-right corners of the plots. Additionally, high mitochondrial gene percentages may be indicative of dying cells, although some cell types may show high mitochondrial gene percentages even when healthy. + +The table at the bottom of the page summarises the number of cells in each sample. If per-sample filtering was performed, the pre- and post-filter cell counts, as well as the percentage of cells filtered per sample, are shown in this table as well. + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false +#| results: asis +#| eval: !expr is_filtered + +markdown_text <- "#### Filtering thresholds + +The samples in the dataset were filtered according to the following thresholds: + +" +``` + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false +#| eval: !expr is_filtered + +filtering_params_cols <- c( + min_ncount = "Min. UMI count", + max_ncount = "Max. UMI count", + min_nfeature = "Min. gene count", + max_nfeature = "Max. gene count", + min_mt_pct = "Min. MT gene proportion (%)", + max_mt_pct = "Max. MT gene proportion (%)", + cells_to_remove = "Cells to remove file path" +) + +filtering_params_non_null <- filtering_params %>% + select(where(~ !all(is.na(.x) | .x == ""))) %>% + filter(!if_all(everything(), ~ is.na(.x) | .x == "")) + +filtering_params_non_null_cols <- filtering_params_cols[colnames(filtering_params_non_null)] + +datatable(filtering_params_non_null, colnames = filtering_params_non_null_cols) ``` :::panel-tabset From 721a3afc974fce48e604b90cb372473221a400d3 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Mon, 20 Apr 2026 17:05:54 +1000 Subject: [PATCH 04/18] update qc filter and cluster reports with text --- assets/qc_cluster.qmd | 28 ++++++++++++++++++++++++++++ assets/qc_filter.qmd | 8 ++++++-- 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/assets/qc_cluster.qmd b/assets/qc_cluster.qmd index b734ad8..ea894e8 100644 --- a/assets/qc_cluster.qmd +++ b/assets/qc_cluster.qmd @@ -56,6 +56,34 @@ for (n in names(feature_count_plot_files)) { } ``` +During the clustering stage, cells are grouped together based on the similarities between their expression profiles. This process is dependent upon a parameter called the "resolution". Smaller resolutions lead to larger clusters, while larger resolutions lead to smaller clusters. It is important to find the right clustering resolution for your dataset so that you capture as much of the information about your sample as possible without creating lots of tiny clusters that have no biological meaning. + +To aid in this process, we have clustered your data at multiple resolutions (see below for all of the resolutions used). We have then used the R package `clustree` to create a graph representation of your clustering results at each of the resolutions. These graphs show a cluster tree where each circle represents a cluster, and each row of clusters represents a given resolution. The rows of clusters are organised from smallest to largest clustering resolution. The arrows represent how cells move between clusters when going from one resolution to the next. The size and opacity of the arrows represents the proportion of cells moving from one cluster to the next. + +A good rule of thumb for selecting a clustering resolution for your data is to look for the resolutions at which the clusters remain relatively stable, i.e. there is little movement of cells between clusters. In the cluster trees this looks like rows of clusters where most of the clusters have a single, solid arrow incoming from the previous resolution. In contrast, unstable resolutions will have clusters with many, faint arrows incoming, representing clusters that are splitting and merging with each other from one resolution to the next. + +Below, you will find cluster trees for each of your samples, plus count vs. feature scatter plots for each of the clusters at each resolution. + +#### Clustering resolutions + +The following clustering resolutions were applied to each of your samples: + +```{r} +#| echo: false +#| results: asis +#| warning: false +#| error: false +#| message: false + +cluster_res_all <- cluster_res[[1]] +cluster_res_all <- as.numeric((cluster_res_all)) +cluster_res_all <- sort(cluster_res_all) +cluster_res_str <- paste("- ", cluster_res_all, collapse = "\n") +cluster_res_str <- paste0(cluster_res_str, "\n") + +cat(cluster_res_str) +``` + :::panel-tabset ```{r} diff --git a/assets/qc_filter.qmd b/assets/qc_filter.qmd index 674eabe..5be0cec 100644 --- a/assets/qc_filter.qmd +++ b/assets/qc_filter.qmd @@ -67,6 +67,7 @@ samplesheet <- read_csv("samplesheets/samplesheet.csv") filtering_params <- samplesheet %>% select(any_of(c( + "sample", "min_ncount", "max_ncount", "min_nfeature", @@ -79,7 +80,7 @@ filtering_params <- samplesheet %>% is_filtered <- !all(is.na(filtering_params) | filtering_params == "") ``` -The following plots show the distributions of UMI counts (called `nCount_RNA`) and unique gene counts (called `nFeature_RNA`) per cell, as well as the mitochondiral gene percentages (`percent.mt`) per cell. These three values are plotted in the "Features vs counts" scatter plots towards the bottom of the page. Ideally, you want to see a roughly linear relationship between these values, up to a plateau point where greater sequencing depth no longer captures additional genes. Poor quality or dying cells may appear as additional clusters of cells towards the bottom-right corners of the plots. Additionally, high mitochondrial gene percentages may be indicative of dying cells, although some cell types may show high mitochondrial gene percentages even when healthy. +The plots on this page show the distributions of UMI counts (called `nCount_RNA`) and unique gene counts (called `nFeature_RNA`) per cell, as well as the mitochondiral gene percentages (`percent.mt`) per cell. These three values are plotted in the "Features vs counts" scatter plots towards the bottom of the page. Ideally, you want to see a roughly linear relationship between these values, up to a plateau point where greater sequencing depth no longer captures additional genes. Poor quality or dying cells may appear as additional clusters of cells towards the bottom-right corners of the plots. Additionally, high mitochondrial gene percentages may be indicative of dying cells, although some cell types may show high mitochondrial gene percentages even when healthy. The table at the bottom of the page summarises the number of cells in each sample. If per-sample filtering was performed, the pre- and post-filter cell counts, as well as the percentage of cells filtered per sample, are shown in this table as well. @@ -96,6 +97,8 @@ markdown_text <- "#### Filtering thresholds The samples in the dataset were filtered according to the following thresholds: " + +cat(markdown_text) ``` ```{r} @@ -106,6 +109,7 @@ The samples in the dataset were filtered according to the following thresholds: #| eval: !expr is_filtered filtering_params_cols <- c( + sample = "Sample name", min_ncount = "Min. UMI count", max_ncount = "Max. UMI count", min_nfeature = "Min. gene count", @@ -121,7 +125,7 @@ filtering_params_non_null <- filtering_params %>% filtering_params_non_null_cols <- filtering_params_cols[colnames(filtering_params_non_null)] -datatable(filtering_params_non_null, colnames = filtering_params_non_null_cols) +datatable(filtering_params_non_null, colnames = unname(filtering_params_non_null_cols)) ``` :::panel-tabset From 64de52ddd8a4224552d2cfe29e5063dc343f119b Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Tue, 21 Apr 2026 13:33:32 +1000 Subject: [PATCH 05/18] add text to doublet and integration reports --- assets/integration_cluster.qmd | 8 ++++++++ assets/integration_qc.qmd | 8 ++++++++ assets/qc_doublets.qmd | 8 ++++++++ 3 files changed, 24 insertions(+) diff --git a/assets/integration_cluster.qmd b/assets/integration_cluster.qmd index 4d41305..51aaf91 100644 --- a/assets/integration_cluster.qmd +++ b/assets/integration_cluster.qmd @@ -41,11 +41,15 @@ integration_res_file <- "integration_resolution.txt" integration_res <- as.numeric(scan(integration_res_file, character())[1]) ``` +This page shows the clustering results for your dataset post-integration. As before with the [per-sample clustering analysis](qc_cluster.qmd), the dataset was clustered at several resolutions (outlined below). A cluster tree plot is shown below to aid you in selecting the best resolution for your dataset. Recall that you want to identify a resolution where the clusters remain relatively stable from one resolution to the next, i.e. most clusters have a single, opaque arrow incoming rather than multiple faint incoming arrows. When you have selected a suitable resolution, re-run the pipeline and pass the parameter `--integrated_resolution `, replacing `` with your chosen resolution. When this parameter is unset, the pipeline choose a default value of `1`, **but you will likely need to change this**. The UMAP and count vs. feature plots below will display the results for your chosen resolution. + #### Dataset name: `{r} cohort_id` #### Integration resolution: `{r} integration_res` ### Integrated clusters UMAP +These UMAP plots show your integrated dataset, colour-coded by either the sample of origin ("orig.ident") or the cluster each cell was assigned to in the chosen integration clustering resolution ("SCT_snn_res.*"). + ```{r} #| results: asis #| echo: false @@ -54,6 +58,8 @@ cat(paste0('integrated clusters umap: ', cohort_ ### Feature-count plot of integrated dataset +This plot shows the count vs. feature plots for each cluster at your chosen clustering resolution. + ```{r} #| results: asis #| echo: false @@ -62,6 +68,8 @@ cat(paste0('<img src= Date: Wed, 22 Apr 2026 16:53:41 +1000 Subject: [PATCH 06/18] update annotation report page, add cluster annotation to report module --- assets/annotation.qmd | 54 +++++++++++++++++++++++++++++++++++++++++++ main.nf | 3 ++- modules/report.nf | 2 ++ 3 files changed, 58 insertions(+), 1 deletion(-) diff --git a/assets/annotation.qmd b/assets/annotation.qmd index ddde05d..6e231a7 100644 --- a/assets/annotation.qmd +++ b/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,10 +151,51 @@ if (is.null(cohort_id)) { } ``` +There are three main annotation options you can choose when running 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. This annotation is automatic for human samples and can be enabled for non-human samples by providing lists of cell cycle marker genes for your species using the `--s_genes` and `--g2m_genes` parameters. + +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, although you can override the database used with the `--annotation_db` parameter. For non-human species, you will need to use this parameter and provide a database file 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. You provide this list of cell types using the `--custom_marker_genes` parameter. The file provided must be a CSV file and contain an annotation label and a list of genes associated with it. Cells that score highly for more than one cell type will be labelled "Ambiguous". + +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". By default, a cell type must make up 2/3 of the cluster for that label to be applied to the cluster as well. You can choose the threshold used with the `--cell_type_proportion_threshold` parameter. + +By default, the annotation category used for cluster annotation is the first in the following list that has been performed: + +- Custom marker gene-based annotation +- Automatic cell type annotation + - For human HPCA annotations, the 'main' level annotation is prioritised over the 'fine' level annotation +- Cell cycle annotation + +The annotation category used for cluster annotation can be overridden using the `--cluster_annotation`. See below for a list of valid values from your dataset. + #### Dataset name: `{r} cohort_id` +### All annotation fields available for cluster annotation + +The following annotation fields are available for cluster annotation: + +```{r} +#| echo: false +#| results: asis +for (ann in annotation_priority) { + cat("- ") + cat(ann) + cat("\n") +} +``` + +The field used for cluster annotation was: **`{r} annotation_used_for_clusters`** + ### All annotation fields available for downstream analyses +The following list contains all of the annotations listed above, as well as `cluster_annotation` and additional metadata fields you provided in the samplesheet. These are all of the fields available to you for performing downstream pseudobulking and subsequent differential expression and functional enrichment analyses. + ```{r} #| echo: false #| results: asis diff --git a/main.nf b/main.nf index 84da383..285beb8 100644 --- a/main.nf +++ b/main.nf @@ -342,7 +342,8 @@ workflow { .map { meta_fields -> [ meta_fields ] } .merge(integrated_resolution) .merge(pseudo_groups) - .map { meta, int_res, pseudo -> [ meta_fields: meta, integration_resolution: int_res, pseudo_groups: pseudo ] } + .merge(cluster_annotation) + .map { meta, int_res, pseudo -> [ meta_fields: meta, integration_resolution: int_res, pseudo_groups: pseudo, cluster_annotation: cluster_annotation ] } // Gather all report templates index_template = channel.fromPath("${projectDir}/assets/index.qmd", checkIfExists: true) diff --git a/modules/report.nf b/modules/report.nf index 1b81482..0ee6952 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -35,6 +35,7 @@ process REPORT { script: def int_res = (metadata.integration_resolution == null) ? '' : metadata.integration_resolution def pseudo_groups = (metadata.pseudo_groups == null) ? '' : metadata.pseudo_groups.tokenize(',').join('\n') + def cluster_annotation = (metadata.cluster_annotation == null) ? '' : metadata.cluster_annotation def meta_fields = (metadata.meta_fields == null) ? '' : metadata.meta_fields.join('\n') """ # Save important metadata to file for reporting @@ -42,6 +43,7 @@ process REPORT { mkdir -p available_annotations echo "${int_res}" > integration_resolution.txt echo -e "${pseudo_groups}" > pseudo_groups.txt + echo -e "${cluster_annotation}" > cluster_annotation.txt echo -e "${meta_fields}" > available_annotations/metadata.txt if [ -d "annotation_clusters" ]; then echo "cluster_annotation" > available_annotations/clusters.txt From 75716074300c06e60954f75f603b10d62d619999 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Wed, 29 Apr 2026 15:38:11 +1000 Subject: [PATCH 07/18] update pseudo report --- assets/analysis_gsea.qmd | 2 +- assets/analysis_pseudo.qmd | 23 ++++++++++++++++++++++- 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/assets/analysis_gsea.qmd b/assets/analysis_gsea.qmd index d59fb88..e19d6fb 100644 --- a/assets/analysis_gsea.qmd +++ b/assets/analysis_gsea.qmd @@ -115,4 +115,4 @@ purrr::walk( ) ``` -::: \ No newline at end of file +::: diff --git a/assets/analysis_pseudo.qmd b/assets/analysis_pseudo.qmd index e96eb6a..0482b59 100644 --- a/assets/analysis_pseudo.qmd +++ b/assets/analysis_pseudo.qmd @@ -37,8 +37,27 @@ pseudo_groups_file <- "pseudo_groups.txt" pseudo_groups <- scan(pseudo_groups_file, character()) ``` +When we run a differential expression analysis on a single cell RNA-seq dataset, we are constructing a statistical model that looks for changes in the mean expression of each gene between case and control samples. These statistical models are built upon several **assumptions** about the underlying data, and if our dataset violates those assumptions, the results can be incorrect. + +One of the major assumptions that these models make is that each data point is an independent, random sample drawn from some population. In the case of a bulk RNA-seq experiment, this assumption is valid: each data point is an RNA sample from an individual person or organism, and those individuals were (hopefully) chosen randomly from their respective populations. In the case of a single-cell RNA-seq experiment though, things get a bit more complicated. + +If we naively pass our scRNA-seq data into a standard differential expression algorithm, it will treat **each cell** as a separate data point and make the **assumption** that each cell is independently sampled from the population. But that isn't the case: cells from the same individual are **not** independent of one another - they originate from the same individual and share the exact same genetic background. This means the cells' expression profiles are correlated with one another, and this will invalidate the independence assumption, bias our results, and ultimately lead to over-inflated test statistics and estimates of significance. This issue is called **pseudoreplication**. + +There are a few ways to deal with pseudoreplication, but the method we have implemented in this pipeline is called **pseudobulking**. Pseudobulking 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. Once complete, the B cell profiles from each sample would be suitably independent of one another and could be compared between cases and controls; the same could be done for the T cell profiles. + +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. + +When the pseudobulking stage of the pipeline runs, it uses the metadata fields provided with the `--pseudo_groups` parameter to decide how to group the cells prior to aggregation. By default, the pipeline will use the `cluster_annotation` field generated in the [annotation stage](annotation.qmd). This means all cells with the same cluster annotation within a sample will be aggregated together. You can provide multiple grouping fields for to this parameter, separated by commas, to refine the groups; for example, if you have a sample of B and T cells and ran both cluster annotation and cell cycle annotation, `--pseudo_groups cluster_annotation,Phase` would result in six pseudobulked expression profiles per cell: S-phase B cells, S-phase T cells, G1-phase B cells, etc. + +See the section titled [All annotation fields available for downstream analyses](annotation.qmd#all-annotation-fields-available-for-downstream-analyses) for a full list of fields relevant to your data that you can provide to `--pseudo_groups`. **Note** that you should include sample-level metadata fields here that are relevant to your differential expression analysis; for example, if you are comparing a particular treatment vs a control, and the relevant metadata field was called `treatment`, you would provide `treatment` to the `--pseudo_groups` parameter, along with any other fields you need. + +The pseudobulked expression profiles will be named after the individual levels within each metadata field. Again, taking the example of B and T cells in a treatment vs control experiment, if the cell labels are `B-cell` and `T-cell` and the labels in the `treatment` field are `drug` and `control`, your pseudobulked groups will be called `B-cell_drug`, `B-cell_control`, `T-cell_drug`, and `T-cell_control`. + #### Dataset name: `{r} cohort_id` -#### Fields used for pseudobulking: + +#### Fields used for pseudobulking + +The following metadata fields were provided to the `--pseudo_groups` parameter and used for aggregation of cells: ```{r} #| echo: false @@ -52,6 +71,8 @@ for (g in pseudo_groups) { ### Pseudobulked groups availablle for comparison +Below is a table summarising all your final pseudobulk groups and the number of samples in your dataset that are represented in each group. You will need to note down these group names to construct a list of your desired differential expression comparisons. **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 From ea4c22680495d141c1d21e79ef9a09adcd418a61 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Thu, 30 Apr 2026 16:16:24 +1000 Subject: [PATCH 08/18] finish adding explainer text --- assets/analysis_de.qmd | 35 ++++++++++++++++++++++++++++------- assets/analysis_gsea.qmd | 20 +++++++++++++------- assets/analysis_ora.qmd | 20 +++++++++++++------- 3 files changed, 54 insertions(+), 21 deletions(-) diff --git a/assets/analysis_de.qmd b/assets/analysis_de.qmd index a806266..b5147c6 100644 --- a/assets/analysis_de.qmd +++ b/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,31 @@ 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") ``` +This page shows the differential expression results for the comparisons you requested on your dataset. 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) + #### Dataset name: `{r} cohort_id` #### All comparisons made: +You provided the following table of comparisons to perform via the `--comparisons` pipeline parameter: + ```{r} -#| echo: false -#| results: asis -for (cmp in all_comparisons) { - cat("- ") - cat(cmp) - cat("\n") -} +datatable(comparisons_df, colnames = c("Reference group", "Test group")) ``` +*Note that the group names above must match the [pseudobulk groups shown in the previous section of the report](analysis_pseudo.qmd#pseudobulked-groups-availablle-for-comparison).* + ### Differential expression results #### 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 shoudl 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 +116,12 @@ cat(paste0('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 (% .$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") }) ``` +Differential expression analyses can generate lots of results that can be quite difficult to interpret by themselves. Functional enrichment analyses are a set of methods designed to take the raw differential expression results and identify patterns of expression that point to altered **biological function**. There are several such approaches, but the two we apply in this pipeline are **Gene Set Enrichment Analysis** (GSEA) and **Over-Representation Analysis** (ORA). + +This page 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. + #### Dataset name: `{r} cohort_id` #### All comparisons made: +One GSEA analysis was conducted for each differential expression comparison that was made. As per the [differential expression page](analysis_de.qmd), these comparisons were: + ```{r} -#| echo: false -#| results: asis -for (cmp in all_comparisons) { - cat("- ") - cat(cmp) - cat("\n") -} +datatable(comparisons_df, colnames = c("Reference group", "Test group")) ``` ### GSEA expression results +Below, you will find plots for each GSEA analysis. These plots show the enrichment score for each **significant** pathway from the analyses; more upregulated pathways have bars pointing to the right and are highlighted in orange; downregulated pathways point to the left and are highlighted in blue. + :::panel-tabset ```{r} diff --git a/assets/analysis_ora.qmd b/assets/analysis_ora.qmd index b227ea3..43a21f9 100644 --- a/assets/analysis_ora.qmd +++ b/assets/analysis_ora.qmd @@ -78,28 +78,34 @@ 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") }) ``` +Differential expression analyses can generate lots of results that can be quite difficult to interpret by themselves. Functional enrichment analyses are a set of methods designed to take the raw differential expression results and identify patterns of expression that point to altered **biological function**. There are several such approaches, but the two we apply in this pipeline are **Gene Set Enrichment Analysis** (GSEA) and **Over-Representation Analysis** (ORA). + +This page 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. + #### Dataset name: `{r} cohort_id` #### All comparisons made: +One GSEA analysis was conducted for each differential expression comparison that was made. As per the [differential expression page](analysis_de.qmd), these comparisons were: + ```{r} -#| echo: false -#| results: asis -for (cmp in all_comparisons) { - cat("- ") - cat(cmp) - cat("\n") -} +datatable(comparisons_df, colnames = c("Reference group", "Test group")) ``` ### ORA expression results +Below, you will find plots for each ORA analysis. These plots show 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; depleted or under-represented pathways point to the left and are highlighted in blue. + :::panel-tabset ```{r} From 477790718880401c634e11e1972e8749886951ae Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Fri, 29 May 2026 12:05:43 +1000 Subject: [PATCH 09/18] WIP: make single-page, collapse cluster plots into tabsets --- assets/index.qmd | 189 ++++++++++++++++++++++++++++++++++++++++- assets/qc_cluster.qmd | 5 +- assets/qc_doublets.qmd | 5 +- 3 files changed, 195 insertions(+), 4 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 3ca7670..14496fc 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -1,11 +1,13 @@ --- title: "scRNAvigator Report" -format: html +format: + html: + embed-resources: true --- Thank you for using scRNAvigator! -This website contains a report of your analysis. Depending on what sections of the pipeline were ran, it will include one or more of the following: +This report contains a summary of your analysis. Depending on what sections of the pipeline were ran, it will include one or more of the following: - Quality control metrics and plots - Doublet detection and filtering results @@ -152,3 +154,186 @@ The following modules were run for this analysis: cat(modules_run_names_str) ``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("qc_filter.qmd") +#| warning: false +#| error: false +#| message: false +cat("# Quality control\n\n") +cat("## Initial filtering\n\n") +``` + +```{r} +#| child: "qc_filter.qmd" +#| eval: !expr file.exists("qc_filter.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("qc_cluster.qmd") +#| warning: false +#| error: false +#| message: false +cat("## Clustering\n\n") +``` + +```{r} +#| child: "qc_cluster.qmd" +#| eval: !expr file.exists("qc_cluster.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("qc_doublets.qmd") +#| warning: false +#| error: false +#| message: false +cat("## Doublet detection and filtering\n\n") +``` + +```{r} +#| child: "qc_doublets.qmd" +#| eval: !expr file.exists("qc_doublets.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("integration_qc.qmd") +#| warning: false +#| error: false +#| message: false +cat("# Integration\n\n") +cat("## Initial quality control\n\n") +``` + +```{r} +#| child: "integration_qc.qmd" +#| eval: !expr file.exists("integration_qc.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("integration_cluster.qmd") +#| warning: false +#| error: false +#| message: false +cat("## Clustering\n\n") +``` + +```{r} +#| child: "integration_cluster.qmd" +#| eval: !expr file.exists("integration_cluster.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("annotation.qmd") +#| warning: false +#| error: false +#| message: false +cat("# Annotation\n\n") +``` + +```{r} +#| child: "annotation.qmd" +#| eval: !expr file.exists("annotation.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("analysis_pseudo.qmd") +#| warning: false +#| error: false +#| message: false +cat("# Analysis\n\n") +cat("## Pseudobulking\n\n") +``` + +```{r} +#| child: "analysis_pseudo.qmd" +#| eval: !expr file.exists("analysis_pseudo.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("analysis_de.qmd") +#| warning: false +#| error: false +#| message: false +cat("## Differential expression analysis\n\n") +``` + +```{r} +#| child: "analysis_de.qmd" +#| eval: !expr file.exists("analysis_de.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("analysis_gsea.qmd") +#| warning: false +#| error: false +#| message: false +cat("## Functional enrichment analysis: GSEA\n\n") +``` + +```{r} +#| child: "analysis_gsea.qmd" +#| eval: !expr file.exists("analysis_gsea.qmd") +#| warning: false +#| error: false +#| message: false +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("analysis_ora.qmd") +#| warning: false +#| error: false +#| message: false +cat("## Functional enrichment analysis: ORA\n\n") +``` + +```{r} +#| child: "analysis_ora.qmd" +#| eval: !expr file.exists("analysis_ora.qmd") +#| warning: false +#| error: false +#| message: false +``` diff --git a/assets/qc_cluster.qmd b/assets/qc_cluster.qmd index ea894e8..102b473 100644 --- a/assets/qc_cluster.qmd +++ b/assets/qc_cluster.qmd @@ -98,11 +98,14 @@ purrr::walk( cat(paste("####", "Clustree plot", "\n\n")) cat(paste0('clustree plot: ', n, '\n\n')) cat(paste("####", "Per-cluster feature-count plots", "\n\n")) + cat(":::panel-tabset\n\n") for (res in sort(names(p_fc_plot_list))) { p_fc_plot <- p_fc_plot_list[[res]] - cat(paste("#####", "Cluster resolution:", res, "\n\n")) + cat(paste("#####", res, "\n\n")) + cat(paste("Cluster resolution:", res, "\n\n")) cat(paste0('Per-cluster feature-count plot: ', n, '; res: ', res, '\n\n')) } + cat(":::\n\n") } ) ``` diff --git a/assets/qc_doublets.qmd b/assets/qc_doublets.qmd index b0b2678..347805b 100644 --- a/assets/qc_doublets.qmd +++ b/assets/qc_doublets.qmd @@ -128,11 +128,14 @@ purrr::walk( cat(paste("####", "Clustree plot", "\n\n")) cat(paste0('clustree plot: ', n, '\n\n')) cat(paste("####", "Per-cluster feature-count plots", "\n\n")) + cat(":::panel-tabset\n\n") for (res in sort(names(p_fc_plot_list))) { p_fc_plot <- p_fc_plot_list[[res]] - cat(paste("#####", "Cluster resolution:", res, "\n\n")) + cat(paste("#####", res, "\n\n")) + cat(paste("Cluster resolution:", res, "\n\n")) cat(paste0('Per-cluster feature-count plot: ', n, '; res: ', res, '\n\n')) } + cat(":::\n\n") } ) ``` From 42e8d8f470c733a8c21013dded3ceff497359be8 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Mon, 1 Jun 2026 12:24:33 +1000 Subject: [PATCH 10/18] add parameter list and executive summary section to report --- assets/analysis_de.qmd | 16 ++-- assets/analysis_gsea.qmd | 8 +- assets/analysis_ora.qmd | 8 +- assets/analysis_pseudo.qmd | 2 +- assets/annotation.qmd | 34 ++++---- assets/index.qmd | 154 +++++++++++++++++++++++++++++---- assets/integration_cluster.qmd | 10 +-- assets/integration_qc.qmd | 6 +- assets/qc_cluster.qmd | 8 +- assets/qc_doublets.qmd | 14 +-- assets/qc_filter.qmd | 8 +- bin/annotate_clusters.R | 6 ++ bin/cluster_integrated.R | 6 ++ modules/report.nf | 11 ++- 14 files changed, 213 insertions(+), 78 deletions(-) diff --git a/assets/analysis_de.qmd b/assets/analysis_de.qmd index b5147c6..0617f01 100644 --- a/assets/analysis_de.qmd +++ b/assets/analysis_de.qmd @@ -90,7 +90,7 @@ This page shows the differential expression results for the comparisons you requ - [Volcano plots](#volcano-plots) - [Differential expression results table](#differential-expression-table) -#### Dataset name: `{r} cohort_id` +##### Dataset name: `r cohort_id` #### All comparisons made: @@ -102,9 +102,9 @@ datatable(comparisons_df, colnames = c("Reference group", "Test group")) *Note that the group names above must match the [pseudobulk groups shown in the previous section of the report](analysis_pseudo.qmd#pseudobulked-groups-availablle-for-comparison).* -### Differential expression results +#### Differential expression results -#### p-value distributions +##### 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 shoudl 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. @@ -114,7 +114,7 @@ The following plots show the distribution of p-values for each comparison made. cat(paste0('p-value distributions: ', cohort_id, '\n\n')) ``` -#### Log-fold-change distributions +##### 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"). @@ -128,7 +128,7 @@ Keep in mind that in many biological contexts, small fold-changes can be very bi cat(paste0('log-fold-change distributions: ', cohort_id, '\n\n')) ``` -#### Volcano plots +##### 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 you have provided a fold-change threshold, 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. @@ -138,13 +138,13 @@ Volcano plots combine both the p-value and fold-change information into one very cat(paste0('volcano plots: ', cohort_id, '\n\n')) ``` -#### Differential expression table +##### 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 +###### Significant results only ```{r} #| echo: false @@ -154,7 +154,7 @@ Finally, below is a table summarising all of the differential expression results datatable(sig_de_table, colnames = de_table_col_names_ordered) ``` -##### All results +###### All results ```{r} #| echo: false diff --git a/assets/analysis_gsea.qmd b/assets/analysis_gsea.qmd index 8764eb5..e9ebc61 100644 --- a/assets/analysis_gsea.qmd +++ b/assets/analysis_gsea.qmd @@ -92,7 +92,7 @@ Differential expression analyses can generate lots of results that can be quite This page 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. -#### Dataset name: `{r} cohort_id` +##### Dataset name: `r cohort_id` #### All comparisons made: @@ -102,7 +102,7 @@ One GSEA analysis was conducted for each differential expression comparison that datatable(comparisons_df, colnames = c("Reference group", "Test group")) ``` -### GSEA expression results +#### GSEA expression results Below, you will find plots for each GSEA analysis. These plots show the enrichment score for each **significant** pathway from the analyses; more upregulated pathways have bars pointing to the right and are highlighted in orange; downregulated pathways point to the left and are highlighted in blue. @@ -114,8 +114,8 @@ 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('gsea results: ', n, '\n\n')) } ) diff --git a/assets/analysis_ora.qmd b/assets/analysis_ora.qmd index 43a21f9..057049a 100644 --- a/assets/analysis_ora.qmd +++ b/assets/analysis_ora.qmd @@ -92,7 +92,7 @@ Differential expression analyses can generate lots of results that can be quite This page 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. -#### Dataset name: `{r} cohort_id` +##### Dataset name: `r cohort_id` #### All comparisons made: @@ -102,7 +102,7 @@ One GSEA analysis was conducted for each differential expression comparison that datatable(comparisons_df, colnames = c("Reference group", "Test group")) ``` -### ORA expression results +#### ORA expression results Below, you will find plots for each ORA analysis. These plots show 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; depleted or under-represented pathways point to the left and are highlighted in blue. @@ -114,8 +114,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('ora results: ', n, '\n\n')) } ) diff --git a/assets/analysis_pseudo.qmd b/assets/analysis_pseudo.qmd index 0482b59..d3ac4d9 100644 --- a/assets/analysis_pseudo.qmd +++ b/assets/analysis_pseudo.qmd @@ -53,7 +53,7 @@ See the section titled [All annotation fields available for downstream analyses] The pseudobulked expression profiles will be named after the individual levels within each metadata field. Again, taking the example of B and T cells in a treatment vs control experiment, if the cell labels are `B-cell` and `T-cell` and the labels in the `treatment` field are `drug` and `control`, your pseudobulked groups will be called `B-cell_drug`, `B-cell_control`, `T-cell_drug`, and `T-cell_control`. -#### Dataset name: `{r} cohort_id` +##### Dataset name: `r cohort_id` #### Fields used for pseudobulking diff --git a/assets/annotation.qmd b/assets/annotation.qmd index 6e231a7..04f1416 100644 --- a/assets/annotation.qmd +++ b/assets/annotation.qmd @@ -174,9 +174,9 @@ By default, the annotation category used for cluster annotation is the first in The annotation category used for cluster annotation can be overridden using the `--cluster_annotation`. See below for a list of valid values from your dataset. -#### Dataset name: `{r} cohort_id` +##### Dataset name: `r cohort_id` -### All annotation fields available for cluster annotation +#### All annotation fields available for cluster annotation The following annotation fields are available for cluster annotation: @@ -190,9 +190,9 @@ for (ann in annotation_priority) { } ``` -The field used for cluster annotation was: **`{r} annotation_used_for_clusters`** +The field used for cluster annotation was: **`r annotation_used_for_clusters`** -### All annotation fields available for downstream analyses +#### All annotation fields available for downstream analyses The following list contains all of the annotations listed above, as well as `cluster_annotation` and additional metadata fields you provided in the samplesheet. These are all of the fields available to you for performing downstream pseudobulking and subsequent differential expression and functional enrichment analyses. @@ -212,7 +212,7 @@ 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('cell cycle umap: ', cohort_id, '\n\n')) ``` @@ -220,10 +220,10 @@ cat(paste0('cell cycle umap: ', cohort_id, '\n\n')) -cat("#### Cell type proportions\n\n") +cat("##### Cell type proportions\n\n") cat(paste0('cell type proportions: ', cohort_id, '\n\n')) ``` @@ -231,7 +231,7 @@ cat(paste0('cell type proportions: ', coho #| results: asis #| echo: false #| eval: !expr db_annotated && file.exists(cell_type_prop_table_hpca_main) -cat(\n\n')) -cat("#### Cell type proportions\n\n") +cat("##### Cell type proportions\n\n") cat(paste0('cell type proportions: ', cohort_id, '\n\n')) ``` @@ -299,10 +299,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('feature plot: ', cell_type, ' - ', cohort_id, '\n\n')) } ``` @@ -311,7 +311,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('cluster cell type assignment umap: ', cohort_id, '\n\n')) ``` diff --git a/assets/index.qmd b/assets/index.qmd index 14496fc..68da15b 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -28,6 +28,10 @@ library(DT) # Get cohort name cohort_name <- scan("cohort_name.txt", character())[1] +# Get parameters TSV file +params_tsv <- "parameters.tsv" +params_df <- read_tsv(params_tsv) + # Read sample sheet samplesheet <- read_csv("samplesheets/samplesheet.csv") @@ -121,11 +125,111 @@ module_names <- list( modules_run_names <- unlist(module_names[unlist(modules_run)]) modules_run_names_str <- paste("- ", modules_run_names, collapse = "\n") modules_run_names_str <- paste0(modules_run_names_str, "\n") + +# Get executive summary info +n_samples <- c(metric = "# Samples", value = nrow(samplesheet_show)) +n_cells_filtered <- NULL +n_cells_per_sample <- NULL +n_clusters <- NULL +cluster_annotations <- NULL +n_de_comparisons <- NULL +n_de_genes <- NULL + +if (file.exists("qc_doublets.qmd")) { + doublet_summary_files <- dir( + "qc_doublets", + pattern = ".*\\.doublet_summary\\.csv$", + recursive = TRUE, + full.names = TRUE + ) + doublet_summaries <- lapply(doublet_summary_files, function(f) { + read_csv(f) + }) + doublet_summary <- bind_rows(doublet_summaries) + n_cells_filtered <- c(metric = "# Cells (filtered)", value = sum(doublet_summary$n_singlets)) +} else if (file.exists("qc_filter.qmd")) { + filter_summary_files <- dir( + "qc_filter", + pattern = ".*\\.hard_filter_summary\\.csv$", + recursive = TRUE, + full.names = TRUE + ) + filter_summaries <- lapply(filter_summary_files, function(f) { + read_csv(f) + }) + filter_summary <- bind_rows(filter_summaries) + n_cells_filtered <- c(metric = "# Cells (filtered)", value = sum(filter_summary$n_cells_post_filtered)) +} + +if (!is.null(n_cells_filtered)) { + n_cells_per_sample <- c(metric = "Avg. cells per sample", value = round(as.integer(n_cells_filtered[["value"]]) / as.integer(n_samples[["value"]]), digits = 2)) +} + +if (file.exists("integration_cluster.qmd")) { + n_clusters_df <- read_csv(paste0("integration_cluster/", cohort_name, ".integrated.clusters.csv")) + n_clusters <- c(metric = "# Clusters", value = nrow(n_clusters_df)) +} + +if (file.exists("annotation.qmd")) { + cluster_annotations <- read_csv(paste0("annotation_clusters/", cohort_name, ".cluster_cell_type_assignments.csv")) %>% + group_by(cell_type) %>% + summarise(n_clusters = n()) +} + +if (file.exists("analysis_de.qmd")) { + comparisons_csv <- "samplesheets/comparisons.csv" + comparisons_df <- read_csv(comparisons_csv) + n_de_comparisons <- c(metric = "# DE comparisons", value = nrow(comparisons_df)) + + sig_de_table <- read_csv(paste0("analysis_de/", cohort_name, ".de.full.csv")) %>% + filter(sig == "sig") + + n_de_genes <- c(metric = "# DE genes", value = length(unique(sig_de_table$Gene.Symbol))) +} + +exec_summary <- bind_rows(n_samples, n_cells_filtered, n_cells_per_sample, n_clusters, n_de_comparisons, n_de_genes) ``` -## Dataset summary +## Analysis summary: `r cohort_name` -### Cohort name: `{r} cohort_name` +The following table summarises the major findings from this analysis: + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false + +datatable( + exec_summary, + colnames = "", + rownames = FALSE +) +``` + +```{r} +#| echo: false +#| results: asis +#| eval: !expr file.exists("annotation.qmd") +#| warning: false +#| error: false +#| message: false +cat("The following annotations were assigned to clusters:\n\n") +``` + +```{r} +#| echo: false +#| eval: !expr (!is.null(cluster_annotations)) +#| warning: false +#| error: false +#| message: false + +datatable( + cluster_annotations, + colnames = c("Cell type", "# Clusters"), + rownames = FALSE +) +``` ### Sample sheet @@ -155,6 +259,26 @@ The following modules were run for this analysis: cat(modules_run_names_str) ``` +### Parameters + +The following parameters were set during the pipeline run: + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false + +datatable( + params_df, + colnames = c("Parameter", "Value") +) +``` + +### Methods + +TODO + ```{r} #| echo: false #| results: asis @@ -162,8 +286,8 @@ cat(modules_run_names_str) #| warning: false #| error: false #| message: false -cat("# Quality control\n\n") -cat("## Initial filtering\n\n") +cat("## Quality control\n\n") +cat("### Initial filtering\n\n") ``` ```{r} @@ -181,7 +305,7 @@ cat("## Initial filtering\n\n") #| warning: false #| error: false #| message: false -cat("## Clustering\n\n") +cat("### Clustering\n\n") ``` ```{r} @@ -199,7 +323,7 @@ cat("## Clustering\n\n") #| warning: false #| error: false #| message: false -cat("## Doublet detection and filtering\n\n") +cat("### Doublet detection and filtering\n\n") ``` ```{r} @@ -217,8 +341,8 @@ cat("## Doublet detection and filtering\n\n") #| warning: false #| error: false #| message: false -cat("# Integration\n\n") -cat("## Initial quality control\n\n") +cat("## Integration\n\n") +cat("### Initial quality control\n\n") ``` ```{r} @@ -236,7 +360,7 @@ cat("## Initial quality control\n\n") #| warning: false #| error: false #| message: false -cat("## Clustering\n\n") +cat("### Clustering\n\n") ``` ```{r} @@ -254,7 +378,7 @@ cat("## Clustering\n\n") #| warning: false #| error: false #| message: false -cat("# Annotation\n\n") +cat("## Annotation\n\n") ``` ```{r} @@ -272,8 +396,8 @@ cat("# Annotation\n\n") #| warning: false #| error: false #| message: false -cat("# Analysis\n\n") -cat("## Pseudobulking\n\n") +cat("## Analysis\n\n") +cat("### Pseudobulking\n\n") ``` ```{r} @@ -291,7 +415,7 @@ cat("## Pseudobulking\n\n") #| warning: false #| error: false #| message: false -cat("## Differential expression analysis\n\n") +cat("### Differential expression analysis\n\n") ``` ```{r} @@ -309,7 +433,7 @@ cat("## Differential expression analysis\n\n") #| warning: false #| error: false #| message: false -cat("## Functional enrichment analysis: GSEA\n\n") +cat("### Functional enrichment analysis: GSEA\n\n") ``` ```{r} @@ -327,7 +451,7 @@ cat("## Functional enrichment analysis: GSEA\n\n") #| warning: false #| error: false #| message: false -cat("## Functional enrichment analysis: ORA\n\n") +cat("### Functional enrichment analysis: ORA\n\n") ``` ```{r} diff --git a/assets/integration_cluster.qmd b/assets/integration_cluster.qmd index 51aaf91..7ffb81c 100644 --- a/assets/integration_cluster.qmd +++ b/assets/integration_cluster.qmd @@ -43,10 +43,10 @@ integration_res <- as.numeric(scan(integration_res_file, character())[1]) This page shows the clustering results for your dataset post-integration. As before with the [per-sample clustering analysis](qc_cluster.qmd), the dataset was clustered at several resolutions (outlined below). A cluster tree plot is shown below to aid you in selecting the best resolution for your dataset. Recall that you want to identify a resolution where the clusters remain relatively stable from one resolution to the next, i.e. most clusters have a single, opaque arrow incoming rather than multiple faint incoming arrows. When you have selected a suitable resolution, re-run the pipeline and pass the parameter `--integrated_resolution `, replacing `` with your chosen resolution. When this parameter is unset, the pipeline choose a default value of `1`, **but you will likely need to change this**. The UMAP and count vs. feature plots below will display the results for your chosen resolution. -#### Dataset name: `{r} cohort_id` -#### Integration resolution: `{r} integration_res` +##### Dataset name: `r cohort_id` +##### Integration resolution: `r integration_res` -### Integrated clusters UMAP +#### Integrated clusters UMAP These UMAP plots show your integrated dataset, colour-coded by either the sample of origin ("orig.ident") or the cluster each cell was assigned to in the chosen integration clustering resolution ("SCT_snn_res.*"). @@ -56,7 +56,7 @@ These UMAP plots show your integrated dataset, colour-coded by either the sample cat(paste0('integrated clusters umap: ', cohort_id, '\n\n')) ``` -### Feature-count plot of integrated dataset +#### Feature-count plot of integrated dataset This plot shows the count vs. feature plots for each cluster at your chosen clustering resolution. @@ -66,7 +66,7 @@ This plot shows the count vs. feature plots for each cluster at your chosen clus cat(paste0('feature-count plot of integrated dataset: ', cohort_id, '\n\n')) ``` -### Clustree plot for integrated dataset +#### Clustree plot for integrated dataset This cluster tree shows how cells move between clusters from one resolution to the next. Use this to identify the ideal clustering resolution for your dataset. Look for resolutions where most clusters have a single, dark arrow incoming, meaning that few cells moved between clusters when moving from the previous resolution. diff --git a/assets/integration_qc.qmd b/assets/integration_qc.qmd index c110cd0..21bc8d8 100644 --- a/assets/integration_qc.qmd +++ b/assets/integration_qc.qmd @@ -45,9 +45,9 @@ We use the anchor-based CCA ("Canonical Correlation Analysis") integration metho Below you will find two UMAP plots. The first plot shows the merged but un-corrected (unintegrated) dataset. Cells are colour-coded by their sample of origin. You may notice several clusters in this plot that appear to originate primarily from a single sample. In the second plot, the integrated, corrected dataset is shown. A good batch correction should show more homogeneous clusters, with cells from all samples mixing together. You should see much fewer clusters originating from a single sample in this plot. -#### Dataset name: `{r} cohort_id` +##### Dataset name: `r cohort_id` -### Merged (unintegrated) dataset UMAP +#### Merged (unintegrated) dataset UMAP ```{r} #| results: asis @@ -55,7 +55,7 @@ Below you will find two UMAP plots. The first plot shows the merged but un-corre cat(paste0('merged dataset umap: ', cohort_id, '\n\n')) ``` -### Integrated dataset UMAP +#### Integrated dataset UMAP ```{r} #| results: asis diff --git a/assets/qc_cluster.qmd b/assets/qc_cluster.qmd index 102b473..6f7f9e4 100644 --- a/assets/qc_cluster.qmd +++ b/assets/qc_cluster.qmd @@ -94,14 +94,14 @@ purrr::walk( function(n) { p_clustree <- clustree_plot_files[[n]] p_fc_plot_list <- feature_count_plot_files[[n]] - cat(paste("###", n, "\n\n")) - cat(paste("####", "Clustree plot", "\n\n")) + cat(paste("####", n, "\n\n")) + cat(paste("#####", "Clustree plot", "\n\n")) cat(paste0('clustree plot: ', n, '\n\n')) - cat(paste("####", "Per-cluster feature-count plots", "\n\n")) + cat(paste("#####", "Per-cluster feature-count plots", "\n\n")) cat(":::panel-tabset\n\n") for (res in sort(names(p_fc_plot_list))) { p_fc_plot <- p_fc_plot_list[[res]] - cat(paste("#####", res, "\n\n")) + cat(paste("######", res, "\n\n")) cat(paste("Cluster resolution:", res, "\n\n")) cat(paste0('Per-cluster feature-count plot: ', n, '; res: ', res, '\n\n')) } diff --git a/assets/qc_doublets.qmd b/assets/qc_doublets.qmd index 347805b..422dc1d 100644 --- a/assets/qc_doublets.qmd +++ b/assets/qc_doublets.qmd @@ -118,20 +118,20 @@ purrr::walk( p_umap <- doublet_umap_plots[[n]] p_clustree <- clustree_plot_files[[n]] p_fc_plot_list <- feature_count_plot_files[[n]] - cat(paste("###", n, "\n\n")) - cat(paste("####", "Doublets per cluster", prim_res_str, "\n\n")) + cat(paste("####", n, "\n\n")) + cat(paste("#####", "Doublets per cluster", prim_res_str, "\n\n")) cat(paste0('doublet per cluster plot: ', n, '; res: ', prim_res, '\n\n')) - cat(paste("####", "Doublet proportions per cluster", prim_res_str, "\n\n")) + cat(paste("#####", "Doublet proportions per cluster", prim_res_str, "\n\n")) cat(paste0('doublet proportions per cluster plot: ', n, '; res: ', prim_res, '\n\n')) - cat(paste("####", "Doublet UMAP", prim_res_str, "\n\n")) + cat(paste("#####", "Doublet UMAP", prim_res_str, "\n\n")) cat(paste0('doublet UMAP: ', n, '; res: ', prim_res, '\n\n')) - cat(paste("####", "Clustree plot", "\n\n")) + cat(paste("#####", "Clustree plot", "\n\n")) cat(paste0('clustree plot: ', n, '\n\n')) - cat(paste("####", "Per-cluster feature-count plots", "\n\n")) + cat(paste("#####", "Per-cluster feature-count plots", "\n\n")) cat(":::panel-tabset\n\n") for (res in sort(names(p_fc_plot_list))) { p_fc_plot <- p_fc_plot_list[[res]] - cat(paste("#####", res, "\n\n")) + cat(paste("######", res, "\n\n")) cat(paste("Cluster resolution:", res, "\n\n")) cat(paste0('Per-cluster feature-count plot: ', n, '; res: ', res, '\n\n')) } diff --git a/assets/qc_filter.qmd b/assets/qc_filter.qmd index 5be0cec..7bbeecd 100644 --- a/assets/qc_filter.qmd +++ b/assets/qc_filter.qmd @@ -140,13 +140,13 @@ purrr::walk( p_log <- cdist_log_graph_files[[n]] f_pre <- pre_filter_plots[[n]] f_post <- post_filter_plots[[n]] - cat(paste("###", n, "\n\n")) - cat(paste("####", "Count distributions", "\n\n")) + cat(paste("####", n, "\n\n")) + cat(paste("#####", "Count distributions", "\n\n")) cat(paste0('count distribution plot: ', n, '\n\n')) cat(paste0('log-count distribution plot: ', n, '\n\n')) - cat(paste("####", "Features vs counts - pre-filter", "\n\n")) + cat(paste("#####", "Features vs counts - pre-filter", "\n\n")) cat(paste0('pre-filter plot: ', n, '\n\n')) - cat(paste("####", "Features vs counts - post-filter", "\n\n")) + cat(paste("#####", "Features vs counts - post-filter", "\n\n")) cat(paste0('post-filter plot: ', n, '\n\n')) } ) diff --git a/bin/annotate_clusters.R b/bin/annotate_clusters.R index 8e0b10c..9a33f53 100755 --- a/bin/annotate_clusters.R +++ b/bin/annotate_clusters.R @@ -27,6 +27,8 @@ species <- annotation_params$value[match("species", annotation_params$param)] manual_annotations_file <- annotation_params$value[match("manual_cluster_annotations", annotation_params$param)] if (!is.na(manual_annotations_file)) { cluster_cell_type_assignments <- read_csv(manual_annotations_file) + cluster_cell_type_assignments %>% + write_csv(paste0(cohort_id, ".cluster_cell_type_assignments.csv")) } else { # Get user-supplied annotation to use or set default cluster_annotation <- annotation_params$value[match("cluster_annotation", annotation_params$param)] @@ -128,5 +130,9 @@ p_cluster_cell_type_assignments_umap <- DimPlot( ) ggsave(paste0("qc_results/", cohort_id, ".umap.cluster_cell_type_assignments.png"), p_cluster_cell_type_assignments_umap) +# Save cluster cell type assignments to qc_results folder as well +cluster_cell_type_assignments %>% + write_csv(paste0("qc_results/", cohort_id, ".cluster_cell_type_assignments.csv")) + # Save pre-processed data to file SaveSeuratRds(integrated, paste0(cohort_id, ".annotated.clusters.rds")) diff --git a/bin/cluster_integrated.R b/bin/cluster_integrated.R index 7f8aa04..d3c14d2 100755 --- a/bin/cluster_integrated.R +++ b/bin/cluster_integrated.R @@ -98,5 +98,11 @@ p_fvc_integrated_clustered <- integrated@meta.data %>% ggtitle(paste0("Integrated Dataset: ", default_res_name)) ggsave(paste0("qc_results/", cohort_id, ".feature_count_plot.integrated.cluster.png"), p_fvc_integrated_clustered) +# Get number of clusters +n_clusters <- table(integrated@meta.data[default_res_name]) +n_clusters <- tibble(data.frame(n_clusters)) +colnames(n_clusters) <- c("cluster", "n_cells") +n_clusters %>% write_csv(paste0("qc_results/", cohort_id, ".integrated.clusters.csv")) + # Save Seurat object to file SaveSeuratRds(integrated, paste0(cohort_id, ".integrated.clustered.rds")) diff --git a/modules/report.nf b/modules/report.nf index 0ee6952..01fc610 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -27,16 +27,17 @@ process REPORT { path samplesheet_csv, stageAs: "samplesheets/samplesheet.csv" path custom_marker_genes_csv, stageAs: "samplesheets/custom_marker_genes.csv" path comparisons_csv, stageAs: "samplesheets/comparisons.csv" + val pipe_params output: - tuple val(cohort_name), path("report"), emit: report - tuple val(cohort_name), path("report.tar"), emit: report_tar + tuple val(cohort_name), path("report.${cohort_name}.html"), emit: report script: def int_res = (metadata.integration_resolution == null) ? '' : metadata.integration_resolution def pseudo_groups = (metadata.pseudo_groups == null) ? '' : metadata.pseudo_groups.tokenize(',').join('\n') def cluster_annotation = (metadata.cluster_annotation == null) ? '' : metadata.cluster_annotation def meta_fields = (metadata.meta_fields == null) ? '' : metadata.meta_fields.join('\n') + def param_tsv = 'parameter\tvalue\n' + pipe_params.collect { p, v -> "${p}\t${v}" }.join('\n') + '\n' """ # Save important metadata to file for reporting echo "${cohort_name}" > cohort_name.txt @@ -48,6 +49,7 @@ process REPORT { if [ -d "annotation_clusters" ]; then echo "cluster_annotation" > available_annotations/clusters.txt fi + echo -e "${param_tsv}" > parameters.tsv # Generate quarto config create_quarto_yml.py @@ -57,9 +59,6 @@ process REPORT { export XDG_CACHE_HOME="\${PWD}/.cache" # Render the quarto website - quarto render . - - # Tar the report directory - tar -cf report.tar report/ + quarto render index.qmd --output report.${cohort_name}.html """ } \ No newline at end of file From 8756f27c48ee8dd368674237b48f3120439d88a0 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Mon, 1 Jun 2026 17:25:06 +1000 Subject: [PATCH 11/18] add params to report call; WIP: add methods section --- assets/index.qmd | 6 +++++- main.nf | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 68da15b..8ff721d 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -277,7 +277,11 @@ datatable( ### Methods -TODO +You can use the following text within your methods sections of any manuscripts using the results of this analysis: + +> #### Bioinformatics and data analysis +> +> ```{r} #| echo: false diff --git a/main.nf b/main.nf index 285beb8..5f80981 100644 --- a/main.nf +++ b/main.nf @@ -433,5 +433,5 @@ workflow { .merge(analysis_ora_results) .merge(available_annotation_files) .merge(report_style) - REPORT(report_input, samplesheet_csv, custom_marker_genes_csv, comparisons_csv) + REPORT(report_input, samplesheet_csv, custom_marker_genes_csv, comparisons_csv, params) } \ No newline at end of file From 7f41423a46a3986d7a4798d7ac67130ea75572a7 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Tue, 2 Jun 2026 11:48:41 +1000 Subject: [PATCH 12/18] WIP: continue adding methods section --- assets/_quarto.yml | 12 ++++ assets/index.qmd | 118 ++++++++++++++++++++++++++++++++++++++- bin/create_quarto_yml.py | 118 --------------------------------------- main.nf | 11 +++- modules/report.nf | 13 +++-- 5 files changed, 145 insertions(+), 127 deletions(-) create mode 100644 assets/_quarto.yml delete mode 100755 bin/create_quarto_yml.py diff --git a/assets/_quarto.yml b/assets/_quarto.yml new file mode 100644 index 0000000..9981f5f --- /dev/null +++ b/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 diff --git a/assets/index.qmd b/assets/index.qmd index 8ff721d..6e49b47 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -188,6 +188,33 @@ if (file.exists("analysis_de.qmd")) { } exec_summary <- bind_rows(n_samples, n_cells_filtered, n_cells_per_sample, n_clusters, n_de_comparisons, n_de_genes) + +# Get information for methods section +# TODO: Get values from environment +pipeline_version <- scan("version.txt", character())[1] +nextflow_version <- scan("nextflow.txt", character())[1] +nextflow_command <- scan("cmd.sh", character(), sep = "\n") +nextflow_command <- paste(nextflow_command, collapse = "\n") +nextflow_command <- paste("```bash", nextflow_command, "```", sep = "\n") +software_versions <- read_csv("software.csv") +r_version <- software_versions$version[match("R", software_versions$tool)] +seurat_version <- software_versions$version[match("Seurat", software_versions$tool)] +tidyverse_ersion <- software_versions$version[match("tidyverse", software_versions$tool)] +clustree_version <- software_versions$version[match("clustree", software_versions$tool)] +doubletfinder_version <- software_versions$version[match("DoubletFinder", software_versions$tool)] +singler_version <- software_versions$version[match("SingleR", software_versions$tool)] +deseq_version <- software_versions$version[match("DESeq2", software_versions$tool)] +webgestaltr_version <- software_versions$version[match("WebGestaltR", software_versions$tool)] +hard_filtered <- TRUE # TODO +mt_qc <- !as.logical(params_df$value[match("no_mt", params_df$parameter)]) +mt_filtered <- TRUE # TODO +cluster_algorithm <- str_to_sentence(params_df$value[match("cluster_method", params_df$parameter)]) +integration_resolution <- params_df$value[match("integrated_resolution", params_df$parameter)] +annotation_species <- params_df$value[match("species", params_df$parameter)] +annotation_db_source <- "the Human Primary Cell Atlas (HPCA)" # TODO +annotation_custom_mad_threshold <- params_df$value[match("custom_annotation_mad_threshold", params_df$parameter)] +cluster_annotation <- "custom marker gene cell type" # TODO +cluster_annotation_threshold <- params_df$value[match("cell_type_proportion_threshold", params_df$parameter)] ``` ## Analysis summary: `r cohort_name` @@ -279,9 +306,94 @@ datatable( You can use the following text within your methods sections of any manuscripts using the results of this analysis: -> #### Bioinformatics and data analysis -> -> +::: {.callout-note appearance="minimal"} + +#### Bioinformatics and data analysis + +Data was processed using scrnavigator-nf `r pipeline_version` (https://github.com/Sydney-Informatics-Hub/scrnavigator-nf), a Nextflow pipeline developed by the Sydney Informatics Hub for analysing single cell RNA sequencing data. The pipeline was executed using version `r nextflow_version` of Nextflow (Di Tommaso, et al. 2017; doi: [10.1038/nbt.3820](http://www.nature.com/nbt/journal/v35/n4/full/nbt.3820.html)) with the following command: + +```{r} +#| echo: false +#| results: asis +#| warning: false +#| error: false +#| message: false +cat(nextflow_command) +``` + +The pipeline is built around the `Seurat` framework (version `r seurat_version`) and the `R` programming language (version `r r_version`). +```{r} +#| echo: false +#| results: asis +#| warning: false +#| error: false +#| message: false +if (file.exists("qc_filter.qmd") && hard_filtered) { + cat("Distributions of UMI and gene counts were inspected to select hard filtering thresholds for removal of low-quality cells.\n") + if (mt_qc && mt_filtered) { + cat("Mitochondrial gene fractions per cell were also inspected to select hard filtering thresholds for removal of potentailly dying cells.\n") + } +} +if (file.exists("qc_cluster.qmd")) { + cat(paste0("Cells were clustered using the ", cluster_algorithm, " algorithm, with clustering resolutions selected using the R package `clustree` (version ", clustree_version, ").\n")) +} + +if (file.exists("qc_doublets.qmd")) { + cat(paste0("Doublets were detected and removed using the R package `DoubletFinder` (version ", doubletfinder_version, ").\n")) +} + +if (file.exists("integration_qc.qmd")) { + cat("Samples were merged into a single dataset and integrated to correct for batch effects with the `Seurat` Canonical Correlation Analysis (CCA) algorithm.\n") +} + +if (file.exists("integration_cluster.qmd")) { + cat(paste0("The integrated dataset was re-clustered using the ", cluster_algorithm, " algorithm, with a clustering resolution of ", integration_resolution, " selected using Clustree.\n")) +} + +if (file.exists("annotation.qmd")) { + cat("\nCells were annotated as follows.\n") + if (dir.exists("annotation_cc")) { + cat("Cell cycle states were annotated using the `Seurat` built-in function `CellCycleScoring`.\n") + } + + if (dir.exists("annotation_db")) { + cat(paste0("Cell types were annotated using the `SingleR` R package (version ", singler_version, ") and a pre-built annotation database (", annotation_db_source, ").\n")) + } + + if (dir.exists("annotation_custom")) { + cat("Cells types were annotated using custom marker gene sets and the `Seurat` built-in function `AddModuleScore`.\n") + cat(paste0("The cell type with the highest ranking module score was used if the difference between the top two module scores was greater than ", annotation_custom_mad_threshold, " median absolute deviation (MAD) of all scores.\n")) + cat("Cells that did not meet this criteria were labelled 'Ambiguous'.\n") + } + + if (dir.exists("annotation_cluster")) { + cat(paste0("Clusters were annotated with the most prevalent ", cluster_annotation, " present within them, if that prevalence exceeded ", cluster_annotation_threshold, ".\n")) + cat("Clusters not meeting this criteria were labelled 'Ambiguous'.\n") + } +} + +if (dir.exists("analysis_pseudo")) { + cat("\nPseudobulking was performed to aggregate expression values from cells within each sample using `Seurat`'s built-in function `AggregateExpression`.\n") +} + +if (dir.exists("analysis_de")) { + cat(paste0("Differential expression analysis was performed on pseudobulked data using `Seurat`'s built-in function `FindMarkers` and the `DESeq2` R package (version ", deseq_version, ").\n")) +} + +if (dir.exists("analysis_gsea") || dir.exists("analysis_ora")) { + cat(paste0("Functional enrichment analysis (FEA) was performed with the `WebGestaltR` R package (version ", webgestaltr_version, "). The following FEA methods were used:\n")) + if (dir.exists("analysis_gsea")) { + cat("Gene Set Enrichment Analysis (GSEA).\n") + } + + if (dir.exists("analysis_ora")) { + cat("Over-Representation Analysis (ORA).\n") + } + +} +``` + +::: ```{r} #| echo: false diff --git a/bin/create_quarto_yml.py b/bin/create_quarto_yml.py deleted file mode 100755 index 0488c54..0000000 --- a/bin/create_quarto_yml.py +++ /dev/null @@ -1,118 +0,0 @@ -#!/usr/bin/env python3 -import yaml -import os - -VALID_QMD_FILES = { - 'Quality Control': [ - ('Filtering', 'qc_filter.qmd'), - ('Clustering', 'qc_cluster.qmd'), - ('Doublet Detection', 'qc_doublets.qmd'), - ], - 'Integration': [ - ('Integration QC', 'integration_qc.qmd'), - ('Integration Clustering', 'integration_cluster.qmd'), - ], - 'Annotation': [ - ('Annotation', 'annotation.qmd'), - ], - 'Analysis': [ - ('Pseudobulking', 'analysis_pseudo.qmd'), - ('Differential Expression', 'analysis_de.qmd'), - ('Gene Set Enrichment Analysis', 'analysis_gsea.qmd'), - ('Over-Representation Analysis', 'analysis_ora.qmd'), - ], -} - -def configure_quarto(navbar_entries: list): - return { - 'project': { - 'type': 'website', - 'output-dir': 'report', - 'render': ['*.qmd'], - }, - 'format': { - 'html': { - 'theme': { - 'light': ['flatly', 'styles.scss'], - }, - 'toc': True, - 'toc-location': 'left', - 'toc-depth': 2, - 'code-link': True, - 'code-fold': True, - 'code-line-numbers': True, - }, - }, - 'website': { - 'title': 'scRNAvigator Report', - 'navbar': { - 'background': 'primary', - 'left': navbar_entries - } - }, - } - -def construct_navbar(qmd_files: list = []): - qmd_basenames = [ - os.path.basename(f) - for f in qmd_files - ] - - # Start with the index page - navbar_list = [ - {'text': 'Home', 'href': 'index.qmd'} - ] - - # Iterate through valid nested QMD files and add them if they exist - i = 1 - for menu in ['Quality Control', 'Integration', 'Annotation', 'Analysis']: - valid_qmd_files = VALID_QMD_FILES.get(menu) - if valid_qmd_files is None: - continue - valid_qmd_files_to_add = [ - (title, f) - for title, f in valid_qmd_files - if f in qmd_basenames - ] - if len(valid_qmd_files_to_add) == 1: - # Only one page to add - add directly to navbar - title, f = valid_qmd_files_to_add[0] - navbar_list.append({ - 'text': f'{i}. {title}', - 'href': f, - }) - elif len(valid_qmd_files_to_add) > 1: - # Multiple pages to add - add a nested menu to navbar - navbar_nested_dict = { - 'text': f'{i}. {menu}', - 'menu': [] - } - j = 1 - for title, f in valid_qmd_files_to_add: - navbar_nested_dict['menu'].append({ - 'text': f'{j}. {title}', - 'href': f, - }) - j += 1 - navbar_list.append(navbar_nested_dict) - else: - continue - i += 1 - - return navbar_list - -if __name__ == '__main__': - all_files = os.listdir() - qmd_files = [ - f for f in all_files - if f.endswith('.qmd') - ] - - navbar_entries = construct_navbar(qmd_files=qmd_files) - - qmd_config = configure_quarto(navbar_entries=navbar_entries) - - config_file = '_quarto.yml' - - with open(config_file, 'w') as f: - yaml.safe_dump(qmd_config, f) \ No newline at end of file diff --git a/main.nf b/main.nf index 5f80981..6f8604b 100644 --- a/main.nf +++ b/main.nf @@ -409,6 +409,7 @@ workflow { .collect() .map { tmp -> [ tmp ] } report_style = channel.fromPath("${projectDir}/assets/styles.scss") + quarto_yaml = channel.fromPath("${projectDir}/assets/_quarto.yml", checkIfExists: true) // Gather CSV inputs custom_marker_genes_csv = !!params.custom_marker_genes ? channel.fromPath(params.custom_marker_genes, checkIfExists: true) : channel.value([]) @@ -433,5 +434,13 @@ workflow { .merge(analysis_ora_results) .merge(available_annotation_files) .merge(report_style) - REPORT(report_input, samplesheet_csv, custom_marker_genes_csv, comparisons_csv, params) + + def env_vars = [ + params:params, + pipe_version:(workflow.manifest.version ?: "0.0.0"), + nxf_version:workflow.nextflow.version, + cmd:workflow.commandLine, + software_versions:[:] + ] + REPORT(report_input, samplesheet_csv, custom_marker_genes_csv, comparisons_csv, quarto_yaml, env_vars) } \ No newline at end of file diff --git a/modules/report.nf b/modules/report.nf index 01fc610..9ac967f 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -27,7 +27,8 @@ process REPORT { path samplesheet_csv, stageAs: "samplesheets/samplesheet.csv" path custom_marker_genes_csv, stageAs: "samplesheets/custom_marker_genes.csv" path comparisons_csv, stageAs: "samplesheets/comparisons.csv" - val pipe_params + path quarto_yaml + val env_vars output: tuple val(cohort_name), path("report.${cohort_name}.html"), emit: report @@ -37,7 +38,8 @@ process REPORT { def pseudo_groups = (metadata.pseudo_groups == null) ? '' : metadata.pseudo_groups.tokenize(',').join('\n') def cluster_annotation = (metadata.cluster_annotation == null) ? '' : metadata.cluster_annotation def meta_fields = (metadata.meta_fields == null) ? '' : metadata.meta_fields.join('\n') - def param_tsv = 'parameter\tvalue\n' + pipe_params.collect { p, v -> "${p}\t${v}" }.join('\n') + '\n' + def param_tsv = 'parameter\tvalue\n' + env_vars.params.collect { p, v -> "${p}\t${v}" }.join('\n') + '\n' + def software_versions_csv = 'tool,version\n' + env_vars.software_versions.collect { p, v -> "${p},${v}" }.join('\n') + '\n' """ # Save important metadata to file for reporting echo "${cohort_name}" > cohort_name.txt @@ -50,9 +52,10 @@ process REPORT { echo "cluster_annotation" > available_annotations/clusters.txt fi echo -e "${param_tsv}" > parameters.tsv - - # Generate quarto config - create_quarto_yml.py + echo -e "${software_versions_csv}" > software.csv + echo "${env_vars.pipe_version}" > version.txt + echo "${env_vars.nxf_version}" > nextflow.txt + echo -e "${env_vars.cmd}" > cmd.sh # Create a cache dir for quarto mkdir -p .cache From ddfa4428795400ac9b79ac796a9035354399e180 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Tue, 2 Jun 2026 12:07:14 +1000 Subject: [PATCH 13/18] WIP: update methods wording - db annotation --- assets/index.qmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/assets/index.qmd b/assets/index.qmd index 6e49b47..2d262b3 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -357,7 +357,7 @@ if (file.exists("annotation.qmd")) { } if (dir.exists("annotation_db")) { - cat(paste0("Cell types were annotated using the `SingleR` R package (version ", singler_version, ") and a pre-built annotation database (", annotation_db_source, ").\n")) + cat(paste0("Cell types were annotated using the `SingleR` R package (version ", singler_version, ") and a pre-built annotation database: ", annotation_db_source, ".\n")) } if (dir.exists("annotation_custom")) { From 30adb9238f8fe3ea1d264999e34d187caf0da695 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Wed, 3 Jun 2026 14:27:50 +1000 Subject: [PATCH 14/18] Added dynamic software versions for final report --- assets/index.qmd | 84 ++++++++++++++++++++++++++++++----- main.nf | 43 +++++++++++++++--- modules/annotate_db.nf | 6 +++ modules/cluster_integrated.nf | 4 ++ modules/de.nf | 4 ++ modules/doublet.nf | 4 ++ modules/gsea.nf | 5 ++- modules/ora.nf | 4 ++ modules/report.nf | 3 +- nextflow.config | 1 + subworkflows/annotate.nf | 2 + subworkflows/integrate.nf | 1 + 12 files changed, 143 insertions(+), 18 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 2d262b3..15ce278 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -197,23 +197,62 @@ nextflow_command <- scan("cmd.sh", character(), sep = "\n") nextflow_command <- paste(nextflow_command, collapse = "\n") nextflow_command <- paste("```bash", nextflow_command, "```", sep = "\n") software_versions <- read_csv("software.csv") -r_version <- software_versions$version[match("R", software_versions$tool)] -seurat_version <- software_versions$version[match("Seurat", software_versions$tool)] -tidyverse_ersion <- software_versions$version[match("tidyverse", software_versions$tool)] +r_version <- paste(version$major, version$minor, sep = ".") +seurat_version <- as.character(packageVersion("Seurat")) +tidyverse_ersion <- as.character(packageVersion("tidyverse")) clustree_version <- software_versions$version[match("clustree", software_versions$tool)] doubletfinder_version <- software_versions$version[match("DoubletFinder", software_versions$tool)] singler_version <- software_versions$version[match("SingleR", software_versions$tool)] deseq_version <- software_versions$version[match("DESeq2", software_versions$tool)] webgestaltr_version <- software_versions$version[match("WebGestaltR", software_versions$tool)] -hard_filtered <- TRUE # TODO +celldex_version <- software_versions$version[match("celldex", software_versions$tool)] +hard_filter_cols <- c("min_ncount", "max_ncount", "min_nfeature", "max_nfeature") +hard_filter_cols <- hard_filter_cols[hard_filter_cols %in% samplesheet] +hard_filtered <- any(!is.na(samplesheet[hard_filter_cols])) mt_qc <- !as.logical(params_df$value[match("no_mt", params_df$parameter)]) -mt_filtered <- TRUE # TODO +mt_filter_cols <- c("min_mt_pct", "max_mt_pct") +mt_filter_cols <- mt_filter_cols[mt_filter_cols %in% samplesheet] +mt_filtered <- any(!is.na(samplesheet[mt_filter_cols])) cluster_algorithm <- str_to_sentence(params_df$value[match("cluster_method", params_df$parameter)]) integration_resolution <- params_df$value[match("integrated_resolution", params_df$parameter)] annotation_species <- params_df$value[match("species", params_df$parameter)] -annotation_db_source <- "the Human Primary Cell Atlas (HPCA)" # TODO +annotation_db_param <- params_df$value[match("annotation_db", params_df$parameter)] +if (is.null(annotation_db_param) || is.na(annotation_db_param) || annotation_db_param == "null" || annotation_db_param == "") { + annotation_db_source <- paste0("the Human Primary Cell Atlas (HPCA) via the R package `celldex` (version ", celldex_version, ")") +} else { + annotation_db_source <- "a pre-built annotation database" + annotation_db_source <- paste0(annotation_db_param) +} annotation_custom_mad_threshold <- params_df$value[match("custom_annotation_mad_threshold", params_df$parameter)] -cluster_annotation <- "custom marker gene cell type" # TODO +available_annotation_files <- dir( + "available_annotations", + pattern = "\\.txt$", + recursive = TRUE, + 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] +cluster_annotation_param <- params_df$value[match("cluster_annotation", params_df$parameter)] +if (is.null(cluster_annotation_param) || is.na(cluster_annotation_param) || cluster_annotation_param == "null" || cluster_annotation_param == "") { + cluster_annotation_param <- annotation_priority[1] +} +cluster_annotation_descriptions <- list( + custom_cell_type = "custom marker gene cell type", + custom_cell_type.max_score = "custom marker gene cell type", + SingleR.annotation = "cell type from the SingleR annotation", + SingleR.hpca_main = "cell type from the SingleR annotation", + SingleR.hpca_fine = "cell type from the SingleR annotation", + Phase = "cell cycle phase" +) +cluster_annotation <- cluster_annotation_descriptions[[cluster_annotation_param]] cluster_annotation_threshold <- params_df$value[match("cell_type_proportion_threshold", params_df$parameter)] ``` @@ -357,16 +396,20 @@ if (file.exists("annotation.qmd")) { } if (dir.exists("annotation_db")) { - cat(paste0("Cell types were annotated using the `SingleR` R package (version ", singler_version, ") and a pre-built annotation database: ", annotation_db_source, ".\n")) + cat(paste0("Cell types were annotated using the `SingleR` R package (version ", singler_version, ") and ", annotation_db_source, ".\n")) } if (dir.exists("annotation_custom")) { cat("Cells types were annotated using custom marker gene sets and the `Seurat` built-in function `AddModuleScore`.\n") - cat(paste0("The cell type with the highest ranking module score was used if the difference between the top two module scores was greater than ", annotation_custom_mad_threshold, " median absolute deviation (MAD) of all scores.\n")) - cat("Cells that did not meet this criteria were labelled 'Ambiguous'.\n") + if (cluster_annotation_param == "custom_cell_type.max_score") { + cat(paste0("The cell type with the highest ranking module score was used.\n")) + } else { + cat(paste0("The cell type with the highest ranking module score was used if the difference between the top two module scores was greater than ", annotation_custom_mad_threshold, " median absolute deviation (MAD) of all scores.\n")) + cat("Cells that did not meet this criteria were labelled 'Ambiguous'.\n") + } } - if (dir.exists("annotation_cluster")) { + if (dir.exists("annotation_clusters")) { cat(paste0("Clusters were annotated with the most prevalent ", cluster_annotation, " present within them, if that prevalence exceeded ", cluster_annotation_threshold, ".\n")) cat("Clusters not meeting this criteria were labelled 'Ambiguous'.\n") } @@ -577,3 +620,22 @@ cat("### Functional enrichment analysis: ORA\n\n") #| error: false #| message: false ``` + +## Software versions + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false + +tibble( + tool = c("R", "Seurat"), + version = c(r_version, seurat_version) +) %>% + bind_rows(software_versions) %>% + datatable( + colnames = c("Tool", "Version"), + rownames = FALSE + ) +``` \ No newline at end of file diff --git a/main.nf b/main.nf index 6f8604b..4bd606e 100644 --- a/main.nf +++ b/main.nf @@ -343,7 +343,7 @@ workflow { .merge(integrated_resolution) .merge(pseudo_groups) .merge(cluster_annotation) - .map { meta, int_res, pseudo -> [ meta_fields: meta, integration_resolution: int_res, pseudo_groups: pseudo, cluster_annotation: cluster_annotation ] } + .map { meta, int_res, pseudo, clust_annot -> [ meta_fields: meta, integration_resolution: int_res, pseudo_groups: pseudo, cluster_annotation: clust_annot ] } // Gather all report templates index_template = channel.fromPath("${projectDir}/assets/index.qmd", checkIfExists: true) @@ -394,7 +394,6 @@ workflow { .merge(analysis_ora_results) .filter { _t, r -> !!r } .map { t, _r -> t } - .map { t, _r -> t } report_templates = index_template .mix(qc_filter_template) .mix(qc_cluster_template) @@ -435,12 +434,46 @@ workflow { .merge(available_annotation_files) .merge(report_style) + // Get software versions for report + clustree_version = INTEGRATE.out.clustree_version + .splitText() + .first() + .map { v -> [ clustree:v ] } + doubletfinder_version = DETECT_DOUBLETS.out.version + .splitText() + .first() + .map { v -> [ DoubletFinder:v ] } + singler_version = ANNOTATE.out.singler_version + .splitText() + .first() + .map { v -> [ SingleR:v ] } + celldex_version = ANNOTATE.out.celldex_version + .splitText() + .first() + .map { v -> [ celldex:v ] } + deseq2_version = DIFFERENTIAL_EXPRESSION.out.version + .splitText() + .first() + .map { v -> [ DESeq2:v ] } + webgestaltr_version = GSEA.out.version + .mix(ORA.out.version) + .first() + .splitText() + .first() + .map { v -> [ WebGestaltR:v ] } + software_versions = clustree_version + .mix(doubletfinder_version) + .mix(singler_version) + .mix(celldex_version) + .mix(deseq2_version) + .mix(webgestaltr_version) + .reduce { acc, v -> acc + v } + def env_vars = [ params:params, pipe_version:(workflow.manifest.version ?: "0.0.0"), nxf_version:workflow.nextflow.version, - cmd:workflow.commandLine, - software_versions:[:] + cmd:workflow.commandLine ] - REPORT(report_input, samplesheet_csv, custom_marker_genes_csv, comparisons_csv, quarto_yaml, env_vars) + REPORT(report_input, samplesheet_csv, custom_marker_genes_csv, comparisons_csv, quarto_yaml, env_vars, software_versions) } \ No newline at end of file diff --git a/modules/annotate_db.nf b/modules/annotate_db.nf index ef9fef9..8deeb50 100644 --- a/modules/annotate_db.nf +++ b/modules/annotate_db.nf @@ -12,6 +12,8 @@ process ANNOTATE_DATABASE { tuple val(cohort_name), path("${cohort_name}.annotated.database.rds"), emit: annotated_rds tuple val(cohort_name), path("available_annotations.txt"), emit: available_annotations tuple val(cohort_name), path("qc_results"), emit: qc_results + path "version.singler.txt", emit: singler_version + path "version.celldex.txt", emit: celldex_version script: def annotation_csv = 'param,value\n' + @@ -27,5 +29,9 @@ process ANNOTATE_DATABASE { export XDG_CACHE_HOME="\${PWD}/.cache" annotate_db.R "${cohort_name}" "${rds_path}" annotation.csv + + # Get R package versions + Rscript -e "cat(as.character(packageVersion('SingleR')), '\\n')" > version.singler.txt + Rscript -e "cat(as.character(packageVersion('celldex')), '\\n')" > version.celldex.txt """ } \ No newline at end of file diff --git a/modules/cluster_integrated.nf b/modules/cluster_integrated.nf index b1e97f6..1d27698 100644 --- a/modules/cluster_integrated.nf +++ b/modules/cluster_integrated.nf @@ -11,6 +11,7 @@ process CLUSTER_INTEGRATED { output: tuple val(cohort_name), path("${cohort_name}.integrated.clustered.rds"), emit: integrated_rds tuple val(cohort_name), path("qc_results"), emit: qc_results + path "version.clustree.txt", emit: version script: def params_csv = 'param,value\n' + @@ -22,5 +23,8 @@ process CLUSTER_INTEGRATED { echo -e "${params_csv}" > params.csv cluster_integrated.R "${cohort_name}" "${rds_path}" params.csv + + # Get clustree version + Rscript -e "cat(as.character(packageVersion('clustree')), '\\n')" > version.clustree.txt """ } \ No newline at end of file diff --git a/modules/de.nf b/modules/de.nf index 93d7454..53fca13 100644 --- a/modules/de.nf +++ b/modules/de.nf @@ -11,9 +11,13 @@ process DIFFERENTIAL_EXPRESSION { output: tuple val(cohort_name), path("${cohort_name}.de.${test_group}.${ref_group}.Rds"), emit: de_rds tuple val(cohort_name), path("${cohort_name}.de.${test_group}.${ref_group}.csv"), val(ref_group), val(test_group), emit: de_csv + path "version.deseq.txt", emit: version script: """ de.R "${cohort_name}" "${rds_path}" "${ref_group}" "${test_group}" + + # Get DESeq2 version + Rscript -e "cat(as.character(packageVersion('DESeq2')), '\\n')" > version.deseq.txt """ } \ No newline at end of file diff --git a/modules/doublet.nf b/modules/doublet.nf index 6560173..c22b09c 100644 --- a/modules/doublet.nf +++ b/modules/doublet.nf @@ -12,6 +12,7 @@ process DETECT_DOUBLETS { tuple val(sample), path("${sample}.doublets_removed.sct_clustered.rds"), emit: doublets_removed_rds tuple val(sample), path("${sample}.doublets_detected.rds"), emit: doublets_marked_rds tuple val(sample), path("qc_results"), emit: qc_results + path "version.doubletfinder.txt", emit: version script: def mr = multiplet_rate == null ? '' : multiplet_rate @@ -30,5 +31,8 @@ process DETECT_DOUBLETS { # Once doublets are removed, re-run SCTransform and clustering sct.R "${sample}" "${sample}.doublets_removed.rds" params.csv mv "${sample}.sct_clustered.rds" "${sample}.doublets_removed.sct_clustered.rds" + + # Get clustree version + Rscript -e "cat(as.character(packageVersion('DoubletFinder')), '\\n')" > version.doubletfinder.txt """ } \ No newline at end of file diff --git a/modules/gsea.nf b/modules/gsea.nf index 511102b..37ac5cf 100644 --- a/modules/gsea.nf +++ b/modules/gsea.nf @@ -13,11 +13,14 @@ process GSEA { tuple val(cohort_name), path("${cohort_name}.gsea.${test_group}.${ref_group}.Rds"), emit: gsea_rds, optional: true tuple val(cohort_name), path("${cohort_name}.gsea.${test_group}.${ref_group}.csv"), val(ref_group), val(test_group), emit: gsea_csv, optional: true tuple val(cohort_name), path("${test_group}_vs_${ref_group}"), val(ref_group), val(test_group), emit: gsea_dir, optional: true - + path "version.webgestaltr.txt", emit: version script: def gsea_db_param = !!gsea_db_file ? gsea_db_file : '' """ gsea.R "${cohort_name}" "${de_rds_path}" "${ref_group}" "${test_group}" "${species}" "${gsea_db_param}" + + # Get DESeq2 version + Rscript -e "cat(as.character(packageVersion('WebGestaltR')), '\\n')" > version.webgestaltr.txt """ } \ No newline at end of file diff --git a/modules/ora.nf b/modules/ora.nf index a014d5d..2eb2634 100644 --- a/modules/ora.nf +++ b/modules/ora.nf @@ -13,10 +13,14 @@ process ORA { tuple val(cohort_name), path("${cohort_name}.ora.${test_group}.${ref_group}.Rds"), emit: ora_rds, optional: true tuple val(cohort_name), path("${cohort_name}.ora.${test_group}.${ref_group}.csv"), val(ref_group), val(test_group), emit: ora_csv, optional: true tuple val(cohort_name), path("${test_group}_vs_${ref_group}"), val(ref_group), val(test_group), emit: ora_dir, optional: true + path "version.webgestaltr.txt", emit: version script: def ora_db_param = !!ora_db_file ? ora_db_file : '' """ ora.R "${cohort_name}" "${de_rds_path}" "${ref_group}" "${test_group}" "${species}" "${ora_db_param}" + + # Get DESeq2 version + Rscript -e "cat(as.character(packageVersion('WebGestaltR')), '\\n')" > version.webgestaltr.txt """ } \ No newline at end of file diff --git a/modules/report.nf b/modules/report.nf index 9ac967f..6264f96 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -29,6 +29,7 @@ process REPORT { path comparisons_csv, stageAs: "samplesheets/comparisons.csv" path quarto_yaml val env_vars + val software_versions output: tuple val(cohort_name), path("report.${cohort_name}.html"), emit: report @@ -39,7 +40,7 @@ process REPORT { def cluster_annotation = (metadata.cluster_annotation == null) ? '' : metadata.cluster_annotation def meta_fields = (metadata.meta_fields == null) ? '' : metadata.meta_fields.join('\n') def param_tsv = 'parameter\tvalue\n' + env_vars.params.collect { p, v -> "${p}\t${v}" }.join('\n') + '\n' - def software_versions_csv = 'tool,version\n' + env_vars.software_versions.collect { p, v -> "${p},${v}" }.join('\n') + '\n' + def software_versions_csv = 'tool,version\n' + software_versions.collect { p, v -> "${p},${v}" }.join('\n') + '\n' """ # Save important metadata to file for reporting echo "${cohort_name}" > cohort_name.txt diff --git a/nextflow.config b/nextflow.config index 8a417e8..3ff64ba 100644 --- a/nextflow.config +++ b/nextflow.config @@ -31,6 +31,7 @@ manifest { ] name = 'scrnavigator-nf' description = 'A Nextflow pipeline for processing single cell RNA sequencing data and performing differential expression and functional enrichment analyses.' + version = '1.0.0' homePage = 'https://github.com/Sydney-Informatics-Hub/scrnavigator-nf' nextflowVersion = '!>=24' } diff --git a/subworkflows/annotate.nf b/subworkflows/annotate.nf index a56aefe..656ec62 100644 --- a/subworkflows/annotate.nf +++ b/subworkflows/annotate.nf @@ -127,4 +127,6 @@ workflow ANNOTATE { custom_annotation_qc_results = ANNOTATE_CUSTOM.out.qc_results cluster_annotation_qc_results = ANNOTATE_CLUSTERS.out.qc_results available_annotations = all_available_annotations + singler_version = ANNOTATE_DATABASE.out.singler_version + celldex_version = ANNOTATE_DATABASE.out.celldex_version } \ No newline at end of file diff --git a/subworkflows/integrate.nf b/subworkflows/integrate.nf index a516a29..fa7e408 100644 --- a/subworkflows/integrate.nf +++ b/subworkflows/integrate.nf @@ -26,4 +26,5 @@ workflow INTEGRATE { qc_results = INTEGRATION.out.qc_results cluster_qc_results = CLUSTER_INTEGRATED.out.qc_results gene_symbols = INTEGRATION.out.gene_symbols + clustree_version = CLUSTER_INTEGRATED.out.version } \ No newline at end of file From 4873c6d8f84662625078be96e0fe9eef3fb2ae3f Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Wed, 3 Jun 2026 15:25:23 +1000 Subject: [PATCH 15/18] adding citations to report methods; WIP: updating wording --- assets/index.qmd | 14 +++++++------- assets/qc_filter.qmd | 8 ++++---- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 15ce278..7c0cd3e 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -349,7 +349,7 @@ You can use the following text within your methods sections of any manuscripts u #### Bioinformatics and data analysis -Data was processed using scrnavigator-nf `r pipeline_version` (https://github.com/Sydney-Informatics-Hub/scrnavigator-nf), a Nextflow pipeline developed by the Sydney Informatics Hub for analysing single cell RNA sequencing data. The pipeline was executed using version `r nextflow_version` of Nextflow (Di Tommaso, et al. 2017; doi: [10.1038/nbt.3820](http://www.nature.com/nbt/journal/v35/n4/full/nbt.3820.html)) with the following command: +Data was processed using scrnavigator-nf `r pipeline_version` (https://github.com/Sydney-Informatics-Hub/scrnavigator-nf), a Nextflow pipeline developed by the Sydney Informatics Hub for analysing single cell RNA sequencing data. The pipeline was executed using version `r nextflow_version` of Nextflow (Di Tommaso, et al., 2017; doi: [10.1038/nbt.3820](http://www.nature.com/nbt/journal/v35/n4/full/nbt.3820.html)) with the following command: ```{r} #| echo: false @@ -360,7 +360,7 @@ Data was processed using scrnavigator-nf `r pipeline_version` (https://github.co cat(nextflow_command) ``` -The pipeline is built around the `Seurat` framework (version `r seurat_version`) and the `R` programming language (version `r r_version`). +The pipeline is built around the `Seurat` framework (version `r seurat_version`; Satija et al., 2015; doi: [10.1038/nbt.3192](https://doi.org/10.1038/nbt.3192)) and the `R` programming language (version `r r_version`). ```{r} #| echo: false #| results: asis @@ -374,11 +374,11 @@ if (file.exists("qc_filter.qmd") && hard_filtered) { } } if (file.exists("qc_cluster.qmd")) { - cat(paste0("Cells were clustered using the ", cluster_algorithm, " algorithm, with clustering resolutions selected using the R package `clustree` (version ", clustree_version, ").\n")) + cat(paste0("Cells were clustered using the ", cluster_algorithm, " algorithm, with clustering resolutions selected using the R package `clustree` (version ", clustree_version, "; Zappia and Oshlack, 2018; doi: [10.1093/gigascience/giy083](https://doi.org/10.1093/gigascience/giy083)).\n")) } if (file.exists("qc_doublets.qmd")) { - cat(paste0("Doublets were detected and removed using the R package `DoubletFinder` (version ", doubletfinder_version, ").\n")) + cat(paste0("Doublets were detected and removed using the R package `DoubletFinder` (version ", doubletfinder_version, "; McGinnis et al., 2019; doi: [10.1016/j.cels.2019.03.003](https://doi.org/10.1016/j.cels.2019.03.003)).\n")) } if (file.exists("integration_qc.qmd")) { @@ -396,7 +396,7 @@ if (file.exists("annotation.qmd")) { } if (dir.exists("annotation_db")) { - cat(paste0("Cell types were annotated using the `SingleR` R package (version ", singler_version, ") and ", annotation_db_source, ".\n")) + cat(paste0("Cell types were annotated using the `SingleR` R package (version ", singler_version, "; Aran et al., 2019; doi: [10.1038/s41590-018-0276-y](https://doi.org/10.1038/s41590-018-0276-y)) and ", annotation_db_source, ".\n")) } if (dir.exists("annotation_custom")) { @@ -420,11 +420,11 @@ if (dir.exists("analysis_pseudo")) { } if (dir.exists("analysis_de")) { - cat(paste0("Differential expression analysis was performed on pseudobulked data using `Seurat`'s built-in function `FindMarkers` and the `DESeq2` R package (version ", deseq_version, ").\n")) + cat(paste0("Differential expression analysis was performed on pseudobulked data using `Seurat`'s built-in function `FindMarkers` and the `DESeq2` R package (version ", deseq_version, "; Love et al., 2014; doi: [10.1186/s13059-014-0550-8](http://dx.doi.org/10.1186/s13059-014-0550-8)).\n")) } if (dir.exists("analysis_gsea") || dir.exists("analysis_ora")) { - cat(paste0("Functional enrichment analysis (FEA) was performed with the `WebGestaltR` R package (version ", webgestaltr_version, "). The following FEA methods were used:\n")) + cat(paste0("Functional enrichment analysis (FEA) was performed with the `WebGestaltR` R package (version ", webgestaltr_version, "; Elizarraras et al., 2024; doi: [10.1093/nar/gkae456](https://doi.org/10.1093/nar/gkae456)). The following FEA methods were used:\n")) if (dir.exists("analysis_gsea")) { cat("Gene Set Enrichment Analysis (GSEA).\n") } diff --git a/assets/qc_filter.qmd b/assets/qc_filter.qmd index 7bbeecd..800b895 100644 --- a/assets/qc_filter.qmd +++ b/assets/qc_filter.qmd @@ -80,10 +80,6 @@ filtering_params <- samplesheet %>% is_filtered <- !all(is.na(filtering_params) | filtering_params == "") ``` -The plots on this page show the distributions of UMI counts (called `nCount_RNA`) and unique gene counts (called `nFeature_RNA`) per cell, as well as the mitochondiral gene percentages (`percent.mt`) per cell. These three values are plotted in the "Features vs counts" scatter plots towards the bottom of the page. Ideally, you want to see a roughly linear relationship between these values, up to a plateau point where greater sequencing depth no longer captures additional genes. Poor quality or dying cells may appear as additional clusters of cells towards the bottom-right corners of the plots. Additionally, high mitochondrial gene percentages may be indicative of dying cells, although some cell types may show high mitochondrial gene percentages even when healthy. - -The table at the bottom of the page summarises the number of cells in each sample. If per-sample filtering was performed, the pre- and post-filter cell counts, as well as the percentage of cells filtered per sample, are shown in this table as well. - ```{r} #| echo: false #| warning: false @@ -128,6 +124,10 @@ filtering_params_non_null_cols <- filtering_params_cols[colnames(filtering_param datatable(filtering_params_non_null, colnames = unname(filtering_params_non_null_cols)) ``` +The following violin plots show the distributions of UMI counts (called `nCount_RNA`) and unique gene counts (called `nFeature_RNA`) per cell, as well as the mitochondiral gene percentages (`percent.mt`) per cell. These three values are plotted against each other in the "Features vs counts" scatter plots further down. Ideally, you want to see a roughly linear relationship between these values, up to a plateau point where greater sequencing depth no longer captures additional genes. Poor quality or dying cells may appear as additional clusters of cells towards the bottom-right corners of the plots. Additionally, high mitochondrial gene percentages may be indicative of dying cells, although be aware that some cell types may show high mitochondrial gene percentages even when healthy. + +Use these plots to identify hard-filtering thresholds to remove poor quality cells. The table below the plots summarises the number of cells in each sample and the number of cells retained after filtering. + :::panel-tabset ```{r} From 6a135fa1170fb10e635e7a1ffc7be874e7cb7d6f Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Wed, 3 Jun 2026 17:09:20 +1000 Subject: [PATCH 16/18] WIP: continue updating wording --- assets/index.qmd | 2 +- assets/integration_cluster.qmd | 5 ++--- assets/integration_qc.qmd | 2 -- assets/qc_cluster.qmd | 41 +++++++++++++++++++++++++++++++++- assets/qc_doublets.qmd | 8 +++---- 5 files changed, 46 insertions(+), 12 deletions(-) diff --git a/assets/index.qmd b/assets/index.qmd index 7c0cd3e..094b12e 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -464,7 +464,7 @@ cat("### Initial filtering\n\n") #| warning: false #| error: false #| message: false -cat("### Clustering\n\n") +cat("### Sample clustering\n\n") ``` ```{r} diff --git a/assets/integration_cluster.qmd b/assets/integration_cluster.qmd index 7ffb81c..5a313a1 100644 --- a/assets/integration_cluster.qmd +++ b/assets/integration_cluster.qmd @@ -41,10 +41,9 @@ integration_res_file <- "integration_resolution.txt" integration_res <- as.numeric(scan(integration_res_file, character())[1]) ``` -This page shows the clustering results for your dataset post-integration. As before with the [per-sample clustering analysis](qc_cluster.qmd), the dataset was clustered at several resolutions (outlined below). A cluster tree plot is shown below to aid you in selecting the best resolution for your dataset. Recall that you want to identify a resolution where the clusters remain relatively stable from one resolution to the next, i.e. most clusters have a single, opaque arrow incoming rather than multiple faint incoming arrows. When you have selected a suitable resolution, re-run the pipeline and pass the parameter `--integrated_resolution `, replacing `` with your chosen resolution. When this parameter is unset, the pipeline choose a default value of `1`, **but you will likely need to change this**. The UMAP and count vs. feature plots below will display the results for your chosen resolution. +This section shows the clustering results for your dataset post-integration. As before with the [per-sample clustering analysis](#sample-clustering), the dataset was clustered at several resolutions (outlined below). A cluster tree plot is shown below to aid you in selecting the best resolution for your dataset. Recall that you want to identify a resolution where the clusters remain relatively stable from one resolution to the next, i.e. most clusters have a single, opaque arrow incoming rather than multiple faint incoming arrows. When you have selected a suitable resolution, re-run the pipeline and pass the parameter `--integrated_resolution `, replacing `` with your chosen resolution. When this parameter is unset, the pipeline choose a default value of `1`, **but you will likely need to change this**. The UMAP and count vs. feature plots below will display the results for your chosen resolution. -##### Dataset name: `r cohort_id` -##### Integration resolution: `r integration_res` +##### Current integration resolution: `r integration_res` #### Integrated clusters UMAP diff --git a/assets/integration_qc.qmd b/assets/integration_qc.qmd index 21bc8d8..f1cc1e3 100644 --- a/assets/integration_qc.qmd +++ b/assets/integration_qc.qmd @@ -45,8 +45,6 @@ We use the anchor-based CCA ("Canonical Correlation Analysis") integration metho Below you will find two UMAP plots. The first plot shows the merged but un-corrected (unintegrated) dataset. Cells are colour-coded by their sample of origin. You may notice several clusters in this plot that appear to originate primarily from a single sample. In the second plot, the integrated, corrected dataset is shown. A good batch correction should show more homogeneous clusters, with cells from all samples mixing together. You should see much fewer clusters originating from a single sample in this plot. -##### Dataset name: `r cohort_id` - #### Merged (unintegrated) dataset UMAP ```{r} diff --git a/assets/qc_cluster.qmd b/assets/qc_cluster.qmd index 6f7f9e4..43e308a 100644 --- a/assets/qc_cluster.qmd +++ b/assets/qc_cluster.qmd @@ -62,7 +62,46 @@ To aid in this process, we have clustered your data at multiple resolutions (see A good rule of thumb for selecting a clustering resolution for your data is to look for the resolutions at which the clusters remain relatively stable, i.e. there is little movement of cells between clusters. In the cluster trees this looks like rows of clusters where most of the clusters have a single, solid arrow incoming from the previous resolution. In contrast, unstable resolutions will have clusters with many, faint arrows incoming, representing clusters that are splitting and merging with each other from one resolution to the next. -Below, you will find cluster trees for each of your samples, plus count vs. feature scatter plots for each of the clusters at each resolution. +Below, you will find cluster trees for each of your samples, plus count vs. feature scatter plots for each of the clusters at each resolution. Use these plots to select a clustering resolution for each of your samples and set these in the samplesheet using the `res` column. + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false + +cluster_res <- read_csv("samplesheets/samplesheet.csv") +if ("res" %in% colnames(cluster_res) && any(!is.na(cluster_res$res))) { + cluster_res <- cluster_res %>% + select(sample, res) +} else { + cluster_res <- NULL +} +``` + +```{r} +#| echo: false +#| eval: !expr (!is.null(cluster_res)) +#| results: asis +#| warning: false +#| error: false +#| message: false + +cat("You have already set the following clustering resolutions for you sample. If you are happy with these values, you can proceed to the next step: [doublet detection and filtering](#doublet-detection-and-filtering).") +``` + +```{r} +#| echo: false +#| eval: !expr (!is.null(cluster_res)) +#| warning: false +#| error: false +#| message: false + +datatable( + cluster_res, + colnames = c("Sample", "Sample clustering resolution") +) +``` #### Clustering resolutions diff --git a/assets/qc_doublets.qmd b/assets/qc_doublets.qmd index 422dc1d..0d2c467 100644 --- a/assets/qc_doublets.qmd +++ b/assets/qc_doublets.qmd @@ -95,13 +95,11 @@ for (n in names(feature_count_plot_files)) { } ``` -Doublets are cases of two cells being captured together during the sample preparation phase, causing both cells to be barcoded with the same DNA barcode. This results in the individual cells' expression profiles being merged together. These doublets (or multiplets in the case of three or more cells being captured together) don't represent real cells and will confound donwstream analyses. As such, they should be removed. +Doublets are cases of two cells being captured together during the sample preparation phase, causing both cells to be barcoded with the same DNA barcode. This results in the individual cells' expression profiles being merged together. These doublets (or multiplets in the case of three or more cells being captured together) don't represent real cells and will confound donwstream analyses. As such, they should be removed. We have used the R package `DoubletFinder` to identify barcodes that are likely doublets and removed them from the dataset. -The doublet detection stage of the workflow uses the R package `DoubletFinder` to compare the expression profile of each barcode against artificial doublets created by averaging the expression profiles of randomly-chosen pairs of barcodes. Barcodes with similar expression profiles to these artificial doublets are considered real doublets and are removed from the dataset. +Below you will find several plots for each sample illustrating the doublets identified. At the top are two bar graphs showing the number and proportion of doublets identified per cluster, respectively. You will also find a UMAP plot of your single cell data with doublets highlighted in red to show which clusters were affected. Below these plots you will find clustering plots similar to [those above](#sample-clustering), including a cluster tree plot and multiple count vs. feature plots for each cluster of cells at each clustering resolution. This clustering analysis has been performed after doublet removal and demonstrate the effect of doublet removal on the clustering results. -Below you will find several plots for each sample illustrating the doublets identified. At the top are two bar graphs showing the number and proportion of doublets identified per cluster, respectively. You will also find a UMAP plot of your single cell data with doublets highlighted in red to show which clusters were affected. Below these plots you will find clustering plots similar to those in the previous [clustering report](qc_cluster.qmd), including a cluster tree plot and multiple count vs. feature plots for each cluster of cells at each clustering resolution. This clustering analysis has been performed after doublet removal and demonstrate the effect of doublet removal on the clustering results. - -Finally, at the bottom of the page, you will find a summary table with the number and proportion of doublets identified per sample. +Finally, below the plots you will find a summary table with the number and proportion of doublets identified per sample. :::panel-tabset From 4989d06add1557bed422ef5e798e787eb57961cb Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Tue, 9 Jun 2026 12:13:28 +1000 Subject: [PATCH 17/18] continue updating explanatory text: up to analysis --- assets/annotation.qmd | 61 +++++++++++++++++++++------------- assets/integration_cluster.qmd | 6 ++-- assets/qc_cluster.qmd | 20 +++++------ assets/qc_filter.qmd | 6 ++-- 4 files changed, 53 insertions(+), 40 deletions(-) diff --git a/assets/annotation.qmd b/assets/annotation.qmd index 04f1416..d082c73 100644 --- a/assets/annotation.qmd +++ b/assets/annotation.qmd @@ -151,50 +151,65 @@ if (is.null(cohort_id)) { } ``` -There are three main annotation options you can choose when running the scRNAvigator pipeline: +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. This annotation is automatic for human samples and can be enabled for non-human samples by providing lists of cell cycle marker genes for your species using the `--s_genes` and `--g2m_genes` parameters. +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, although you can override the database used with the `--annotation_db` parameter. For non-human species, you will need to use this parameter and provide a database file in order to perform automatic cell type annotation. +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. You provide this list of cell types using the `--custom_marker_genes` parameter. The file provided must be a CSV file and contain an annotation label and a list of genes associated with it. Cells that score highly for more than one cell type will be labelled "Ambiguous". +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". By default, a cell type must make up 2/3 of the cluster for that label to be applied to the cluster as well. You can choose the threshold used with the `--cell_type_proportion_threshold` parameter. - -By default, the annotation category used for cluster annotation is the first in the following list that has been performed: - -- Custom marker gene-based annotation -- Automatic cell type annotation - - For human HPCA annotations, the 'main' level annotation is prioritised over the 'fine' level annotation -- Cell cycle annotation - -The annotation category used for cluster annotation can be overridden using the `--cluster_annotation`. See below for a list of valid values from your dataset. - -##### Dataset name: `r cohort_id` +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 -The following annotation fields are available for cluster annotation: +Your cells were annotated with the following annotation fields: ```{r} #| echo: false -#| results: asis +#| 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) { - cat("- ") - cat(ann) - cat("\n") + all_annotations_df <- all_annotations_df %>% + bind_rows(tibble( + annotation_field = ann, + description = annotation_descriptions[[ann]] + )) } ``` -The field used for cluster annotation was: **`r annotation_used_for_clusters`** +```{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 analyses -The following list contains all of the annotations listed above, as well as `cluster_annotation` and additional metadata fields you provided in the samplesheet. These are all of the fields available to you for performing downstream pseudobulking and subsequent differential expression and functional enrichment 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 diff --git a/assets/integration_cluster.qmd b/assets/integration_cluster.qmd index 5a313a1..23ba76f 100644 --- a/assets/integration_cluster.qmd +++ b/assets/integration_cluster.qmd @@ -41,9 +41,7 @@ integration_res_file <- "integration_resolution.txt" integration_res <- as.numeric(scan(integration_res_file, character())[1]) ``` -This section shows the clustering results for your dataset post-integration. As before with the [per-sample clustering analysis](#sample-clustering), the dataset was clustered at several resolutions (outlined below). A cluster tree plot is shown below to aid you in selecting the best resolution for your dataset. Recall that you want to identify a resolution where the clusters remain relatively stable from one resolution to the next, i.e. most clusters have a single, opaque arrow incoming rather than multiple faint incoming arrows. When you have selected a suitable resolution, re-run the pipeline and pass the parameter `--integrated_resolution `, replacing `` with your chosen resolution. When this parameter is unset, the pipeline choose a default value of `1`, **but you will likely need to change this**. The UMAP and count vs. feature plots below will display the results for your chosen resolution. - -##### Current integration resolution: `r integration_res` +This section shows the clustering results for your dataset post-integration. As before with the [per-sample clustering analysis](#sample-clustering), the dataset was clustered at several resolutions (outlined below). A cluster tree plot is shown below to aid you in selecting the best resolution for your dataset. Recall that we want to identify a resolution where the clusters remain relatively stable from one resolution to the next, i.e. most clusters have a single, opaque arrow incoming rather than multiple faint incoming arrows. By default, the pipeline will choose a value of `1`, **but often this will need to be altered to suit the dataset**. For your dataset, a value of **`r integration_res`** was selected. The UMAP and count vs. feature plots below will display the clustering results for your chosen resolution. #### Integrated clusters UMAP @@ -67,7 +65,7 @@ cat(paste0('feature-count plot of integrated #### Clustree plot for integrated dataset -This cluster tree shows how cells move between clusters from one resolution to the next. Use this to identify the ideal clustering resolution for your dataset. Look for resolutions where most clusters have a single, dark arrow incoming, meaning that few cells moved between clusters when moving from the previous resolution. +This cluster tree shows how cells move between clusters from one resolution to the next. Stable resolutions are where most clusters have a single, dark incoming arrow, meaning that few cells moved between clusters when moving from the previous resolution. ```{r} #| results: asis diff --git a/assets/qc_cluster.qmd b/assets/qc_cluster.qmd index 43e308a..6b5a82e 100644 --- a/assets/qc_cluster.qmd +++ b/assets/qc_cluster.qmd @@ -62,43 +62,41 @@ To aid in this process, we have clustered your data at multiple resolutions (see A good rule of thumb for selecting a clustering resolution for your data is to look for the resolutions at which the clusters remain relatively stable, i.e. there is little movement of cells between clusters. In the cluster trees this looks like rows of clusters where most of the clusters have a single, solid arrow incoming from the previous resolution. In contrast, unstable resolutions will have clusters with many, faint arrows incoming, representing clusters that are splitting and merging with each other from one resolution to the next. -Below, you will find cluster trees for each of your samples, plus count vs. feature scatter plots for each of the clusters at each resolution. Use these plots to select a clustering resolution for each of your samples and set these in the samplesheet using the `res` column. - ```{r} #| echo: false #| warning: false #| error: false #| message: false -cluster_res <- read_csv(% +cluster_res_df <- read_csv("samplesheets/samplesheet.csv") +if ("res" %in% colnames(cluster_res_df) && any(!is.na(cluster_res_df$res))) { + cluster_res_df <- cluster_res_df %>% select(sample, res) } else { - cluster_res <- NULL + cluster_res_df <- NULL } ``` ```{r} #| echo: false -#| eval: !expr (!is.null(cluster_res)) +#| eval: !expr (!is.null(cluster_res_df)) #| results: asis #| warning: false #| error: false #| message: false -cat("You have already set the following clustering resolutions for you sample. If you are happy with these values, you can proceed to the next step: [doublet detection and filtering](#doublet-detection-and-filtering).") +cat("The following clustering resolutions have been set for each sample:") ``` ```{r} #| echo: false -#| eval: !expr (!is.null(cluster_res)) +#| eval: !expr (!is.null(cluster_res_df)) #| warning: false #| error: false #| message: false datatable( - cluster_res, + cluster_res_df, colnames = c("Sample", "Sample clustering resolution") ) ``` @@ -123,6 +121,8 @@ cluster_res_str <- paste0(cluster_res_str, "\n") cat(cluster_res_str) ``` +Below, you will find cluster trees for each of your samples, plus feature vs. count scatter plots for each of the clusters at each resolution. + :::panel-tabset ```{r} diff --git a/assets/qc_filter.qmd b/assets/qc_filter.qmd index 800b895..658164b 100644 --- a/assets/qc_filter.qmd +++ b/assets/qc_filter.qmd @@ -124,9 +124,7 @@ filtering_params_non_null_cols <- filtering_params_cols[colnames(filtering_param datatable(filtering_params_non_null, colnames = unname(filtering_params_non_null_cols)) ``` -The following violin plots show the distributions of UMI counts (called `nCount_RNA`) and unique gene counts (called `nFeature_RNA`) per cell, as well as the mitochondiral gene percentages (`percent.mt`) per cell. These three values are plotted against each other in the "Features vs counts" scatter plots further down. Ideally, you want to see a roughly linear relationship between these values, up to a plateau point where greater sequencing depth no longer captures additional genes. Poor quality or dying cells may appear as additional clusters of cells towards the bottom-right corners of the plots. Additionally, high mitochondrial gene percentages may be indicative of dying cells, although be aware that some cell types may show high mitochondrial gene percentages even when healthy. - -Use these plots to identify hard-filtering thresholds to remove poor quality cells. The table below the plots summarises the number of cells in each sample and the number of cells retained after filtering. +The following violin plots show the distributions of UMI counts (called `nCount_RNA`) and unique gene counts (called `nFeature_RNA`) per cell, as well as the mitochondiral gene percentages (`percent.mt`) per cell. These three values are plotted against each other in the "Features vs counts" scatter plots further down. These plots can be used together to identify hard-filtering thresholds for removing poor quality cells. Ideally, you want to see a roughly linear relationship between these values, up to a plateau point where greater sequencing depth no longer captures additional genes. Poor quality or dying cells may appear as additional clusters of cells towards the bottom-right corners of the plots. Additionally, high mitochondrial gene percentages may be indicative of dying cells, although be aware that some cell types may show high mitochondrial gene percentages even when healthy. :::panel-tabset @@ -156,6 +154,8 @@ purrr::walk( #### Filtering summary +The following table summarises the number of cells in each sample and the number of cells retained after filtering. + ```{r} #| echo: false datatable( From d9fc639c8f9362d7a1d44e51d88858cb19cd0f57 Mon Sep 17 00:00:00 2001 From: mgeaghan Date: Tue, 9 Jun 2026 13:54:35 +1000 Subject: [PATCH 18/18] Finish updating explanatory text --- assets/analysis_de.qmd | 40 ++++++++++++++------------------------ assets/analysis_gsea.qmd | 18 ++--------------- assets/analysis_ora.qmd | 18 ++--------------- assets/analysis_pseudo.qmd | 22 ++++----------------- assets/annotation.qmd | 2 +- assets/index.qmd | 19 ++++++++++++++++-- 6 files changed, 41 insertions(+), 78 deletions(-) diff --git a/assets/analysis_de.qmd b/assets/analysis_de.qmd index 0617f01..1f9c8ef 100644 --- a/assets/analysis_de.qmd +++ b/assets/analysis_de.qmd @@ -83,30 +83,22 @@ 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") ``` -This page shows the differential expression results for the comparisons you requested on your dataset. 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) - -##### Dataset name: `r cohort_id` - -#### All comparisons made: - -You provided the following table of comparisons to perform via the `--comparisons` pipeline parameter: +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} datatable(comparisons_df, colnames = c("Reference group", "Test group")) ``` -*Note that the group names above must match the [pseudobulk groups shown in the previous section of the report](analysis_pseudo.qmd#pseudobulked-groups-availablle-for-comparison).* +The results are presented in four sections: -#### Differential expression results +- [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 +#### 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 shoudl 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. +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 @@ -114,13 +106,11 @@ The following plots show the distribution of p-values for each comparison made. cat(paste0('p-value distributions: ', cohort_id, '\n\n')) ``` -##### Log-fold-change distributions +#### 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"). -You can use these plots to determine whether you want to apply an additional fold-change threshold to your data. Such a threshold can help to remove noise from very lowly perturbed genes that may not represent a significant biological change. You can apply a fold-change threshold with the `--fc_threshold` parameter. **Note** that this parameter takes a **positive, non-log-transformed value**; i.e., to filter out all results with less than a 1.5x change in either direction (equivalent to ±0.58 in these log-scaled plots), you would provide `--fc_threshold 1.5`. - -Keep in mind that in many biological contexts, small fold-changes can be very biologically significant, so only use this parameter if you are sure you want to discard low fold-change results. +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 @@ -128,9 +118,9 @@ Keep in mind that in many biological contexts, small fold-changes can be very bi cat(paste0('log-fold-change distributions: ', cohort_id, '\n\n')) ``` -##### Volcano plots +#### 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 you have provided a fold-change threshold, 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. +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 @@ -138,13 +128,13 @@ Volcano plots combine both the p-value and fold-change information into one very cat(paste0('volcano plots: ', cohort_id, '\n\n')) ``` -##### Differential expression table +#### 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 +##### Significant results only ```{r} #| echo: false @@ -154,7 +144,7 @@ Finally, below is a table summarising all of the differential expression results datatable(sig_de_table, colnames = de_table_col_names_ordered) ``` -###### All results +##### All results ```{r} #| echo: false diff --git a/assets/analysis_gsea.qmd b/assets/analysis_gsea.qmd index e9ebc61..099a560 100644 --- a/assets/analysis_gsea.qmd +++ b/assets/analysis_gsea.qmd @@ -88,23 +88,9 @@ gsea_plots <- sapply(all_comparisons_str, function(cmp) { }) ``` -Differential expression analyses can generate lots of results that can be quite difficult to interpret by themselves. Functional enrichment analyses are a set of methods designed to take the raw differential expression results and identify patterns of expression that point to altered **biological function**. There are several such approaches, but the two we apply in this pipeline are **Gene Set Enrichment Analysis** (GSEA) and **Over-Representation Analysis** (ORA). +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. -This page 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. - -##### Dataset name: `r cohort_id` - -#### All comparisons made: - -One GSEA analysis was conducted for each differential expression comparison that was made. As per the [differential expression page](analysis_de.qmd), these comparisons were: - -```{r} -datatable(comparisons_df, colnames = c("Reference group", "Test group")) -``` - -#### GSEA expression results - -Below, you will find plots for each GSEA analysis. These plots show the enrichment score for each **significant** pathway from the analyses; more upregulated pathways have bars pointing to the right and are highlighted in orange; downregulated pathways point to the left and are highlighted in blue. +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 diff --git a/assets/analysis_ora.qmd b/assets/analysis_ora.qmd index 057049a..4f255da 100644 --- a/assets/analysis_ora.qmd +++ b/assets/analysis_ora.qmd @@ -88,23 +88,9 @@ ora_plots <- sapply(all_comparisons_str, function(cmp) { }) ``` -Differential expression analyses can generate lots of results that can be quite difficult to interpret by themselves. Functional enrichment analyses are a set of methods designed to take the raw differential expression results and identify patterns of expression that point to altered **biological function**. There are several such approaches, but the two we apply in this pipeline are **Gene Set Enrichment Analysis** (GSEA) and **Over-Representation Analysis** (ORA). +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. -This page 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. - -##### Dataset name: `r cohort_id` - -#### All comparisons made: - -One GSEA analysis was conducted for each differential expression comparison that was made. As per the [differential expression page](analysis_de.qmd), these comparisons were: - -```{r} -datatable(comparisons_df, colnames = c("Reference group", "Test group")) -``` - -#### ORA expression results - -Below, you will find plots for each ORA analysis. These plots show 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; depleted or under-represented pathways point to the left and are highlighted in blue. +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 diff --git a/assets/analysis_pseudo.qmd b/assets/analysis_pseudo.qmd index d3ac4d9..0c69bae 100644 --- a/assets/analysis_pseudo.qmd +++ b/assets/analysis_pseudo.qmd @@ -37,27 +37,13 @@ pseudo_groups_file <- "pseudo_groups.txt" pseudo_groups <- scan(pseudo_groups_file, character()) ``` -When we run a differential expression analysis on a single cell RNA-seq dataset, we are constructing a statistical model that looks for changes in the mean expression of each gene between case and control samples. These statistical models are built upon several **assumptions** about the underlying data, and if our dataset violates those assumptions, the results can be incorrect. +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. -One of the major assumptions that these models make is that each data point is an independent, random sample drawn from some population. In the case of a bulk RNA-seq experiment, this assumption is valid: each data point is an RNA sample from an individual person or organism, and those individuals were (hopefully) chosen randomly from their respective populations. In the case of a single-cell RNA-seq experiment though, things get a bit more complicated. - -If we naively pass our scRNA-seq data into a standard differential expression algorithm, it will treat **each cell** as a separate data point and make the **assumption** that each cell is independently sampled from the population. But that isn't the case: cells from the same individual are **not** independent of one another - they originate from the same individual and share the exact same genetic background. This means the cells' expression profiles are correlated with one another, and this will invalidate the independence assumption, bias our results, and ultimately lead to over-inflated test statistics and estimates of significance. This issue is called **pseudoreplication**. - -There are a few ways to deal with pseudoreplication, but the method we have implemented in this pipeline is called **pseudobulking**. Pseudobulking 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. Once complete, the B cell profiles from each sample would be suitably independent of one another and could be compared between cases and controls; the same could be done for the T cell profiles. - -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. - -When the pseudobulking stage of the pipeline runs, it uses the metadata fields provided with the `--pseudo_groups` parameter to decide how to group the cells prior to aggregation. By default, the pipeline will use the `cluster_annotation` field generated in the [annotation stage](annotation.qmd). This means all cells with the same cluster annotation within a sample will be aggregated together. You can provide multiple grouping fields for to this parameter, separated by commas, to refine the groups; for example, if you have a sample of B and T cells and ran both cluster annotation and cell cycle annotation, `--pseudo_groups cluster_annotation,Phase` would result in six pseudobulked expression profiles per cell: S-phase B cells, S-phase T cells, G1-phase B cells, etc. - -See the section titled [All annotation fields available for downstream analyses](annotation.qmd#all-annotation-fields-available-for-downstream-analyses) for a full list of fields relevant to your data that you can provide to `--pseudo_groups`. **Note** that you should include sample-level metadata fields here that are relevant to your differential expression analysis; for example, if you are comparing a particular treatment vs a control, and the relevant metadata field was called `treatment`, you would provide `treatment` to the `--pseudo_groups` parameter, along with any other fields you need. - -The pseudobulked expression profiles will be named after the individual levels within each metadata field. Again, taking the example of B and T cells in a treatment vs control experiment, if the cell labels are `B-cell` and `T-cell` and the labels in the `treatment` field are `drug` and `control`, your pseudobulked groups will be called `B-cell_drug`, `B-cell_control`, `T-cell_drug`, and `T-cell_control`. - -##### Dataset name: `r cohort_id` +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 were provided to the `--pseudo_groups` parameter and used for aggregation of cells: +The following [metadata fields](#all-annotation-fields-available-for-downstream-analysis) were used for aggregation of cells: ```{r} #| echo: false @@ -71,7 +57,7 @@ for (g in pseudo_groups) { ### Pseudobulked groups availablle for comparison -Below is a table summarising all your final pseudobulk groups and the number of samples in your dataset that are represented in each group. You will need to note down these group names to construct a list of your desired differential expression comparisons. **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. +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 diff --git a/assets/annotation.qmd b/assets/annotation.qmd index d082c73..27156e8 100644 --- a/assets/annotation.qmd +++ b/assets/annotation.qmd @@ -207,7 +207,7 @@ datatable( Cell **clusters** were annotated with the labels from the **`r annotation_used_for_clusters`** annotation field. -#### All annotation fields available for downstream analyses +#### All annotation fields available for downstream analysis 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. diff --git a/assets/index.qmd b/assets/index.qmd index 094b12e..c7be38b 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -585,6 +585,21 @@ cat("### Differential expression analysis\n\n") #| message: false ``` +```{r} +#| echo: false +#| results: asis +#| eval: !expr (file.exists("analysis_gsea.qmd") || file.exists("analysis_ora.qmd")) +#| warning: false +#| error: false +#| message: false + +cat("### Functional enrichment analysis (FEA)\n\n") + +cat("Differential expression analyses can generate lots of results that can be quite difficult to interpret by themselves.\n") +cat("Functional enrichment analyses are a set of methods designed to take the raw differential expression results and identify patterns of expression that point to altered **biological function**.\n") +cat("There are several such approaches, but the two we apply in this pipeline are **Gene Set Enrichment Analysis** (GSEA) and **Over-Representation Analysis** (ORA).\n\n") +``` + ```{r} #| echo: false #| results: asis @@ -592,7 +607,7 @@ cat("### Differential expression analysis\n\n") #| warning: false #| error: false #| message: false -cat("### Functional enrichment analysis: GSEA\n\n") +cat("#### Gene Set Enrichment Analysis (GSEA)\n\n") ``` ```{r} @@ -610,7 +625,7 @@ cat("### Functional enrichment analysis: GSEA\n\n") #| warning: false #| error: false #| message: false -cat("### Functional enrichment analysis: ORA\n\n") +cat("#### Over-Representation Analysis (ORA)\n\n") ``` ```{r}