### Running ichorCNA R Script

Source: https://rdrr.io/github/broadinstitute/ichorCNA/f/scripts/README.md

Example of how to launch the ichorCNA R script from the command line with various parameters.

```APIDOC
## POST /runIchorCNA

### Description
Executes the ichorCNA R script to perform copy number alteration analysis on WIG files.

### Method
POST

### Endpoint
/runIchorCNA

### Parameters
#### Query Parameters
- **libdir** (string) - Required - Script library path
- **datadir** (string) - Required - Reference wig dir path
- **id** (string) - Required - Patient ID
- **WIG** (string) - Required - Path to tumor WIG file.
- **NORMWIG** (string) - Optional - Path to normal WIG file.
- **normalPanel** (string) - Optional - Median corrected depth from panel of normals
- **exons.bed** (string) - Optional - Path to bed file containing exon regions.
- **normal** (string) - Optional - Initial normal contamination (e.g., "c(0.5,0.6,0.7,0.8,0.9)")
- **scStates** (string) - Optional - Subclonal states to consider (e.g., "c(1,3)")
- **ploidy** (string) - Optional - Initial tumour ploidy (e.g., "c(2,3)")
- **maxCN** (integer) - Optional - Total clonal CN states (e.g., 5)
- **estimateNormal** (boolean) - Optional - Estimate normal contamination (e.g., True)
- **estimateScPrevalence** (boolean) - Optional - Estimate subclonal prevalence (e.g., True)
- **estimatePloidy** (boolean) - Optional - Estimate tumour ploidy (e.g., True)
- **chrs** (string) - Optional - Chromosomes to consider (e.g., "c(1:22, \"X\")")
- **chrTrain** (string) - Optional - Chromosomes for training (e.g., "c(1:22)")
- **centromere** (string) - Required - Path to centromere file (e.g., "../../inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt")
- **txnE** (float) - Optional - Transactional E value threshold (e.g., 0.9999)
- **txnStrength** (integer) - Optional - Transactional strength (e.g., 10000)
- **plotFileType** (string) - Optional - Type of plot file to generate (e.g., "png")
- **plotYLim** (string) - Optional - Y-axis limits for plots (e.g., "c(-2,4)")
- **outDir** (string) - Required - Output directory for results

### Request Example
```json
{
  "libdir": "../../R/",
  "datadir": "../../inst/extdata/",
  "id": "tumor_sample1",
  "WIG": "results/readDepth/tumor_sample1.bin1000000.wig",
  "ploidy": "c(2,3)",
  "normal": "c(0.5,0.6,0.7,0.8,0.9)",
  "maxCN": 5,
  "includeHOMD": false,
  "chrs": "c(1:22, \"X\")",
  "chrTrain": "c(1:22)",
  "estimateNormal": true,
  "estimatePloidy": true,
  "estimateScPrevalence": true,
  "scStates": "c(1,3)",
  "centromere": "../../inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt",
  "exons.bed": "None",
  "txnE": 0.9999,
  "txnStrength": 10000,
  "plotFileType": "png",
  "plotYLim": "c(-2,4)",
  "outDir": "results/ichorCNA/tumor_sample1/"
}
```

### Response
#### Success Response (200)
- **status** (string) - Indicates the success of the operation.
- **message** (string) - A message detailing the outcome of the script execution.

#### Response Example
```json
{
  "status": "success",
  "message": "ichorCNA analysis completed successfully. Results saved to results/ichorCNA/tumor_sample1/"
}
```
```

--------------------------------

### Install ichorCNA Package

Source: https://rdrr.io/github/broadinstitute/ichorCNA

Install the latest version of the ichorCNA package from GitHub using the remotes package. Ensure remotes is installed first.

```r
install.packages("remotes")
remotes::install_github("broadinstitute/ichorCNA")
```

--------------------------------

### Get Default Parameters

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/segmentation.R

Initializes the default parameters for the HMM model, including copy number states and prior distributions.

