SVCFit

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).

Resources

Open example and benchmark data on Mendeley Data.
Controlled-access study data: European Genome-phenome Archive accession EGAD00001001343.
Prostate mixture scripts.

Installation

SVCFit is public. Installation does not require a GitHub account or personal access token. Install the current release from GitHub with remotes:

if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}
remotes::install_github("KarchinLab/SVCFit")

library(SVCFit)
packageVersion("SVCFit")

Building the source package and vignette

Ordinary installation does not require building the vignette. Developers and users who install with build_vignettes = TRUE need the Quarto CLI and the R quarto package. With Quarto on PATH, the repository launcher uses it directly:

Rscript -e 'install.packages("quarto")'
tools/with_quarto.sh quarto check
tools/with_quarto.sh R CMD build .

For an isolated Conda installation, the launcher also supports this setup:

mamba create -y -p "$HOME/.local/share/svcfit-quarto" -c conda-forge quarto
tools/with_quarto.sh quarto check

Set SVCFIT_QUARTO_PREFIX if that environment is stored elsewhere. To install from GitHub and build the vignette, use remotes::install_github("KarchinLab/SVCFit", build_vignettes = TRUE) after Quarto is available.

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		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		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]. It retains FACETS total copy number, minor copy number, and segment cellular fraction when annotating overlapping CNVs. The SVCF equations use total/minor copy state and the SNP-derived allele copy ratio; FACETS cf.em is retained as cncf for annotation and audit rather than used directly by calc_svcf().

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 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)

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])

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])

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])

# 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:

# 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)

# 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])

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

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

# 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:

Load input data — load_data()
Process BND events — proc_bnd()
Parse SV metadata — parse_sv_info()
Parse heterozygous SNPs — parse_het_snps(), parse_snp_on_sv()

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:

Assign SV IDs to SNPs — assign_svids()
Summarizes phasing + zygosity — sum_sv_info()
Assign CNV to SV — assign_cnv()
Annotate overlapping CNV — annotate_cnv(), parse_snp_on_sv()

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-like data frame.

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	Noise buffer around the sign-based SV/CNV ordering criterion; `ss1 <= thresh` selects the alternate branch, and deletions always use it.
`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	Rows with `sv_ref = 0` that independent evidence confirms are genuine clonal events; reviewed rows are recovered at the SVCF boundary of 1.

Output: An annotated VCF-like data frame with additional fields for VAF, the read-count normalization, carrier-copy multiplicity, and SVCF. The alternate overlapping-CNV estimate is retained as ss2_raw; ss2_constraint_status records whether it was in range or constrained to 0 or 1. On hemizygous chromosomes the output also carries pl (local normal ploidy), svcf_status, sv_cnv_order, and svcf_is_bound.

VAF: variant allele frequency
Rbar: class-specific read-count normalization; for the isolated diploid-duplication derivation, 2(A+B)/B estimates sample-average copy load under the read model
r: carrier-cell total copy number from FACETS under the duplication model; r-2 is the added-copy/junction multiplicity
SVCF: structural variant cellular fraction.

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

This step builds a tumor-evolution tree from SV clusters obtained with the Dirichlet-process Gaussian mixture model (DP-GMM). The current interface is designed for paired longitudinal samples.

data_dir may point directly to the directory containing <sample>.bed files. For compatibility with existing workflow output, SVCFit also recognizes SVCFit_output/ and COMBAT/SVCFit_output/ beneath that directory.

cluster_result <- cluster_data(
  pair_path = pair_path,
  pur_path = pur_path,
  data_dir = data_dir,
  pair_num = 1
)
clones <- cluster_result[[3]]

tree_result <- 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	—	Directory containing per-sample SVCF BED files, or a run root containing `SVCFit_output/` or `COMBAT/SVCFit_output/`.
`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	CCF values below this threshold are set to zero before clustering.

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

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:

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.

5.2 attach clonal assignment to output

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:

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. FACETS supplies the autosomal copy-number inputs used by the standard estimator. On a hemizygous chromosome, allele-specific estimates may be unreliable because there are no heterozygous germline SNPs to provide an allelic-imbalance signal. Supply cn_bar from an independent read-depth analysis. Do not derive it from the same SV’s SVCFit estimate, because 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:

