diff --git a/DESCRIPTION b/DESCRIPTION index a580015..d62c2a2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: ichorCNA Type: Package Title: A tool for estimating the fraction of tumor in ultra-low-pass whole genome sequencing (ULP-WGS) of cell-free DNA -Version: 0.3.4 -Date: 2020-08-17 +Version: 0.6.1 +Date: 2025-12-04 Author: Gavin Ha, Maintainer: Gavin Ha Depends: R (>= 3.6.0) diff --git a/scripts/runIchorCNA.R b/scripts/runIchorCNA.R old mode 100644 new mode 100755 index bf1be7e..107be6f --- a/scripts/runIchorCNA.R +++ b/scripts/runIchorCNA.R @@ -58,11 +58,50 @@ option_list <- list( make_option(c("--txnStrength"), type="numeric", default=1e7, help = "Transition pseudo-counts. Exponent should be the same as the number of decimal places of --txnE. Default: [%default]"), make_option(c("--multSampleTxnStrength"), type="numeric", default=1, help="Strength of same state transition between multiple samples. Default: [%default]"), make_option(c("--plotFileType"), type="character", default="pdf", help = "File format for output plots. Default: [%default]"), - make_option(c("--plotYLim"), type="character", default="c(-2,2)", help = "ylim to use for chromosome plots. Default: [%default]"), + make_option(c("--plotYLim"), type="character", default="c(-2,2)", help = "ylim to use for chromosome plots. Default: [%default]"), make_option(c("--outDir"), type="character", default="./", help = "Output Directory. Default: [%default]"), make_option(c("--libdir"), type = "character", default=NULL, help = "Script library path. Usually exclude this argument unless custom modifications have been made to the ichorCNA R package code and the user would like to source those R files. Default: [%default]"), make_option(c("--cores"), type="numeric", default = 1, help = "Number of cores to use for EM. Default: [%default]") ) + +# Function to parse character vectors +parse_numeric_vector <- function(str_value) { + if (is.null(str_value) || str_value == "NULL" || str_value == "None") { + return(NULL) + } + # Remove "c(" and ")" if present + str_value <- gsub("^c\\(", "", str_value) + str_value <- gsub("\\)$", "", str_value) + # Handle ranges like "1:22" or "c(1:22,\"X\")" + if (grepl(":", str_value)) { + parts <- strsplit(str_value, ",")[[1]] + result <- c() + for (part in parts) { + part <- trimws(part) + if (grepl(":", part)) { + range_parts <- as.numeric(strsplit(part, ":")[[1]]) + result <- c(result, seq(range_parts[1], range_parts[2])) + } else { + # Remove quotes if present + part <- gsub("\"", "", part) + result <- c(result, part) + } + } + return(result) + } + # Split by comma and clean up + values <- strsplit(str_value, ",")[[1]] + values <- trimws(values) + values <- gsub("\"", "", values) + # Convert to numeric, otherwise keep as character + numeric_values <- suppressWarnings(as.numeric(values)) + if (all(!is.na(numeric_values))) { + return(numeric_values) + } else { + return(values) + } +} + parseobj <- OptionParser(option_list=option_list) opt <- parse_args(parseobj) print(opt) @@ -88,12 +127,20 @@ exons.bed <- opt$exons.bed # "0" if none specified centromere <- opt$centromere minMapScore <- opt$minMapScore flankLength <- opt$rmCentromereFlankLength -normal <- eval(parse(text = opt$normal)) -normal.init <- eval(parse(text = opt$normal.init)) -scStates <- eval(parse(text = opt$scStates)) + +# Parse vector parameters without eval +normal <- parse_numeric_vector(opt$normal) +normal.init <- parse_numeric_vector(opt$normal.init) +scStates <- parse_numeric_vector(opt$scStates) +lambda <- parse_numeric_vector(opt$lambda) +ploidy <- parse_numeric_vector(opt$ploidy) +plotYLim <- parse_numeric_vector(opt$plotYLim) +chrs <- as.character(parse_numeric_vector(opt$chrs)) +chrTrain <- as.character(parse_numeric_vector(opt$chrTrain)) +chrNormalize <- as.character(parse_numeric_vector(opt$chrNormalize)) + subclone.penalty <- opt$scPenalty likModel <- opt$likModel -lambda <- eval(parse(text = opt$lambda)) lambdaScaleHyperParam <- opt$lambdaScaleHyperParam estimateNormal <- opt$estimateNormal estimatePloidy <- opt$estimatePloidy @@ -102,7 +149,6 @@ maxFracCNASubclone <- opt$maxFracCNASubclone maxFracGenomeSubclone <- opt$maxFracGenomeSubclone minSegmentBins <- opt$minSegmentBins altFracThreshold <- opt$altFracThreshold -ploidy <- eval(parse(text = opt$ploidy)) coverage <- opt$coverage maxCN <- opt$maxCN txnE <- opt$txnE @@ -116,13 +162,10 @@ normal2IgnoreSC <- opt$normal2IgnoreSC outDir <- opt$outDir libdir <- opt$libdir plotFileType <- opt$plotFileType -plotYLim <- eval(parse(text=opt$plotYLim)) outImage <- paste0(outDir,"/", patientID,".RData") genomeBuild <- opt$genomeBuild genomeStyle <- opt$genomeStyle -chrs <- as.character(eval(parse(text = opt$chrs))) -chrTrain <- as.character(eval(parse(text=opt$chrTrain))); -chrNormalize <- as.character(eval(parse(text=opt$chrNormalize))); + seqlevelsStyle(chrs) <- genomeStyle seqlevelsStyle(chrNormalize) <- genomeStyle seqlevelsStyle(chrTrain) <- genomeStyle @@ -130,13 +173,13 @@ cores <- opt$cores ## load ichorCNA library or source R scripts if (!is.null(libdir) && libdir != "None"){ - source(paste0(libdir,"/R/utils.R")) - source(paste0(libdir,"/R/segmentation.R")) - source(paste0(libdir,"/R/EM.R")) - source(paste0(libdir,"/R/output.R")) - source(paste0(libdir,"/R/plotting.R")) + source(paste0(libdir,"/R/utils.R")) + source(paste0(libdir,"/R/segmentation.R")) + source(paste0(libdir,"/R/EM.R")) + source(paste0(libdir,"/R/output.R")) + source(paste0(libdir,"/R/plotting.R")) } else { - library(ichorCNA) + library(ichorCNA) } ## load seqinfo @@ -158,12 +201,12 @@ if (is.null(exons.bed) || exons.bed == "None" || exons.bed == "NULL"){ ## load PoN if (is.null(normal_panel) || normal_panel == "None" || normal_panel == "NULL"){ - normal_panel <- NULL + normal_panel <- NULL } if (is.null(centromere) || centromere == "None" || centromere == "NULL"){ # no centromere file provided - centromere <- system.file("extdata", "GRCh37.p13_centromere_UCSC-gapTable.txt", - package = "ichorCNA") + centromere <- system.file("extdata", "GRCh37.p13_centromere_UCSC-gapTable.txt", + package = "ichorCNA") } centromere <- read.delim(centromere,header=T,stringsAsFactors=F,sep="\t") save.image(outImage) @@ -172,7 +215,7 @@ save.image(outImage) message("Reading GC and mappability files") gc <- wigToGRanges(gcWig) if (is.null(gc)){ - stop("GC wig file not provided but is required") + stop("GC wig file not provided but is required") } map <- wigToGRanges(mapWig) if (is.null(map)){ @@ -201,10 +244,10 @@ for (i in 1:numSamples) { tumour_reads <- wigToGRanges(wigFiles[i,2]) message("Correcting Tumour") counts[[id]] <- loadReadCountsFromWig(tumour_reads, chrs = chrs, gc = gc, map = map, repTime = repTime, - centromere = centromere, flankLength = flankLength, - targetedSequences = targetedSequences, chrXMedianForMale = chrXMedianForMale, - genomeStyle = genomeStyle, fracReadsInChrYForMale = fracReadsInChrYForMale, - chrNormalize = chrNormalize, mapScoreThres = minMapScore) + centromere = centromere, flankLength = flankLength, + targetedSequences = targetedSequences, chrXMedianForMale = chrXMedianForMale, + genomeStyle = genomeStyle, fracReadsInChrYForMale = fracReadsInChrYForMale, + chrNormalize = chrNormalize, mapScoreThres = minMapScore) gender <- counts[[id]]$gender if ((!is.null(sex) && sex != "None" && sex != "NULL") && gender$gender != sex ){ #compare with user-defined sex @@ -214,48 +257,48 @@ for (i in 1:numSamples) { ## load in normal file if provided if (!is.null(normal_file) && normal_file != "None" && normal_file != "NULL"){ - message("Loading normal file:", normal_file) - normal_reads <- wigToGRanges(normal_file) - message("Correcting Normal") - counts.normal <- loadReadCountsFromWig(normal_reads, chrs=chrs, gc=gc, map=map, repTime = repTime, - centromere=centromere, flankLength = flankLength, targetedSequences=targetedSequences, - genomeStyle = genomeStyle, chrNormalize = chrNormalize, mapScoreThres = minMapScore) - normal_copy <- counts.normal$counts #as(counts$counts, "GRanges") - counts[[id]]$counts$cor.gc.normal <- counts.normal$counts$cor.gc - counts[[id]]$counts$cor.map.normal <- counts.normal$counts$cor.map - counts[[id]]$counts$cor.rep.normal <- counts.normal$counts$cor.rep - gender.normal <- counts[[id]]$gender + message("Loading normal file:", normal_file) + normal_reads <- wigToGRanges(normal_file) + message("Correcting Normal") + counts.normal <- loadReadCountsFromWig(normal_reads, chrs=chrs, gc=gc, map=map, repTime = repTime, + centromere=centromere, flankLength = flankLength, targetedSequences=targetedSequences, + genomeStyle = genomeStyle, chrNormalize = chrNormalize, mapScoreThres = minMapScore) + normal_copy <- counts.normal$counts #as(counts$counts, "GRanges") + counts[[id]]$counts$cor.gc.normal <- counts.normal$counts$cor.gc + counts[[id]]$counts$cor.map.normal <- counts.normal$counts$cor.map + counts[[id]]$counts$cor.rep.normal <- counts.normal$counts$cor.rep + gender.normal <- counts[[id]]$gender }else{ - normal_copy <- NULL + normal_copy <- NULL } ### DETERMINE GENDER ### ## if normal file not given, use chrY, else use chrX message("Determining gender...", appendLF = FALSE) gender.mismatch <- FALSE if (!is.null(normal_copy)){ - if (gender$gender != gender.normal$gender){ #use tumour # use normal if given - # check if normal is same gender as tumour - gender.mismatch <- TRUE - } + if (gender$gender != gender.normal$gender){ #use tumour # use normal if given + # check if normal is same gender as tumour + gender.mismatch <- TRUE + } } message("Gender ", gender$gender) - + ## NORMALIZE GENOME-WIDE BY MATCHED NORMAL OR NORMAL PANEL (MEDIAN) ## tumour_copy[[id]] <- counts[[id]]$counts tumour_copy[[id]] <- normalizeByPanelOrMatchedNormal(tumour_copy[[id]], chrs = chrs, - normal_panel = normal_panel, normal_copy = normal_copy, - gender = gender$gender, normalizeMaleX = normalizeMaleX) - - ### OUTPUT FILE ### - ### PUTTING TOGETHER THE COLUMNS IN THE OUTPUT ### - outMat <- as.data.frame(tumour_copy[[id]]) - #outMat <- outMat[,c(1,2,3,12)] - outMat <- outMat[,c("seqnames","start","end","copy")] - colnames(outMat) <- c("chr","start","end","log2_TNratio_corrected") - outFile <- paste0(outDir,"/",id,".correctedDepth.txt") - message(paste("Outputting to:", outFile)) - write.table(outMat, file=outFile, row.names= FALSE, col.names = TRUE, quote = FALSE, sep="\t") - + normal_panel = normal_panel, normal_copy = normal_copy, + gender = gender$gender, normalizeMaleX = normalizeMaleX) + + ### OUTPUT FILE ### + ### PUTTING TOGETHER THE COLUMNS IN THE OUTPUT ### + outMat <- as.data.frame(tumour_copy[[id]]) + #outMat <- outMat[,c(1,2,3,12)] + outMat <- outMat[,c("seqnames","start","end","copy")] + colnames(outMat) <- c("chr","start","end","log2_TNratio_corrected") + outFile <- paste0(outDir,"/",id,".correctedDepth.txt") + message(paste("Outputting to:", outFile)) + write.table(outMat, file=outFile, row.names= FALSE, col.names = TRUE, quote = FALSE, sep="\t") + } ## end of for each sample chrInd <- as.character(seqnames(tumour_copy[[1]])) %in% chrTrain @@ -279,8 +322,8 @@ ptmTotalSolutions <- proc.time() # start total timer results <- list() numCombinations <- (length(normal) * length(ploidy)) ^ S loglik <- as.data.frame(matrix(NA, nrow = numCombinations, ncol = 7, - dimnames = list(c(), c("n_0", "phi_0", "n_est", "phi_est", - "Frac_genome_subclonal", "Frac_CNA_subclonal", "loglik")))) + dimnames = list(c(), c("n_0", "phi_0", "n_est", "phi_est", + "Frac_genome_subclonal", "Frac_CNA_subclonal", "loglik")))) fracGenomeSub <- as.data.frame(matrix(NA, nrow = numCombinations, ncol = S)) fracAltSub <- as.data.frame(matrix(NA, nrow = numCombinations, ncol = S)) # prepare normal/ploidy restarts for different solutions @@ -313,10 +356,10 @@ for (i in 1:length(ploidy)){ logR <- as.data.frame(lapply(tumour_copy, function(x) { x$copy })) # NEED TO EXCLUDE CHR X # param <- getDefaultParameters(logR[valid & chrInd, , drop=F], n_0 = n, maxCN = maxCN, - includeHOMD = includeHOMD, - ct.sc=scStates, normal2IgnoreSC = normal2IgnoreSC, ploidy_0 = floor(p), - e=txnE, e.subclone = subclone.penalty, e.sameState = multSampleTxnStrength, - strength=txnStrength, likModel = likModel) + includeHOMD = includeHOMD, + ct.sc=scStates, normal2IgnoreSC = normal2IgnoreSC, ploidy_0 = floor(p), + e=txnE, e.subclone = subclone.penalty, e.sameState = multSampleTxnStrength, + strength=txnStrength, likModel = likModel) ############################################ ######## CUSTOM PARAMETER SETTINGS ######### @@ -330,59 +373,59 @@ for (i in 1:length(ploidy)){ # param$betaLambda[param$ct.sc.status] <- param$betaLambda[param$ct.sc.status] /2 # param$alphaVar[param$ct.sc.status, 1] <- param$alphaVar[param$ct.sc.status, 1] *2 # } - ############################################# - ################ RUN HMM #################### - ############################################# + ############################################# + ################ RUN HMM #################### + ############################################# hmmResults.cor <- HMMsegment(tumour_copy, valid, dataType = "copy", param = param, chrTrain = chrTrain, maxiter = 20, estimateNormal = estimateNormal, estimatePloidy = estimatePloidy, estimatePrecision = TRUE, estimateVar = TRUE, estimateSubclone = estimateScPrevalence, estimateTransition = TRUE, estimateInitDist = TRUE, likChangeConvergence = 1e-4, verbose = TRUE) - + for (s in 1:numSamples){ - iter <- hmmResults.cor$results$iter - id <- names(hmmResults.cor$cna)[s] - + iter <- hmmResults.cor$results$iter + id <- names(hmmResults.cor$cna)[s] + ## correct integer copy number based on estimated purity and ploidy - correctedResults <- correctIntegerCN(cn = hmmResults.cor$cna[[s]], - segs = hmmResults.cor$results$segs[[s]], - purity = 1 - hmmResults.cor$results$n[s, iter], ploidy = hmmResults.cor$results$phi[s, iter], - cellPrev = 1 - hmmResults.cor$results$sp[s, iter], - maxCNtoCorrect.autosomes = maxCN, maxCNtoCorrect.X = maxCN, minPurityToCorrect = 0.05, - gender = gender$gender, chrs = chrs, correctHOMD = includeHOMD) - hmmResults.cor$results$segs[[s]] <- correctedResults$segs - hmmResults.cor$cna[[s]] <- correctedResults$cn - ## convert full diploid solution (of chrs to train) to have 1.0 normal or 0.0 purity - ## check if there is an altered segment that has at least a minimum # of bins - segsS <- hmmResults.cor$results$segs[[s]] - segsS <- segsS[segsS$chr %in% chrTrain, ] - segAltInd <- which(segsS$event != "NEUT") - maxBinLength = -Inf - if (length(segAltInd) > 0){ - maxInd <- which.max(segsS$end[segAltInd] - segsS$start[segAltInd] + 1) - maxSegRD <- GRanges(seqnames=segsS$chr[segAltInd[maxInd]], - ranges=IRanges(start=segsS$start[segAltInd[maxInd]], end=segsS$end[segAltInd[maxInd]])) - hits <- findOverlaps(query=maxSegRD, subject=tumour_copy[[s]][valid, ]) - maxBinLength <- length(subjectHits(hits)) - } - ## check if there are proportion of total bins altered - # if segment size smaller than minSegmentBins, but altFrac > altFracThreshold, then still estimate TF - cnaS <- hmmResults.cor$cna[[s]] - altInd <- cnaS[cnaS$chr %in% chrTrain, "event"] == "NEUT" - altFrac <- sum(!altInd, na.rm=TRUE) / length(altInd) - if ((maxBinLength <= minSegmentBins) & (altFrac <= altFracThreshold)){ - hmmResults.cor$results$n[s, iter] <- 1.0 - } + correctedResults <- correctIntegerCN(cn = hmmResults.cor$cna[[s]], + segs = hmmResults.cor$results$segs[[s]], + purity = 1 - hmmResults.cor$results$n[s, iter], ploidy = hmmResults.cor$results$phi[s, iter], + cellPrev = 1 - hmmResults.cor$results$sp[s, iter], + maxCNtoCorrect.autosomes = maxCN, maxCNtoCorrect.X = maxCN, minPurityToCorrect = 0.05, + gender = gender$gender, chrs = chrs, correctHOMD = includeHOMD) + hmmResults.cor$results$segs[[s]] <- correctedResults$segs + hmmResults.cor$cna[[s]] <- correctedResults$cn + ## convert full diploid solution (of chrs to train) to have 1.0 normal or 0.0 purity + ## check if there is an altered segment that has at least a minimum # of bins + segsS <- hmmResults.cor$results$segs[[s]] + segsS <- segsS[segsS$chr %in% chrTrain, ] + segAltInd <- which(segsS$event != "NEUT") + maxBinLength = -Inf + if (length(segAltInd) > 0){ + maxInd <- which.max(segsS$end[segAltInd] - segsS$start[segAltInd] + 1) + maxSegRD <- GRanges(seqnames=segsS$chr[segAltInd[maxInd]], + ranges=IRanges(start=segsS$start[segAltInd[maxInd]], end=segsS$end[segAltInd[maxInd]])) + hits <- findOverlaps(query=maxSegRD, subject=tumour_copy[[s]][valid, ]) + maxBinLength <- length(subjectHits(hits)) + } + ## check if there are proportion of total bins altered + # if segment size smaller than minSegmentBins, but altFrac > altFracThreshold, then still estimate TF + cnaS <- hmmResults.cor$cna[[s]] + altInd <- cnaS[cnaS$chr %in% chrTrain, "event"] == "NEUT" + altFrac <- sum(!altInd, na.rm=TRUE) / length(altInd) + if ((maxBinLength <= minSegmentBins) & (altFrac <= altFracThreshold)){ + hmmResults.cor$results$n[s, iter] <- 1.0 + } ## plot solution ## mainName[[s]] <- c(mainName[[s]], paste0("Solution: ", counter, ", Sample: ", id, ", n: ", n[s], ", p: ", p[s], - ", log likelihood: ", signif(hmmResults.cor$results$loglik[hmmResults.cor$results$iter], digits = 4))) + ", log likelihood: ", signif(hmmResults.cor$results$loglik[hmmResults.cor$results$iter], digits = 4))) ## plot individual samples if (numSamples == 1){ # if only single-sample analysis outPlotFile <- paste0(outDir, "/", id, "/", id, "_genomeWide_", "n", n[s], "-p", p[s]) plotGWSolution(hmmResults.cor, s=s, outPlotFile=outPlotFile, plotFileType=plotFileType, - logR.column = "logR", call.column = "Corrected_Call", - plotYLim=plotYLim, estimateScPrevalence=estimateScPrevalence, seqinfo=seqinfo, main=mainName[[s]][counter]) + logR.column = "logR", call.column = "Corrected_Call", + plotYLim=plotYLim, estimateScPrevalence=estimateScPrevalence, seqinfo=seqinfo, main=mainName[[s]][counter]) } } iter <- hmmResults.cor$results$iter @@ -416,22 +459,22 @@ save.image(outImage) fracAltSub[is.nan(fracAltSub$V1), ] <- 0 if (estimateScPrevalence){ ## sort but excluding solutions with too large % subclonal if (numSamples > 1){ - fracInd <- which(rowSums(fracAltSub <= maxFracCNASubclone) == S & - rowSums(fracGenomeSub <= maxFracGenomeSubclone) == S) + fracInd <- which(rowSums(fracAltSub <= maxFracCNASubclone) == S & + rowSums(fracGenomeSub <= maxFracGenomeSubclone) == S) }else{ - fracInd <- which(fracAltSub <= maxFracCNASubclone & - fracGenomeSub <= maxFracGenomeSubclone) + fracInd <- which(fracAltSub <= maxFracCNASubclone & + fracGenomeSub <= maxFracGenomeSubclone) } - if (length(fracInd) > 0){ ## if there is a solution satisfying % subclonal + if (length(fracInd) > 0){ ## if there is a solution satisfying % subclonal # sort by highest loglik and then by lowest fracCNAsubclonal (if ties) - ind <- fracInd[order(-loglik[fracInd, "loglik"], loglik[fracInd, "Frac_CNA_subclonal"])] + ind <- fracInd[order(-loglik[fracInd, "loglik"], loglik[fracInd, "Frac_CNA_subclonal"])] # sort by highest loglik for solutions that don't pass fracCNAsubclonal threshold fracInd.exclude <- setdiff(1:nrow(loglik), ind) ind.excludeSC <- fracInd.exclude[order(-loglik[fracInd.exclude, "loglik"])] ind <- c(ind, ind.excludeSC) - }else{ # otherwise just take largest likelihood - ind <- order(as.numeric(loglik[, "loglik"]), decreasing=TRUE) - } + }else{ # otherwise just take largest likelihood + ind <- order(as.numeric(loglik[, "loglik"]), decreasing=TRUE) + } }else{#sort by likelihood only ind <- order(as.numeric(loglik[, "loglik"]), decreasing=TRUE) } @@ -447,16 +490,16 @@ for (s in 1:numSamples){ turnDevOff <- FALSE turnDevOn <- FALSE if (i == 1){ - turnDevOn <- TRUE + turnDevOn <- TRUE } if (i == length(ind)){ - turnDevOff <- TRUE + turnDevOff <- TRUE } plotGWSolution(hmmResults.cor, s=s, outPlotFile=outPlotFile, plotFileType="pdf", - logR.column = "logR", call.column = "Corrected_Call", - plotYLim=plotYLim, estimateScPrevalence=estimateScPrevalence, - seqinfo = seqinfo, - turnDevOn = turnDevOn, turnDevOff = turnDevOff, main=mainName[[s]][ind[i]]) + logR.column = "logR", call.column = "Corrected_Call", + plotYLim=plotYLim, estimateScPrevalence=estimateScPrevalence, + seqinfo = seqinfo, + turnDevOn = turnDevOn, turnDevOff = turnDevOff, main=mainName[[s]][ind[i]]) } } # multisample, combined plots @@ -465,7 +508,7 @@ if (numSamples > 1){ if (plotFileType == "png"){ outPlotFile <- paste0(outPlotFile, ".png") png(outPlotFile,width=20,height=numSamples*6,units="in",res=300) - }else{ + }else{ outPlotFile <- paste0(outPlotFile, ".pdf") pdf(outPlotFile,width=20,height=numSamples*6) } @@ -474,11 +517,11 @@ if (numSamples > 1){ par(mfrow=c(numSamples, 1)) for (s in 1:numSamples){ plotGWSolution(hmmResults.cor, s=s, outPlotFile=outPlotFile, plotFileType="pdf", - logR.column = "logR", call.column = "Corrected_Call", - plotYLim=plotYLim, estimateScPrevalence=estimateScPrevalence, - seqinfo = seqinfo, spacing = 4, - cex.text = 1.25, cex=0.5, - turnDevOn = FALSE, turnDevOff = FALSE, main=mainName[[s]][ind[i]]) + logR.column = "logR", call.column = "Corrected_Call", + plotYLim=plotYLim, estimateScPrevalence=estimateScPrevalence, + seqinfo = seqinfo, spacing = 4, + cex.text = 1.25, cex=0.5, + turnDevOn = FALSE, turnDevOff = FALSE, main=mainName[[s]][ind[i]]) } } dev.off() @@ -493,7 +536,7 @@ hmmResults.cor$results$chrXMedian <- gender$chrXMedian hmmResults.cor$results$coverage <- coverage outputHMM(cna = hmmResults.cor$cna, segs = hmmResults.cor$results$segs, - results = hmmResults.cor$results, patientID = patientID, outDir=outDir) + results = hmmResults.cor$results, patientID = patientID, outDir=outDir) outFile <- paste0(outDir, "/", patientID, ".params.txt") outputParametersToFile(hmmResults.cor, file = outFile) @@ -501,5 +544,4 @@ outputParametersToFile(hmmResults.cor, file = outFile) plotSolutions(hmmResults.cor, tumour_copy, chrs, outDir, counts, numSamples=numSamples, logR.column = "logR", call.column = "event", likModel = likModel, plotFileType=plotFileType, plotYLim=plotYLim, seqinfo = seqinfo, - estimateScPrevalence=estimateScPrevalence, maxCN=maxCN) - + estimateScPrevalence=estimateScPrevalence, maxCN=maxCN) \ No newline at end of file diff --git a/scripts/snakemake/config/config.yaml b/scripts/snakemake/config/config.yaml index d19122d..02d2aaa 100644 --- a/scripts/snakemake/config/config.yaml +++ b/scripts/snakemake/config/config.yaml @@ -58,4 +58,4 @@ ichorCNA_txnE: 0.9999 # lower (e.g. 100) leads to higher sensitivity and more segments ichorCNA_txnStrength: 10000 ichorCNA_plotFileType: pdf -ichorCNA_plotYlim: c(-2,4) +ichorCNA_plotYlim: c(-2,4) \ No newline at end of file diff --git a/scripts/snakemake/config/samples.yaml b/scripts/snakemake/config/samples.yaml index c5f3d43..e7f40ef 100644 --- a/scripts/snakemake/config/samples.yaml +++ b/scripts/snakemake/config/samples.yaml @@ -1,2 +1,2 @@ samples: - tumor_sample1: /path/to/bam/tumor.bam + tumor_sample1: /path/to/bam/tumor.bam \ No newline at end of file