```R
getDefaultParameters <- function(x, maxCN = 5, ct.sc = NULL, ploidy = 2, e = 0.9999999, e.sameState = 10, strength = 10000000, includeHOMD = FALSE){
  if (includeHOMD){
    ct <- 0:maxCN
  }else{
    ct <- 1:maxCN
  }
	param <- list(
		strength = strength, e = e,
		ct = c(ct, ct.sc),
		ct.sc.status = c(rep(FALSE, length(ct)), rep(TRUE, length(ct.sc))),
		phi_0 = 2, alphaPhi = 4, betaPhi = 1.5,
```

--------------------------------

### Initialize R Environment and Variables

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Loads required libraries and maps parsed command-line options to local R variables for downstream processing.

```R
library(HMMcopy)
library(GenomicRanges)
library(GenomeInfoDb)
options(stringsAsFactors=FALSE)
options(bitmapType='cairo')

patientID <- opt$id
tumour_file <- opt$WIG
normal_file <- opt$NORMWIG
gcWig <- opt$gcWig
mapWig <- opt$mapWig
normal_panel <- opt$normalPanel
exons.bed <- opt$exons.bed  # "0" if none specified
centromere <- opt$centromere
minMapScore <- opt$minMapScore
flankLength <- opt$rmCentromereFlankLength
normal <- eval(parse(text = opt$normal))
scStates <- eval(parse(text = opt$scStates))
lambda <- eval(parse(text = opt$lambda))
lambdaScaleHyperParam <- opt$lambdaScaleHyperParam
estimateNormal <- opt$estimateNormal
estimatePloidy <- opt$estimatePloidy
estimateScPrevalence <- opt$estimateScPrevalence
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
txnStrength <- opt$txnStrength
normalizeMaleX <- as.logical(opt$normalizeMaleX)
includeHOMD <- as.logical(opt$includeHOMD)
minTumFracToCorrect <- opt$minTumFracToCorrect
fracReadsInChrYForMale <- opt$fracReadsInChrYForMale
chrXMedianForMale <- -0.1
outDir <- opt$outDir
libdir <- opt$libdir
plotFileType <- opt$plotFileType
plotYLim <- eval(parse(text=opt$plotYLim))
gender <- NULL
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
```

--------------------------------

### Initialize Results for HMM Analysis

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Sets up a data frame to store results from different normal and ploidy solutions for the HMM analysis. Initializes timers and result storage.

```R
ptmTotalSolutions <- proc.time() # start total timer
results <- list()
loglik <- as.data.frame(matrix(NA, nrow = length(normal) * length(ploidy), ncol = 7, 
                 dimnames = list(c(), c("init", "n_est", "phi_est", "BIC", 
                 										"Frac_genome_subclonal", "Frac_CNA_subclonal", "loglik"))))
counter <- 1
compNames <- rep(NA, nrow(loglik))
mainName <- rep(NA, length(normal) * length(ploidy))
```

--------------------------------

### Get Genome-Wide Positions

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/plotting.R

Transforms chromosome-specific positions to genome-wide coordinates. It calculates chromosome breakpoints and adds them to the positions. Optionally uses seqinfo for accurate length calculations.

```R
getGenomeWidePositions <- function(chrs, posns, seqinfo = NULL) {
    # create genome coordinate scaffold
    positions <- as.numeric(posns)
    chrs <- as.character(chrs)
    chrsNum <- unique(chrs)
    chrBkpt <- rep(0, length(chrsNum) + 1)
    prevChrPos <- 0
    i <- 1
    if (length(chrsNum) > 1){
		for (i in 2:length(chrsNum)) {
			chrInd <- which(chrs == chrsNum[i])
			if (!is.null(seqinfo)){
				prevChrPos <- seqlengths(seqinfo)[i-1] + prevChrPos
			}else{
				prevChrPos <- positions[chrInd[1] - 1]
			}
			chrBkpt[i] = prevChrPos
			positions[chrInd] = positions[chrInd] + prevChrPos
		}
	}
    chrBkpt[i + 1] <- positions[length(positions)]
    return(list(posns = positions, chrBkpt = chrBkpt))
}
```

--------------------------------

### Define and Parse Command-Line Options

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Defines the available command-line arguments for the ichorCNA script and parses them into an object.

