README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# SVCFit <img src="inst/extdata/svcfit_logo.png" align="left" width="60"/>

SVCFit is a fast and scalable computational tool designed to estimate the Structural Variant Cellular Fraction (SVCF) of inversions, deletions, tandem duplications, and translocations. Developed for the R environment, SVCFit integrates structural variant (SV) calls with Copy Number Variation (CNV) and Single Nucleotide Polymorphism (SNP) data to provide accurate cellular fraction estimates.

SVCFit also handles **hemizygous (single-copy) chromosomes** — for example chrX and chrY in a male subject — where heterozygous germline SNPs do not exist and the standard allele-copy-ratio correction cannot be estimated. On these chromosomes SVCFit switches to a ploidy-aware estimator that uses a read-depth–derived mean copy number (`cn_bar`) in place of SNP-based phasing. See [Hemizygous chromosomes (chrX / chrY)](#hemizygous-chromosomes-chrx--chry).

**Resources**

-   Open access data: It is available on mendeley (doi: 10.17632/2nhhdjx225.6)

-   Protected Data: Available via European Genome-phenome Archive (EGAD00001001343).

-   Prostate mixture scripts: [GitHub Repository](https://github.com/mcmero/SVclone_Rmarkdown/blob/master/make_insilico_mixtures.sh)

## Installation

SVCFit is hosted on GitHub. You can install it directly within R using the `remotes` package.

**Note:** Installation requires a GitHub Personal Access Token (PAT) because the repository is hosted on GitHub.\\.

``` r
if (!requireNamespace("remotes", quietly = TRUE))
    install.packages("remotes")

# 1. Setup GitHub Credentials (if not already configured)
if (!requireNamespace("usethis", quietly = TRUE))
    install.packages("usethis")

# Create a token in your browser
usethis::create_github_token() 

# Store the token (paste when prompted)
credentials::set_github_pat()

# 2. Install SVCFit
remotes::install_github("KarchinLab/SVCFit", build_vignettes = TRUE, dependencies = TRUE)
```

## Input Requirements

### 1. Structural variants

SVCFit accepts standard Variant Call Format (VCF) files. By default, the parser is optimized for VCFs produced by the SVTyper package [1].

| CHROM | POS | ID | REF | ALT | QUAL | FILTER | INFO | FORMAT | normal | tumor |
|----|----|----|----|----|----|----|----|----|----|----|
| chr1 | 1000 | INV:6:0:1:0:0:0 | T | <INV> | 100 | PASS | END=1500;SVTYPE=INV;SVLEN=500;... | GT:PR:SR:... | 0/1:76,0:70,0:... | 0/1:76,0:70,0:... |
| chr2 | 5000 | DEL:7:0:1:0:0:0 | G | <DEL> | 100 | PASS | END=5300;SVTYPE=DEL;SVLEN=300;... | GT:PR:SR:... | 0/1:76,0:70,0:... | 0/1:76,0:70,0:... |

Required **INFO** fields:

-   `SVTYPE` (e.g. INV, DEL, DUP, BND)
-   `END`

### 2. Copy number variants

SVCFit currently utilizes copy number calls from the FACETS package [2]. The tool uses allele-specific copy number, total copy number, and the cellular fraction (cncf) to annotate SVs with overlapping CNVs. No modification is required for standard FACETS output.

### 3. Heterozygous SNP near SV

SVCFit requires heterozygous SNP calls to phase overlapping CNVs. These can be generated using GATK4 [3] HaplotypeCaller and filtered using bcftools [4].

For computational efficiency, VCF can be filtered to include only heterozygous SNPs within 500bp of SV breakpoints.

The following code is not included in SVCFit and should be run separately.

```         
# 1. Call SNPs using GATK
gatk --java-options "-Xmx4g" HaplotypeCaller \
    -R $ref \
    -I $normal_bam \
    -O $snp_dir/SNP.vcf.gz

# 2. Filter for heterozygous SNPs using bcftools
bcftools view -v snps -g het -Oz -o $snp_dir/het_snp.vcf.gz $snp_dir/SNP.vcf.gz
tabix -p vcf $snp_dir/het_snp.vcf.gz
```

### 4. Heterozygous SNP on SV supporting reads

To infer SV phasing, SVCFit specifically examines heterozygous SNPs found on reads that support the structural variant. This is done following the steps: 1) Extract SV-supporting reads from the BAM file using samtools [4]. 2) Generate a read pileup using bcftools restricted to the heterozygous SNP positions identified in step 3.

