Skip to content

Commit 7326af5

Browse files
d4straubd4straub
committed
add perfomance metrics such as recall and precision
1 parent e651b99 commit 7326af5

4 files changed

Lines changed: 266 additions & 4 deletions

File tree

‎bin/compare_performance.r‎

Lines changed: 201 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,201 @@
1+
#!/usr/bin/env Rscript
2+
3+
# compare_sequences.r
4+
5+
# Get params and files from the command line
6+
args <- commandArgs(trailingOnly=TRUE)
7+
tag <- args[1] # tag for each sample
8+
alignmentFILE <- args[2] # alignment file; expected columns: query, target, mismatches_final
9+
obsabundFILE <- args[3] # observed abundance*
10+
expabundFILE <- args[4] # expected abundance*
11+
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.
12+
mismatch_threshold <- as.numeric(args[6])
13+
# 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
14+
15+
# function to produce statistics
16+
get_stats <- function(i_exp,i_obs,df,sample) {
17+
18+
# if none are expected, skip!
19+
if( length(i_exp)==0 ) {
20+
print( paste("No expected seq in sample",sample, "- skipping") ); next
21+
} else {
22+
print( paste(length(i_exp),"expected seq in sample",sample) )
23+
}
24+
25+
# stats
26+
TP <- intersect(i_exp, i_obs)
27+
FN <- setdiff(i_exp, i_obs)
28+
FP <- setdiff(i_obs, i_exp)
29+
if ( length(TP) > 0 ) {
30+
Fone <- ( 2 * length(TP)/length(i_obs) * length(TP)/length(i_exp) ) / ( length(TP)/length(i_obs) + length(TP)/length(i_exp) )
31+
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) )
32+
} else {
33+
Fone <- 0
34+
Fbeta <- 0
35+
}
36+
37+
# save
38+
types <- c(
39+
"detected",
40+
"expected",
41+
"TP",
42+
"FN",
43+
"FP",
44+
"recall",
45+
"precision",
46+
"F1",
47+
"Fbeta",
48+
"fdr",
49+
"jaccard",
50+
"TPs_exp",
51+
"FNs_exp",
52+
"FPs_obs"
53+
)
54+
values <- c(
55+
length(i_obs),
56+
length(i_exp),
57+
length(TP),
58+
length(FN),
59+
length(FP),
60+
length(TP)/length(i_exp),
61+
length(TP)/length(i_obs),
62+
Fone,
63+
Fbeta,
64+
length(FP)/(length(i_obs)),
65+
length(TP)/(length(i_obs)+length(FN)),
66+
paste(head(TP, n=100), collapse=','),
67+
paste(head(FN, n=100), collapse=','),
68+
paste(head(FP, n=100), collapse=',')
69+
)
70+
ids <- rep(sample, length(types))
71+
72+
df_append <- data.frame(
73+
sample=ids,
74+
type= types,
75+
value= values,
76+
stringsAsFactors=FALSE)
77+
df <- rbind( df, df_append )
78+
79+
return(df)
80+
}
81+
82+
# PREPARE INPUT
83+
84+
# Read sequence alignment table
85+
alignment = read.table( alignmentFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
86+
alignment <- alignment[alignment$target !="*" & alignment$mismatch_final <= mismatch_threshold,]
87+
88+
# Read pipeline's ASV abundance table (either from QIIME2 [skipping first line] or from previous steps)
89+
if (grepl("^# Constructed from biom file", readLines(obsabundFILE, n = 1))) {
90+
observed = read.table( obsabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE, comment.char = "", skip=1 )
91+
} else {
92+
observed = read.table( obsabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE, comment.char = "" )
93+
}
94+
colnames(observed)[1] <- "ID"
95+
96+
# select available samples
97+
observed_samples <- colnames(observed)[2:ncol(observed)]
98+
print(paste( "Observed samples:", paste(observed_samples,collapse=",")))
99+
SAMPLES <- observed_samples
100+
101+
# Read expected abundance table
102+
if( file.exists(expabundFILE) ) {
103+
exp = read.table( expabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
104+
colnames(exp)[1] <- "ID"
105+
# extract samples to analyse
106+
exp_samples <- colnames(exp)[2:ncol(exp)]
107+
print(paste( "Expected samples:", paste( exp_samples ,collapse=",")))
108+
SAMPLES <- exp_samples[exp_samples %in% observed_samples]
109+
print(paste( "Investigate samples:", paste( SAMPLES ,collapse=",")))
110+
}
111+
112+
# check if there are any samples to analyse
113+
if (length(SAMPLES) == 0) {
114+
stop("ERROR - Found no samples to investigate")
115+
}
116+
117+
# ANALYSE
118+
119+
# Initialize dataframe
120+
df <- data.frame(
121+
sample=character(),
122+
type=character(),
123+
value=character(),
124+
stringsAsFactors=FALSE)
125+
126+
# for each sample
127+
for (sample in SAMPLES) {
128+
129+
# keep only non-empty columns and sample abundances
130+
keep_cols <- c("ID",sample)
131+
s_exp <- subset(exp, select = keep_cols)
132+
s_obs <- subset(observed, select = keep_cols)
133+
134+
# keep only detected (abundance > 0)
135+
s_exp = s_exp[s_exp[,2] > 0,]
136+
s_obs = s_obs[s_obs[,2] > 0,]
137+
138+
# keep only alignments where both, query and target are present
139+
s_alignment <- alignment[alignment$query %in% s_obs$ID & alignment$target %in% s_exp$ID,]
140+
141+
# translate s_obs$ID to alignment$target if available
142+
obsIDs <- c()
143+
for (obsID in s_obs$ID){
144+
i_alignment <- s_alignment[s_alignment$query == obsID,]
145+
# collapse several targets if multiple matches
146+
if( nrow(i_alignment) > 0 ){
147+
target <- paste(i_alignment$target, collapse="-")
148+
} else {
149+
target <- obsID
150+
}
151+
obsIDs <- c(obsIDs, target)
152+
}
153+
154+
# aggregate s_exp$ID if several alignment$query match
155+
expIDs <- s_exp$ID[!s_exp$ID %in% s_alignment$target]
156+
i_alignment <- aggregate(target ~ query, s_alignment, paste, collapse = "-")
157+
expIDs <- c(expIDs, i_alignment$target)
158+
159+
# calculate stats
160+
df <- get_stats( unique(expIDs), unique(obsIDs), df, sample)
161+
}
162+
163+
# Write detailed output
164+
outfile <- "performance_per-sample.tsv"
165+
df$tag <- rep(tag, nrow(df))
166+
print(paste("write",outfile))
167+
write.table(df, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
168+
169+
170+
# Make & write summary output
171+
outfile <- "performance_summary.tsv"
172+
df_sum <- data.frame(type=character(), median=numeric(), mean=numeric(), min=numeric(), max=numeric(), count=numeric(), stringsAsFactors=FALSE)
173+
print("Warnings are expected at that point:")
174+
for (type in unique(df$type)) {
175+
df_subset <- df[df$type == type, ]
176+
MEDIAN <- median( as.numeric(df_subset$value) )
177+
MEAN <- mean( as.numeric(df_subset$value) )
178+
MIN <- min( as.numeric(df_subset$value) )
179+
MAX <- max( as.numeric(df_subset$value) )
180+
COUNT <- length( as.numeric(df_subset$value) )
181+
df_sum[nrow(df_sum) + 1,] <- list(type, MEDIAN, MEAN, MIN, MAX, COUNT)
182+
}
183+
df_sum$tag <- rep(tag, nrow(df_sum))
184+
print(paste("write",outfile))
185+
write.table(df_sum, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
186+
187+
# Plot Values
188+
df_subset <- subset(df, type %in% c("recall","precision","F1","Fbeta","fdr","jaccard") )
189+
df_subset$value <- as.numeric(df_subset$value)
190+
191+
outfile <- "performance_boxplot.svg"
192+
print(paste("write",outfile))
193+
svg(outfile, height = 8, width = 10)
194+
boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1))
195+
invisible(dev.off())
196+
197+
outfile <- "performance_boxplot.png"
198+
print(paste("write",outfile))
199+
png(outfile, height = 400, width = 500)
200+
boxplot(value~type, data=df_subset, xlab="Type", ylab="Value", ylim = c(0, 1))
201+
invisible(dev.off())

‎conf/modules.config‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1241,6 +1241,23 @@ process {
12411241
]
12421242
}
12431243

1244+
withName: COMPARE_PERFORMANCE {
1245+
ext.mismatch_threshold = 0 // thereshold for comparing matches (0=exact match)
1246+
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.
1247+
publishDir = [
1248+
[
1249+
path: { "${params.outdir}/comparison/${meta.id}" },
1250+
mode: params.publish_dir_mode,
1251+
pattern: "*{performance_summary.tsv,.log,.svg,.png}"
1252+
],
1253+
[
1254+
path: { "${params.outdir}/comparison/${meta.id}/per-sample" },
1255+
mode: params.publish_dir_mode,
1256+
pattern: "*{_per-sample.tsv}"
1257+
]
1258+
]
1259+
}
1260+
12441261
withName: 'MULTIQC' {
12451262
ext.args = { params.multiqc_title ? "--title \"$params.multiqc_title\"" : '' }
12461263
publishDir = [

‎modules/local/compare_performance.nf‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
process COMPARE_PERFORMANCE {
2+
tag "${meta.id}"
3+
label 'process_single'
4+
5+
conda "conda-forge::r-base=4.2.1"
6+
container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
7+
'https://depot.galaxyproject.org/singularity/r-base:4.2.1' :
8+
'biocontainers/r-base:4.2.1' }"
9+
10+
input:
11+
tuple val(meta), path(alignment) // alignment file with "query,target,mismatch_final"
12+
path(observed_abundances) // ASV table from within pipeline
13+
path(expected_abundances) // ASV table from params.expected_abundances
14+
15+
output:
16+
path("*.svg") , emit: svg
17+
path("*.png") , emit: png
18+
path("performance_summary.tsv") , emit: performance
19+
path("performance_per-sample.tsv") , emit: performance_per_sample
20+
path("performance.log") , emit: log
21+
path("Warnings.txt") , emit: warnings
22+
path "versions.yml" , emit: versions, topic: versions
23+
24+
script:
25+
def fbeta = task.ext.fbeta ?: 2
26+
def mismatch_threshold = task.ext.mismatch_threshold ?: 0
27+
"""
28+
compare_performance.r \\
29+
"${meta.id}" \\
30+
"$alignment" \\
31+
"$observed_abundances" \\
32+
"$expected_abundances" \\
33+
"$fbeta" \\
34+
"$mismatch_threshold" \\
35+
>"performance.log"
36+
37+
# isolate warnings (when grep find no match, exit code is 1, "|| true" fixes that)
38+
grep "WARN -" "performance.log" | sed 's/^\\[1\\] //g' | sed 's/"//g' > Warnings.txt || true
39+
40+
cat <<-END_VERSIONS > versions.yml
41+
"${task.process}":
42+
R: \$(R --version | sed -n 1p | sed 's/R version //g' | sed 's/\\s.*\$//')
43+
END_VERSIONS
44+
"""
45+
}

‎subworkflows/local/comparison_wf.nf‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
include { VSEARCH_USEARCHGLOBAL as VSEARCH_USEARCHGLOBAL_BM } from '../../modules/nf-core/vsearch/usearchglobal/main'
22
include { COMPARE_SEQUENCES } from '../../modules/local/compare_sequences'
3+
include { COMPARE_PERFORMANCE } from '../../modules/local/compare_performance'
34

45
workflow COMPARISON_WF {
56
take:
@@ -28,10 +29,8 @@ workflow COMPARISON_WF {
2829
if( it.countLines() > 0 ) { log.warn "about comparing sequences\n\n" + it.splitText().join("") }
2930
}
3031

31-
// (2) input: val_md5sum_version, "${val_md5sum_version}_nucleotide_differences.tsv", ch_detected_abundances, ch_expected_abundances
32-
// -> compare exact matches per sample: "${val_md5sum_version}_nucleotide_differences.tsv" = 0 & prevalence ch_detected_abundances vs ch_expected_abundances
33-
// -> per sample count FN/FP/... & calculate F1/...
34-
//COMPARISON_SEQUENCES ( val_md5sum_version, COMPARE_SEQUENCES.out.matches, ch_detected_abundances, ch_expected_abundances )
32+
// Calculate performance metrics per sample such as precision, recall, F1 score
33+
COMPARE_PERFORMANCE ( COMPARE_SEQUENCES.out.matches.map { it = [ [id: val_md5sum_version], file(it) ] }, ch_observed_abundances, ch_expected_abundances )
3534

3635
// (3) input: val_md5sum_version, "${val_md5sum_version}_matches.txt", ch_detected_abundances, ch_expected_abundances
3736
// -> compare abundance of exact matches

0 commit comments

Comments
 (0)