From 1b75298dde6dbf768809468d65e140cc8affffcb Mon Sep 17 00:00:00 2001 From: Dalil Taher Date: Mon, 22 Jun 2026 13:44:58 +0200 Subject: [PATCH 1/2] Add input validation (min binsize 30k) and skip cells with missing VCFs --- R/getAS_CNA.R | 63 ++++++++++++++++++++++++++++++++----------- R/run_sc_sequencing.R | 36 ++++++++++++++++--------- 2 files changed, 71 insertions(+), 28 deletions(-) diff --git a/R/getAS_CNA.R b/R/getAS_CNA.R index dbd6713..7f1ceb9 100644 --- a/R/getAS_CNA.R +++ b/R/getAS_CNA.R @@ -71,10 +71,29 @@ getAS_CNA <- function(res, readPhases <- function(phasing_paths) { - phasing <- lapply(phasing_paths,function(x) as.data.frame(data.table::fread(x))) + phasing <- lapply(phasing_paths, function(x) { + tryCatch({ + cmd_str <- if(grepl("\\.gz$", x, ignore.case=TRUE)) paste0("zgrep -v '^##' ", x) else paste0("grep -v '^##' ", x) + df <- as.data.frame(data.table::fread(cmd=cmd_str)) + if(ncol(df) > 0 && grepl("CHROM", colnames(df)[1], ignore.case=TRUE)) colnames(df)[1] <- "#CHROM" + df + }, error = function(e) { + tryCatch({ + as.data.frame(data.table::fread(x, skip="#CHROM")) + }, error = function(e2) { + as.data.frame(data.table::fread(x)) + }) + }) + }) phases <- lapply(phasing,function(x) { - x <- x[grep("0\\|1|1\\|0",x[, 10]), ] + if(ncol(x) >= 10 && !("REF" %in% colnames(x))) { + colnames(x)[4] <- "REF" + colnames(x)[5] <- "ALT" + } + if(ncol(x) < 10) return(list(chr=integer(0), pos=integer(0), phases1=integer(0), phases2=integer(0))) + x <- x[grep("0\\|1|1\\|0",x[, 10]), , drop=FALSE] + if(nrow(x) == 0) return(list(chr=integer(0), pos=integer(0), phases1=integer(0), phases2=integer(0))) phase <- gsub("(.*)\\|(.*)","\\1",x[,10]) phases1 <- x[,"REF"] phases1[phase=="1"] <- x[phase=="1","ALT"] @@ -426,17 +445,26 @@ getAS_CNA <- function(res, res$allProfiles_AS <- parallel::mclapply(1:length(res$allTracks.processed), function(x) { cat(".") + profile_to_use <- if(any(grepl("refitted",names(res)))) res$allProfiles.refitted.auto[[x]] else res$allProfiles[[x]] + if(is.null(profile_to_use) || inherits(profile_to_use, "try-error")) return(NULL) + + ac_p <- if(length(list_ac_counts_paths)==1) list_ac_counts_paths[[1]] else list_ac_counts_paths[[x]] + ph_p <- if(length(path_to_phases)==1) path_to_phases[[1]] else path_to_phases[[x]] + + if(!all(file.exists(ac_p)) || (!is.null(ph_p) && !all(file.exists(ph_p)))) { + warning(paste("Missing allele count or phasing files for cell", x, "- skipping AS CNA for this cell.")) + return(NULL) + } + getAS_CNA_sample(track=res$allTracks.processed[[x]], - profile=if(any(grepl("refitted",names(res)))) res$allProfiles.refitted.auto[[x]] else res$allProfiles[[x]], - ac_counts_paths=list_ac_counts_paths[[x]], - phases=phases, - purity=if(any(grepl("refitted",names(res)))) res$allSolutions.refitted.auto[[x]]$purity - else res$allSolutions[[x]]$purity, - ploidy=if(any(grepl("refitted",names(res)))) res$allSolutions.refitted.auto[[x]]$ploidy - else res$allSolutions[[x]]$ploidy, + profile=profile_to_use, + ac_counts_paths=ac_p, purs=purs[[x]], ploidies=ploidies[[x]], - path_to_phases=if(length(path_to_phases)>1) path_to_phases[[x]] else NULL, + purity=if(any(grepl("refitted",names(res)))) res$allProfiles.refitted.auto[[x]]$purity else res$allSolutions[[x]]$purity, + ploidy=if(any(grepl("refitted",names(res)))) res$allProfiles.refitted.auto[[x]]$ploidy else res$allSolutions[[x]]$ploidy, + phases=if(!is.null(phases)) phases, + path_to_phases=ph_p, steps=steps, betabinom=betabinom) },mc.cores=mc.cores) @@ -445,25 +473,30 @@ getAS_CNA <- function(res, print("Estimate best overdispersion parameters per cell") res$allProfiles_AS <- lapply(res$allProfiles_AS, function(x) { + if(is.null(x) || inherits(x, "try-error")) return(NULL) tmp <- get_best_overdispersion_profile(x$nprof.fixed) x$nprof.fixed <- tmp$prof x$overdispersion_best_fit <- tmp$best_rho x }) - odps <- sapply(res$allProfiles_AS,function(x) x$overdispersion_best_fit) - print(paste("Overdispersion quantiles:",quantile(odps,probs=seq(0,1,.1)),collapse=" ")) + odps <- sapply(res$allProfiles_AS,function(x) if(!is.null(x)) x$overdispersion_best_fit else NA) + odps <- odps[!is.na(odps)] + if(length(odps)>0) print(paste("Overdispersion quantiles:",paste(quantile(odps,probs=seq(0,1,.1),na.rm=TRUE),collapse=" "))) } print("## write to disk and plot Allele-specific Profiles") pdf(paste0(outdir,"/all_as_cna_profiles_",projectname,".pdf"),width=15,height=5) tnull <- lapply(1:length(res$allProfiles_AS), function(x) { + if(is.null(res$allProfiles_AS[[x]]) || inherits(res$allProfiles_AS[[x]], "try-error")) return(NULL) try({ plot_AS_profile(res$allProfiles_AS[[x]]$nprof.fixed) title(paste0(names(res$allTracks)[x]," - bam",x) ,cex=.5) }) - write.table(res$allProfiles_AS[[x]]$nprof.fixed, - sep="\t",col.names=T,row.names=F,quote=F, - file=paste0(outdir,"/as_cna_profile_",names(res$allTracks)[x],"_bam",x,".txt")) + try({ + write.table(res$allProfiles_AS[[x]]$nprof.fixed, + sep="\t",col.names=T,row.names=F,quote=F, + file=paste0(outdir,"/as_cna_profile_",names(res$allTracks)[x],"_bam",x,".txt")) + }) if(betabinom) { try({ diff --git a/R/run_sc_sequencing.R b/R/run_sc_sequencing.R index 4cef645..e34c3b4 100644 --- a/R/run_sc_sequencing.R +++ b/R/run_sc_sequencing.R @@ -16,7 +16,7 @@ run_sc_sequencing <- function(tumour_bams, barcodes_10x=NULL, normal_bams=NULL, outdir="./", - is_pdf=F, + is_pdf=FALSE, probs_filters=.1, path_to_phases=NULL, list_ac_counts_paths=NULL, @@ -33,6 +33,14 @@ run_sc_sequencing <- function(tumour_bams, betabinom=FALSE) { checkArguments_scs(c(as.list(environment()))) + + # --- Strict input validation guards --- + if(is.null(tumour_bams)) stop("tumour_bams cannot be NULL.") + if(!all(file.exists(tumour_bams))) stop("One or more files in tumour_bams do not exist.") + if(!is.null(normal_bams) && !all(file.exists(normal_bams))) stop("One or more files in normal_bams do not exist.") + build <- match.arg(build, choices = c("hg19", "hg38", "mm39")) + # -------------------------------------- + suppressPackageStartupMessages(require(parallel)) suppressPackageStartupMessages(require(Rsamtools)) suppressPackageStartupMessages(require(Biostrings)) @@ -54,6 +62,7 @@ run_sc_sequencing <- function(tumour_bams, if(binsize<30000) { print("Current minimum bin size is 30000 - resetting to 30000") + binsize <- 30000 } if(is.null(res)) res <- list() @@ -140,15 +149,15 @@ run_sc_sequencing <- function(tumour_bams, } if(is.null(barcodes_10x)) { - if(!is.null(normal_bams[1]) & is.null(res$nlCTS.normal)) + if(!is.null(normal_bams[1]) && is.null(res$nlCTS.normal)) { print("## get all tracks from normal bams") timetoread_normals <- system.time(res$lCTS.normal <- mclapply(normal_bams,function(bamfile) { lCTS.normal <- lapply(paste0(chrstring_bam,allchr), function(chr) getCoverageTrack(bamPath=bamfile, chr=chr, - lSe[[chr]]$starts, - lSe[[chr]]$ends, + res$lSe[[chr]]$starts, + res$lSe[[chr]]$ends, mapqFilter=30)) list(lCTS.normal=lCTS.normal, nlCTS.normal=treatTrack(lCTS=lCTS.normal, @@ -165,7 +174,7 @@ run_sc_sequencing <- function(tumour_bams, res$lNormals <- NULL res$timetoread_normals <- NULL } - if(any(names(res)=="allTracks") & res$binsize!=binsize) + if(any(names(res)=="allTracks") && res$binsize!=binsize) { print("## adjust Tracks for bin size ") res$timetoread_tumours <- system.time(res$allTracks <- mclapply(names(res$allTracks),function(bamfile) @@ -217,7 +226,7 @@ run_sc_sequencing <- function(tumour_bams, } else { - if(!any("allTracks.processed"%in%names(res)) | res$binsize!=binsize) + if(!any("allTracks.processed"%in%names(res)) || res$binsize!=binsize) { print("## smooth Tracks") res$timetoprocessed <- system.time(res$allTracks.processed <- mclapply(1:length(res$allTracks), function(x) @@ -253,17 +262,18 @@ run_sc_sequencing <- function(tumour_bams, purs = purs[[x]], ploidies = ploidies[[x]], maxTumourPhi=maxtumourpsi, - ismale=if(sex[x]=="male") T else F, - isPON=res$isPON),silent=F) + ismale=if(sex[x]=="male") TRUE else FALSE, + isPON=res$isPON),silent=FALSE) },mc.cores=MC.CORES)) print("## get Fitted CNA Profiles") res$allProfiles <- mclapply(1:length(res$allTracks.processed), function(x) { + if(inherits(res$allSols[[x]], "try-error")) return(NULL) try(getProfile(fitProfile(res$allTracks.processed[[x]], purity=res$allSols[[x]]$purity, ploidy=res$allSols[[x]]$ploidy, - ismale=if(sex[x]=="male") T else F), - CHRS=allchr),silent=F) + ismale=if(sex[x]=="male") TRUE else FALSE), + CHRS=allchr),silent=FALSE) },mc.cores=MC.CORES) names(res$allProfiles) <- names(res$allSols) <- names(res$allTracks) print("## compile Results") @@ -308,7 +318,7 @@ run_sc_sequencing <- function(tumour_bams, outdir=outdir, projectname=projectname) } - if(!is.null(list_ac_counts_paths) & !is.null(path_to_phases)) + if(!is.null(list_ac_counts_paths) && !is.null(path_to_phases)) { print("## get Allele-specific CNA") res <- getAS_CNA(res, @@ -323,13 +333,13 @@ run_sc_sequencing <- function(tumour_bams, mc.cores=MC.CORES) names(res$allProfiles_AS) <- names(res$allProfiles) } - if(smooth_sc & any(grepl("_AS",names(res)))) + if(smooth_sc && any(grepl("_AS",names(res)))) { print("## smooth Across Single-Cells/Nulcei") res <- getAS_CNA_smoothed(res, mc.cores=MC.CORES) } - if(smooth_sc & !any(grepl("_AS",names(res)))) + if(smooth_sc && !any(grepl("_AS",names(res)))) print("Warning: Smoothing is only possible for allele-specific copy numbers") res } From 3ff6a0a0adca7f1a1a2c6c7b0749ec3361f79928 Mon Sep 17 00:00:00 2001 From: galder-max Date: Thu, 30 Jul 2026 11:24:05 +0200 Subject: [PATCH 2/2] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- R/run_sc_sequencing.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/run_sc_sequencing.R b/R/run_sc_sequencing.R index e34c3b4..15ca07f 100644 --- a/R/run_sc_sequencing.R +++ b/R/run_sc_sequencing.R @@ -35,9 +35,9 @@ run_sc_sequencing <- function(tumour_bams, checkArguments_scs(c(as.list(environment()))) # --- Strict input validation guards --- - if(is.null(tumour_bams)) stop("tumour_bams cannot be NULL.") + if(is.null(tumour_bams) || length(tumour_bams) == 0) stop("tumour_bams must contain at least one BAM path.") if(!all(file.exists(tumour_bams))) stop("One or more files in tumour_bams do not exist.") - if(!is.null(normal_bams) && !all(file.exists(normal_bams))) stop("One or more files in normal_bams do not exist.") + if(!is.null(normal_bams) && length(normal_bams) > 0 && !all(file.exists(normal_bams))) stop("One or more files in normal_bams do not exist.") build <- match.arg(build, choices = c("hg19", "hg38", "mm39")) # --------------------------------------