```R
  make_option(c("--fracReadsInChrYForMale"), type="numeric", default=0.001, help = "Threshold for fraction of reads in chrY to assign as male. Default: [%default]"),
  make_option(c("--includeHOMD"), type="logical", default=FALSE, help="If FALSE, then exclude HOMD state. Useful when using large bins (e.g. 1Mb). Default: [%default]"),
  make_option(c("--txnE"), type="numeric", default=0.9999999, help = "Self-transition probability. Increase to decrease number of segments. Default: [%default]"),
  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("--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("--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]")
)
parseobj <- OptionParser(option_list=option_list)
opt <- parse_args(parseobj)
print(opt)
options(scipen=0, stringsAsFactors=F)
```

--------------------------------

### Initialize ichorcna Parameters

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/segmentation.R

Initializes parameters for the ichorcna model, including hyperparameters for copy number states, precision, and transition matrices. Handles both single and multiple sample inputs.

```R
initialize_params <- function(x, param, strength, e, e.sameState, includeHOMD, K){
	# initialize hyperparameters for precision using observed data ##
	if (!is.null(dim(x))){
    param$numberSamples <- ncol(x)
    betaLambdaVal <- ((apply(x, 2, sd, na.rm = TRUE) / sqrt(length(param$ct))) ^ 2)   
	}else{
	  param$numberSamples <- 1
	  betaLambdaVal <- ((sd(x, na.rm = TRUE) / sqrt(length(param$ct))) ^ 2)
	}
	param$betaLambda <- matrix(betaLambdaVal, ncol = param$numberSamples, nrow = length(param$ct), byrow = TRUE)
  param$alphaLambda <- rep(param$alphaLambda, K)
  
	# increase prior precision for -1, 0, 1 copies at ploidy
	S <- param$numberSamples
	logR.var <- 1 / ((apply(x, 2, sd, na.rm = TRUE) / sqrt(length(param$ct))) ^ 2)
	if (!is.null(dim(x))){
		param$lambda <- matrix(logR.var, nrow=K, ncol=S, byrow=T, dimnames=list(c(),colnames(x)))
	}else{
    param$lambda <- matrix(logR.var, length(param$ct))
    param$lambda[param$ct %in% c(2)] <- logR.var 
    param$lambda[param$ct %in% c(1,3)] <- logR.var 
    param$lambda[param$ct >= 4] <- logR.var / 5
    param$lambda[param$ct == max(param$ct)] <- logR.var / 15
    param$lambda[param$ct.sc.status] <- logR.var / 10
  }
  # define joint copy number states #
  param$jointCNstates <- expand.grid(rep(list(param$ct), S))
  param$jointSCstatus <- expand.grid(rep(list(param$ct.sc.status), S))
  colnames(param$jointCNstates) <- paste0("Sample.", 1:param$numberSamples)
  colnames(param$jointSCstatus) <- paste0("Sample.", 1:param$numberSamples)
  
	# Initialize transition matrix to the prior
	xn <- getTransitionMatrix(K ^ S, e, strength)
  ## set higher transition probs for same CN states across samples ##
  cnStateDiff <- apply(param$jointCNstates, 1, function(x){ (abs(max(x) - min(x)))})
  if (e.sameState > 0 & S > 1){
		xn$A[, cnStateDiff == 0] <- txn$A[, cnStateDiff == 0] * e.sameState * K 
		xn$A[, cnStateDiff >= 3] <- txn$A[, cnStateDiff >=3]  / e.sameState / K
	}
  for (i in 1:nrow(txn$A)){
    for (j in 1:ncol(txn$A)){
      if (i == j){
        txn$A[i, j] <- e
      }
    }
  }
  txn$A <- normalize(txn$A)
	param$A <- txn$A
	param$dirPrior <- txn$A * strength[1] 
  param$A[, param$ct.sc.status] <- param$A[, param$ct.sc.status] / 10
  param$A <- normalize(param$A)
  param$dirPrior[, param$ct.sc.status] <- param$dirPrior[, param$ct.sc.status] / 10
  
  if (includeHOMD){
    K <- length(param$ct)
    param$A[1, 2:K] <- param$A[1, 2:K] * 1e-5; param$A[2:K, 1] <- param$A[2:K, 1] * 1e-5;
    param$A[1, 1] <- param$A[1, 1] * 1e-5
    param$A <- normalize(param$A); param$dirPrior <- param$A * param$strength
  }

  param$kappa <- rep(75, K ^ S)
  param$kappa[cnStateDiff == 0] <- param$kappa[cnStateDiff == 0] + 125
	param$kappa[cnStateDiff >=3] <- param$kappa[cnStateDiff >=3] - 50
	param$kappa[which(rowSums(param$jointCNstates==2) == S)] <- 800
  
  return(param)
}
```

