Skip to content

Add input validation (min binsize 30k) and skip cells with missing VCFs #35

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
tdalil wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Add input validation (min binsize 30k) and skip cells with missing VCFs #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1b75298
Add input validation (min binsize 30k) and skip cells with missing VCFs
tdalil Jun 22, 2026
3ff6a0a
Potential fix for pull request finding
galder-max Jul 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 48 additions & 15 deletions R/getAS_CNA.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
Expand Down Expand Up @@ -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)
Expand All @@ -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({
Expand Down
36 changes: 23 additions & 13 deletions R/run_sc_sequencing.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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) || 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) && 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"))

@galder-max galder-max Jul 30, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

build different from hg19, hg38 and mm39 should be accepted

# --------------------------------------

suppressPackageStartupMessages(require(parallel))
suppressPackageStartupMessages(require(Rsamtools))
suppressPackageStartupMessages(require(Biostrings))
Expand All @@ -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()
Expand Down Expand Up @@ -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,
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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,
Expand All @@ -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
}