JaBbA (Junction Balance Analysis)

 _____         ___    _      _____ 
(___  )       (  _`\ ( )    (  _  )
    | |   _ _ | (_) )| |_   | (_) |
 _  | | /'_` )|  _ <'| '_`\ |  _  |
( )_| |( (_| || (_) )| |_) )| | | |
`\___/'`\__,_)(____/'(_,__/'(_) (_)

(Junction     Balance     Analysis)

JaBbA builds a genome graph based on junctions and read depth from whole genome sequencing, inferring optimal copy numbers for both vertices (DNA segments) and edges (bonds between segments). It can be used for discovering various patterns of structural variations.

If you use JaBbA in your work, please cite: Distinct Classes of Complex Structural Variation Uncovered across Thousands of Cancer Genome Graphs

Table of contents

• Table of contents
• Installation
• Usage
• Output
• Best practice
• FAQ
• The design of JaBbA
• Attributions
• Funding sources
• Fun fact

Installation

Dockerized installation

A Docker container with the latest stable release of JaBbA v1.1 and its dependencies installed can be found here.

:warning: Warning:

JaBbA depends on CPLEX MIP Optimizer to work. Because CPLEX is a proprietary software, it isn't included in the image and needs to be installed by the user. To add CPLEX:

1. Download CPLEX (Linux x86-64). (You may need to use the HTTP method.)
2. Run the add_mip_solver.sh script and enter the requested information.

Alternatively, you can follow these steps to install CPLEX manually:

docker pull mskilab/jabba:latest
docker run -it --rm --platform linux/amd64 mskilab/jabba:latest

3. Copy CPLEX binary into the container: docker cp /PATH/TO/DOWNLOADED_CPLEX.bin CONTAINER_ID:/opt/cplex_studio
4. Install CPLEX: /opt/cplex_studio/DOWNLOADED_CPLEX.bin (If you get a Permission denied error, run chmod 777 /PATH/TO/DOWNLOADED_CPLEX.bin before copying it into the container.)
5. When prompted for an installation path, type /opt/cplex. This is what the CPLEX_DIR environmental variable is set to.
6. Save changes to a new image for future use: Exit container (type exit or press Ctrl-D) Run docker commit CONTAINER_ID NEW_IMAGE_ID

Installation from GitHub

1. Install IBM ILOG CPLEX Studio. The software is proprietary, but can be obtained for free under IBM's academic initiative.

   Note: as of v1.1, JaBbA now supports the use of Gurobi for MIP optimization as an alternative to CPLEX. Please see below for installation instructions.

2. Set CPLEX_DIR variable to your CPLEX Studio installation directory

export CPLEX_DIR=/path/to/your/copy/of/CPLEX_Studio/

NOTE: if CPLEX_DIR is set correctly then $CPLEX_DIR/cplex/include and $CPLEX_DIR/cplex/lib should both exist.

3. Install JaBbA

devtools::install_github('mskilab-org/JaBbA')

4. For convenience, add jba executable to your PATH

$ JABBA_PATH=$(Rscript -e 'cat(paste0(installed.packages()["JaBbA", "LibPath"], "/JaBbA/extdata/"))')
$ export PATH=${PATH}:${JABBA_PATH}
$ jba ## to see usage

5. Test run jba executable on provided toy data

$ jba ${JABBA_PATH}/junctions.vcf ${JABBA_PATH}/coverage.txt 

:warning: WARNING : JaBbA uses some functions from gGnome package as a wrapper that needsCPLEX functionality to be added to gGnome. If the dependency package gGnome was installed before installing JaBbA or if CPLEX functionality was not added while installing gGnome or if you get an error while running JaBbA as shown below:

Error in names(x) <- c("xopt", "obj", "status", "extra", "epgap") :
  'names' attribute [5] must be the same length as the vector [1]
Error in .C("Rcplex_free", PACKAGE = "gGnome") :
  "Rcplex_free" not available for .C() for package "gGnome"

This means gGnome was installed without CPLEX functionality. Currently, JaBbA checks if gGnome has CPLEX functionality added while installing and loading the JaBbA package, if not, it will reinstall gGnome to add the CPLEX functionality, so it should take some time for applying this optimization during installation.

Running JaBbA with Gurobi

As of v1.1 JaBbA can now be run with Gurobi, which offers a free academic license. To use Gurobi with JaBbA, please follow the following steps for installation:

1. Install Gurobi and set up a license.
2. Install the R interface
3. Set up your GUROBI_HOME environment variable:

export GUROBI_HOME=[your path to gurobi]

4. Install JaBbA

devtools::install_github('mskilab-org/JaBbA')

5. For convenience, add jba executable to your PATH

$ JABBA_PATH=$(Rscript -e 'cat(paste0(installed.packages()["JaBbA", "LibPath"], "/JaBbA/extdata/"))')
$ export PATH=${PATH}:${JABBA_PATH}
$ jba ## to see usage

6. Note that to run JaBbA with Gurobi you will need to use the ---gurobi flag:

$ jba ${JABBA_PATH}/junctions.vcf ${JABBA_PATH}/coverage.txt --gurobi TRUE

Usage

 _____         ___    _      _____ 
(___  )       (  _`\ ( )    (  _  )
    | |   _ _ | (_) )| |_   | (_) |
 _  | | /'_` )|  _ <'| '_`\ |  _  |
( )_| |( (_| || (_) )| |_) )| | | |
`\___/'`\__,_)(____/'(_,__/'(_) (_)

(Junction     Balance     Analysis)

Usage: jba JUNCTIONS COVERAGE [options]
	JUNCTIONS can be BND style vcf, bedpe, rds of GrangesList
 	COVERAGE is a .wig, .bw, .bedgraph, .bed., .rds of a granges, or .tsv  .csv /.txt  file that is coercible to a GRanges
	use --field=FIELD argument so specify which column to use if specific meta field of a multi-column table


Options:
	--j.supp=J.SUPP
		supplement junctions in the same format as junctions

	--blacklist.junctions=BLACKLIST.JUNCTIONS
		rearrangement junctions to be excluded from consideration

	--whitelist.junctions=WHITELIST.JUNCTIONS
		rearrangement junctions to be forced to be incorporated

	--geno
		whether the junction has genotype information

	--indel=INDEL
		character of the decision to 'exclude' or 'include' small(< min.nbins * coverage bin width) isolated INDEL-like events into the model. Default NULL, do nothing.

	--cfield=CFIELD
		junction confidence meta data field in ra

	--tfield=TFIELD
		tier confidence meta data field in ra. tiers are 1 = must use, 2 = may use, 3 = use only in iteration>1 if near loose end. Default 'tier'.

	--iterate=ITERATE
		the number of extra re-iterations allowed, to rescue lower confidence junctions that are near loose end. Default 0. This requires junctions to be tiered via a metadata field tfield.

	--window=WINDOW
		window size in bp within which to look for lower confidence junctions. Default 1000.

	--nudgebalanced=NUDGEBALANCED
		whether to attempt to add a small incentive for chains of quasi-reciprocal junctions.

	--edgenudge=EDGENUDGE
		numeric hyper-parameter of how much to nudge or reward aberrant junction incorporation. Default 0.1 (should be several orders of magnitude lower than average 1/sd on individual segments), a nonzero value encourages incorporation of perfectly balanced rearrangements which would be equivalently optimal with 0 copies or more copies.

	--strict
		if used will only include junctions that exactly overlap segs

	--allin
		if TRUE will put all tiers in the first round of iteration

	--field=FIELD
		name of the metadata column of coverage that contains the data. Default 'ratio' (coverage ratio between tumor and normal). If using dryclean, it is 'foreground'.

	-s SEG, --seg=SEG
		Path to .rds file of GRanges object of intervals corresponding to initial segmentation (required)

	--maxna=MAXNA
		Any node with more NA than this fraction will be ignored

	--blacklist.coverage=BLACKLIST.COVERAGE
		Path to .rds, BED, TXT, containing the blacklisted regions of the reference genome

	--nseg=NSEG
		Path to .rds file of GRanges object of intervals corresponding to normal tissue copy number, needs to have $cn field

	--hets=HETS
		Path to tab delimited hets file output of pileup with fields seqnames, start, end, alt.count.t, ref.count.t, alt.count.n, ref.count.n

	--ploidy=PLOIDY
		Ploidy guess, can be a length 2 range

	--purity=PURITY
		Purity guess, can be a length 2 range

	--ppmethod=PPMETHOD
		select from 'ppgrid', 'ppurple', and 'sequenza' to infer purity and ploidy if not both given. Default, 'sequenza'.

	--cnsignif=CNSIGNIF
		alpha value for CBS

	--slack=SLACK
		Slack penalty to apply per loose end

	--linear
		if TRUE will use L1 loose end penalty

	-t TILIM, --tilim=TILIM
		Time limit for JaBbA MIP

	--epgap=EPGAP
		threshold for calling convergence

	-o OUTDIR, --outdir=OUTDIR
		Directory to dump output into

	-n NAME, --name=NAME
		Sample / Individual name

	--cores=CORES
		Number of cores for JaBBa MIP
		
	--lp
		Run as LP
	
	--ism
		Include infinite sites constraints
	
	--gurobi
		Use Gurobi for MIP optimization instead of CPLEX

	-v, --verbose
		verbose output

	-h, --help
		Show this help message and exitsage: jba JUNCTIONS COVERAGE [options]
        JUNCTIONS can be BND style vcf, bedpe, rds of GrangesList
        COVERAGE is a .wig