--------------------------------

### Initialize Plotting Colors and Parameters

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/plotting.R

Sets up color palettes and alpha values for plotting copy number alterations, including options for subclonal data. Alpha values are converted to hexadecimal format.

```R
alphaVal <- ceiling(alphaVal * 255); class(alphaVal) = "hexmode"
  alphaSubcloneVal <- ceiling(alphaVal / 2 * 255); class(alphaVal) = "hexmode"
  cnCol <- c("#00FF00","#006400","#0000FF","#8B0000",rep("#FF0000", 26))
  subcloneCol <- c("#00FF00")
  cnCol <- paste(cnCol,alphaVal,sep="")
  names(cnCol) <- c("HOMD","HETD","NEUT","GAIN","AMP","HLAMP",paste0(rep("HLAMP", 8), 2:25))
#  segCol <- cnCol
#  ## add in colors for subclone if param provided
#  if (!is.null(param)){
#    ind <- ((which.max(param$ct) + 1) : length(param$ct)) + 1
#    cnCol[ind] <- paste0(cnCol[ind], alphaSubcloneVal / 2)
#    segCol[ind] <- "#00FF00"

```

--------------------------------

### Prepare Data for HMM and Save Environment

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Identifies valid chromosomal positions for training and saves the current R environment. This step is crucial before running the Hidden Markov Model (HMM).

```R
chrInd <- as.character(seqnames(tumour_copy[[1]])) %in% chrTrain
## get positions that are valid
valid <- tumour_copy[[1]]$valid
if (length(tumour_copy) >= 2) {
  for (i in 2:length(tumour_copy)){
    valid <- valid & tumour_copy[[i]]$valid 
  }
}
save.image(outImage)
```

--------------------------------

### ichorCNA R Script Help - Usage

Source: https://rdrr.io/github/broadinstitute/ichorCNA/f/scripts/README.md

Displays the usage information for the runIchorCNA.R script, outlining the available command-line options and their purpose.

```bash
>Rscript runIchorCNA.R --help

Usage: runIchorCNA.R [options] 


Options:
	-l LIBDIR, --libdir=LIBDIR
		Script library path

	--datadir=DATADIR
		Reference wig dir path

	-t WIG, --WIG=WIG
		Path to tumor WIG file.

	--NORMWIG=NORMWIG
		Path to normal WIG file.

	--normalPanel=NORMALPANEL
		Median corrected depth from panel of normals

	-e EXONS.BED, --exons.bed=EXONS.BED
		Path to bed file containing exon regions.

	--id=ID
		Patient ID.


	-n NORMAL, --normal=NORMAL
		Initial normal contamination

	--scStates=SCSTATES
		Subclonal states to consider

	-p PLOIDY, --ploidy=PLOIDY
		Initial tumour ploidy

	-m MAXCN, --maxCN=MAXCN
		Total clonal CN states

	--estimateNormal=ESTIMATENORMAL
		Estimate normal.

	--estimateScPrevalence=ESTIMATESCPREVALENCE
		Estimate subclonal prevalence.

	--estimatePloidy=ESTIMATEPLOIDY
		Estimate tumour ploidy.

 (plus more arguments)


```

--------------------------------

### ichorCNA Script Help

Source: https://rdrr.io/github/broadinstitute/ichorCNA/f/scripts/README.md

Displays the usage and available options for the ichorCNA R script.

