diff --git a/.github/actions/nf-test/action.yml b/.github/actions/nf-test/action.yml index ad686e8..5e66435 100644 --- a/.github/actions/nf-test/action.yml +++ b/.github/actions/nf-test/action.yml @@ -19,6 +19,16 @@ inputs: runs: using: "composite" steps: + # The pipeline pulls ~17 container images, three of them heavy bioconductor/R stacks, which + # overruns the runner disk part-way through (docker exit 125, "no space left on device"). + # Reclaim the preinstalled toolchain first: the RunsOn default image is ubuntu24-full-x64, so the + # same Android/.NET/Haskell bloat as a GitHub-hosted runner is present and removable. + - name: Free disk space + uses: jlumbroso/free-disk-space@54081f138730dfa15788a46383842cd2f914a1be # v1.3.1 + with: + tool-cache: false # later steps need the toolchain + large-packages: false # slowest sweep, smallest reliable gain + - name: Setup Nextflow uses: nf-core/setup-nextflow@b4ec1bc7c16a94435159de94a05253542fddf6ef # v3 with: diff --git a/conf/modules.config b/conf/modules.config index 0dd3406..7fb9b4b 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -43,6 +43,9 @@ process { // --- Reference-derived, merged into intermediate_files --- withName: 'BWA_INDEX' { + cpus = 1 + memory = 1.GB + time = 30.min publishDir = [ path: "${params.outdir}/intermediate_files", mode: 'copy', diff --git a/modules/local/dmsanalysis/error_correction_false_doubles/templates/error_correction_false_doubles.R b/modules/local/dmsanalysis/error_correction_false_doubles/templates/error_correction_false_doubles.R index 9282ca3..4828a4c 100644 --- a/modules/local/dmsanalysis/error_correction_false_doubles/templates/error_correction_false_doubles.R +++ b/modules/local/dmsanalysis/error_correction_false_doubles/templates/error_correction_false_doubles.R @@ -21,7 +21,7 @@ suppressMessages({ ## ------------------------------------------------------------------------------------------------- seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw, input_count_path_processed, output_file_path, seq_error_rate_path, codon_window){ - ## load data (nucleotide-level counts from GATK) + ## load data (nucleotide-level counts) # WT sequence wt.seq <- readDNAStringSet(wt_path) @@ -57,11 +57,53 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw all.false.doubles <- all.false.doubles[which(all.false.doubles[,"varying_bases"] >= 3 & all.false.doubles[,"varying_bases"] < 5 & all.false.doubles[,"varying_codons"] == 2),,drop = F] all.false.doubles\$codon_dist <- vapply(regmatches(all.false.doubles[,"codon_mut"], gregexpr("\\\\d+(?=:)", all.false.doubles[,"codon_mut"], perl = TRUE)), function(z) abs(diff(as.integer(z))), integer(1)) + ### additional filters: + + #### remove outlier double codon mutants with very high counts + all.false.doubles <- all.false.doubles[which(all.false.doubles\$counts <= quantile(all.false.doubles\$counts, c(0.995))),] + + #### only keep 2+1 nt and 3+1 nt double codon mutants + all.false.doubles.nt <- all.false.doubles\$codon_mut + all.false.doubles.nt <- str_split_fixed(all.false.doubles.nt, ", ", 2) + all.false.doubles.nt[,1] <- str_split_fixed(all.false.doubles.nt[,1], ":", 2)[,2] + all.false.doubles.nt[,2] <- str_split_fixed(all.false.doubles.nt[,2], ":", 2)[,2] + classify_double_codon_variant <- function(x) { + parts <- strsplit(x, ">", fixed = TRUE) + n_changes <- vapply(parts, function(p) { + if (length(p) != 2L || nchar(p[1]) != nchar(p[2])) { + return(NA_integer_) + } + sum(strsplit(p[1], "")[[1]] != strsplit(p[2], "")[[1]]) + }, integer(1)) + labels <- c(1, 2, 3) + ifelse(n_changes %in% 1:3, + labels[n_changes], + ifelse(n_changes == 0L, "no change", NA_character_)) + } + all.false.doubles.nt[,1] <- classify_double_codon_variant(all.false.doubles.nt[,1]) + all.false.doubles.nt[,2] <- classify_double_codon_variant(all.false.doubles.nt[,2]) + if(length(which(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 2) != 0)){ + all.false.doubles <- all.false.doubles[-which(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 2),] + all.false.doubles.nt <- all.false.doubles.nt[-which(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 2),] + } + if(length(which(all.false.doubles.nt[,1] == 1 & all.false.doubles.nt[,2] == 1) != 0)){ + all.false.doubles <- all.false.doubles[-which(all.false.doubles.nt[,1] == 1 & all.false.doubles.nt[,2] == 1),] + all.false.doubles.nt <- all.false.doubles.nt[-which(all.false.doubles.nt[,1] == 1 & all.false.doubles.nt[,2] == 1),] + } + if(length(which(c(all.false.doubles.nt[,1] == 3 & all.false.doubles.nt[,2] == 2) | c(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 3)) != 0)){ + all.false.doubles <- all.false.doubles[-which(c(all.false.doubles.nt[,1] == 3 & all.false.doubles.nt[,2] == 2) | c(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 3)),] + all.false.doubles.nt <- all.false.doubles.nt[-which(c(all.false.doubles.nt[,1] == 3 & all.false.doubles.nt[,2] == 2) | c(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 3)),] + } + if(length(which(is.na(all.false.doubles.nt[,1]) == T | is.na(all.false.doubles.nt[,2]) == T) != 0)){ + all.false.doubles <- all.false.doubles[-which(is.na(all.false.doubles.nt[,1]) == T | is.na(all.false.doubles.nt[,2]) == T),] + all.false.doubles.nt <- all.false.doubles.nt[-which(is.na(all.false.doubles.nt[,1]) == T | is.na(all.false.doubles.nt[,2]) == T),] + } + all.false.doubles.codons <- str_split_fixed(all.false.doubles\$codon_mut, ", ", 2) all.false.doubles.codons[,1] <- as.integer(sub(":.*", "", all.false.doubles.codons[,1])) all.false.doubles.codons[,2] <- as.integer(sub(":.*", "", all.false.doubles.codons[,2])) - median.high.conf.coverage.per.pos <- rep(NA, max(as.numeric(str_split_fixed(input.counts.processed\$codon_mut, ":", 2)[,1]))) - names(median.high.conf.coverage.per.pos) <- 1:max(as.numeric(str_split_fixed(input.counts.processed\$codon_mut, ":", 2)[,1])) + median.high.conf.coverage.per.pos <- rep(NA, max.codon) + names(median.high.conf.coverage.per.pos) <- 1:max.codon for(i in 1:length(median.high.conf.coverage.per.pos)){ median.high.conf.coverage.per.pos[i] <- median(input.counts.processed[which(names(median.high.conf.coverage.per.pos)[i] == str_split_fixed(input.counts.processed\$codon_mut, ":", 2)[,1]),"cov"]) } @@ -86,11 +128,16 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw return(invisible(NULL)) } + ## keep the full distance range before dropping rows, so the lookup below stays long enough + .max_dist <- max(c(all.false.doubles\$codon_dist, codon_window), na.rm = TRUE) + ## fit on the finite rows only - lm() errors on the -Inf that log() gives a zero-coverage double + all.false.doubles <- all.false.doubles[.usable, , drop = FALSE] + fit <- lm(log(all.false.doubles\$perc_cov_to_max) ~ all.false.doubles\$codon_dist + I(all.false.doubles\$codon_dist^2)) all.false.doubles\$pred_cov_perc <- exp(predict(fit)) ### convert this into a look-up table: "% max. possible coverage" depending on codon-codon distance - dual.codon.coverage.estimate <- matrix(NA, nrow = max(all.false.doubles\$codon_dist), ncol = 2) + dual.codon.coverage.estimate <- matrix(NA, nrow = .max_dist, ncol = 2) colnames(dual.codon.coverage.estimate) <- c("Codon distance", "Estimate % max. possible coverage") dual.codon.coverage.estimate[,"Codon distance"] <- 1:nrow(dual.codon.coverage.estimate) ## deepmutscan fix (Maxi-approved deviation from the verbatim script): fill the coverage-decay @@ -104,8 +151,8 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw .cf <- coef(fit) dual.codon.coverage.estimate[,"Estimate % max. possible coverage"] <- exp(.cf[1] + .cf[2] * .dccd + .cf[3] * .dccd^2) - ## Process the GATK file, only look at single nucleotide variants - cat("Sequencing error correction of GATK counts...\\n") + ## Process the count file, only look at single nucleotide variants + cat("Sequencing error correction of raw counts...\\n") for (i in grep("[,]", input.counts.processed[,"base_mut"], invert = T)){ ### algorithm: @@ -114,7 +161,7 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw ### -> if there are zero false double mutants: skip / set correction factor to zero ### -> if there are any false double mutants: continue ### 3.) generate a table for ALL possible 2/3nt variants in the selected window: - ### true high-conf. variant count | true high-conf. variant coverage | false double double variant count (often 0 or 1) | false double double variant coverage + ### true high-conf. variant count | true high-conf. variant coverage | false double variant count (often 0 or 1) | false double variant coverage ### 4.) Maximum likelihood estimation (MLE) over all events ### 5.) Apply correction factor @@ -123,26 +170,13 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw seq.error.rate[tmp.single, "1nt counts/coverage"] <- input.counts.processed[i,"counts_per_cov"] ## 2.) look for 2/3nt hits in up to 40 codons (120 bp) upstream or downstream - tmp.false.doubles <- input.counts.raw[grep(paste0("(?= 3 & tmp.false.doubles[,"varying_bases"] < 5 & tmp.false.doubles[,"varying_codons"] == 2),,drop = F] - if(nrow(tmp.false.doubles) == 0){ - next - } + tmp.false.doubles <- all.false.doubles[grep(paste0("(?= tmp.min.from.ref & tmp.double.codon.n1 <= tmp.max.from.ref) & (tmp.double.codon.n2 >= tmp.min.from.ref & tmp.double.codon.n2 <= tmp.max.from.ref) - tmp.false.doubles <- tmp.false.doubles[which(tmp.double.codon.within.range == T),,drop = F] - if(nrow(tmp.false.doubles) == 0){ - next - } + tmp.false.doubles <- tmp.false.doubles[which(tmp.false.doubles\$codon_dist < codon_window),,drop = F] ## 3.) generate a table for ALL possible 2/3nt variants in the selected window tmp.summary.table.entries <- input.counts.processed @@ -156,14 +190,14 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw colnames(tmp.summary.table) <- c("codon distance", "true high-conf. variant count", "true high-conf. variant coverage", - "false double double variant count", - "false double double variant coverage") + "false double variant count", + "false double variant coverage") rownames(tmp.summary.table) <- paste0(c(tmp.summary.table.entries\$codon_mut),", ", c(tmp.summary.table.entries\$base_mut)) tmp.summary.table[,"codon distance"] <- abs(as.numeric(str_split_fixed(tmp.summary.table.entries\$codon_mut, ":", 2)[,1]) - tmp.ref.codon) tmp.summary.table[,"true high-conf. variant count"] <- tmp.summary.table.entries\$counts tmp.summary.table[,"true high-conf. variant coverage"] <- tmp.summary.table.entries\$cov - tmp.summary.table[,"false double double variant count"] <- 0 + tmp.summary.table[,"false double variant count"] <- 0 ## input the actual hits tmp.double.codon <- str_split_fixed(tmp.false.doubles\$codon_mut, ", ", 2) @@ -182,8 +216,8 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw } } - tmp.summary.table[tmp.match,"false double double variant count"] <- tmp.false.doubles[,"counts"] - tmp.summary.table[tmp.match,"false double double variant coverage"] <- tmp.false.doubles[,"cov"] + tmp.summary.table[tmp.match,"false double variant count"] <- tmp.false.doubles[,"counts"] + tmp.summary.table[tmp.match,"false double variant coverage"] <- tmp.false.doubles[,"cov"] tmp.summary.table <- as.data.frame(tmp.summary.table) ## 4.) Maximum likelihood estimation (MLE) over all events, iterating over codon positions @@ -197,20 +231,22 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw tmp.name <- rownames(tmp.summary.table.pos)[j] tmp.summary.table.pos[j,"true high-conf. variant count"] <- sum(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"true high-conf. variant count"]) tmp.summary.table.pos[j,"true high-conf. variant coverage"] <- max(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"true high-conf. variant coverage"]) - tmp.summary.table.pos[j,"false double double variant count"] <- sum(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double double variant count"]) - tmp.summary.table.pos[j,"false double double variant coverage"] <- max(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double double variant coverage"], na.rm = T) + tmp.summary.table.pos[j,"false double variant count"] <- sum(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double variant count"]) + tmp.summary.table.pos[j,"false double variant coverage"] <- max(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double variant coverage"], na.rm = T) ## infer a good proxy for the false double variant coverage with 0 counts (based estimated distance-dependent % coverage decay observed across all observed double mutants) - if(tmp.summary.table.pos[j,"false double double variant coverage"] == "-Inf"){ - tmp.summary.table.pos[j,"false double double variant coverage"] <- c(dual.codon.coverage.estimate[abs(tmp.summary.table.pos[j,"codon distance"]),"Estimate % max. possible coverage"] / 100) * tmp.summary.table.pos[j,"true high-conf. variant coverage"] - tmp.summary.table.pos[j,"false double double variant coverage"] <- round(tmp.summary.table.pos[j,"false double double variant coverage"]) + if(tmp.summary.table.pos[j,"false double variant coverage"] == "-Inf"){ + tmp.summary.table.pos[j,"false double variant coverage"] <- c(dual.codon.coverage.estimate[abs(tmp.summary.table.pos[j,"codon distance"]),"Estimate % max. possible coverage"] / 100) * tmp.summary.table.pos[j,"true high-conf. variant coverage"] + tmp.summary.table.pos[j,"false double variant coverage"] <- round(tmp.summary.table.pos[j,"false double variant coverage"]) } } tmp.summary.table.pos <- as.data.frame(tmp.summary.table.pos) - e_MLE <- sum(tmp.summary.table.pos\$`false double double variant count`) / - sum(tmp.summary.table.pos\$`true high-conf. variant count` * c(tmp.summary.table.pos\$`false double double variant coverage` / tmp.summary.table.pos\$`true high-conf. variant coverage`)) + e_MLE <- sum(tmp.summary.table.pos\$`false double variant count`) / + sum(tmp.summary.table.pos\$`true high-conf. variant count` * c(tmp.summary.table.pos\$`false double variant coverage` / tmp.summary.table.pos\$`true high-conf. variant coverage`)) seq.error.rate[tmp.single,"1nt false counts/coverage"] <- e_MLE - input.counts.processed[i,"counts_per_cov_corrected"] <- input.counts.processed[i,"counts_per_cov"] - e_MLE + ## clamped here too - the guard below only tests the rounded count, so a subtraction landing in + ## (-0.5, 0] rounds to 0, skips the guard, and leaves a negative frequency behind + input.counts.processed[i,"counts_per_cov_corrected"] <- max(0, input.counts.processed[i,"counts_per_cov"] - e_MLE) ## 5.) Apply correction factor input.counts.processed[i,"counts_corrected"] <- input.counts.processed[i,"counts"] - c(e_MLE * input.counts.processed[i,"cov"]) @@ -235,7 +271,7 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw ## ------------------------------------------------------------------------------------------------- seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, input_count_path_processed, output_file_path, seq_error_rate_path, codon_window){ - ## load data (nucleotide-level counts from GATK) + ## load data (nucleotide-level counts) # WT sequence wt.seq <- readDNAStringSet(wt_path) @@ -271,6 +307,48 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, all.false.doubles <- all.false.doubles[which(all.false.doubles[,"varying_bases"] >= 3 & all.false.doubles[,"varying_bases"] < 5 & all.false.doubles[,"varying_codons"] == 2),,drop = F] all.false.doubles\$codon_dist <- vapply(regmatches(all.false.doubles[,"codon_mut"], gregexpr("\\\\d+(?=:)", all.false.doubles[,"codon_mut"], perl = TRUE)), function(z) abs(diff(as.integer(z))), integer(1)) + ### additional filters: + + #### remove outlier double codon mutants with very high counts + all.false.doubles <- all.false.doubles[which(all.false.doubles\$counts <= quantile(all.false.doubles\$counts, c(0.995))),] + + #### only keep 2+1 nt and 3+1 nt double codon mutants + all.false.doubles.nt <- all.false.doubles\$codon_mut + all.false.doubles.nt <- str_split_fixed(all.false.doubles.nt, ", ", 2) + all.false.doubles.nt[,1] <- str_split_fixed(all.false.doubles.nt[,1], ":", 2)[,2] + all.false.doubles.nt[,2] <- str_split_fixed(all.false.doubles.nt[,2], ":", 2)[,2] + classify_double_codon_variant <- function(x) { + parts <- strsplit(x, ">", fixed = TRUE) + n_changes <- vapply(parts, function(p) { + if (length(p) != 2L || nchar(p[1]) != nchar(p[2])) { + return(NA_integer_) + } + sum(strsplit(p[1], "")[[1]] != strsplit(p[2], "")[[1]]) + }, integer(1)) + labels <- c(1, 2, 3) + ifelse(n_changes %in% 1:3, + labels[n_changes], + ifelse(n_changes == 0L, "no change", NA_character_)) + } + all.false.doubles.nt[,1] <- classify_double_codon_variant(all.false.doubles.nt[,1]) + all.false.doubles.nt[,2] <- classify_double_codon_variant(all.false.doubles.nt[,2]) + if(length(which(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 2) != 0)){ + all.false.doubles <- all.false.doubles[-which(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 2),] + all.false.doubles.nt <- all.false.doubles.nt[-which(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 2),] + } + if(length(which(all.false.doubles.nt[,1] == 1 & all.false.doubles.nt[,2] == 1) != 0)){ + all.false.doubles <- all.false.doubles[-which(all.false.doubles.nt[,1] == 1 & all.false.doubles.nt[,2] == 1),] + all.false.doubles.nt <- all.false.doubles.nt[-which(all.false.doubles.nt[,1] == 1 & all.false.doubles.nt[,2] == 1),] + } + if(length(which(c(all.false.doubles.nt[,1] == 3 & all.false.doubles.nt[,2] == 2) | c(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 3)) != 0)){ + all.false.doubles <- all.false.doubles[-which(c(all.false.doubles.nt[,1] == 3 & all.false.doubles.nt[,2] == 2) | c(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 3)),] + all.false.doubles.nt <- all.false.doubles.nt[-which(c(all.false.doubles.nt[,1] == 3 & all.false.doubles.nt[,2] == 2) | c(all.false.doubles.nt[,1] == 2 & all.false.doubles.nt[,2] == 3)),] + } + if(length(which(is.na(all.false.doubles.nt[,1]) == T | is.na(all.false.doubles.nt[,2]) == T) != 0)){ + all.false.doubles <- all.false.doubles[-which(is.na(all.false.doubles.nt[,1]) == T | is.na(all.false.doubles.nt[,2]) == T),] + all.false.doubles.nt <- all.false.doubles.nt[-which(is.na(all.false.doubles.nt[,1]) == T | is.na(all.false.doubles.nt[,2]) == T),] + } + all.false.doubles.codons <- str_split_fixed(all.false.doubles\$codon_mut, ", ", 2) all.false.doubles.codons[,1] <- as.integer(sub(":.*", "", all.false.doubles.codons[,1])) all.false.doubles.codons[,2] <- as.integer(sub(":.*", "", all.false.doubles.codons[,2])) @@ -297,11 +375,16 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, return(invisible(NULL)) } + ## keep the full distance range before dropping rows, so the lookup below stays long enough + .max_dist <- max(c(all.false.doubles\$codon_dist, codon_window), na.rm = TRUE) + ## fit on the finite rows only - lm() errors on the -Inf that log() gives a zero-coverage double + all.false.doubles <- all.false.doubles[.usable, , drop = FALSE] + fit <- lm(log(all.false.doubles\$perc_cov_to_max) ~ all.false.doubles\$codon_dist + I(all.false.doubles\$codon_dist^2)) all.false.doubles\$pred_cov_perc <- exp(predict(fit)) ### convert this into a look-up table: "% max. possible coverage" depending on codon-codon distance - dual.codon.coverage.estimate <- matrix(NA, nrow = max(all.false.doubles\$codon_dist), ncol = 2) + dual.codon.coverage.estimate <- matrix(NA, nrow = .max_dist, ncol = 2) colnames(dual.codon.coverage.estimate) <- c("Codon distance", "Estimate % max. possible coverage") dual.codon.coverage.estimate[,"Codon distance"] <- 1:nrow(dual.codon.coverage.estimate) ## deepmutscan fix (Maxi-approved deviation from the verbatim script): fill the coverage-decay @@ -315,8 +398,8 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, .cf <- coef(fit) dual.codon.coverage.estimate[,"Estimate % max. possible coverage"] <- exp(.cf[1] + .cf[2] * .dccd + .cf[3] * .dccd^2) - ## Process the GATK file, only look at single nucleotide variants - cat("Sequencing error correction of GATK counts...\\n") + ## Process the count file, only look at single nucleotide variants + cat("Sequencing error correction of raw counts...\\n") out.summary <- matrix(NA, ncol = 7, nrow = length(grep("[,]", input.counts.processed[,"base_mut"], invert = T))) colnames(out.summary) <- c("SNV", "Category", "False_double_variant_count", "False_double_variant_coverage", @@ -340,7 +423,7 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, ### -> if there are zero false double mutants: skip / set correction factor to zero ### -> if there are any false double mutants: continue ### 3.) generate a table for ALL possible 2/3nt variants in the selected window: - ### true high-conf. variant count | true high-conf. variant coverage | false double double variant count (often 0 or 1) | false double double variant coverage + ### true high-conf. variant count | true high-conf. variant coverage | false double variant count (often 0 or 1) | false double variant coverage ### 4.) Iterate of all positions ### 5.) Calculate the empirical Bayes ### 6.) Apply correction factor @@ -350,26 +433,11 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, seq.error.rate[tmp.single, "1nt counts/coverage"] <- input.counts.processed[i,"counts_per_cov"] ## 2.) look for 2/3nt hits in up to 40 codons (120 bp) upstream or downstream - tmp.false.doubles <- input.counts.raw[grep(paste0("(?= 3 & tmp.false.doubles[,"varying_bases"] < 5 & tmp.false.doubles[,"varying_codons"] == 2),,drop = F] - if(nrow(tmp.false.doubles) == 0){ - next - } - - ## set and enforce distance threshold + tmp.false.doubles <- all.false.doubles[grep(paste0("(?= tmp.min.from.ref & tmp.double.codon.n1 <= tmp.max.from.ref) & (tmp.double.codon.n2 >= tmp.min.from.ref & tmp.double.codon.n2 <= tmp.max.from.ref) - tmp.false.doubles <- tmp.false.doubles[which(tmp.double.codon.within.range == T),,drop = F] - if(nrow(tmp.false.doubles) == 0){ - next - } + tmp.false.doubles <- tmp.false.doubles[which(tmp.false.doubles\$codon_dist < codon_window),,drop = F] ## 3.) generate a table for ALL possible 2/3nt variants in the selected window tmp.summary.table.entries <- input.counts.processed @@ -383,14 +451,14 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, colnames(tmp.summary.table) <- c("codon distance", "true high-conf. variant count", "true high-conf. variant coverage", - "false double double variant count", - "false double double variant coverage") + "false double variant count", + "false double variant coverage") rownames(tmp.summary.table) <- paste0(c(tmp.summary.table.entries\$codon_mut),", ", c(tmp.summary.table.entries\$base_mut)) tmp.summary.table[,"codon distance"] <- abs(as.numeric(str_split_fixed(tmp.summary.table.entries\$codon_mut, ":", 2)[,1]) - tmp.ref.codon) tmp.summary.table[,"true high-conf. variant count"] <- tmp.summary.table.entries\$counts tmp.summary.table[,"true high-conf. variant coverage"] <- tmp.summary.table.entries\$cov - tmp.summary.table[,"false double double variant count"] <- 0 + tmp.summary.table[,"false double variant count"] <- 0 ## input the actual hits tmp.double.codon <- str_split_fixed(tmp.false.doubles\$codon_mut, ", ", 2) @@ -409,8 +477,8 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, } } - tmp.summary.table[tmp.match,"false double double variant count"] <- tmp.false.doubles[,"counts"] - tmp.summary.table[tmp.match,"false double double variant coverage"] <- tmp.false.doubles[,"cov"] + tmp.summary.table[tmp.match,"false double variant count"] <- tmp.false.doubles[,"counts"] + tmp.summary.table[tmp.match,"false double variant coverage"] <- tmp.false.doubles[,"cov"] tmp.summary.table <- as.data.frame(tmp.summary.table) ## 4.) Iterate over codon positions @@ -424,28 +492,30 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, tmp.name <- rownames(tmp.summary.table.pos)[j] tmp.summary.table.pos[j,"true high-conf. variant count"] <- sum(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"true high-conf. variant count"]) tmp.summary.table.pos[j,"true high-conf. variant coverage"] <- max(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"true high-conf. variant coverage"]) - tmp.summary.table.pos[j,"false double double variant count"] <- sum(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double double variant count"]) - tmp.summary.table.pos[j,"false double double variant coverage"] <- max(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double double variant coverage"], na.rm = T) + tmp.summary.table.pos[j,"false double variant count"] <- sum(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double variant count"]) + tmp.summary.table.pos[j,"false double variant coverage"] <- max(tmp.summary.table[grep(paste0("^",tmp.name), rownames(tmp.summary.table)),"false double variant coverage"], na.rm = T) ## infer a good proxy for the false double variant coverage with 0 counts (based estimated distance-dependent % coverage decay observed across all observed double mutants) - if(tmp.summary.table.pos[j,"false double double variant coverage"] == "-Inf"){ - tmp.summary.table.pos[j,"false double double variant coverage"] <- c(dual.codon.coverage.estimate[abs(tmp.summary.table.pos[j,"codon distance"]),"Estimate % max. possible coverage"] / 100) * tmp.summary.table.pos[j,"true high-conf. variant coverage"] - tmp.summary.table.pos[j,"false double double variant coverage"] <- round(tmp.summary.table.pos[j,"false double double variant coverage"]) + if(tmp.summary.table.pos[j,"false double variant coverage"] == "-Inf"){ + tmp.summary.table.pos[j,"false double variant coverage"] <- c(dual.codon.coverage.estimate[abs(tmp.summary.table.pos[j,"codon distance"]),"Estimate % max. possible coverage"] / 100) * tmp.summary.table.pos[j,"true high-conf. variant coverage"] + tmp.summary.table.pos[j,"false double variant coverage"] <- round(tmp.summary.table.pos[j,"false double variant coverage"]) } } tmp.summary.table.pos <- as.data.frame(tmp.summary.table.pos) ### fill the master table - out.summary[match(tmp.single,out.summary[,"SNV"]),"False_double_variant_count"] <- sum(tmp.summary.table.pos\$`false double double variant count`) - out.summary[match(tmp.single,out.summary[,"SNV"]),"False_double_variant_coverage"] <- sum(tmp.summary.table.pos\$`false double double variant coverage`) + out.summary[match(tmp.single,out.summary[,"SNV"]),"False_double_variant_count"] <- sum(tmp.summary.table.pos\$`false double variant count`) + out.summary[match(tmp.single,out.summary[,"SNV"]),"False_double_variant_coverage"] <- sum(tmp.summary.table.pos\$`false double variant coverage`) out.summary[match(tmp.single,out.summary[,"SNV"]),"True_high_conf_variant_count"] <- sum(tmp.summary.table.pos\$`true high-conf. variant count`) out.summary[match(tmp.single,out.summary[,"SNV"]),"True_high_conf_variant_coverage"] <- sum(tmp.summary.table.pos\$`true high-conf. variant coverage`) - out.summary[match(tmp.single,out.summary[,"SNV"]),"Exposure_m"] <- sum(tmp.summary.table.pos\$`true high-conf. variant count` * c(tmp.summary.table.pos\$`false double double variant coverage` / tmp.summary.table.pos\$`true high-conf. variant coverage`)) + out.summary[match(tmp.single,out.summary[,"SNV"]),"Exposure_m"] <- sum(tmp.summary.table.pos\$`true high-conf. variant count` * c(tmp.summary.table.pos\$`false double variant coverage` / tmp.summary.table.pos\$`true high-conf. variant coverage`)) } ## 5.) Calculate the empirical-Bayes sequencing error estimates - out.summary <- out.summary[-which(is.na(out.summary[,"False_double_variant_count"]) == T),] + if(length(which(is.na(out.summary[,"False_double_variant_count"]) == T)) != 0){ + out.summary <- out.summary[-which(is.na(out.summary[,"False_double_variant_count"]) == T),] + } out.summary <- as.data.frame(out.summary) class(out.summary\$False_double_variant_count) <- "numeric" class(out.summary\$False_double_variant_coverage) <- "numeric" @@ -543,14 +613,10 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, ##### posterior mean = shrunk EB estimate out.summary\$e_EB <- out.summary\$post_shape / out.summary\$post_rate - ##### other posterior intervals - # out.summary\$e_EB_upper_1SD <- qgamma(0.68, shape = out.summary\$post_shape, rate = out.summary\$post_rate) - # out.summary\$e_EB_upper_2SD <- qgamma(0.95, shape = out.summary\$post_shape, rate = out.summary\$post_rate) - # out.summary\$e_EB_upper_3SD <- qgamma(0.997, shape = out.summary\$post_shape, rate = out.summary\$post_rate) - ## 6.) Apply correction factors systematically (use upper 95% percentile to correct, a.k.a. "conservative") seq.error.rate[match(out.summary\$SNV, rownames(seq.error.rate)),"1nt false counts/coverage"] <- out.summary\$e_EB - input.counts.processed[match(out.summary\$SNV, input.counts.processed\$base_mut),"counts_per_cov_corrected"] <- input.counts.processed[match(out.summary\$SNV, input.counts.processed\$base_mut),"counts_per_cov"] - out.summary\$e_EB + ## clamped here too - the guard below only tests the rounded count (see the MLE branch) + input.counts.processed[match(out.summary\$SNV, input.counts.processed\$base_mut),"counts_per_cov_corrected"] <- pmax(0, input.counts.processed[match(out.summary\$SNV, input.counts.processed\$base_mut),"counts_per_cov"] - out.summary\$e_EB) input.counts.processed[match(out.summary\$SNV, input.counts.processed\$base_mut),"counts_corrected"] <- input.counts.processed[match(out.summary\$SNV, input.counts.processed\$base_mut),"counts"] - c(out.summary\$e_EB * input.counts.processed[match(out.summary\$SNV, input.counts.processed\$base_mut),"cov"]) ## round to nearest integer, do not allow for negative counts diff --git a/modules/local/dmsanalysis/process_variant_counts/templates/process_variant_counts.R b/modules/local/dmsanalysis/process_variant_counts/templates/process_variant_counts.R index e1871b3..43b0093 100644 --- a/modules/local/dmsanalysis/process_variant_counts/templates/process_variant_counts.R +++ b/modules/local/dmsanalysis/process_variant_counts/templates/process_variant_counts.R @@ -82,6 +82,7 @@ filter_counts_by_codon_library <- function(counts_file_path, codon_library_path, } # Filter the variant-count table + counts_table\$codon_mut <- toupper(counts_table\$codon_mut) filtered_counts <- counts_table %>% filter(varying_codons == 1) %>% # Keep rows with single-codon mutations rowwise() %>% @@ -190,6 +191,7 @@ complete_prefiltered_counts <- function(possible_nnk_path, prefiltered_counts_pa mutate(codon_mut = paste0(Codon_Number, ":", wt_codon, ">", Variant)) # Merge both dataframes based on the codon_mut column (full join to include all) + prefiltered_counts\$codon_mut <- toupper(prefiltered_counts\$codon_mut) merged_data <- full_join(prefiltered_counts, possible_nnk, by = "codon_mut") # Fill missing values in counts_per_cov and counts with 0.0000001 diff --git a/nextflow.config b/nextflow.config index beb79d6..d134d79 100644 --- a/nextflow.config +++ b/nextflow.config @@ -11,8 +11,8 @@ params { input = null - min_counts = 10 - base_qual = 30 + min_counts = 1 + base_qual = 40 min_flank = 2 error_correction = 'false_doubles' false_doubles_method = 'mle' diff --git a/nextflow_schema.json b/nextflow_schema.json index dfe43e8..7716c0f 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -50,13 +50,13 @@ "type": "integer", "description": "minimum counts for variant to be recognized. All variants below min_counts will be set to 0", "minimum": 1, - "default": 10 + "default": 1 }, "base_qual": { "type": "integer", "description": "Minimum base quality (Phred) for a mutation to be counted by the variant counter.", "minimum": 0, - "default": 30 + "default": 40 }, "min_flank": { "type": "integer", diff --git a/tests/default.nf.test.snap b/tests/default.nf.test.snap index c8ecb88..40e04f9 100644 --- a/tests/default.nf.test.snap +++ b/tests/default.nf.test.snap @@ -342,26 +342,26 @@ "GID1A.sa:md5,1f697dcb3fedd2f3f29642598bdd795b", "possible_mutations.csv:md5,9ce59fe91b726b1d88587a13fa899567", "seq_error_rate.csv:md5,350a018c981eca2506ac6291a6e9e161", - "variantCounts_for_heatmaps.csv:md5,842bd5c00760d594e2a450aec929e4aa", - "variantCounts_for_heatmaps_error_corrected.csv:md5,842bd5c00760d594e2a450aec929e4aa", + "variantCounts_for_heatmaps.csv:md5,8794f9fc716405d1891d24af6b636a2f", + "variantCounts_for_heatmaps_error_corrected.csv:md5,8794f9fc716405d1891d24af6b636a2f", "seq_error_rate.csv:md5,350a018c981eca2506ac6291a6e9e161", - "variantCounts_for_heatmaps.csv:md5,44069d5abca76323ccaed3c182447104", - "variantCounts_for_heatmaps_error_corrected.csv:md5,44069d5abca76323ccaed3c182447104", + "variantCounts_for_heatmaps.csv:md5,63c358cb04210b9dda2416ab3352494a", + "variantCounts_for_heatmaps_error_corrected.csv:md5,63c358cb04210b9dda2416ab3352494a", "seq_error_rate.csv:md5,350a018c981eca2506ac6291a6e9e161", - "variantCounts_for_heatmaps.csv:md5,3a9ce900a567cb663bc71033d3771d8a", - "variantCounts_for_heatmaps_error_corrected.csv:md5,3a9ce900a567cb663bc71033d3771d8a", + "variantCounts_for_heatmaps.csv:md5,33c6c3af1825b9ea7fbc949a18bb38aa", + "variantCounts_for_heatmaps_error_corrected.csv:md5,33c6c3af1825b9ea7fbc949a18bb38aa", "seq_error_rate.csv:md5,350a018c981eca2506ac6291a6e9e161", - "variantCounts_for_heatmaps.csv:md5,f263448ec0d7ec68f6e23aa050b961e6", - "variantCounts_for_heatmaps_error_corrected.csv:md5,f263448ec0d7ec68f6e23aa050b961e6", + "variantCounts_for_heatmaps.csv:md5,5277b202849d27c0b1773cacd8aae416", + "variantCounts_for_heatmaps_error_corrected.csv:md5,5277b202849d27c0b1773cacd8aae416", "variant_counts_columns.tsv:md5,c94801255a553ce5ff4a00d3170f6768", "variant_counts_columns.tsv:md5,c94801255a553ce5ff4a00d3170f6768", "variant_counts_columns.tsv:md5,c94801255a553ce5ff4a00d3170f6768", "variant_counts_columns.tsv:md5,c94801255a553ce5ff4a00d3170f6768", - "seqdepth_curve.csv:md5,eceecf3eff21b17b1d818c6ccefa72d9", - "seqdepth_curve.csv:md5,04649f2dd97ae3f7f4d7a94efb299d78", - "seqdepth_curve.csv:md5,5bde4b2a2b0cce57f4696913f4cdf8e7", - "seqdepth_curve.csv:md5,0deb002959055376128c81bcfba8e0e2", - "error_correction_report.html:md5,45455144c768473a8242381380deeba8", + "seqdepth_curve.csv:md5,b3d9041a2812b56a813bb69f61c7163a", + "seqdepth_curve.csv:md5,47bc7210d898e07966ac9a7adfec9ebf", + "seqdepth_curve.csv:md5,148be1fe719cf8194e4539cba7a7361a", + "seqdepth_curve.csv:md5,96bbb8a5e431acfa97a2e3960cce1e32", + "error_correction_report.html:md5,802b1cca71ffb34758810cea748c73ea", "fastqc-status-check-heatmap.txt:md5,58a9580c11af4210b1af65753383e120", "fastqc_adapter_content_plot.txt:md5,2990ad34b6fa29151dd4b9a38828afbe", "fastqc_overrepresented_sequences_plot.txt:md5,c703a9df0ea3a5c1c5f515a1fd88bbe8", @@ -377,11 +377,11 @@ "multiqc_general_stats.txt:md5,7c2fd44ded06835128a208d3365aa23b" ], { - "edger_rowCount": 1231, - "limma_rowCount": 1231 + "edger_rowCount": 1230, + "limma_rowCount": 1230 } ], - "timestamp": "2026-07-22T20:08:15.261389", + "timestamp": "2026-08-25T00:14:40.583305", "meta": { "nf-test": "0.9.4", "nextflow": "26.04.1"