The following code is not included in SVCFit and should be run separately.

```         
# Extract SV supporting reads
samtools view -f 1 -F 2 -b $tumor_bam > $snp_dir/sup_$samp_name.bam
samtools index $snp_dir/sup_$samp_name.bam

# Generate pileup at known heterozygous sites
# Note: pos_$samp_name.bed should contain the positions from het_snp.vcf.gz above
bcftools mpileup -f $ref -a DP,AD -A \
    -R $snp_dir/pos_$samp_name.bed \
    $snp_dir/sup_$samp_name.bam -Ov > $snp_dir/het_on_sv_$samp_name.vcf
```

## Example Preprocessing Pipeline

The commands below provide a complete worked example for generating SVCFit inputs from raw BAM files. Each tumor-versus-normal biopsy pair is processed as a separate sample. Tool versions used in the manuscript are listed in the [Tool Versions](#tool-versions) table.

The expected input layout:

```         
sample/
├── tumor.bam        (+ .bai)
├── normal.bam       (+ .bai)
├── ref.fa           (+ .fai, .dict)
└── intervals.bed    (optional)
```

### 1. Trim and align reads (if starting from FASTQ)

``` bash
trim_galore --paired tumor_R1.fq.gz tumor_R2.fq.gz -o trimmed/
bwa mem -t 8 ref.fa trimmed/tumor_R1_val_1.fq.gz trimmed/tumor_R2_val_2.fq.gz \
  | samtools sort -@ 4 -o tumor.bam -
samtools index tumor.bam
# repeat for the matched normal
```

### 2. Mark duplicates and recalibrate (GATK4[3])

``` bash
gatk MarkDuplicates -I tumor.bam -O tumor.md.bam -M tumor.metrics.txt
gatk BaseRecalibrator -I tumor.md.bam -R ref.fa --known-sites known.vcf.gz -O tumor.bqsr.table
gatk ApplyBQSR        -I tumor.md.bam -R ref.fa --bqsr-recal-file tumor.bqsr.table -O tumor.recal.bam
# repeat for the matched normal
```

### 3. Somatic SV calling (Manta[5])

``` bash
configManta.py --tumorBam tumor.recal.bam --normalBam normal.recal.bam \
               --referenceFasta ref.fa --runDir manta_run/
manta_run/runWorkflow.py
# output: manta_run/results/variants/somaticSV.vcf.gz
```

### 4. SV genotyping (SVtyper[1])

``` bash
# Ensure CIPOS and CIEND INFO fields are present; SVtyper requires them.
svtyper -B tumor.recal.bam -i tumor.vcf -o tumor.gt.vcf
```

The output VCF has `AO` (SV-supporting reads) and `RO` (reference reads) per breakpoint, read by SVCFit as `BPC` and `BEC` after scaling by mean read depth.

For multi-caller consensus (recommended for clinical samples), merge with SURVIVOR:

``` bash
# vcf_list.txt contains one column where each row is the path to VCF from a SV caller 
SURVIVOR merge vcf_list.txt 500 1 1 1 0 30 consensus.vcf
```

### 5. Germline heterozygous SNP detection (GATK4 + bcftools)

``` bash
# 1. Call SNPs using GATK
gatk --java-options "-Xmx4g" HaplotypeCaller \
    -R $ref \
    -I $normal_bam \
    -O $snp_dir/SNP.vcf.gz

# 2. Filter for heterozygous SNPs using bcftools (can be filtered for computation efficiency)
bcftools view -v snps -g het -Oz -o $snp_dir/het_snp.vcf.gz $snp_dir/SNP.vcf.gz
tabix -p vcf $snp_dir/het_snp.vcf.gz
```

### 6. Allele-specific copy-number profile (FACETS[2])

``` bash
snp-pileup -g -q15 -Q20 -P100 -r25,0 normal.het.vcf.gz tumor.snp.csv \
           normal.recal.bam tumor.recal.bam
```

``` r
library(facets)
rcmat = readSnpMatrix(SNP_file)
xx=preProcSample(rcmat,ndepth=20, gbuild="hg38", cval=50)
oo = procSample(xx, cval=150, dipLogR=NULL)
ooo= procSample(xx, cval=500, dipLogR=oo$dipLogR)
fit = emcncf(ooo)
write.table(fit$cncf, file = "tumor.facets.tsv", sep = "\t", quote = FALSE, row.names = FALSE)
```

### 7. Heterozygous SNPs pileup on SV supporting reads

``` bash
# 1. Extract SV supporting reads
samtools view -f 1 -F 2 -b $tumor_bam > $snp_dir/sup_$samp_name.bam
samtools index $snp_dir/sup_$samp_name.bam

# 2. Generate pileup at known heterozygous sites that's on SV supporting reads
# Note: pos_$samp_name.bed should contain the positions from het_snp.vcf.gz above
bcftools mpileup -f $ref -a DP,AD -A \
    -R $snp_dir/pos_$samp_name.bed \
    $snp_dir/sup_$samp_name.bam -Ov > $snp_dir/het_on_sv_$samp_name.vcf
```

### Caller-specific quirks

-   **Manta[5]**: writes `INFO/CIPOS` and `INFO/CIEND` natively — no reformatting needed.
-   **Delly[6]**: writes `INFO/CIPOS` but not `INFO/CIEND` — append `CIEND=-50,50` (or the per-call confidence interval if available).
-   **GRIDSS[7]**: uses paired-end / split-read counts in dedicated INFO fields; convert to LUMPY-style `INFO/MATEID`, `INFO/CIPOS`, and `INFO/CIEND` before SVtyper.
-   **Manta + Delly + GRIDSS consensus via SURVIVOR[8]**: SURVIVOR drops `CIPOS`/`CIEND` from some merged records — backfill them before genotyping.

## Usage Workflow

The SVCFit pipeline consists of three main steps: **Extraction**, **Characterization**, and **Calculation**.

### 1. Extract SV and SNP Information — `extract_info()`

Load and preprocess input VCF and CNV files. This step parses metadata, processes breakends (BND), and handles heterozygous SNPs. `extract_info()` internally performs:

1.  **Load input data** — `load_data()`
2.  **Process BND events** — `proc_bnd()`
3.  **Parse SV metadata** — `parse_sv_info()`
4.  **Parse heterozygous SNPs** — `parse_het_snps()`, `parse_snp_on_sv()`

``` r
info <- extract_info(
  p_het = "path/to/het_snps.vcf",
  p_onsv = "path/to/snps_on_sv.vcf",
  p_sv = "path/to/structural_variants.vcf",
  p_cnv = "path/to/cnv_file.txt",
  chr_lst = NULL,
  flank_del = 50, 
  QUAL_thresh = 100, 
  min_alt = 2, 
  tumor_only = FALSE
)
```

#### Function Arguments

#### `extract_info()`

| Argument | Type | Default | Description |
|----|----|----|----|
| `p_het` | Character | — | Path to VCF of heterozygous SNPs. |
| `p_onsv` | Character | — | Path to VCF of SNPs overlapping SV-supporting reads. |
| `p_sv` | Character | — | Path to SV VCF. |
| `p_cnv` | Character | — | Path to CNV file. |
| `chr_lst` | Character | NULL | Chromosomes to include. |
| `flank_del` | numeric | 50 | Max distance to consider deletion overlapping a BND. |
| `QUAL_thresh` | numeric | 100 | Minimum QUAL score. |
| `min_alt` | numeric | 2 | Minimum alternative reads. |
| `tum_only` | Logical | — | Whether SVs come from tumor-only calling. |

**Output:** A list of data frames containing parsed SV + SNP information.

### 2. Annotate SVs Using CNV and SNP Information — `characterize_sv()`

This step integrates CNV and heterozygous SNPs to infer phasing, zygosity, and overlapping CNV. `characterize_sv()` internally performs:

1.  **Assign SV IDs to SNPs** — `assign_svids()`
2.  **Summarizes phasing + zygosity** — `sum_sv_info()`
3.  **Assign CNV to SV** — `assign_cnv()`
4.  **Annotate overlapping CNV** — `annotate_cnv()`, `parse_snp_on_sv()`

``` r
sv_char <- characterize_sv(
  sv_phase = info$sv_phase, 
  sv_info = info$sv_info, 
  cnv = info$cnv,
  flank_snp = 500,
  flank_cnv = 1000
)
```

#### Function Arguments

#### `characterize_sv()`

| Argument    | Type       | Default | Description                       |
|-------------|------------|---------|-----------------------------------|
| `sv_phase`  | data.frame | —       | Phasing/zygosity from SNPs.       |
| `sv_info`   | data.frame | —       | Parsed SV metadata.               |
| `cnv`       | data.frame | —       | CNV data.                         |
| `flank_snp` | numeric    | 500     | Max assignment distance for SNPs. |
| `flank_cnv` | numeric    | 1000    | Max assignment distance for CNVs. |

### 3. Calculate SVCF for Structural Variants — `calc_svcf()`

This step computes the **Structural Variant Cellular Fraction (SVCF)**. and returns an annotated VCF file in data.frame format.

``` r

svcf_out <- calc_svcf(
  anno_sv_cnv = sv_char$anno_sv_cnv,
  sv_info     = sv_char$sv_info,
  thresh      = 0.1,
  samp        = "SampleID",
  exper       = "ExperimentID"
)
```

#### Function Arguments

| Argument      | Type       | Default | Description                            |
|---------------|------------|---------|----------------------------------------|
| `anno_sv_cnv` | data.frame | —       | CNV-annotated SVs.                     |
| `sv_info`     | data.frame | —       | Parsed SV info.                        |
| `thresh`      | numeric    | 0.1     | Threshold for SV-before-CNV inference. |
| `samp`        | character  | —       | Sample name.                           |
| `exper`       | character  | —       | Experiment name.                       |
| `hemizygous_chr` | character | NULL | Chromosomes single-copy in the germline (e.g. `c("chrX","chrY")`). `NULL` = diploid-only behavior. |
| `hemi_cn_bar` | data.frame / numeric | NULL | Read-depth mean copy number (`cn_bar`) for hemizygous rows — a `data.frame(CHROM, POS, cn_bar)` or numeric vector. Without it (or for an unmatched row) hemizygous **duplications** are left unresolved; deletions, inversions, and translocations fall back to `SVCF = VAF`. |
| `hemi_bg_cn` | data.frame / numeric | NULL | Flanking (background) copy number for the CNV-first deletion form; optional. |
| `hemi_dup_r` | numeric | 2 | Copies in carrier cells for a hemizygous tandem duplication; SVCFs reported as upper bounds. |
| `zero_ref_allowlist` | data.frame | NULL | Hemizygous rows with `sv_ref = 0` that BAM evidence confirms are genuine clonal losses (SVCF = VAF = 1). |

**Output:** An annotated VCF-like data frame with additional fields for VAF, Rbar, r, and SVCF. On hemizygous chromosomes it also carries `pl` (local normal ploidy), `svcf_status`, `sv_cnv_order`, and `svcf_is_bound`.

1.  VAF: variant allele frequency
2.  Rbar: average break interval count in a sample
3.  r: inferred integer copy number of break intervals
4.  SVCF: structural variant cellular fraction.

### 4. Build tumor evolution tree — `build_tree()`

This step build the tumor evolutionary tree based on SV clusters obtained from Dirichlet process Gaussian Mixture Model (DP-GMM).Currently, this step is optimized for two sample longitudinal data.

``` r
output <- cluster_data(
  pair_path,
  pur_path,
  data_dir,
  pair_num = 1)
clone2=output[[3]]

build_tree(
  clones,
  lineage_precedence_thresh=0.2, 
  sum_filter_thresh=0.2)
```

#### Function Arguments

cluster_data()

| Argument | Type | Default | Description |
|----|----|----|----|
| `pair_path` | character | — | Path to a tab-separated file with columns for 'pre_BAT sample' and 'on_BAT sample'. |
| `pur_path` | character | — | Path to a tab-separated file with columns for 'sample' and 'purity'. |
| `data_dir` | character | — | Path to the directory containing per-sample SVCF output files. |
| `Kmax` | numeric | 10 | Maximum number of clusters for DP-GMM. |
| `n_steps` | numeric | 100 | Number of DP-GMM iterations. |
| `thr_min_w` | numeric | 0.01 | Minimum cluster weight threshold. |
| `random_state` | integer | 0 | Random seed for reproducibility. |
| `concentration` | numeric | 1 | Dirichlet concentration parameter. |
| `min_n` | numeric | 5 | Minimum cluster size to retain. |
| `min_dist` | numeric | 0.2 | Minimum distance for merging nearby clusters. |
| `pair_num` | numeric | 1 | The identifier (index or ID) for the specific sample pair (patient) being analyzed. |
| `pairs` | integer vector | NULL | Subset of pair IDs to process; defaults to all pairs. |
| `exclude_pairs` | integer vector | integer(0) | Pair IDs to exclude from analysis. |
| `deduplicate` | Logical | TRUE | Whether to deduplicate events before clustering. |
| `ccf_floor` | numeric | 0.1 | Minimum CCF value before flooring. |

build_tree()

| Argument | Type | Default | Description |
|----|----|----|----|
| `clones` | data.frame | — | SV clustering result. |
| `lineage_precedence_thresh` | numeric | 0.2 | Maximum violation of lineage precedence rule. |
| `sum_filter_thresh` | numeric | 0.2 | Maximum violation of sum condition rule. |
| `linear_penalty` | numeric | 0 | Penalty applied to linear (chain) topologies during tree scoring. |

**Output:** A tumor evolutionary tree rooted at the germline (G). Node numbers correspond to SV cluster numbers. The branching depicts the chronological occurrence of SV clusters.

### 5. Simulation & Benchmarking

SVCFit includes utility functions for processing simulation data from VISOR and attaching "ground truth" labels to structural variants for benchmarking.

5.1 read clonal assignment

``` r
truth <- load_truth(
  truth_path = "path/to/truth_beds", 
  overlap = FALSE
  )
```

This function has 1 arguments:

| Argument | Type | Default | Description |
|----|----|----|----|
| `truth_path` | Character | N/A | Path to BED files storing true structural variant information with clonal assignment. Each BED file should be named like `"c1.bed, c2.bed"`, etc for non-overlapping simulations and `"c11.bed, c22.bed"`, etc for overlapping simulations. Structural variants should be saved in separate BED files if they belong to different (sub)clones. |
| `overlap` | Logical | FALSE | Whether the simulation has SV-CNV overlap. |

The file path should follow this structure:

``` r
root/
├── true_clone/
│   ├── c1.bed/
│   ├── c2.bed/
│   ├── c3.bed/
│   └── .../
```

Parent nodes should always have lower number in name than its children (i.e. c1.bed instead of c3.bed) and all child node bed file should conatin its ancestors mutations.

![](inst/extdata/tree.png){width="487"}

5.2 attach clonal assignment to output

``` r
svcf_truth <- attach_truth(svcf_out, truth)
```

This function has 2 arguments:

| Variable | Type | Default | Description |
|----|----|----|----|
| `svcf_out` | DataFrame | N/A | The output from `calc_svcf` |
| `truth` | DataFrame | N/A | Stores the clone assignment for each structural variant designed in a simulation. |

This appends the known clonal assignment to the calculated SVCF output for performance evaluation.

## Hemizygous chromosomes (chrX / chrY)

On a diploid autosome SVCFit uses heterozygous germline SNPs to estimate the allele copy ratio and phase each SV against overlapping CNVs. A **hemizygous** chromosome — a male X or Y, or any chromosome that is single-copy in the germline — has no heterozygous SNPs, so that route is undefined and the diploid conversion factor of 2 no longer applies. Left uncorrected, a clonal hemizygous SV is estimated at roughly twice its true cellular fraction and is then silently dropped downstream.

SVCFit handles these chromosomes with a ploidy-aware estimator. The key quantity is `cn_bar`, the **mean copy number of the locus per cell**, measured from read depth (`cn_bar = R * psi_sample / 2`, where `R` is the tumor/normal depth ratio and `psi_sample` the autosomal ploidy scaling). Total alleles per cell is always `cn_bar`, and with `VAF = BPC / (BPC + BEC)` the two orderings of an SV relative to an overlapping CNV are:

```         
H1  (SV precedes CNV):  SVCF = cn_bar * VAF - (cn_bar - 1)
H2  (CNV precedes SV):  SVCF = cn_bar * VAF
```

Both reduce to `SVCF = VAF` at `cn_bar = 1` (copy-neutral). The ordering is decided by the sign of the H1 form — no integer copy number or CNV cellular fraction is needed. A tandem duplication is its own copy-number change and uses `SVCF = (cn_bar - 1) / (r - 1)`; because `r` is not identifiable from a single locus, these are reported as **upper bounds**.

### Enabling it

Pass the single-copy chromosomes to `run_svcfit()` (or `calc_svcf()`), together with a per-SV `cn_bar` table from a depth segmentation of the chromosome:

``` r
result <- run_svcfit(
  p_het = p_het, p_onsv = p_onsv, p_sv = p_sv, p_cnv = p_cnv,
  samp = "SampleID", exper = "ExperimentID",
  hemizygous_chr = c("chrX", "chrY"),   # single-copy in this subject's germline
  hemi_cn_bar    = cn_bar_table         # data.frame(CHROM, POS, cn_bar), from read depth
)
```

With `hemi_cn_bar = NULL`, hemizygous SVs are treated as copy-neutral and resolve to `SVCF = VAF` — except tandem duplications, which need `cn_bar` and are left unresolved (`final_svcf = NA`, flagged `hemizygous_dup_needs_cn_bar` in `svcf_status`). A deletion or inversion on a genuinely copy-altered segment is *not* detectable without `cn_bar` on a hemizygous chromosome (there are no heterozygous SNPs to flag it), so it too falls back to `SVCF = VAF` rather than being singled out. Autosomes are unaffected: with `hemizygous_chr = NULL` every result is identical to the diploid path.

### Obtaining `cn_bar` (chrX/chrY depth segmentation with DNAcopy)

SVCFit does not compute `cn_bar` itself — it consumes it. On the autosomes copy number comes from FACETS, but **FACETS cannot fit a hemizygous chromosome**: a male X has no heterozygous germline SNPs to supply the allelic-imbalance signal FACETS relies on, so it returns a single whole-chromosome segment with an inflated, unusable total copy number (`tcn.em` of 3–9 in our cohort, contradicted by read depth). The containing SV's own SVCFit estimate cannot be used either — it is hemizygous and copy-altered, so that would be circular.

`cn_bar` is therefore measured directly from **read depth**, segmented with the Bioconductor **DNAcopy** package (circular binary segmentation, CBS). The recipe:

1.  **Bin** chrX read depth (e.g. 100 kb bins) in both the tumor and the matched-normal BAM (`samtools depth`).
2.  **Normalize** each BAM by its own autosomal 2-copy baseline (the median depth of a few known-diploid autosomal regions), then take the tumor/normal ratio. Because the matched normal is *also single-copy* on chrX, this ratio is the tumor's mean copy number per cell relative to the germline.
3.  **Scale by `psi_sample / 2`** to recover absolute copy number, where `psi_sample = purity * psi_tumor + (1 - purity) * 2` is the sample mean autosomal ploidy (from the FACETS autosomal fit, which *is* trustworthy). This step is not optional: across our cohort `psi_sample / 2` ran from 0.87 to 1.72.
4.  **Segment** `log2(cn_bar)` per bin with DNAcopy and report each segment's mean as `cn_bar = 2^seg.mean`.

```
library(DNAcopy)
# cn_bar_bin : per-bin (tumor/normal depth ratio) * (psi_sample / 2)
# pos        : bin start positions on the hemizygous contig
cna <- CNA(log2(cn_bar_bin), rep("chrX", length(pos)), pos,
           data.type = "logratio", sampleid = "chrX")
seg <- segment(smooth.CNA(cna), alpha = 0.01, min.width = 2)$output
seg_cn_bar <- 2^seg$seg.mean          # mean copy number per cell, per segment
```

`cn_bar` is the **final** product and is never rounded to an integer — the corrected hemizygous forms take it directly (there is no integer copy number `c` or CNV cellular fraction `f_CNV` to solve for; read depth alone cannot separate them on a hemizygous locus). A segment at true germline copy number reads near `1.0`; if none does, either the chromosome is wholly altered or `psi_sample` is wrong.

DNAcopy's own defaults split a flat chrX into dozens of noise segments on real data; the `alpha`, `min.width`, and `undo.SD` knobs control that and should be tuned against a known copy-neutral control before use.

#### The format SVCFit actually needs — DNAcopy is replaceable

**SVCFit never sees DNAcopy or its segment objects.** It consumes a plain per-SV table, so any method that yields a mean copy number per cell — a different segmenter (CBS, PSCBS, HMMcopy, GATK ModelSegments, CNVkit, Battenberg/ASCAT on the autosome-analogous signal), or even a single tumor/normal depth ratio per breakpoint — works as long as the output is coerced into this shape:

`hemi_cn_bar` — `data.frame` with **exactly these columns**:

| Column   | Type      | Meaning                                                                 |
|----------|-----------|-------------------------------------------------------------------------|
| `CHROM`  | character | Contig name, matching the SV's `CHROM` exactly (e.g. `"chrX"`, or `"X"`). |
| `POS`    | integer   | The SV breakpoint position, matching the SV's `POS` **exactly**.        |
| `cn_bar` | numeric   | Mean copy number of that locus per cell (not rounded to an integer).    |

The lookup is an exact `(CHROM, POS)` match against the SV rows — **one row per SV, keyed on the breakpoint**, not per segment. That is the one manual step: if your source produces segment *intervals* (`chrom, start, end, cn_bar`), assign each SV the `cn_bar` of the segment spanning its `POS` before passing the table. An SV whose breakpoint has no matching row gets no depth correction — a tandem duplication is then left unresolved (`final_svcf = NA`, flagged `hemizygous_dup_needs_cn_bar`), while a deletion, inversion, or translocation falls back to the copy-neutral `SVCF = VAF`. A bare numeric vector is also accepted, but it is recycled across **all** SV rows in internal order (not just the hemizygous ones), so in practice only a single value — the same `cn_bar` for every hemizygous SV — is reliable; use the keyed data frame otherwise.

```
# minimal hemi_cn_bar — however you produced cn_bar, this is all SVCFit reads
hemi_cn_bar <- data.frame(
  CHROM  = c("chrX",     "chrX"),
  POS    = c(31200000L,  67500000L),   # exact SV breakpoint positions
  cn_bar = c(1.02,       1.54)         # mean copies/cell at each locus
)
```

The optional `hemi_bg_cn` (flanking/background copy number for the CNV-first deletion form) uses the same convention with columns `CHROM, POS, bg_cn`. DNAcopy is only needed for this hemizygous depth step — it is **not** a dependency of the autosomal pipeline.

### Estimator helpers

These building blocks are exported so a `cn_bar` and read counts can be scored directly, without the full VCF pipeline:

| Function | Purpose |
|----|----|
| `local_ploidy(chrom, hemizygous_chr)` | Local normal ploidy (1 on a hemizygous chromosome, 2 otherwise). |
| `classify_cn(cna, minor, pl)` | Ploidy-aware DUP / norm / DEL classification (reduces to the diploid tests at `pl = 2`). |
| `resolve_hemizygous_svcf(bpc, bec, cn_bar)` | SVCF for an SV inside a CNV made by another event; picks H1/H2 by the sign rule. |
| `hemizygous_dup_svcf(cn_bar, r = 2)` | SVCF for a hemizygous tandem duplication (upper bound). |
| `hemizygous_del_svcf(bpc, bec, cn_bar)` | SVCF for a hemizygous deletion, with a depth-vs-read consistency check. |
| `svcf_status(pl, cn_type, sv_ref)` | Per-row status label so exclusions are countable, not silent. |

``` r
resolve_hemizygous_svcf(bpc = 4, bec = 1, cn_bar = 2.0)
#>   svcf      ordering status  h1  h2
#> 1  0.6 sv_before_cnv     ok 0.6 1.6

hemizygous_dup_svcf(cn_bar = 1.5, r = 2)
#>   svcf is_upper_bound status
#> 1  0.5           TRUE     ok
```

## Tutorial

``` r
library(SVCFit)
vignette("SVCFit_guide", package = "SVCFit")
```

## Frequently Asked Questions

**Can I use a different SV caller?** Yes. Any caller that produces a per-SV breakpoint VCF compatible with SVtyper will work. Differences across callers in breakpoint detection sensitivity, split-read versus discordant-pair definitions, and quality filtering may produce different REF/ALT counts and therefore different SVCF estimates from the same data — see Discussion in the manuscript.

**Can I skip FACETS and use Battenberg / ASCAT instead?** Yes, as long as you provide per-segment total copy number and gain/loss classification in the same TSV format.

**Do I need a matched normal?** For the COMBAT analysis we used matched normals throughout. SVCFit can run on tumor-only when a matched normal is unavailable, but tumor purity must be supplied externally and ASCN inference becomes less reliable without germline-heterozygous SNP calls; this is treated as an unsupported configuration in the current release.

**What about complex SVs (chromothripsis, BFB, chromoplexy)?** The closed-form SVCF estimators cover deletions, tandem duplications, inversions, and the three classes of translocations. Multi-breakpoint complex SVs are not yet explicitly modeled — the per-breakpoint estimates are still produced, but their interpretation as a single cellular fraction is approximate. Future releases will add structure-aware handling.

## Tool Versions {#tool-versions}

Tool versions used in the manuscript:

| Step                       | Tool                                  | Version  |
|----------------------------|---------------------------------------|----------|
| Read trimming              | trim_galore                           | v0.6.1   |
| Alignment                  | bwa mem                               | v0.7.19  |
| MarkDuplicates / BQSR      | GATK4                                 | v4.6.2.0 |
| Somatic SV (single-caller) | Manta                                 | v1.6.0   |
| Somatic SV (multi-caller)  | Manta + Delly v1.5.0 + GRIDSS v2.13.2 | —        |
| Multi-caller merge         | SURVIVOR                              | v1.0.7   |
| SV genotyping              | SVtyper                               | v0.7.1   |
| Germline SNP calling       | GATK4 HaplotypeCaller                 | v4.6.2.0 |
| Het-SNP filter             | bcftools                              | v1.20    |
| Copy number                | FACETS                                | v0.6.2   |
| Copy number (hemizygous chrX/chrY) | DNAcopy (circular binary segmentation) | Bioconductor |
| Breakpoint read filter     | samtools                              | v1.21    |
| SNP pileup                 | bcftools mpileup                      | v1.20    |

## Reference

1.  Chiang, C. *et al.* SpeedSeq: ultra-fast personal genome analysis and interpretation. *Nat Methods* **12**, 966–968 (2015).
2.  Shen, R. & Seshan, V. E. FACETS: allele-specific copy number and clonal heterogeneity analysis tool for high-throughput DNA sequencing. *Nucleic Acids Res* **44**, e131–e131 (2016).
3.  Van Der Auwera, G. A. & O’Connor, B. D. *Genomics in the Cloud*. (O’Reilly Media, 2020).
4.  Danecek, P. *et al.* Twelve years of SAMtools and BCFtools. *GigaScience* **10**, giab008 (2021).
5.  Chen, X. *et al.* Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications. *Bioinformatics* **32**, 1220–1222 (2016).
6.  Rausch, T. *et al.* DELLY: structural variant discovery by integrated paired-end and split-read analysis. *Bioinformatics* **28**, i333–i339 (2012).
7.  Cameron, D. L. *et al.* GRIDSS2: comprehensive characterisation of somatic structural variation using single breakend variants and structural variant phasing. *Genome Biol* **22**, 202 (2021).
8.  Jeffares, D. C. *et al.* Transient structural variations have strong effects on quantitative traits and reproductive isolation in fission yeast. *Nat. Commun.* **8**, 14061 (2017).