Skip to content

Commit e522174

Browse files
tom-brekketom-brekke
authored
Merge branch 'dev' into GloSED_2
2 parents 68c013d + fd3b565 commit e522174

27 files changed

Lines changed: 792 additions & 102 deletions

‎CHANGELOG.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
77

88
### `Added`
99

10-
- [[#1022](https://github.com/nf-core/ampliseq/issues/1022)] - added new GloSED database for classification of Eukaryotic ITS. Not support for SH numbers as of yet. (by @tom-brekke)
10+
- [#1022](https://github.com/nf-core/ampliseq/issues/1022) - added new GloSED database for classification of Eukaryotic ITS. Not support for SH numbers as of yet. (by @tom-brekke)
11+
- [#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)
1112

1213
### `Changed`
1314

‎bin/compare_performance.r‎

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
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+
"observed",
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+
100+
# Read expected abundance table
101+
exp = read.table( expabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
102+
colnames(exp)[1] <- "ID"
103+
# extract samples to analyse
104+
exp_samples <- colnames(exp)[2:ncol(exp)]
105+
print(paste( "Expected samples:", paste( exp_samples ,collapse=",")))
106+
SAMPLES <- exp_samples[exp_samples %in% observed_samples]
107+
print(paste( "Investigate samples:", paste( SAMPLES ,collapse=",")))
108+
109+
# check if there are any samples to analyse
110+
if (length(SAMPLES) == 0) {
111+
stop("ERROR - Found no samples to investigate")
112+
}
113+
114+
# ANALYSE
115+
116+
# Initialize dataframe
117+
df <- data.frame(
118+
sample=character(),
119+
type=character(),
120+
value=character(),
121+
stringsAsFactors=FALSE)
122+
123+
# for each sample
124+
for (sample in SAMPLES) {
125+
126+
# keep only non-empty columns and sample abundances
127+
keep_cols <- c("ID",sample)
128+
s_exp <- subset(exp, select = keep_cols)
129+
s_obs <- subset(observed, select = keep_cols)
130+
131+
# keep only observed (abundance > 0)
132+
s_exp = s_exp[s_exp[,2] > 0,]
133+
s_obs = s_obs[s_obs[,2] > 0,]
134+
135+
# keep only alignments where both, query and target are present
136+
s_alignment <- alignment[alignment$query %in% s_obs$ID & alignment$target %in% s_exp$ID,]
137+
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+
}
150+
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)
155+
156+
# calculate stats
157+
df <- get_stats( unique(expIDs), unique(obsIDs), df, sample)
158+
}
159+
160+
# Write detailed output
161+
outfile <- "performance_per-sample.tsv"
162+
df$tag <- rep(tag, nrow(df))
163+
print(paste("write",outfile))
164+
write.table(df, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
165+
166+
167+
# Make & write summary output
168+
outfile <- "performance_summary.tsv"
169+
df_sum <- data.frame(type=character(), median=numeric(), mean=numeric(), min=numeric(), max=numeric(), count=numeric(), stringsAsFactors=FALSE)
170+
print("Warnings are expected at that point:")
171+
for (type in unique(df$type)) {
172+
df_subset <- df[df$type == type, ]
173+
MEDIAN <- median( as.numeric(df_subset$value) )
174+
MEAN <- mean( as.numeric(df_subset$value) )
175+
MIN <- min( as.numeric(df_subset$value) )
176+
MAX <- max( as.numeric(df_subset$value) )
177+
COUNT <- length( as.numeric(df_subset$value) )
178+
df_sum[nrow(df_sum) + 1,] <- list(type, MEDIAN, MEAN, MIN, MAX, COUNT)
179+
}
180+
df_sum$tag <- rep(tag, nrow(df_sum))
181+
print(paste("write",outfile))
182+
write.table(df_sum, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
183+
184+
# Plot Values
185+
df_subset <- subset(df, type %in% c("recall","precision","F1","Fbeta","fdr","jaccard") )
186+
df_subset$value <- as.numeric(df_subset$value)
187+
188+
outfile <- "performance_boxplot.svg"
189+
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))
192+
invisible(dev.off())
193+
194+
outfile <- "performance_boxplot.png"
195+
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))
198+
invisible(dev.off())

‎bin/compare_sequences.r‎

Lines changed: 169 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
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+
alignmentFILE <- args[1] # expected columns: "query+target+ql+tl+qilo+qihi+tilo+tihi+gaps+mism+qstrand", see https://torognes.github.io/vsearch/misc/vsearch-userfields.7.html
8+
obsabundFILE <- args[2] # observed abundance*
9+
expabundFILE <- args[3] # expected abundance*
10+
prefix <- args[4] # prefix string for output files
11+
similarity_threshold <- as.numeric(args[5]) # similarity threshold for alignmentFILE
12+
query_or_target <- args[6] # use mismatches to query or to target?
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+
# Input
16+
if (file.exists(expabundFILE)) {
17+
print(paste("Compare",obsabundFILE,"to",expabundFILE))
18+
} else {
19+
print(paste("Compare",obsabundFILE,"to all sequences"))
20+
}
21+
22+
# PREPARE
23+
24+
# Read sequence alignment table
25+
alignment = read.table( alignmentFILE, header = FALSE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
26+
colnames(alignment) <- c("query","target","ql","tl","qilo","qihi","tilo","tihi","gaps","mism","qstrand")
27+
28+
# Read pipeline's ASV abundance table (either from QIIME2 [skipping first line] or from previous steps)
29+
if (grepl("^# Constructed from biom file", readLines(obsabundFILE, n = 1))) {
30+
observed = read.table( obsabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE, comment.char = "", skip=1 )
31+
} else {
32+
observed = read.table( obsabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE, comment.char = "" )
33+
}
34+
colnames(observed)[1] <- "ID"
35+
36+
# select available samples
37+
observed_samples <- colnames(observed)[2:ncol(observed)]
38+
print(paste( "Observed samples:", paste(observed_samples,collapse=",")))
39+
SAMPLES <- observed_samples
40+
41+
# Read expected abundance table if it exists
42+
if (file.exists(expabundFILE)) {
43+
exp = read.table( expabundFILE, header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE, strip.white = TRUE)
44+
colnames(exp)[1] <- "ID"
45+
# extract samples to analyse
46+
exp_samples <- colnames(exp)[2:ncol(exp)]
47+
print(paste( "Expected samples:", paste( exp_samples ,collapse=",")))
48+
SAMPLES <- exp_samples[exp_samples %in% observed_samples]
49+
print(paste( "Investigate samples:", paste( SAMPLES ,collapse=",")))
50+
# warn if samples are missing due to incomplete overlap
51+
if ( !all(observed_samples %in% SAMPLES) ) {
52+
print(paste("WARN - sample(s) not expected (file:",expabundFILE,") and omitted:",paste(observed_samples[!observed_samples %in% SAMPLES], collapse="," )))
53+
}
54+
if ( !all(exp_samples %in% SAMPLES) ) {
55+
print(paste("WARN - sample(s) not observed and omitted:",paste(exp_samples[!exp_samples %in% SAMPLES], collapse="," )))
56+
}
57+
# warn about missing sequences
58+
alignment_target <- alignment$target[alignment$target != "*"]
59+
if( !all(alignment_target %in% exp$ID) ) {
60+
print(paste("WARN - Some expected sequences in alignment results are not in",expabundFILE,". The following are missing:",paste( unique(alignment_target[!alignment_target %in% exp$ID]),collapse=",")))
61+
}
62+
if( !all(exp$ID %in% alignment$target) ) {
63+
print(paste("WARN - Some expected sequences of",expabundFILE,"are not in the alignment results. The following are missing:",paste( unique(exp$ID[!exp$ID %in% alignment$target]),collapse=",")))
64+
}
65+
}
66+
67+
# check if there are any samples to analyse
68+
if (length(SAMPLES) == 0) {
69+
stop("ERROR - Found no samples to investigate")
70+
}
71+
72+
# Investigate alignment, required columns: query,target,ql,tl,qilo,qihi,tilo,tihi,gaps,mism
73+
alignment$qterminalgaps <- alignment$ql - (abs(alignment$qihi-alignment$qilo)+1)
74+
alignment$qmism <- alignment$gaps + alignment$qterminalgaps + alignment$mism
75+
alignment$tterminalgaps <- alignment$tl - (abs(alignment$tihi-alignment$tilo)+1)
76+
alignment$tmism <- alignment$gaps + alignment$tterminalgaps + alignment$mism
77+
if( query_or_target=="query" ) {
78+
alignment$mismatch_final <- alignment$qmism
79+
} else if( query_or_target=="target" ) {
80+
alignment$mismatch_final <- alignment$tmism
81+
} else if( query_or_target=="alignment" ) {
82+
alignment$mismatch_final <- alignment$gaps + alignment$mism
83+
} else {
84+
stop( paste("ERROR -",query_or_target,"is not valid (valid: query,target,alignment)") )
85+
}
86+
87+
# order for reproducibility
88+
alignment <- alignment[order(alignment$query, alignment$target),]
89+
90+
# write file
91+
outfile <- paste0(prefix,"nucleotide-differences.tsv")
92+
print(paste("write",outfile))
93+
write.table(alignment, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
94+
95+
# ANALYSE
96+
97+
# select matches above threshold
98+
matches_above_threshold <- alignment[alignment$target !="*",]
99+
100+
for (sample in SAMPLES) {
101+
if( file.exists(expabundFILE) ) {
102+
# filter expected sequences in that sample (matches_above_threshold$target) with expected abundances (abundance > 0) of exp$ID
103+
keep_cols <- c("ID",sample)
104+
print(paste("Found",nrow(exp),"expected sequences overall in",expabundFILE))
105+
s_exp <- subset(exp, select = keep_cols)
106+
s_exp = s_exp[s_exp[,2] > 0,]
107+
print(paste("Found",nrow(s_exp),"expected sequences in sample", sample,"in",expabundFILE))
108+
s_matches <- matches_above_threshold[matches_above_threshold$target %in% s_exp$ID,]
109+
} else {
110+
s_matches <- matches_above_threshold
111+
}
112+
113+
# sample will be skipped if s_matches has no data
114+
if( nrow(s_matches)==0 ) {
115+
print(paste("WARN - Skipping sample",sample,"because found no matches to expected sequences"))
116+
next
117+
} else {
118+
print(paste("Found",length(unique(s_matches$query)),"to be expected sequences of",length(unique(matches_above_threshold$query)),"in sample",sample))
119+
print(paste("Found",length(unique(s_matches$target)),"matches of",length(unique(matches_above_threshold$target)),"total in sample",sample))
120+
}
121+
122+
# select best match (sort by mismatches and retain only first unique entry)
123+
s_matches <- s_matches[order(s_matches$mismatch_final, as.numeric(s_matches$mismatch_final)), ]
124+
s_matches <- s_matches[!duplicated(s_matches$query), ]
125+
print(paste("Found",length(unique(s_matches$target)),"best matches of",length(unique(matches_above_threshold$target)),"total in sample",sample))
126+
127+
# filter for ASVs in that sample
128+
keep_cols <- c("ID",sample)
129+
s_observed <- subset(observed, select = keep_cols)
130+
# keep only observed entries in that sample (abundance > 0)
131+
s_observed = s_observed[s_observed[,2] > 0,]
132+
print(paste("Found",nrow(s_observed),"of",nrow(observed),"observed sequences in sample",sample))
133+
134+
# filter alignment result by observed
135+
s_matches <- s_matches[s_matches$query %in% s_observed$ID,]
136+
print(paste("Found",nrow(s_matches),"observed sequences with match in sample",sample))
137+
s_below_threshold <- s_observed[!s_observed$ID %in% s_matches$query,]
138+
print(paste("Found",nrow(s_below_threshold),"observed sequences without match in sample",sample))
139+
140+
# make barplot
141+
df <- as.data.frame( table(s_matches$mismatch_final) )
142+
df <- rbind(df, data.frame(Var1 = paste0(">=",(1-similarity_threshold)*100,"%"), Freq = length(unique(s_below_threshold$ID))))
143+
colnames(df) <- c("Mismatches to expected sequences","Number of sequences")
144+
145+
outfile <- paste0(prefix,sample,"_nucleotide-distance_barplot.svg")
146+
print(paste("write",outfile))
147+
svg(outfile, height = 4, width = 5)
148+
barplot(df[,2],
149+
names.arg=df[,1],
150+
main=paste("Sample",sample,"with",length(unique(s_observed$ID)),"observed sequences\n(",length(unique(s_matches$query)),"sequences above",similarity_threshold*100,"% similarity)"),
151+
xlab="Mismatches to expected sequences",
152+
ylab="Number of sequences")
153+
invisible(dev.off())
154+
155+
outfile <- paste0(prefix,sample,"_nucleotide-distance_barplot.png")
156+
print(paste("write",outfile))
157+
png(outfile, height = 400, width = 500)
158+
barplot(df[,2],
159+
names.arg=df[,1],
160+
main=paste("Sample",sample,"with",length(unique(s_observed$ID)),"observed sequences\n(",length(unique(s_matches$query)),"sequences above",similarity_threshold*100,"% similarity)"),
161+
xlab="Mismatches to expected sequences",
162+
ylab="Number of sequences")
163+
invisible(dev.off())
164+
165+
# write file
166+
outfile <- paste0(prefix,sample,"_nucleotide-differences_per-sample.tsv")
167+
print(paste("write",outfile))
168+
write.table(df, file = outfile, row.names = FALSE, col.names = TRUE, quote = FALSE, na = '', sep="\t")
169+
}

0 commit comments

Comments
 (0)