Bin chrX read depth (e.g. 100 kb bins) in both the tumor and the matched-normal BAM (samtools depth).
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.
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. This scaling is required whenever the sample mean autosomal ploidy differs from two.
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.

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

The installed vignette is available when SVCFit was installed with build_vignettes = TRUE:

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

Compare two SVCFit runs

The package includes a command-line utility for event-level comparison of two SVCFit result sets. Each input may be a tab-delimited table, an RDS data frame, or a run directory containing per-sample BED files:

COMPARE_SCRIPT=$(Rscript -e \
  'cat(system.file("scripts", "compare_svcf_runs.R", package = "SVCFit"))')
Rscript "$COMPARE_SCRIPT" old_run_or_table new_run_or_table comparison_output

The utility requires unique event keys and writes the joined event table, overall and grouped summaries, source-file checksums, and a Markdown report. It never modifies either input.

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? A matched normal is recommended. For a tumor-only SV VCF, set tum_only = TRUE; SVCFit still requires compatible heterozygous-SNP and copy-number inputs. Tumor purity is not used to calculate SVCF. It is required only when converting SVCF to cancer cell fraction (CCF) for downstream clustering.

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 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

Chiang, C. et al. SpeedSeq: ultra-fast personal genome analysis and interpretation. Nat Methods 12, 966–968 (2015).
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).
Van Der Auwera, G. A. & O’Connor, B. D. Genomics in the Cloud. (O’Reilly Media, 2020).
Danecek, P. et al. Twelve years of SAMtools and BCFtools. GigaScience 10, giab008 (2021).
Chen, X. et al. Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications. Bioinformatics 32, 1220–1222 (2016).
Rausch, T. et al. DELLY: structural variant discovery by integrated paired-end and split-read analysis. Bioinformatics 28, i333–i339 (2012).
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).
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).

Repository scope and example data

This repository contains the reusable SVCFit R package, API documentation, usage guide, tests, and the de-identified example data used by the vignette. The example and plotting data are bundled under inst/extdata and available through system.file() after installation; no external data path is needed.

Analysis pipelines, benchmark drivers, cluster submission scripts, and figure-generation tools belong in the separate svcfit_workflows repository. Large research datasets and generated results remain outside both repositories.

Name		Name	Last commit message	Last commit date
R		R
inst		inst
man		man
tests		tests
tools		tools
vignettes		vignettes
.Rbuildignore		.Rbuildignore
.gitignore		.gitignore
DESCRIPTION		DESCRIPTION
LICENSE		LICENSE
LICENSE.md		LICENSE.md
NAMESPACE		NAMESPACE
README.Rmd		README.Rmd
README.md		README.md
SVCFit.Rproj		SVCFit.Rproj

Latest commit

History

Folders and files

Repository files navigation

SVCFit

Installation

Building the source package and vignette

Input Requirements

1. Structural variants

2. Copy number variants

3. Heterozygous SNP near SV

4. Heterozygous SNP on SV supporting reads

Example Preprocessing Pipeline

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

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

3. Somatic SV calling (Manta[5])

4. SV genotyping (SVtyper[1])

5. Germline heterozygous SNP detection (GATK4 + bcftools)

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

7. Heterozygous SNPs pileup on SV supporting reads

Caller-specific quirks

Usage Workflow

1. Extract SV and SNP Information — extract_info()

Function Arguments

extract_info()

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

Function Arguments

characterize_sv()

3. Calculate SVCF for Structural Variants — calc_svcf()

Function Arguments

4. Build tumor evolution tree — build_tree()

Function Arguments

5. Simulation & Benchmarking

Hemizygous chromosomes (chrX / chrY)

Enabling it

Obtaining cn_bar (chrX/chrY depth segmentation with DNAcopy)

The format SVCFit actually needs — DNAcopy is replaceable

Estimator helpers

Tutorial

Compare two SVCFit runs

Frequently Asked Questions

Tool Versions

Reference

Repository scope and example data

About

Resources

Stars

Watchers

Forks

Releases

Packages

Contributors

Languages

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

`extract_info()`

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

`characterize_sv()`

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

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

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