```APIDOC
## GET /runIchorCNA/help

### Description
Retrieves the help information for the ichorCNA R script, listing all available command-line options.

### Method
GET

### Endpoint
/runIchorCNA/help

### Parameters
No parameters are required for this endpoint.

### Response
#### Success Response (200)
- **usage** (string) - The usage string for the script.
- **options** (object) - An object containing detailed descriptions of each command-line option.

#### Response Example
```json
{
  "usage": "Usage: runIchorCNA.R [options]",
  "options": {
    "-l, --libdir": {
      "description": "Script library path",
      "argument": "LIBDIR"
    },
    "--datadir": {
      "description": "Reference wig dir path",
      "argument": "DATADIR"
    },
    "-t, --WIG": {
      "description": "Path to tumor WIG file.",
      "argument": "WIG"
    },
    "--NORMWIG": {
      "description": "Path to normal WIG file.",
      "argument": "NORMWIG"
    },
    "--normalPanel": {
      "description": "Median corrected depth from panel of normals",
      "argument": "NORMALPANEL"
    },
    "-e, --exons.bed": {
      "description": "Path to bed file containing exon regions.",
      "argument": "EXONS.BED"
    },
    "--id": {
      "description": "Patient ID.",
      "argument": "ID"
    },
    "-n, --normal": {
      "description": "Initial normal contamination",
      "argument": "NORMAL"
    },
    "--scStates": {
      "description": "Subclonal states to consider",
      "argument": "SCSTATES"
    },
    "-p, --ploidy": {
      "description": "Initial tumour ploidy",
      "argument": "PLOIDY"
    },
    "-m, --maxCN": {
      "description": "Total clonal CN states",
      "argument": "MAXCN"
    },
    "--estimateNormal": {
      "description": "Estimate normal.",
      "argument": "ESTIMATENORMAL"
    },
    "--estimateScPrevalence": {
      "description": "Estimate subclonal prevalence.",
      "argument": "ESTIMATESCPREVALENCE"
    },
    "--estimatePloidy": {
      "description": "Estimate tumour ploidy.",
      "argument": "ESTIMATEPLOIDY"
    }
  }
}
```
```

--------------------------------

### Export HMM Parameters to File

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/output.R

Writes sample-specific HMM parameters such as tumor fraction, ploidy, and subclone fractions to a file.

```R
outputParametersToFile <- function(hmmResults, file){
  S <- hmmResults$results$param$numberSamples
  x <- hmmResults$results
  i <- x$iter
  fc <- file(file, "w+")
  outMat <- as.data.frame(cbind(`Tumor Fraction` = signif(1 - x$n[, i], digits = 4), Ploidy = signif(x$phi[, i], digits = 4)))
  outMat <- cbind(Sample = names(hmmResults$cna), outMat)
  rownames(outMat) <- names(hmmResults$cna)
  write.table(outMat, file = fc, col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")
  write.table("\n", file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
  for (s in 1:S){
    id <- names(hmmResults$cna)[s]
    write.table(id, file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
    write.table(paste0("Gender:\t", x$gender), file = fc, col.names = FALSE, 
                row.names = FALSE, quote = FALSE, sep = "\t")
    write.table(paste0("Tumor Fraction:\t", signif(1 - x$n[s, i], digits = 4)), file = fc, col.names = FALSE, 
                row.names = FALSE, quote = FALSE, sep = "\t")
    ploidy <- (1 - x$n[i]) * x$phi[i] + x$n[i] * 2
    write.table(paste0("Ploidy:\t", signif(x$phi[s, i], digits = 4)), file = fc, col.names = FALSE, 
                row.names = FALSE, quote = FALSE, sep = "\t")
   	subcloneGenomeFrac <- sum(hmmResults.cor$cna[[s]]$subclone.status) / nrow(hmmResults.cor$cna[[s]])
   	subcloneCNAFrac <- sum(hmmResults.cor$cna[[s]]$subclone.status) / sum(hmmResults.cor$cna[[s]]$copy.num != 2)
   	scFrac <- 1 - x$sp[s, i]
   	if (subcloneGenomeFrac == 0){
   		scFrac <- NA
   	}else{
   		scFrac <- signif(scFrac, digits = 4)
   	}
    #if (sum(x$param$ct.sc.status) != 0){    	
      write.table(paste0("Subclone Fraction:\t", scFrac), file = fc, col.names = FALSE, 
                  row.names = FALSE, quote = FALSE, sep = "\t")
      write.table(paste0("Fraction Genome Subclonal:\t", signif(subcloneGenomeFrac, digits = 4)), file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
      write.table(paste0("Fraction CNA Subclonal:\t", signif(subcloneCNAFrac, digits = 4)), file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
    #}
    if (is.null(x$coverage)){
    	coverage <- NA
    }else{
    	coverage <- signif(coverage, digits = 4)
    }
		write.table(paste0("Coverage:\t", coverage), file = fc, col.names = FALSE, 
								row.names = FALSE, quote = FALSE, sep = "\t")
	
    if (!is.null(x$chrYCov)){
```

