run some abundance metrics without FP and FN

Author: d4straub

bin/compare_performance.r

@@ -14,14 +14,12 @@ mismatch_threshold <- as.numeric(args[6])
 
 # 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)
@@ -33,7 +31,6 @@ get_stats <- function(i_exp,i_obs,df,sample) {
 		Fone <- 0
 		Fbeta <- 0
 	}
-
 	# save
 	types <- c(
 		"observed",
@@ -67,32 +64,21 @@ get_stats <- function(i_exp,i_obs,df,sample) {
 		paste(head(FN, n=100), collapse=','),
 		paste(head(FP, n=100), collapse=',')
 	)
+	# prepare output
 	ids <- rep(sample, length(types))
-
 	df_append <- data.frame(
 		sample=ids,
 		type= types,
 		value= values,
 		stringsAsFactors=FALSE)
 	df <- rbind( df, df_append )
-
 	return(df)
 }
 
-# function to produce abundance statistics
-get_stats_abundance <- function(expected_abund,observed_abund,df,sample) {
-
-	# if none are expected, skip!
-	if( length(expected_abund)==0 ) {
-		print( paste("No expected abundances in sample",sample, "- skipping") ); next
-	} else {
-		print( paste(length(expected_abund),"expected seq in sample",sample) )
-	}
-
+# function to produce distances based on complete abundance statistics (incl. FN & FP)
+get_stats_distances <- function(expected_abund,observed_abund,df,sample) {
+	if( length(expected_abund)==0 ) { next }
 	# stats
-	pearson_cor <- cor.test(expected_abund, observed_abund, method = "pearson")
-	spearman_rho <- cor.test(expected_abund, observed_abund, method = "spearman")
-	percent_dev <- abs((observed_abund - expected_abund) / expected_abund) * 100
 	jensen_shannon <- function(p, q) {
 		# KL divergence function (with pseudocount to handle zeros)
 		kl_div <- function(x, y, pseudocount = 1e-10) {
@@ -102,46 +88,65 @@ get_stats_abundance <- function(expected_abund,observed_abund,df,sample) {
 			y <- y / sum(y)  # normalize
 			sum(x * log(x / y))
 		}
-
 		# Midpoint distribution
 		m <- (p + q) / 2
-
 		# Jensen-Shannon is the square root of the average of two KL divergences
 		sqrt(0.5 * kl_div(p, m) + 0.5 * kl_div(q, m))
 	}
+	# save
+	types <- c(
+		"bray-curtis",
+		"hellinger",
+		"jensen-shannon"
+	)
+	values <- c(
+		1 - (2 * sum(pmin(observed_abund, expected_abund))) / (sum(observed_abund) + sum(expected_abund)), # Bray-Curtis Dissimilarity
+		sqrt(sum((sqrt(observed_abund) - sqrt(expected_abund))^2)), # Hellinger Distance
+		jensen_shannon(observed_abund, expected_abund) # Jensen-Shannon Divergence
+	)
+	# prepare output
+	ids <- rep(sample, length(types))
+	df_append <- data.frame(
+		sample=ids,
+		type=types,
+		value=values,
+		stringsAsFactors=FALSE)
+	df <- rbind( df, df_append )
+	return(df)
+}
 
+# function to produce abundance statistics based on filtered abundance statistics (excl. FN & FP)
+get_stats_abundance <- function(expected_abund,observed_abund,df,sample) {
+	if( length(expected_abund)==0 ) { next }
+	# stats
+	pearson_cor <- cor.test(expected_abund, observed_abund, method = "pearson")
+	spearman_rho <- cor.test(expected_abund, observed_abund, method = "spearman")
+	percent_dev <- abs((observed_abund - expected_abund) / expected_abund) * 100
 	# save
 	types <- c(
 		"pearson_cor",
 		"spearman_rho",
 		"deviation",
 		"mae",
 		"rmse",
-		"ps",
-		"bray-curtis",
-		"hellinger",
-		"jensen-shannon"
+		"ps"
 	)
 	values <- c(
 		as.list(pearson_cor$estimate)[[1]], # Pearson's Correlation (cor)
 		as.list(spearman_rho$estimate)[[1]], # Spearman's Rank Correlation (rho)
 		median(percent_dev, na.rm = TRUE), # Median Percent Abundance Deviation
 		mean(abs(observed_abund - expected_abund), na.rm = TRUE), # Mean Absolute Error (MAE)
 		sqrt(mean((observed_abund - expected_abund)^2, na.rm = TRUE)), # Root Mean Square Error (RMSE)
-		sum(pmin(observed_abund, expected_abund)) / sum(expected_abund), # Proportional Similarity (PS)
-		1 - (2 * sum(pmin(observed_abund, expected_abund))) / (sum(observed_abund) + sum(expected_abund)), # Bray-Curtis Dissimilarity
-		sqrt(sum((sqrt(observed_abund) - sqrt(expected_abund))^2)), # Hellinger Distance
-		jensen_shannon(observed_abund, expected_abund) # Jensen-Shannon Divergence
+		sum(pmin(observed_abund, expected_abund)) / sum(expected_abund) # Proportional Similarity (PS)
 	)
+	# prepare output
 	ids <- rep(sample, length(types))
-
 	df_append <- data.frame(
 		sample=ids,
-		type= types,
-		value= values,
+		type=types,
+		value=values,
 		stringsAsFactors=FALSE)
 	df <- rbind( df, df_append )
-
 	return(df)
 }
 
@@ -258,8 +263,12 @@ for (sample in SAMPLES) {
 	obsIDs <- data$ID[!is.na(data$query)]
 	df <- get_stats( expIDs, obsIDs, df, sample)
 
-	# calculate abundance stats, this includes also non-matching sequences
-	df <- get_stats_abundance( data$expected_abund, data$observed_abund, df, sample)
+	# Calculate distance metrics, this *includes* non-matching sequences (FP & FN)
+	df <- get_stats_distances( data$expected_abund, data$observed_abund, df, sample)
+
+	# Calculate relationship and abundance error, this *excludes* non-matching sequences (FP & FN)
+	data_matches <- data[data$observed_abund >0 & data$expected_abund >0,] # only on matched observed to expected
+	df <- get_stats_abundance( data_matches$expected_abund, data_matches$observed_abund, df, sample)
 
 	# (4) PLOTS - pairwise comparisons
 
@@ -282,7 +291,6 @@ for (sample in SAMPLES) {
 	invisible(dev.off())
 
 	# Scatter plot: Observed vs. Expected Abundance (log-log)
-	data_matches <- data[data$observed_abund >0 & data$expected_abund >0,] # only on matched observed to expected
 	outfile <- paste0(sample,"_scatter_loglog")
 	print(paste("write",outfile))
 	svg(paste0(outfile,".svg"), width = 8, height = 6)

docs/output.md

@@ -681,7 +681,7 @@ When expected sequences are supplied, the following files will be produced:
 
 When additionally expected abundances are available, additional performance metrics will be generated.
 
-(1) Metrics base on presence/absence:
+#### Metrics base on presence/absence:
 
 - `observed`: Number of observed sequences
 - `expected`: Number of expected sequences
@@ -698,24 +698,31 @@ When additionally expected abundances are available, additional performance metr
 - `FNs_exp`: FN IDs (max 100) corresponding to expected sequences
 - `FPs_obs`: FP IDs (max 100) corresponding to observed sequences
 
-(2) Metrics based on abundances:
+> [!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 and vice versa. 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.
+
+#### Metrics based on abundances:
+
+(1) based on filtered abundance tables, **excluding** non-matching sequences (FP & FN)
+
+These metrics would be inflated, dominated or skewed by false positives and false negatives.
 
 - `pearson_cor`: Pearson's Correlation (cor)
 - `spearman_rho`: Spearman's Rank Correlation (rho)
 - `deviation`: Median Percent Abundance Deviation
 - `mae`: Mean Absolute Error (MAE)
 - `rmse`: Root Mean Square Error (RMSE)
 - `ps`: Proportional Similarity (PS)
+
+(2) based on complete abundance tables, **including** non-matching sequences (FP & FN)
+
 - `bray-curtis`: Bray-Curtis Dissimilarity
 - `hellinger`: Hellinger Distance
 - `jensen-shannon`: Jensen-Shannon Divergence
 
-> [!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">

docs/usage.md

@@ -343,10 +343,12 @@ To investigate the quality of data generation and/or data analysis, analysis out
 The observed sequences will be aligned globally (using `VSEASRCH --usearch_global`) to the expected sequences (`--expected_sequences`).
 Depending on the region to analyse (`--expected_sequences_region`) the mismatches and gaps within the alignment will be summarized with or without terminal gaps.
 The nucleotide differences will be evaluated for each observed sequence to its best match.
+
 Expected abundances per sequence (`--expected_abundances`) enable sample specific presence/absence metrics and abundance-based comparisons.
-For those, observed and expected sequences will be aggregated by their analysed region: when one observed sequence matches exactely to several expected sequences, the expected sequence IDs (and abundances) will be concatenated and vice versa. That means, for example, if observed sequences are shorter than expected sequences and the analysed region is "query", each expected sequence (ID) that matches the same observed sequence will be aggregated.
+For those, observed and expected sequences will be aggregated by their analysed region: when one observed sequence matches exactely to several expected sequences, the expected sequence IDs will be concatenated and vice versa.
+That means, for example, if observed sequences are shorter than expected sequences and the analysed region is "query", each expected sequence (ID) that matches the same observed sequence will be aggregated.
 
-The aggregation of obsevered and expected IDs and abundances based on perfect matches:
+The aggregation of obsevered and expected IDs and abundances based on perfect matches
 
 | obsID | expID | exp_abund | obs_abund |
 | ----- | ----- | --------- | --------- |
@@ -373,7 +375,7 @@ will be transformed to:
 | h     |       | 0         | 1         |
 |       | 4     | 1         | 0         |
 
-This aggregation will not work properly when many sequences are not observed, e.g. in the above axample observed ID "a" links expected IDs "1" and "2", which would not have been aggregated if "a" would not have been observed. This would inflate expected sequences. Therefore, optimal sequence and abundance input are tailored towards the actual sequenced region and de-duplicated.
+This aggregation will not work properly when many sequences are not observed, e.g. in the above example observed ID "a" links expected IDs "1" and "2", which would not have been aggregated if "a" would not have been observed. This would inflate expected sequences. Therefore, optimal sequence and abundance input are tailored towards the actual sequenced region and de-duplicated.
 
 ### Differential abundance analysis