Error in rho[, 1] : subscript out of bounds

[Dries-B]

Hello,

Running the following code:

ichorCNArepo=~/FM_seq_001/ichorCNA/ichorCNA

Rscript ../runIchorCNA.R --id S21 \

        --WIG results/readDepth/S21.bin50000.wig \

        --ploidy "c(2,3)" \

        --normal "c(0.5,0.6,0.7,0.8,0.9)" \

        --maxCN 5 \

        --gcWig $ichorCNArepo/inst/extdata/gc_hg19_50kb.wig \

  --mapWig $ichorCNArepo/inst/extdata/map_hg19_50kb.wig \

  --centromere $ichorCNArepo/inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt \

  --normalPanel $ichorCNArepo/inst/extdata/HD_ULP_PoN_1Mb_median_normAutosome_mapScoreFiltered_median.rds \

  --includeHOMD False \

  --chrs "c(1:22, \"X\")" \

  --chrTrain "c(1:22)" \

  --estimateNormal True \

  --estimatePloidy True \

  --estimateScPrevalence True \

  --scStates "c(1,3)" \

  --txnE 0.9999 \

  --txnStrength 10000 \

  --outDir results/ichorCNA

Resulted in the following error:

Sorting by decreasing chromosome size
Correcting Tumour
Removed 1495 bins near centromeres.
Applying filter on data...
Correcting for GC bias...
Correcting for mappability bias...
Filtering low uniqueness regions with mappability score < 0.75
Removed 65 bins near centromeres.
Determining gender...Gender
Outputting to: results/ichorCNA/S21/S21.correctedDepth.txt
Warning messages:
1: In dir.create(paste0(outDir, "/", id, "/"), recursive = TRUE) :
  'results/ichorCNA/S21/S21' already exists
2: In regularize.values(x, y, ties, missing(ties)) :
  collapsing to unique 'x' values
runEM: Initialization
runEM iter1: Expectation
runEM iter1: Maximization

Error in rho[, 1] : subscript out of bounds

Calls: HMMsegment -> runEM -> estimateParamsMap
No traceback available
Total ULP-WGS HMM Runtime: 0.012 min.
Error in results[[ind[i]]] : subscript out of bounds
8: fn(par, ...)
7: (function (par)
   fn(par, ...))(0.5)
6: optim(n_prev, fn = completeLikelihoodFun, pType = rep("n", S),
       n = n_prev, sp = sp_prev, phi = phi_prev, lambda = lambda_prev,
       piG = pi_hat, A = A_hat, params = params, D = D, rho = rho,
       Zcounts = Zcounts, estimateNormal = estimateNormal, estimatePloidy = estimatePloidy,
       estimatePrecision = estimatePrecision, estimateInitDist = estimateInitDist,
       estimateTransition = estimateTransition, estimateSubclone = estimateSubclone,
       method = "L-BFGS-B", lower = intervalNormal[1], upper = intervalNormal[2],
       control = list(trace = 0, fnscale = -1))
5: withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning"))
4: suppressWarnings(estNorm <- optim(n_prev, fn = completeLikelihoodFun,
       pType = rep("n", S), n = n_prev, sp = sp_prev, phi = phi_prev,
       lambda = lambda_prev, piG = pi_hat, A = A_hat, params = params,
       D = D, rho = rho, Zcounts = Zcounts, estimateNormal = estimateNormal,
       estimatePloidy = estimatePloidy, estimatePrecision = estimatePrecision,
       estimateInitDist = estimateInitDist, estimateTransition = estimateTransition,
       estimateSubclone = estimateSubclone, method = "L-BFGS-B",
       lower = intervalNormal[1], upper = intervalNormal[2], control = list(trace = 0,
           fnscale = -1)))
3: estimateParamsMap(copy[chrTrain, ], n[, i - 1], sp[, i - 1],
       phi[, i - 1], lambdas[, , i - 1], piG[, i - 1], A, param,
       rho[, chrTrain], Zcounts, estimateNormal = estimateNormal,
       estimatePloidy = estimatePloidy, estimatePrecision = estimatePrecision,
       estimateTransition = estimateTransition, estimateInitDist = estimateInitDist,
       estimateSubclone = estimateSubclone)
2: runEM(dataMat, chr, chrInd, param, maxiter, verbose, estimateNormal = estimateNormal,
       estimatePloidy = estimatePloidy, estimateSubclone = estimateSubclone,
       estimatePrecision = estimatePrecision, estimateTransition = estimateTransition,
       estimateInitDist = estimateInitDist)
1: HMMsegment(tumour_copy, valid, dataType = "copy", param = param,
       chrTrain = chrTrain, maxiter = 50, estimateNormal = estimateNormal,
       estimatePloidy = estimatePloidy, estimateSubclone = estimateScPrevalence,
       verbose = TRUE)
Error in results[[ind[1]]] : subscript out of bounds
No traceback available
Error in hmmResults.cor$results$loglik <- as.data.frame(loglik) :
  object 'hmmResults.cor' not found
No traceback available
Error in hmmResults.cor$results$gender <- gender$gender :
  object 'hmmResults.cor' not found
No traceback available
Error in hmmResults.cor$results$chrYCov <- gender$chrYCovRatio :
  object 'hmmResults.cor' not found
No traceback available
Error in hmmResults.cor$results$chrXMedian <- gender$chrXMedian :
  object 'hmmResults.cor' not found
No traceback available
Error in hmmResults.cor$results$coverage <- coverage :
  object 'hmmResults.cor' not found
No traceback available
Error in outputHMM(cna = hmmResults.cor$cna, segs = hmmResults.cor$results$segs,  :
  object 'hmmResults.cor' not found
No traceback available
Error in outputParametersToFile(hmmResults.cor, file = outFile) :
  object 'hmmResults.cor' not found
1: outputHMM(cna = hmmResults.cor$cna, segs = hmmResults.cor$results$segs,
       results = hmmResults.cor$results, patientID = patientID,
       outDir = outDir)
Error in plotSolutions(hmmResults.cor, tumour_copy, chrs, outDir, numSamples = numSamples,  :
  object 'hmmResults.cor' not found
1: outputParametersToFile(hmmResults.cor, file = outFile)

This excludes earlier slurping and parsing, and includes traceback() information.

The same error occurs for 30 different samples when run through Snakemake. It does not occur for the test sample MBC_315. The latter shows that the software seems to work correctly.

Do you know what could be causing this? Just like in #39, it occurs during optimization, but in all 30 samples.

[Dries-B]

Some updates on my trouble-shooting:

1. The WIG file S21.bin50000.wig looks normal compared to the test sample MBC_315.

2. In runIchorCNA_debug.log I share extra debugging information, such as the variables at time of error.

[gavinha]

Hi @Dries-B

Are you able to share one of the .RData files so that I can debug it for you? If so, please send me an email at gha@fredhutch.org

Thanks, Gavin

[Dries-B]

Hi Gavin,

Thank you for your reply! I shared a file with you, which should enable you to replicate the error.

By the way, one of my colleagues suggested that the error might arise because of some bins having log(2) ratios which are outliers. What do you think about that?

Cheers, Dries

[gavinha]

Hi @Dries-B

Please try running again but removing this argument --normalPanel $ichorCNArepo/inst/extdata/HD_ULP_PoN_1Mb_median_normAutosome_mapScoreFiltered_median.rds

Gavin

[Dries-B]

Hi Gavin,

Indeed, I had overlooked the 1Mb in that filename. The script now ran without errors.

Thank you for your help!