--------------------------------

### Run ichorCNA R Script from Command Line

Source: https://rdrr.io/github/broadinstitute/ichorCNA/f/scripts/README.md

Execute the ichorCNA R script with specified WIG files and analysis parameters. Ensure all paths and arguments are correctly set for your data.

```bash
Rscript ../runIchorCNA.R --libdir ../../R/ --datadir ../../inst/extdata/ \
  --id tumor_sample1 --WIG results/readDepth/tumor_sample1.bin1000000.wig \
  --ploidy "c(2,3)" --normal "c(0.5,0.6,0.7,0.8,0.9)" --maxCN 5 \
  --includeHOMD False --chrs "c(1:22, \"X\")" --chrTrain "c(1:22)" \
  --estimateNormal True --estimatePloidy True --estimateScPrevalence True \
  --scStates "c(1,3)" --centromere ../../inst/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt \
  --exons.bed None --txnE 0.9999 --txnStrength 10000 --plotFileType png --plotYLim "c(-2,4)" \
  --outDir results/ichorCNA/tumor_sample1/ > logs/ichorCNA/tumor_sample1.log 2> logs/ichorCNA/tumor_sample1.log

```

--------------------------------

### Estimate Model Parameters using MAP

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/EM.R

Estimates model parameters (n, sp, phi, lambda, pi, A) using Maximum A Posteriori (MAP) estimation. Allows selective estimation of parameters based on boolean flags.

```R
estimateParamsMap <- function(D, n_prev, sp_prev, phi_prev, lambda_prev, pi_prev, A_prev, 
                              params, rho, Zcounts, 
                              estimateNormal = TRUE, estimatePloidy = TRUE, 
                              estimatePrecision = TRUE, estimateInitDist = TRUE, 
                              estimateTransition = TRUE, estimateSubclone = TRUE, 
                              verbose = TRUE) {
  KS <- nrow(rho)
  K <- length(params$ct)
  T <- ncol(rho)
  S <- length(n_prev)
  #dimen <- ncol(D)
  intervalNormal <- c(1e-6, 1 - 1e-6)
  intervalSubclone <- c(1e-6, 1 - 1e-6)
  intervalPhi <- c(.Machine$double.eps, 10)
  intervalLambda <- c(1e-5, 1e4)
  
  # initialize params to be estimated
  n_hat <- n_prev
  sp_hat <- sp_prev
  phi_hat <- phi_prev
  lambda_hat <- lambda_prev
  pi_hat <- pi_prev
  A_hat <- A_prev
  
  # Update transition matrix A
  if (estimateTransition){
    for (k in 1:KS) {
      A_hat[k, ] <- Zcounts[k, ] + params$dirPrior[k, ]
      A_hat[k, ] <- normalize(A_hat[k, ])
    }
  }
  # map estimate for pi (initial state dist)
  if (estimateInitDist){
    pi_hat <- estimateMixWeightsParamMap(rho, params$kappa)
  }
  
  # map estimate for normal 
  if (estimateNormal){
    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))
    )
    n_hat <- estNorm$par
  }
  if (estimateSubclone){
    suppressWarnings(
      estSubclone <- optim(sp_prev, fn = completeLikelihoodFun, pType = rep("sp", S), n = n_hat, 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))
    )
    sp_hat <- estSubclone$par
  }
  if (estimatePloidy){
    suppressWarnings(
      estPhi <- optim(phi_prev, fn = completeLikelihoodFun, pType = rep("phi", length(phi_prev)), n = n_hat, sp = sp_hat, 
                      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",

```

--------------------------------

### getSeqInfo

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/utils.R

Retrieves sequence information for a specified genome build.

```APIDOC
## getSeqInfo

### Description
Retrieves sequence information (Seqinfo) for a given genome build and style.

### Parameters
#### Arguments
- **genomeBuild** (string) - Optional - The genome build version (e.g., "hg19").
- **genomeStyle** (string) - Optional - The sequence naming style (e.g., "NCBI").
```

