Skip to content

Commit 63e4f06

Browse files
d4straubd4straub
committed
add abundance metrics
1 parent fd3b565 commit 63e4f06

4 files changed

Lines changed: 199 additions & 40 deletions

File tree

‎bin/compare_performance.r‎

Lines changed: 194 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,73 @@ get_stats <- function(i_exp,i_obs,df,sample) {
7979
return(df)
8080
}
8181

82-
# PREPARE INPUT
82+
# function to produce abundance statistics
83+
get_stats_abundance <- function(expected_abund,observed_abund,df,sample) {
84+
85+
# if none are expected, skip!
86+
if( length(expected_abund)==0 ) {
87+
print( paste("No expected abundances in sample",sample, "- skipping") ); next
88+
} else {
89+
print( paste(length(expected_abund),"expected seq in sample",sample) )
90+
}
91+
92+
# stats
93+
pearson_cor <- cor.test(expected_abund, observed_abund, method = "pearson")
94+
spearman_rho <- cor.test(expected_abund, observed_abund, method = "spearman")
95+
percent_dev <- abs((observed_abund - expected_abund) / expected_abund) * 100
96+
jensen_shannon <- function(p, q) {
97+
# KL divergence function (with pseudocount to handle zeros)
98+
kl_div <- function(x, y, pseudocount = 1e-10) {
99+
x <- x + pseudocount
100+
y <- y + pseudocount
101+
x <- x / sum(x) # normalize
102+
y <- y / sum(y) # normalize
103+
sum(x * log(x / y))
104+
}
105+
106+
# Midpoint distribution
107+
m <- (p + q) / 2
108+
109+
# Jensen-Shannon is the square root of the average of two KL divergences
110+
sqrt(0.5 * kl_div(p, m) + 0.5 * kl_div(q, m))
111+
}
112+
113+
# save
114+
types <- c(
115+
"pearson_cor",
116+
"spearman_rho",
117+
"deviation",
118+
"mae",
119+
"rmse",
120+
"ps",
121+
"bray-curtis",
122+
"hellinger",
123+
"jensen-shannon"
124+
)
125+
values <- c(
126+
as.list(pearson_cor$estimate)[[1]], # Spearman's Rank Correlation (rho)
127+
as.list(spearman_rho$estimate)[[1]], # Spearman's Rank Correlation (rho)
128+
median(percent_dev, na.rm = TRUE), # Median Percent Abundance Deviation
129+
mean(abs(observed_abund - expected_abund), na.rm = TRUE), # Mean Absolute Error (MAE)
130+
sqrt(mean((observed_abund - expected_abund)^2, na.rm = TRUE)), # Root Mean Square Error (RMSE)
131+
sum(pmin(observed_abund, expected_abund)) / sum(expected_abund), # Proportional Similarity (PS)
132+
1 - (2 * sum(pmin(observed_abund, expected_abund))) / (sum(observed_abund) + sum(expected_abund)), # Bray-Curtis Dissimilarity
133+
sqrt(sum((sqrt(observed_abund) - sqrt(expected_abund))^2)), # Hellinger Distance
134+
jensen_shannon(observed_abund, expected_abund) # Jensen-Shannon Divergence
135+
)
136+
ids <- rep(sample, length(types))
137+
138+
df_append <- data.frame(
139+
sample=ids,
140+
type= types,
141+
value= values,
142+
stringsAsFactors=FALSE)
143+
df <- rbind( df, df_append )
144+
145+
return(df)
146+
}
147+
148+
# (A) PREPARE INPUT
83149

