Comparison of observed sequences to expected, the second

CHANGELOG.md

@@ -7,7 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### `Added`
 
-- [#1021](https://github.com/nf-core/ampliseq/pull/1021) - Comparison of observed ASVs against expected sequences or their abundance is now available with `--expected_*` parameters (by @d4straub, reviewed by @erikrikarddaniel)
+- [#1021](https://github.com/nf-core/ampliseq/pull/1021),[#1023](https://github.com/nf-core/ampliseq/pull/1023) - Comparison of observed ASVs against expected sequences or their abundance is now available with `--expected_*` parameters (by @d4straub, reviewed by @erikrikarddaniel)
 
 ### `Changed`
 

bin/compare_performance.r

@@ -0,0 +1,198 @@
+#!/usr/bin/env Rscript
+
+# compare_sequences.r
+
+# Get params and files from the command line
+args            <- commandArgs(trailingOnly=TRUE)
+tag             <- args[1] # tag for each sample
+alignmentFILE   <- args[2] # alignment file; expected columns: query, target, mismatches_final
+obsabundFILE    <- args[3] # observed abundance*
+expabundFILE    <- args[4] # expected abundance*
+fbeta           <- as.numeric(args[5]) # Fbeta weight, default 2; 1=precision and recall equally weighted (=F1 score), 2=weighs recall higher than precision, 0.5=weighs recall lower than precision.
+mismatch_threshold <- as.numeric(args[6])
+# tab-separated file with header: first column with sequences name, following one or many columns (=samples) with numeric values (=abundance), only presence (>0)/absence are used here
+
+# function to produce statistics
+get_stats <- function(i_exp,i_obs,df,sample) {
+
+	# if none are expected, skip!
+	if( length(i_exp)==0 ) {
+		print( paste("No expected seq in sample",sample, "- skipping") ); next
+	} else {
+		print( paste(length(i_exp),"expected seq in sample",sample) )
+	}
+
+	# stats
+	TP <- intersect(i_exp, i_obs)
+	FN <- setdiff(i_exp, i_obs)
+	FP <- setdiff(i_obs, i_exp)
+	if ( length(TP) > 0 ) {
+		Fone <- ( 2 * length(TP)/length(i_obs) * length(TP)/length(i_exp) ) / ( length(TP)/length(i_obs) + length(TP)/length(i_exp) )
+		Fbeta <- ( (1+fbeta^2) * length(TP)/length(i_obs) * length(TP)/length(i_exp) ) / ( (fbeta^2) * length(TP)/length(i_obs) + length(TP)/length(i_exp) )
+	} else {
+		Fone <- 0
+		Fbeta <- 0
+	}
+
+	# save
+	types <- c(
+		"observed",
+		"expected",
+		"TP",
+		"FN",
+		"FP",
+		"recall",
+		"precision",
+		"F1",
+		"Fbeta",
+		"fdr",
+		"jaccard",
+		"TPs_exp",
+		"FNs_exp",
+		"FPs_obs"
+    )
+	values <- c(
+		length(i_obs),
+		length(i_exp),
+		length(TP),
+		length(FN),
+		length(FP),
+		length(TP)/length(i_exp),
+		length(TP)/length(i_obs),
+		Fone,
+		Fbeta,
+		length(FP)/(length(i_obs)),
+		length(TP)/(length(i_obs)+length(FN)),
+		paste(head(TP, n=100), collapse=','),
+		paste(head(FN, n=100), collapse=','),
+		paste(head(FP, n=100), collapse=',')
+	)
+	ids <- rep(sample, length(types))
+
+	df_append <- data.frame(
+		sample=ids,
+		type= types,
+		value= values,
+		stringsAsFactors=FALSE)
+	df <- rbind( df, df_append )
+
+	return(df)
+}
+
+# PREPARE INPUT
+
+# Read sequence alignment table
+alignment = read.table( alignmentFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
+alignment <- alignment[alignment$target !="*" & alignment$mismatch_final <= mismatch_threshold,]
+
+# Read pipeline's ASV abundance table (either from QIIME2 [skipping first line] or from previous steps)
+if (grepl("^# Constructed from biom file", readLines(obsabundFILE, n = 1))) {
+	observed = read.table( obsabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE, comment.char = "", skip=1 )
+} else {
+	observed = read.table( obsabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE, comment.char = "" )
+}
+colnames(observed)[1] <- "ID"
+
+# select available samples
+observed_samples <- colnames(observed)[2:ncol(observed)]
+print(paste( "Observed samples:", paste(observed_samples,collapse=",")))
+
+# Read expected abundance table
+exp = read.table( expabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
+colnames(exp)[1] <- "ID"
+# extract samples to analyse
+exp_samples <- colnames(exp)[2:ncol(exp)]
+print(paste( "Expected samples:", paste( exp_samples ,collapse=",")))
+SAMPLES <- exp_samples[exp_samples %in% observed_samples]
+print(paste( "Investigate samples:", paste( SAMPLES ,collapse=",")))
+
+# check if there are any samples to analyse
+if (length(SAMPLES) == 0) {
+	stop("ERROR - Found no samples to investigate")
+}
+
+# ANALYSE
+
+# Initialize dataframe
+df <- data.frame(
+	sample=character(),
+	type=character(),
+	value=character(),
+	stringsAsFactors=FALSE)
+
+# for each sample
+for (sample in SAMPLES) {
+
+	# keep only non-empty columns and sample abundances
+	keep_cols <- c("ID",sample)
+	s_exp <- subset(exp, select = keep_cols)
+	s_obs <- subset(observed, select = keep_cols)
+
+	# keep only observed (abundance > 0)
+	s_exp = s_exp[s_exp[,2] > 0,]
+	s_obs = s_obs[s_obs[,2] > 0,]
+
+	# keep only alignments where both, query and target are present
+	s_alignment <- alignment[alignment$query %in% s_obs$ID & alignment$target %in% s_exp$ID,]
+
+	# translate s_obs$ID to alignment$target if available
+	obsIDs <- c()
+	for (obsID in s_obs$ID){
+		i_alignment <- s_alignment[s_alignment$query == obsID,]
+		# collapse several targets if multiple matches
+		if( nrow(i_alignment) > 0 ){
+			target <- paste(i_alignment$target, collapse="-")
+		} else {
+			target <- obsID
+		}
+		obsIDs <- c(obsIDs, target)
+	}
+
+	# aggregate s_exp$ID if several alignment$query match
+	expIDs <- s_exp$ID[!s_exp$ID %in% s_alignment$target]
+	i_alignment <- aggregate(target ~ query, s_alignment, paste, collapse = "-")
+	expIDs <- c(expIDs, i_alignment$target)
+
+	# calculate stats
+	df <- get_stats( unique(expIDs), unique(obsIDs), df, sample)
+}
+
+# Write detailed output
+outfile <- "performance_per-sample.tsv"
+df$tag <- rep(tag, nrow(df))
+print(paste("write",outfile))
+write.table(df, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
+
+
+# Make & write summary output
+outfile <- "performance_summary.tsv"
+df_sum <- data.frame(type=character(), median=numeric(), mean=numeric(), min=numeric(), max=numeric(), count=numeric(), stringsAsFactors=FALSE)
+print("Warnings are expected at that point:")
+for (type in unique(df$type)) {
+	df_subset <- df[df$type == type, ]
+	MEDIAN <- median( as.numeric(df_subset$value) )
+	MEAN <- mean( as.numeric(df_subset$value) )
+	MIN <- min( as.numeric(df_subset$value) )
+	MAX <- max( as.numeric(df_subset$value) )
+	COUNT <- length( as.numeric(df_subset$value) )
+	df_sum[nrow(df_sum) + 1,] <- list(type, MEDIAN, MEAN, MIN, MAX, COUNT)
+}
+df_sum$tag <- rep(tag, nrow(df_sum))
+print(paste("write",outfile))
+write.table(df_sum, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
+
+# Plot Values
+df_subset <- subset(df, type %in% c("recall","precision","F1","Fbeta","fdr","jaccard") )
+df_subset$value <- as.numeric(df_subset$value)
+
+outfile <- "performance_boxplot.svg"
+print(paste("write",outfile))
+svg(outfile, height = 8, width = 10)
+boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1))
+invisible(dev.off())
+
+outfile <- "performance_boxplot.png"
+print(paste("write",outfile))
+png(outfile, height = 400, width = 500)
+boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1))
+invisible(dev.off())

conf/modules.config

@@ -1241,6 +1241,23 @@ process {
         ]
     }
 
+    withName: COMPARE_PERFORMANCE {
+        ext.mismatch_threshold = 0 // thereshold for comparing matches (0=exact match)
+        ext.fbeta = 2 // Fbeta weight, default 2; 1=precision and recall equally weighted (=F1 score), 2=weighs recall higher than precision, 0.5=weighs recall lower than precision.
+        publishDir = [
+            [
+                path: { "${params.outdir}/comparison/${meta.id}" },
+                mode: params.publish_dir_mode,
+                pattern: "*{performance_summary.tsv,.log,.svg,.png}"
+            ],
+            [
+                path: { "${params.outdir}/comparison/${meta.id}/per-sample" },
+                mode: params.publish_dir_mode,
+                pattern: "*{_per-sample.tsv}"
+            ]
+        ]
+    }
+
     withName: 'MULTIQC' {
         ext.args   = { params.multiqc_title ? "--title \"$params.multiqc_title\"" : '' }
         publishDir = [

docs/output.md

@@ -662,21 +662,61 @@ On request (`--ancombc2`), ANCOM-BC2 is applied to each suitable or specified me
 
 Comparison evaluates observed results against expected outcomes, typically for samples with known composition such as mock communities, to assess the data and analysis. Steps to evaluate the produced ASVs are implemented in the pipeline.
 
+When expected sequences are supplied, the following files will be produced:
+
 <details markdown="1">
 <summary>Output files</summary>
 
 - `comparison/<parameter md5sum>_<pipeline version>`
   - `vsearch_usearchglobal.tsv`: VSEARCH --usearch_global output for ASV to sequence comparisons.
   - `md5sum_version.txt`: Contains parameter md5sum and pipeline version that is also the folder name, additionally the time stamp used in files in `pipeline_info/`.
   - `nucleotide-differences.log`: Log file for comparing matches of ASVs to expected sequences.
-  - `nucleotide-differences.tsv`: Tabl-separated table based on VSEARCH results comparing matches of ASVs to expected sequences.
+  - `nucleotide-differences.tsv`: Tab-separated table based on VSEARCH results comparing matches of ASVs to expected sequences.
 - `comparison/<parameter md5sum>_<pipeline version>/per-sample`
   - `<sample>_nucleotide-differences_per-sample.tsv`: Number of sequences and mismatches to expected sequences.
   - `<sample>_nucleotide-distance_barplot.png`: Barplot of number of sequences versus number of mismatches in png format.
   - `<sample>_nucleotide-distance_barplot.svg`: Barplot of number of sequences versus number of mismatches in svg format.
 
 </details>
 
+When additionally expected abundances are available, additional performance metrics will be generated:
+
+- `observed`: Number of observed sequences
+- `expected`: Number of expected sequences
+- `TP`: Number of true positive sequences (observed and expected)
+- `FN`: Number of false negative sequences (not observed but expected)
+- `FP`: Number of false positive sequences (observed but not expected)
+- `recall`: Recall of expected sequences = TP/(TP+FN)
+- `precision`: Precision = TP/(TP+FP)
+- `F1`: F1 score, the harmonic mean of the precision and recall = 2TP/(2TP+FP+FN)
+- `Fbeta`: Fbeta score, the weighted F1 score. By default beta=2, i.e. weighs recall twice higher than precision
+- `fdr`: False discovery rate = FP/(TP+FP)
+- `jaccard`: Jaccard index = TP/(TP+FP+FN)
+- `TPs_exp`: TP IDs (max 100) corresponding to expected sequences
+- `FNs_exp`: FN IDs (max 100) corresponding to expected sequences
+- `FPs_obs`: FP IDs (max 100) corresponding to observed sequences
+
+> [!NOTE]
+> The F1 score is unreliable with strongly unbalanced data.
+
+> [!WARNING]
+> If the supplied expected sequences are not unique in the region of the alignment with the observed sequences, alignment matches are used to aggregate identical expected sequences. In case there isnt an exact match to a set of identical expected sequences, those will not be aggregated and inflate the number of expected sequences.
+
+The following additional files will be produced:
+
+<details markdown="1">
+<summary>Output files</summary>
+
+- `comparison/<parameter md5sum>_<pipeline version>`
+  - `performance_summary.tsv`: Tab-separated table with aggregated performance metrics.
+  - `performance.log`: Log file, complementary to `nucleotide-differences.log`.
+  - `performance_boxplot.png`: Boxplot of number of aggregated performance metrics in png format.
+  - `performance_boxplot.svg`: Boxplot of number of aggregated performance metrics in svg format.
+- `comparison/<parameter md5sum>_<pipeline version>/per-sample`
+  - `performance_per-sample.tsv`: Tab-separated table with performance metrics per sample, in long format.
+
+</details>
+
 ### PICRUSt2
 
 PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) is a software for predicting functional abundances based only on marker gene sequences. On demand (`--picrust`), Enzyme Classification numbers (EC), KEGG orthologs (KO) and MetaCyc ontology predictions will be made for each sample.

modules/local/compare_performance.nf

@@ -0,0 +1,45 @@
+process COMPARE_PERFORMANCE {
+    tag "${meta.id}"
+    label 'process_single'
+
+    conda "conda-forge::r-base=4.2.1"
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/r-base:4.2.1' :
+        'biocontainers/r-base:4.2.1' }"
+
+    input:
+    tuple val(meta), path(alignment) // alignment file with "query,target,mismatch_final"
+    path(observed_abundances)        // ASV table from within pipeline
+    path(expected_abundances)        // ASV table from params.expected_abundances
+
+    output:
+    path("*.svg")                            , emit: svg
+    path("*.png")                            , emit: png
+    path("performance_summary.tsv")          , emit: performance
+    path("performance_per-sample.tsv")       , emit: performance_per_sample
+    path("performance.log")                  , emit: log
+    path("Warnings.txt")                     , emit: warnings
+    path "versions.yml"                      , emit: versions, topic: versions
+
+    script:
+    def fbeta = task.ext.fbeta ?: 2
+    def mismatch_threshold = task.ext.mismatch_threshold ?: 0
+    """
+    compare_performance.r \\
+        "${meta.id}" \\
+        "$alignment" \\
+        "$observed_abundances" \\
+        "$expected_abundances" \\
+        "$fbeta" \\
+        "$mismatch_threshold" \\
+        >"performance.log"
+
+    # isolate warnings (when grep find no match, exit code is 1, "|| true" fixes that)
+    grep "WARN -" "performance.log" | sed 's/^\\[1\\] //g' | sed 's/"//g' > Warnings.txt || true
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        R: \$(R --version | sed -n 1p | sed 's/R version //g' | sed 's/\\s.*\$//')
+    END_VERSIONS
+    """
+}

nextflow_schema.json

@@ -978,7 +978,7 @@
                     "mimetype": "text/tsv",
                     "pattern": "^\\S+\\.(fasta|fas|fna|fa|ffn)$",
                     "description": "Path to file of expected sequences in fasta format",
-                    "help_text": "Supplying only sequences will compare ASVs of each sample to the given sequence set. Consider other parameters for more comprehensive comparisons."
+                    "help_text": "Ideally the expected sequences should be all unique within the alignment region with the observed sequences, because the alignment matches are used to aggregate identical expected sequences and in case there isnt an exact match to a set of identical expected sequences, those will not be aggregated and inflate the number of expected unique sequences. Supplying only sequences will compare ASVs of each sample to the given sequence set. Consider other parameters for more comprehensive comparisons."
                 },
                 "expected_sequences_region": {
                     "type": "string",

subworkflows/local/comparison_wf.nf

@@ -1,5 +1,6 @@
 include { VSEARCH_USEARCHGLOBAL as VSEARCH_USEARCHGLOBAL_BM } from '../../modules/nf-core/vsearch/usearchglobal/main'
 include { COMPARE_SEQUENCES       } from '../../modules/local/compare_sequences'
+include { COMPARE_PERFORMANCE     } from '../../modules/local/compare_performance'
 
 workflow COMPARISON_WF {
     take:
@@ -28,10 +29,8 @@ workflow COMPARISON_WF {
             if( it.countLines() > 0 ) { log.warn "about comparing sequences\n\n" + it.splitText().join("") }
         }
 
-    // (2) input: val_md5sum_version, "${val_md5sum_version}_nucleotide_differences.tsv", ch_detected_abundances, ch_expected_abundances
-    //     -> compare exact matches per sample: "${val_md5sum_version}_nucleotide_differences.tsv" = 0 & prevalence ch_detected_abundances vs ch_expected_abundances
-    //     -> per sample count FN/FP/... & calculate F1/...
-    //COMPARISON_SEQUENCES ( val_md5sum_version, COMPARE_SEQUENCES.out.matches, ch_detected_abundances, ch_expected_abundances )
+    // Calculate performance metrics per sample such as precision, recall, F1 score
+    COMPARE_PERFORMANCE ( COMPARE_SEQUENCES.out.matches.map { it = [ [id: val_md5sum_version], file(it) ] }, ch_observed_abundances, ch_expected_abundances )
 
     // (3) input: val_md5sum_version, "${val_md5sum_version}_matches.txt", ch_detected_abundances, ch_expected_abundances
     //     -> compare abundance of exact matches