--------------------------------

### Measure and Report Runtime

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Calculates and prints the total elapsed time for HMM solution processing.

```R
## get total time for all solutions ##
elapsedTimeSolutions <- proc.time() - ptmTotalSolutions
message("Total ULP-WGS HMM Runtime: ", format(elapsedTimeSolutions[3] / 60, digits = 2), " min.")
```

--------------------------------

### Estimate Model Parameters via Optimization

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/EM.R

Uses the optim function to estimate ploidy (phi) and precision (lambda) parameters based on provided likelihood functions.

```R
                      lower = intervalPhi[1], upper = intervalPhi[2],
                      control = list(trace = 0, fnscale = -1))
    )
    phi_hat <- estPhi$par
  }
  if (estimatePrecision){
    suppressWarnings(
      estLambda <- optim(c(lambda_prev), fn = completeLikelihoodFun, pType = rep("lambda", K*S), n = n_hat, sp = sp_hat,
                         phi = phi_hat, 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 = intervalLambda[1],
                         control = list(trace = 0, fnscale = -1))
    )
    lambda_hat <- matrix(estLambda$par, ncol = S, byrow = FALSE)
  }
  
  #}
  #}
  
  #map estimate for phi (ploidy) 
  # if (estimatePloidy) {
  #    suppressWarnings(
  #   estPhi <- optim(phi_prev, fn = phiLikelihoodFun, n = n_prev, lambda = lambda_prev, 
  #								params = params, D = D, rho = rhoS, method = "L-BFGS-B",
  #                 control = list(fnscale = -1, ndeps = 1e-5), lower = intervalPhi[1])#, upper = intervalPhi[2])
  #  )
  #  phi_hat <- estPhi$par
  #}
  
  # map estimate for lambda (Student's t precision)
  # if (estimatePrecision){
  #for (s in 1:S){
  #  paramsS <- params
  #  paramsS$betaLambda <- paramsS$betaLambda[s]
  #  for (k in 1:K){
  #		  suppressWarnings(
  #	  estLambda <- optim(lambda_prev[k, s], fn = lambdaLikelihoodFun, n = n_prev[s], phi = phi_prev, 
  #										params = paramsS, D = D[, s], rho = rhoS[[s]][k, , drop = FALSE], k = k,
  #										control = list(fnscale = -1), method = "L-BFGS-B",
  #										lower = intervalLambda[1])#, upper = intervalLambda[2])
  
  
  #lambda_hat <- matrix(estLambda$par, ncol = S, byrow = TRUE)
  #  }
  #}
  #}
  #lambda_hat <- estimatePrecisionParamMap(lambda_prev, n_prev, phi_prev, params, D, rho)
  
  
  #rm(a, b, c, d, e)
  gc(verbose = FALSE, reset = TRUE)
  #if (verbose == TRUE) {
  #  message("n=", format(n_hat, digits = 2), ", phi=", format(phi_hat, digits = 4), ", lambda=", paste0(format(lambda_hat, digits=2), collapse = " ", sep = ","))
  #}
  return(list(n = n_hat, sp = sp_hat, phi = phi_hat, lambda = lambda_hat, piG = pi_hat, A = A_hat, F = estLambda$value))
}
```

--------------------------------

### Calculate Prior Probabilities

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/EM.R

Calculates log-prior probabilities for various model parameters including transitions, precision, normal contamination, subclone fraction, ploidy, and initial distribution.

