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/analysis_de.qmd b/assets/analysis_de.qmd index a806266..1f9c8ef 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,23 @@ p_dist_plot_file <- paste0("analysis_de/results/", cohort_id, ".p_val_dist.png") log_fc_dist_plot_file <- paste0("analysis_de/results/", cohort_id, ".log_fc_dist.png") ``` -#### Dataset name: `{r} cohort_id` - -#### All comparisons made: +This section shows the differential expression (DE) results for the comparisons made for your dataset. Each DE analysis is a simple case vs. control comparison of two [pseudobulk groups](#pseudobulked-groups-availablle-for-comparison) from your dataset. The following table lists the DE comparisons that were made: ```{r} -#| echo: false -#| results: asis -for (cmp in all_comparisons) { - cat("- ") - cat(cmp) - cat("\n") -} +datatable(comparisons_df, colnames = c("Reference group", "Test group")) ``` -### Differential expression results +The results are presented in four sections: + +- [p-value distribution results](#p-value-distributions) +- [Log-fold change distribution results](#log-fold-change-distributions) +- [Volcano plots](#volcano-plots) +- [Differential expression results table](#differential-expression-table) #### p-value distributions +The following plots show the distribution of p-values for each comparison made. These can be useful to determine the quality of the results: a well-powered experiment with true positive results should show an enrichment of low p-values (i.e. towards the left side of each plot). Uniformly distributed p-values indicate no true positive results or low power to detect true positives. + ```{r} #| results: asis #| echo: false @@ -105,6 +108,10 @@ cat(paste0('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") }) ``` -#### Dataset name: `{r} cohort_id` - -#### All comparisons made: - -```{r} -#| echo: false -#| results: asis -for (cmp in all_comparisons) { - cat("- ") - cat(cmp) - cat("\n") -} -``` +This section presents the results from the GSEA analysis conducted on your dataset. GSEA works by first ranking the entire set of differential expression results by their fold-change values; this includes **both** significant **and** non-significant genes. The method then looks at a database of **biological pathways** and the **gene sets** associated with each. For each gene set, GSEA assesses whether these genes appear clustered together within your ranked dataset, more so than would be expected by chance. The degree of clustering towards the upregulated or downregulated ends of your dataset is used to calculate an **enrichment score** for that gene set. -### GSEA expression results +One GSEA analysis was conducted for [each differential expression comparison](#differential-expression-analysis) that was made. Below, you will find plots showing the enrichment score for each **significant** pathway from each analysis: more upregulated pathways have bars pointing to the right and are highlighted in orange, while downregulated pathways point to the left and are highlighted in blue. :::panel-tabset @@ -108,11 +100,11 @@ purrr::walk( sort(all_comparisons), function(n) { p <- gsea_plots[[n]] - cat(paste("####", n, "\n\n")) - cat(paste("#####", "Gene Set Enrichment Analysis Results:", n)) + cat(paste("#####", n, "\n\n")) + cat(paste("######", "Gene Set Enrichment Analysis Results:", n)) cat(paste0('gsea results: ', n, '\n\n')) } ) ``` -::: \ No newline at end of file +::: diff --git a/assets/analysis_ora.qmd b/assets/analysis_ora.qmd index b227ea3..4f255da 100644 --- a/assets/analysis_ora.qmd +++ b/assets/analysis_ora.qmd @@ -78,27 +78,19 @@ all_comparisons_str <- unique_comparisons %>% .$comparison names(all_comparisons_str) <- all_comparisons +# Get comparisons CSV +comparisons_csv <- "samplesheets/comparisons.csv" +comparisons_df <- read_csv(comparisons_csv) + # Get ORA plots ora_plots <- sapply(all_comparisons_str, function(cmp) { paste0("analysis_ora/plots/", cohort_id, ".ora.", cmp, ".png") }) ``` -#### Dataset name: `{r} cohort_id` - -#### All comparisons made: - -```{r} -#| echo: false -#| results: asis -for (cmp in all_comparisons) { - cat("- ") - cat(cmp) - cat("\n") -} -``` +This section presents the results from the ORA analysis conducted on your dataset. ORA works by taking the set of your **significantly differentially expressed** genes and comparing that against each set of genes mapped to **biological pathways** in a curated database. For each pathway, ORA looks at the proportion of those genes that appear in your list of significant genes. If that proportion is higher than what would be expected by chance, the gene set is considered **over-represented** or **enriched** in your dataset. In contrast, some pathways may have very low representation in your dataset, in which case they can be considered **under-represented** or **depleted**. The ORA method generates a list of **enrichment ratios** for each pathway. -### ORA expression results +One ORA analysis was conducted for [each differential expression comparison](#differential-expression-analysis) that was made. Below, you will find plots showing the log2-transformed *enrichment ratio* for each **significant** pathway from the analysis: more enriched pathways have bars pointing to the right and are highlighted in orange, while depleted or under-represented pathways point to the left and are highlighted in blue. :::panel-tabset @@ -108,8 +100,8 @@ purrr::walk( sort(all_comparisons), function(n) { p <- ora_plots[[n]] - cat(paste("####", n, "\n\n")) - cat(paste("#####", "Over-Representation Analysis Results:", n)) + cat(paste("#####", n, "\n\n")) + cat(paste("######", "Over-Representation Analysis Results:", n)) cat(paste0('ora results: ', n, '\n\n')) } ) diff --git a/assets/analysis_pseudo.qmd b/assets/analysis_pseudo.qmd index e96eb6a..0c69bae 100644 --- a/assets/analysis_pseudo.qmd +++ b/assets/analysis_pseudo.qmd @@ -37,8 +37,13 @@ pseudo_groups_file <- "pseudo_groups.txt" pseudo_groups <- scan(pseudo_groups_file, character()) ``` -#### Dataset name: `{r} cohort_id` -#### Fields used for pseudobulking: +With the dataset filtered, integrated, and annotated, the next step is to perform pseudobulking. This is the process of aggregating the RNA counts across many cells within a sample into a single data point. For example, we may have a sample with 5000 B cells and 4000 T cells. During pseudobulking, we would aggregate (i.e. sum together) the expression profiles of all B cells and all T cells, leaving us with one pseudobulked B cell expression profile and one pseudobulked T cell expression profile for that sample. We would then repeat this procedure for every sample. + +The process of pseudobulking effectively performs a bulk RNA-seq experiment for each cell type - hence the name "pseudobulk". But, instead of relying on cell sorting methods to isolate each cell type from one another, we can use our single cell resolution data and cell annotations. The main advantage of pseudobulking is that we can use well-established statistical techniques for assessing differentially expressed genes that were developed for bulk RNA-seq and are more robust than newer, more complex scRNA-seq techniques. See the [single-cell best practices guide](https://www.sc-best-practices.org/conditions/differential_gene_expression.html) for more information on the benefits of pseudobulking. + +#### Fields used for pseudobulking + +The following [metadata fields](#all-annotation-fields-available-for-downstream-analysis) were used for aggregation of cells: ```{r} #| echo: false @@ -52,6 +57,8 @@ for (g in pseudo_groups) { ### Pseudobulked groups availablle for comparison +Below is a table summarising all of your final pseudobulk groups and the number of samples in your dataset that are represented in each group. These pseudobulk group names are now available for comparison during differential expression analysis. **Note** that for differential expression analysis, you will ideally have *at least three* samples for each of your case and control groups, if not more. + ```{r} #| echo: false #| warning: false diff --git a/assets/annotation.qmd b/assets/annotation.qmd index ddde05d..27156e8 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,9 +151,65 @@ if (is.null(cohort_id)) { } ``` -#### Dataset name: `{r} cohort_id` +There are three main annotation options supported by the scRNAvigator pipeline: + +- Cell cycle annotation +- Automatic cell type annotation +- Custom marker gene-based annotation + +Cell cycle annotation is useful if you are studying rapidly-dividing cells such as tumor samples and wish to capture information about which stage of the cell cycle each cell is in. + +Automatic cell type annotation is performed by comparing your data against curated and pre-annotated public datasets and assigning annotations to your cells based on their similarities to the curated datasets. For human samples, the Human Primary Cell Atlas (HPCA) database is used for this purpose by default, while for non-human species, an external database needs to be provided in order to perform automatic cell type annotation. + +Custom marker gene-based annotation takes a list of cell types and associated marker genes and runs a scoring algorithm to determine the closest matching cell type for each cell. Custom cell type annotations creates two annotation fields in the sample metadata: `custom_cell_type` and `custom_cell_type.max_score`. The `custom_cell_type.max_score` annotation simply picks the custom cell type that scored the highest. The `custom_cell_type` annotation, however, also takes into account how much the highest score differs from the second-highest score; when these are too close together for a given cell, the cell will be labelled "Ambiguous". Of these two annotation fields, we recommend using `custom_cell_type` for downstream analysis due to it being more conservative. + +After annotation with one or more of the above methods, your cell **clusters** will be annotated as well. This takes one of the annotation categories from above and determines whether that cluster predominantly consists of a single cell type or label. If so, that label is applied to the whole cluster. If there is too much heterogeneity within a cluster, it is labelled "Ambiguous". + +#### All annotation fields available for cluster annotation + +Your cells were annotated with the following annotation fields: + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false +all_annotations_df <- tibble( + annotation_field = c(), + description = c() +) +annotation_descriptions <- list( + custom_cell_type = "Custom marker gene-based cell type annotation (highest score significantly higher than second-highest score)", + custom_cell_type.max_score = "Custom marker gene-based cell type annotation (highest score only)", + SingleR.annotation = "SingleR annotation with custom cell annotation database", + SingleR.hpca_main = "SingleR annotation with HPCA database ('main' cell type category)", + SingleR.hpca_fine = "SingleR annotation with HPCA database ('fine' cell type category)", + Phase = "Cell cycle phase" +) +for (ann in annotation_priority) { + all_annotations_df <- all_annotations_df %>% + bind_rows(tibble( + annotation_field = ann, + description = annotation_descriptions[[ann]] + )) +} +``` + +```{r} +#| echo: false +#| eval: !expr (dim(all_annotations_df)[1] > 0) + +datatable( + all_annotations_df, + colnames = c("Annotation field", "Description") +) +``` + +Cell **clusters** were annotated with the labels from the **`r annotation_used_for_clusters`** annotation field. + +#### All annotation fields available for downstream analysis -### All annotation fields available for downstream analyses +The following list contains all of the annotations listed above, as well as the `cluster_annotation` field and additional metadata fields that were provided in the samplesheet. These are all of the fields available for performing downstream pseudobulking and subsequent differential expression and functional enrichment analyses. ```{r} #| echo: false @@ -158,7 +227,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')) ``` @@ -166,10 +235,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')) ``` @@ -177,7 +246,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')) ``` @@ -245,10 +314,10 @@ datatable(custom_cell_type_props_table) #| results: asis #| echo: false #| eval: !expr custom_annotated -cat("#### Custom cell type feature plots\n\n") +cat("##### Custom cell type feature plots\n\n") for (cell_type in names(custom_cell_type_feature_plots)) { p <- custom_cell_type_feature_plots[[cell_type]] - cat(paste0("##### ", cell_type, "\n\n")) + cat(paste0("###### ", cell_type, "\n\n")) cat(paste0('feature plot: ', cell_type, ' - ', cohort_id, '\n\n')) } ``` @@ -257,7 +326,7 @@ for (cell_type in names(custom_cell_type_feature_plots)) { #| results: asis #| echo: false #| eval: !expr clusters_annotated -cat("### Cluster cell type annotation\n\n") +cat("#### Cluster cell type annotation\n\n") cat(paste0('cluster cell type assignment umap: ', cohort_id, '\n\n')) ``` diff --git a/assets/index.qmd b/assets/index.qmd index 66ead11..c7be38b 100644 --- a/assets/index.qmd +++ b/assets/index.qmd @@ -1,4 +1,656 @@ --- title: "scRNAvigator Report" -format: html +format: + html: + embed-resources: true --- + +Thank you for using scRNAvigator! + +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 +- Dataset integration results +- Cell annotation results +- 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] + +# Get parameters TSV file +params_tsv <- "parameters.tsv" +params_df <- read_tsv(params_tsv) + +# 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") + +# 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) + +# 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 <- 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)] +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_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_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)] +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)] +``` + +## Analysis summary: `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 + +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) +``` + +### 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 + +You can use the following text within your methods sections of any manuscripts using the results of this 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`; 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 +#| 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, "; 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, "; 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")) { + 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, "; 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")) { + cat("Cells types were annotated using custom marker gene sets and the `Seurat` built-in function `AddModuleScore`.\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_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") + } +} + +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, "; 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, "; 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") + } + + if (dir.exists("analysis_ora")) { + cat("Over-Representation Analysis (ORA).\n") + } + +} +``` + +::: + +```{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("### Sample 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") || 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 +#| eval: !expr file.exists("analysis_gsea.qmd") +#| warning: false +#| error: false +#| message: false +cat("#### Gene Set 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("#### Over-Representation Analysis (ORA)\n\n") +``` + +```{r} +#| child: "analysis_ora.qmd" +#| eval: !expr file.exists("analysis_ora.qmd") +#| warning: false +#| 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/assets/integration_cluster.qmd b/assets/integration_cluster.qmd index 4d41305..23ba76f 100644 --- a/assets/integration_cluster.qmd +++ b/assets/integration_cluster.qmd @@ -41,10 +41,11 @@ integration_res_file <- "integration_resolution.txt" integration_res <- as.numeric(scan(integration_res_file, character())[1]) ``` -#### Dataset name: `{r} cohort_id` -#### 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 +#### 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 @@ -52,7 +53,9 @@ integration_res <- as.numeric(scan(integration_res_file, character())[1]) 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. ```{r} #| results: asis @@ -60,7 +63,9 @@ cat(paste0('integrated clusters umap: ', cohort_ cat(paste0('<img src=\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. 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/integration_qc.qmd b/assets/integration_qc.qmd index 8f7b700..f1cc1e3 100644 --- a/assets/integration_qc.qmd +++ b/assets/integration_qc.qmd @@ -37,9 +37,15 @@ integrated_umap_file <- dir( cohort_id <- gsub("^integration_qc\\/(.*)\\.umap\\.merged\\.png$", "\\1", merged_umap_file) ``` -#### Dataset name: `{r} cohort_id` +Following initial quality control of your data, the next phase of the scRNAvigator workflow is to conduct dataset integration. Integration is a vital step where all of your samples are merged into a single dataset and undergo batch effect correction. -### Merged (unintegrated) dataset UMAP +Batch effects arise in single cell data due to slight differences in how each sample is handled, processed, or sequenced. These small differences are enough to cause cells to cluster together based on their sample of origin rather than by biologically significant differences, such as cell type or experimental condition. This erroneous clustering can impact downstream analysis and so must be corrected. + +We use the anchor-based CCA ("Canonical Correlation Analysis") integration method implemented in the `Seurat` package to perform this batch correction. Briefly, this method looks for shared patterns in the correlation between cells across each of the samples and uses this to estimate and correct for the batch effects in the data. + +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. + +#### Merged (unintegrated) dataset UMAP ```{r} #| results: asis @@ -47,7 +53,7 @@ cohort_id <- gsub("^integration_qc\\/(.*)\\.umap\\.merged\\.png$", "\\1", merged 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 b734ad8..6b5a82e 100644 --- a/assets/qc_cluster.qmd +++ b/assets/qc_cluster.qmd @@ -56,6 +56,73 @@ 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. + +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false + +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_df <- NULL +} +``` + +```{r} +#| echo: false +#| eval: !expr (!is.null(cluster_res_df)) +#| results: asis +#| warning: false +#| error: false +#| message: false + +cat("The following clustering resolutions have been set for each sample:") +``` + +```{r} +#| echo: false +#| eval: !expr (!is.null(cluster_res_df)) +#| warning: false +#| error: false +#| message: false + +datatable( + cluster_res_df, + colnames = c("Sample", "Sample clustering 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) +``` + +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} @@ -66,15 +133,18 @@ 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("#####", "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 1cd8a61..0d2c467 100644 --- a/assets/qc_doublets.qmd +++ b/assets/qc_doublets.qmd @@ -95,6 +95,12 @@ 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. We have used the R package `DoubletFinder` to identify barcodes that are likely doublets and removed them 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. + +Finally, below the plots you will find a summary table with the number and proportion of doublets identified per sample. + :::panel-tabset ```{r} @@ -110,21 +116,24 @@ 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("#####", "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_filter.qmd b/assets/qc_filter.qmd index 7beba75..658164b 100644 --- a/assets/qc_filter.qmd +++ b/assets/qc_filter.qmd @@ -61,8 +61,71 @@ 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( + "sample", + "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 == "") +``` + +```{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: + +" + +cat(markdown_text) ``` +```{r} +#| echo: false +#| warning: false +#| error: false +#| message: false +#| 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", + 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 = 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. 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 ```{r} @@ -75,13 +138,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')) } ) @@ -91,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( 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/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 e224ae6..4bd606e 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) + samplesheet_csv = channel.fromPath(params.input, checkIfExists: true) + samplesheet = samplesheet_csv .splitCsv( header: true ) .map { row -> { def sample = row.sample @@ -341,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, 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) @@ -392,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) @@ -407,6 +408,11 @@ 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([]) + comparisons_csv = !!params.comparisons ? channel.fromPath(params.comparisons, checkIfExists: true) : channel.value([]) // Run the report module report_input = cohort_id @@ -427,5 +433,47 @@ workflow { .merge(analysis_ora_results) .merge(available_annotation_files) .merge(report_style) - REPORT(report_input) + + // 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 + ] + 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 29a4521..6264f96 100644 --- a/modules/report.nf +++ b/modules/report.nf @@ -24,36 +24,45 @@ process REPORT { path(analysis_ora_results, stageAs: "analysis_ora/*"), path(available_annotation_files, stageAs: "available_annotations/??.txt"), path(report_style) + 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" + path quarto_yaml + val env_vars + val software_versions 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' + env_vars.params.collect { p, v -> "${p}\t${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 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 fi - - # Generate quarto config - create_quarto_yml.py + echo -e "${param_tsv}" > parameters.tsv + 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 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 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