84150
# Read sequence alignment table
85151
alignment = read.table( alignmentFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
@@ -111,7 +177,7 @@ if (length(SAMPLES) == 0) {
111177
stop("ERROR - Found no samples to investigate")
112178
}
113179

114-
# ANALYSE
180+
# (B) ANALYSE
115181

116182
# Initialize dataframe
117183
df <- data.frame(
@@ -123,47 +189,143 @@ df <- data.frame(
123189
# for each sample
124190
for (sample in SAMPLES) {
125191

192+
# (1) FILTER - for presence/absence in that sample
126193
# keep only non-empty columns and sample abundances
127194
keep_cols <- c("ID",sample)
128195
s_exp <- subset(exp, select = keep_cols)
129196
s_obs <- subset(observed, select = keep_cols)
130-
131197
# keep only observed (abundance > 0)
132198
s_exp = s_exp[s_exp[,2] > 0,]
133199
s_obs = s_obs[s_obs[,2] > 0,]
134-
135200
# keep only alignments where both, query and target are present
136201
s_alignment <- alignment[alignment$query %in% s_obs$ID & alignment$target %in% s_exp$ID,]
137202

138-
# translate s_obs$ID to alignment$target if available
139-
obsIDs <- c()
140-
for (obsID in s_obs$ID){
141-
i_alignment <- s_alignment[s_alignment$query == obsID,]
142-
# collapse several targets if multiple matches
143-
if( nrow(i_alignment) > 0 ){
144-
target <- paste(i_alignment$target, collapse="-")
145-
} else {
146-
target <- obsID
147-
}
148-
obsIDs <- c(obsIDs, target)
149-
}
203+
# (2) AGGREGATE - combine IDs and abundance if multiple query or targets match each other
204+
colnames(s_exp) <- c("ID","expected_abund")
205+
colnames(s_obs) <- c("ID","observed_abund")
206+
s_alignment <- subset(s_alignment, select = c("query","target"))
207+
merged <- merge(s_alignment, s_exp, by.x="target", by.y="ID", all=TRUE)
208+
merged <- merge(merged, s_obs, by.x="query", by.y="ID", all=TRUE)
209+
merged <- unique(merged)
210+
# retain exp without query & obs without target
211+
nomatch_exp <- merged[is.na(merged$query),]
212+
nomatch_obs <- merged[is.na(merged$target),]
213+
# aggregate query target by query (exp without query is lost)
214+
merged <- merged[order(merged$target), ]
215+
merged_exp <- aggregate(target ~ query, merged, paste, collapse = "-")
216+
merged_exp_abund <- aggregate(expected_abund ~ query, merged, sum)
217+
data_exp <- merge(merged_exp, merged_exp_abund, by.x="query", by.y="query", all=TRUE)
218+
data_exp <- unique( merge(data_exp, s_obs, by.x="query", by.y="ID", all=TRUE) )
219+
# aggregate query by target (obs without target is lost)
220+
data_exp <- data_exp[order(data_exp$query), ]
221+
merged_obs <- aggregate(query ~ target, data_exp, paste, collapse = "-")
222+
merged_obs_abund <- aggregate(observed_abund ~ target, data_exp, sum)
223+
data <- merge(merged_obs, merged_obs_abund, by="target", all=TRUE)
224+
# add exp abundance
225+
data_exp_abund <- subset(data_exp, select=c("target","expected_abund"))
226+
data_exp_abund <- unique(data_exp_abund[!is.na(data_exp_abund$target),])
227+
data <- merge(data, data_exp_abund, by="target", all=TRUE)
228+
# re-add exp without query & obs without target
229+
data <- rbind(data, nomatch_obs)
230+
data <- rbind(data, nomatch_exp)
231+
# add ID list
232+
ID_list <- ifelse(is.na(data$target), data$query, data$target)
233+
data$ID <- factor(ID_list, levels = ID_list)
234+
data$expected_abund[is.na(data$expected_abund)] <- 0
235+
data$observed_abund[is.na(data$observed_abund)] <- 0
236+
# make relative abundances
237+
data$expected_abund <- data$expected_abund/sum(data$expected_abund,na.rm=TRUE)
238+
data$observed_abund <- data$observed_abund/sum(data$observed_abund,na.rm=TRUE)
239+
# Sort by expected abundance
240+
data <- data[order(data$expected_abund, decreasing = TRUE), ]
241+
# write table
242+
outfile <- paste0(sample,"_abundances.tsv")
243+
print(paste("write",outfile))
244+
write.table(data, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
245+
246+
# (3) STATISTICS
247+
248+
# calculate absence/presence stats
249+
expIDs <- data$target[!is.na(data$target)]
250+
obsIDs <- data$ID[!is.na(data$query)]
251+
df <- get_stats( expIDs, obsIDs, df, sample)
252+
253+
# calculate abundance stats, this includes also non-matching sequences
254+
df <- get_stats_abundance( data$expected_abund, data$observed_abund, df, sample)
255+
256+
# (4) PLOTS - pairwise comparisons
150257

151-
# aggregate s_exp$ID if several alignment$query match
152-
expIDs <- s_exp$ID[!s_exp$ID %in% s_alignment$target]
153-
i_alignment <- aggregate(target ~ query, s_alignment, paste, collapse = "-")
154-
expIDs <- c(expIDs, i_alignment$target)
258+
# Side-by-Side Bar Plots
259+
outfile <- paste0(sample,"_abundance_barplot.svg")
260+
print(paste("write",outfile))
261+
svg(outfile, height = 8, width = 10)
262+
barplot(
263+
rbind(data$expected_abund, data$observed_abund),
264+
beside = TRUE,
265+
names.arg = data$ID,
266+
col = c("skyblue", "salmon"),
267+
legend.text = c("Expected", "Observed"),
268+
ylab = "Abundance",
269+
xlab = "Taxon",
270+
main = "Side-by-Side Bar Plot: Expected vs. Observed Abundance",
271+
cex.names = 0.7,
272+
las = 2
273+
)
274+
invisible(dev.off())
275+
276+
# Scatter plot: Observed vs. Expected Abundance (log-log)
277+
outfile <- paste0(sample,"_scatter_loglog")
278+
print(paste("write",outfile))
279+
svg(paste0(outfile,".svg"), width = 8, height = 6)
280+
plot(
281+
log10(data$expected_abund),
282+
log10(data$observed_abund),
283+
xlab = "log10(Expected Abundance)",
284+
ylab = "log10(Observed Abundance)",
285+
main = "Scatter Plot: Observed vs. Expected Abundance (log-log)",
286+
pch = 19,
287+
col = "blue"
288+
)
289+
abline(a = 0, b = 1, col = "red", lty = 2)
290+
invisible(dev.off())
155291

156-
# calculate stats
157-
df <- get_stats( unique(expIDs), unique(obsIDs), df, sample)
292+
# Rank Abundance Curves
293+
outfile <- paste0(sample,"_rank_abundance_curve")
294+
print(paste("write",outfile))
295+
svg(paste0(outfile,".svg"), width = 8, height = 6)
296+
plot(
297+
1:nrow(data),
298+
sort(data$expected_abund, decreasing = TRUE),
299+
type = "l",
300+
col = "blue",
301+
lwd = 2,
302+
xlab = "Rank",
303+
ylab = "Abundance",
304+
main = "Rank Abundance Curves",
305+
ylim = c(0, max(data$expected_abund, data$observed_abund))
306+
)
307+
lines(
308+
1:nrow(data),
309+
sort(data$observed_abund, decreasing = TRUE),
310+
col = "red",
311+
lwd = 2
312+
)
313+
legend(
314+
"topright",
315+
legend = c("Expected", "Observed"),
316+
col = c("blue", "red"),
317+
lwd = 2
318+
)
319+
invisible(dev.off())
158320
}
159321

322+
# (5) TABLES - metrics overall
323+
160324
# Write detailed output
161325
outfile <- "performance_per-sample.tsv"
162326
df$tag <- rep(tag, nrow(df))
163327
print(paste("write",outfile))
164328
write.table(df, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
165-
166-
167329
# Make & write summary output
168330
outfile <- "performance_summary.tsv"
169331
df_sum <- data.frame(type=character(), median=numeric(), mean=numeric(), min=numeric(), max=numeric(), count=numeric(), stringsAsFactors=FALSE)
@@ -181,18 +343,18 @@ df_sum$tag <- rep(tag, nrow(df_sum))
181343
print(paste("write",outfile))
182344
write.table(df_sum, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
183345

184-
# Plot Values
185-
df_subset <- subset(df, type %in% c("recall","precision","F1","Fbeta","fdr","jaccard") )
346+
# (6) PLOT - selected metrics, overall
347+
df_subset <- subset(df, type %in% c("recall","precision","F1","Fbeta","fdr","jaccard","pearson_cor","spearman_rho","mae","rmse","ps","bray-curtis","hellinger","jensen-shannon") )
186348
df_subset$value <- as.numeric(df_subset$value)
187-
188349
outfile <- "performance_boxplot.svg"
189350
print(paste("write",outfile))
190-
svg(outfile, height = 8, width = 10)
191-
boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1))
351+
svg(outfile, height = 8, width = 12)
352+
par(mar=c(8, 4, 4, 2))
353+
boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1), las=2)
192354
invisible(dev.off())
193-
194355
outfile <- "performance_boxplot.png"
195356
print(paste("write",outfile))
196-
png(outfile, height = 400, width = 500)
197-
boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1))
357+
png(outfile, height = 400, width = 600)
358+
par(mar=c(8, 4, 4, 2)) #First value (8): Bottom margin (increase this if labels are still cut off); Second value (4): Left margin; Third value (4): Top margin; Fourth value (2): Right margin
359+
boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1), las=2)
198360
invisible(dev.off())