```R
priorProbs <- function(n, sp, phi, lambda, piG, A, params, 
                       estimateNormal = TRUE, estimatePloidy = TRUE,
                       estimatePrecision = TRUE, estimateTransition = TRUE,
                       estimateInitDist = TRUE, estimateSubclone = TRUE){
  S <- params$numberSamples
  K <- length(params$ct)
  KS <- K ^ S
  ## prior terms ##
  priorA <- 0
  if (estimateTransition){
    for (ks in 1:KS){
      priorA <- priorA + dirichletpdflog(A[ks, ], params$dirPrior[ks, ])
    }
  }
  priorLambda <- 0
  if (estimatePrecision){
    for (s in 1:S){
      for (k in 1:K){
        priorLambda <- priorLambda + 
          gammapdflog(as.data.frame(lambda)[k, s], params$alphaLambda[k], params$betaLambda[k, s])
      }
    } 
  }
  priorLambda <- as.numeric(priorLambda)
  priorN <- 0
  if (estimateNormal){
    priorN <- sum(betapdflog(n, params$alphaN, params$betaN))
  }
  priorSP <- 0
  if (estimateNormal){
    priorSP <- sum(betapdflog(sp, params$alphaSp, params$betaSp))
  }
  priorPhi <- 0
  if (estimatePloidy){
    priorPhi <- sum(gammapdflog(phi, params$alphaPhi, params$betaPhi))
  }
  priorPi <- 0
  if (estimateInitDist){
    priorPi <- dirichletpdflog(piG, params$kappa)
  }
  prior <- priorA + priorLambda + priorN + priorSP + priorPhi + priorPi  
  return(list(prior = prior, priorA = priorA, priorLambda = priorLambda, 
              priorN = priorN, priorSP = priorSP, priorPhi = priorPhi, priorPi = priorPi))
}
```

--------------------------------

### Load ichorCNA and Prepare Data

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Loads the ichorCNA package or custom scripts, retrieves sequence information, and processes input files.

```R
## 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"))
} else {
    library(ichorCNA)
}

## load seqinfo 
seqinfo <- getSeqInfo(genomeBuild, genomeStyle)

if (substr(tumour_file,nchar(tumour_file)-2,nchar(tumour_file)) == "wig") {
  wigFiles <- data.frame(cbind(patientID, tumour_file))
} else {
  wigFiles <- read.delim(tumour_file, header=F, as.is=T)
}

## FILTER BY EXONS IF PROVIDED ##
## add gc and map to GRanges object ##
if (is.null(exons.bed) || exons.bed == "None" || exons.bed == "NULL"){
  targetedSequences <- NULL
}else{
  targetedSequences <- read.delim(exons.bed, header=T, sep="\t")  
}

## load PoN
if (is.null(normal_panel) || normal_panel == "None" || 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 <- read.delim(centromere,header=T,stringsAsFactors=F,sep="\t")
save.image(outImage)
```

--------------------------------

### Select Optimal HMM Solution

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Filters and sorts HMM solutions based on log-likelihood and subclonal prevalence constraints.

```R
### SELECT SOLUTION WITH LARGEST LIKELIHOOD ###
loglik <- loglik[!is.na(loglik$init), ]
if (estimateScPrevalence){ ## sort but excluding solutions with too large % subclonal 
	fracInd <- which(loglik[, "Frac_CNA_subclonal"] <= maxFracCNASubclone & 
						 		   loglik[, "Frac_genome_subclonal"] <= maxFracGenomeSubclone)
	if (length(fracInd) > 0){ ## if there is a solution satisfying % subclonal
		ind <- fracInd[order(loglik[fracInd, "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) 
}
```

--------------------------------

### Save R Workspace

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/scripts/runIchorCNA.R

Saves the current R environment to a file.

```R
### SAVE R IMAGE ###
save.image(outImage)
#save(tumour_copy, results, loglik, file=paste0(outDir,"/",id,".RData"))
```

--------------------------------

### Retrieve Sequence Information

Source: https://rdrr.io/github/broadinstitute/ichorCNA/src/R/utils.R

Fetches sequence information for a specified genome build and style, falling back to generic Seqinfo if the BSgenome package is unavailable.

```R
getSeqInfo <- function(genomeBuild = "hg19", genomeStyle = "NCBI"){
	bsg <- paste0("BSgenome.Hsapiens.UCSC.", genomeBuild)
	if (!require(bsg, character.only=TRUE, quietly=TRUE, warn.conflicts=FALSE)) {
		seqinfo <- Seqinfo(genome=genomeBuild)
	} else {
		seqinfo <- seqinfo(get(bsg))
	}
	seqlevelsStyle(seqinfo) <- genomeStyle
	seqinfo <- keepSeqlevels(seqinfo, value = chrs)
	#seqinfo <- cbind(seqnames = seqnames(seqinfo), as.data.frame(seqinfo))
	return(seqinfo)	
}
```