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..f04d508 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) @@ -56,12 +56,54 @@ seq_error_correct_by_false_doubles_MLE <- function(wt_path, input_count_path_raw all.false.doubles <- input.counts.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"]) } @@ -104,8 +146,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 +156,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 +165,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 +185,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 +211,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,18 +226,18 @@ 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 @@ -235,7 +264,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) @@ -270,6 +299,48 @@ seq_error_correct_by_false_doubles_EB <- function(wt_path, input_count_path_raw, all.false.doubles <- input.counts.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])) @@ -315,8 +386,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 +411,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 +421,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 +439,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 +465,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 +480,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,11 +601,6 @@ 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 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'