‎conf/modules.config‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1248,12 +1248,12 @@ process {
12481248
[
12491249
path: { "${params.outdir}/comparison/${meta.id}" },
12501250
mode: params.publish_dir_mode,
1251-
pattern: "*{performance_summary.tsv,.log,.svg,.png}"
1251+
pattern: "*{performance_summary.tsv,.log,performance_boxplot.svg,performance_boxplot.png,performance_per-sample.tsv}"
12521252
],
12531253
[
12541254
path: { "${params.outdir}/comparison/${meta.id}/per-sample" },
12551255
mode: params.publish_dir_mode,
1256-
pattern: "*{_per-sample.tsv}"
1256+
pattern: "*{rank_abundance_curve.svg,_scatter_loglog.svg,_abundance_barplot.svg,_abundances.tsv}"
12571257
]
12581258
]
12591259
}

‎modules/local/compare_performance.nf‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ process COMPARE_PERFORMANCE {
1717
path("*.png") , emit: png
1818
path("performance_summary.tsv") , emit: performance
1919
path("performance_per-sample.tsv") , emit: performance_per_sample
20+
path("*_abundances.tsv") , emit: abundance_table_per_sample
2021
path("performance.log") , emit: log
2122
path("Warnings.txt") , emit: warnings
2223
path "versions.yml" , emit: versions, topic: versions

‎subworkflows/local/comparison_wf.nf‎

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -29,14 +29,10 @@ workflow COMPARISON_WF {
2929
if( it.countLines() > 0 ) { log.warn "about comparing sequences\n\n" + it.splitText().join("") }
3030
}
3131

32-
// Calculate performance metrics per sample such as precision, recall, F1 score
32+
// Calculate absence/presence performance metrics per sample such as precision, recall, F1 score
33+
// Calculate abundance performance metrics per sample such as spearman's rho, RMSE, Bray-Curtis distance
3334
COMPARE_PERFORMANCE ( COMPARE_SEQUENCES.out.matches.map { it = [ [id: val_md5sum_version], file(it) ] }, ch_observed_abundances, ch_expected_abundances )
3435

35-
// (3) input: val_md5sum_version, "${val_md5sum_version}_matches.txt", ch_detected_abundances, ch_expected_abundances
36-
// -> compare abundance of exact matches
37-
// -> per sample % total deviation from expected, potentially normalized somehow
38-
//COMPARISON_ABUNDANCES ( val_md5sum_version, "${val_md5sum_version}_exact_matches.txt", ch_detected_abundances, ch_expected_abundances )
39-
4036
emit:
4137
mismatch_barplot_png = COMPARE_SEQUENCES.out.png.collect()
4238
mismatch_barplot_svg = COMPARE_SEQUENCES.out.svg.collect()

0 commit comments

Comments
 (0)