A Java tool for gene-level rare-variant burden analysis in case-control cohorts.
Developed in the Professor Polychronakos Laboratory, McGill University (2025).
Author: Milad Eidi
Gene Burden Analyzer tests whether rare, potentially damaging variants accumulate in cases more than expected by chance. For each gene it runs:
- A Dominant model — counts samples carrying at least one rare qualifying alt allele
- A Recessive model — counts samples carrying two or more qualifying alt alleles (homozygous or compound-heterozygous)
Expected carrier rates are derived from population allele frequencies (gnomAD-style control VCF) using closed-form probability formulas, and tested with a one-sided binomial exact test. Benjamini–Hochberg FDR correction is applied independently within each model.
The tool ships with a full graphical user interface (launch the JAR with no arguments) as well as a traditional command-line interface for batch/HPC use.
| Feature | Details |
|---|---|
| Dominant model | Carrier = ≥ 1 qualifying alt allele; expected rate = 1 − ∏(1 − AF)² |
| Recessive model | Carrier = ≥ 2 qualifying alt alleles; expected rate via haplotype convolution |
| Statistical test | One-sided binomial exact test + Clopper–Pearson 95% CI on carrier rate |
| FDR correction | Benjamini–Hochberg, independently per model |
| Pathogenicity filter | REVEL score threshold (default 0.65) OR gene-level CADD/MSC 95th percentile |
| Population AF filter | Max allele frequency cutoff in gnomAD population (default 0.005) |
| Joint AF | Expected rates use case+control joint AF by default for accuracy |
| Callable-sites filter | BED file restricts analysis to callable genomic positions |
| Variant types | Configurable list of mutation types (nonsynonymous, frameshift, stop, etc.) |
| Gene annotation | Built-in 71-column annotation database bundled in the JAR |
| HPO prioritization | HPO term → gene mapping; matching genes appear at the top of results |
| Visualizations | QQ plot (−log₁₀ p with λ) and top-genes bar chart — saved as PNG |
| GUI | Swing application with 5 tabs: Parameters, Results, QQ Plot, Top Genes, Log |
| CLI | Full headless mode for servers and pipelines |
| Requirement | Version |
|---|---|
| Java | 17 or later |
| Maven | 3.6+ (only needed to build from source) |
No other installation is required. All dependencies are either bundled in lib/ (for IDE use) or included in dist/GeneBurden.jar at runtime.
The compiled JAR is available at dist/GeneBurden.jar in this repository.
# GUI mode — opens automatically if no arguments are given
java -jar dist/GeneBurden.jar
# CLI mode
java -jar dist/GeneBurden.jar \
--cases cases.vcf.gz \
--control gnomad_exomes_r2.1.1_sites.vcf.bgz \
--population NFE \
--output-dir ./resultsgit clone https://github.com/MiladEidi/GeneBurdenAnalyzer.git
cd GeneBurdenAnalyzer
mvn package
# Output: dist/GeneBurden.jarLaunch the GUI by running the JAR with no arguments:
java -jar dist/GeneBurden.jarThe application opens a window with five tabs:
| Tab | Purpose |
|---|---|
| Parameters | Fill in all settings with file-browser buttons; click Run Analysis |
| Results | Sortable/filterable gene table; colour-coded by significance; Export TSV |
| QQ Plot | Genomic-inflation QQ plot displayed inline |
| Top Genes | Bar chart of the 20 most significant genes |
| Log | Real-time analysis log (stdout/stderr) |
The analysis runs in a background thread so the window stays responsive. Progress is visible in the Log tab.
| Argument | Description |
|---|---|
--cases <file> |
Case VCF/VCF.GZ with sample genotypes |
--control <file> |
Population VCF (e.g. gnomAD) with AC/AN INFO fields |
--population <CODE> |
gnomAD population code: NFE, AFR, AMR, EAS, SAS, FIN, ASJ |
--output-dir <path> |
Directory where output files are written (created if absent) |
| Argument | Default | Description |
|---|---|---|
--af-threshold <float> |
0.005 |
Maximum population allele frequency |
--revel-threshold <float> |
0.65 |
Minimum REVEL pathogenicity score |
--revel-field <name> |
REVEL |
VCF INFO field for REVEL score |
--use-cadd |
off | Enable gene-based CADD/MSC 95th-percentile filter |
--cadd-field <name> |
CADD |
VCF INFO field for CADD score |
--callable-bed <file> |
none | BED file of callable genomic sites (0-based half-open) |
--dominant-only |
off | Run dominant model only |
--recessive-only |
off | Run recessive model only |
--no-joint-af |
off | Use control-only AF for expected rates |
--mutation-types <list> |
see below | Comma-separated allowed variant types |
--variant-gene-field <name> |
Gene.refGene |
VCF INFO field for gene name |
--variant-type-field <name> |
ExonicFunc.refGene |
VCF INFO field for variant type |
--hpo-id <HP:XXXXXXX,...> |
none | Comma-separated HPO IDs to prioritise |
Default mutation types:
nonsynonymous_SNV, nonsynonymous SNV,
frameshift_deletion, frameshift_insertion,
frameshift deletion, frameshift insertion,
stopgain, stoploss, startgain, startloss
java -jar dist/GeneBurden.jar \
--cases cohort_exomes.vcf.gz \
--control gnomad_exomes_NFE.vcf.bgz \
--population NFE \
--output-dir ./burden_results \
--af-threshold 0.005 \
--revel-threshold 0.65 \
--use-cadd \
--callable-bed callable_sites_hg38.bed \
--hpo-id HP:0001250,HP:0002353 \
--dominant-onlyThree files are written to --output-dir, prefixed with GeneBurdenAnalysis_REVEL<score>_AFMax<af>:
One row per gene × model combination, sorted by BH-adjusted p-value.
HPO-matched genes appear at the top.
| Column | Description |
|---|---|
Gene |
HGNC gene symbol |
Model |
Dominant or Recessive |
Observed_Carriers |
Case samples carrying qualifying variants |
Total_Samples |
Total genotyped case samples |
Expected_Rate |
Population-derived expected carrier rate |
Het_Carriers |
Heterozygous carriers (dominant het / compound-het) |
Hom_CompHet_Carriers |
Homozygous carriers |
BinomialP |
One-sided binomial exact p-value |
BinomialAdjP |
Benjamini–Hochberg adjusted p-value |
OE_Ratio |
Observed / expected carrier rate |
CarrierRate_Lower95CI |
Clopper–Pearson 95% lower confidence bound |
CarrierRate_Upper95CI |
Clopper–Pearson 95% upper confidence bound |
| (71 annotation columns) | From the built-in gene annotation database |
One row per gene listing all qualifying variants and which case samples carry them.
Gene Variants(Carriers)
BRCA2 chr13-32316527-A-T(sample1,sample4), chr13-32340300-C-G(sample2)
Genomic-inflation QQ plot — expected vs observed −log₁₀(p-values) with:
- Separate series per model (blue = Dominant, red = Recessive)
- Reference diagonal y = x
- Genomic inflation factor λ annotated for each model
Bar chart of the 20 most significant genes (by BH-adjusted p), coloured by model, with a dashed FDR 5% significance line.
A variant is included if it passes all of:
- Mutation type is in the allowed set
- Population AF ≤
--af-thresholdin the specified population - Pathogenicity: REVEL ≥ threshold OR CADD ≥ gene-level MSC cutoff (if CADD enabled)
- Position is callable according to the BED file (if provided)
Variants missing both REVEL and CADD scores when score filtering is active are excluded.
Dominant model — probability a sample carries ≥1 alt allele across all qualifying variants:
rate = 1 − ∏ᵢ (1 − AFᵢ)²
Recessive model — probability a sample has ≥2 alt alleles (compound-het or homozygous) via per-haplotype convolution:
h0 = 1.0, h1 = 0.0
for each variant with allele frequency AF:
h0_new = h0 × (1 − AF)
h1_new = h0 × AF + h1 × (1 − AF)
h0, h1 = h0_new, h1_new
rate = 1 − h0² − 2·h0·h1
By default, expected rates use the joint case+control AF:
AF = (case_AC + control_AC) / (case_AN + control_AN)
This gives a more accurate estimate when case enrichment is mild. Disable with --no-joint-af to use control-only AF.
For each gene and model, a one-sided upper-tail binomial exact test is performed:
P(X ≥ observed_carriers | n = total_samples, p = expected_rate)
Benjamini–Hochberg FDR correction is applied independently within each model (dominant and recessive), so the two tests do not inflate each other's correction.
GeneBurdenAnalyzer/
├── src/cumulativegeneburden/
│ ├── CumulativeGeneBurden.java # Entry point, CLI parsing, GUI launcher
│ ├── AnalysisConfig.java # Immutable config (builder pattern)
│ ├── BurdenAnalysisEngine.java # Full pipeline orchestration
│ ├── filter/
│ │ ├── VariantFilter.java # Per-variant, per-gene filter logic
│ │ └── CallableBedFilter.java # BED interval filter (binary search)
│ ├── stats/
│ │ └── BurdenStatistics.java # Fisher, binomial, CI, BH correction
│ ├── model/
│ │ ├── GeneResult.java # Gene-level result DTO
│ │ ├── VariantInfo.java # Per-variant data container
│ │ └── VariantCaseInfo.java # Variant → carrier sample mapping
│ ├── io/
│ │ └── ResultWriter.java # TSV output with HPO prioritisation
│ ├── visualization/
│ │ ├── QQPlotter.java # Genomic-inflation QQ plot (PNG)
│ │ └── GeneSummaryPlotter.java # Top-genes bar chart (PNG)
│ └── gui/
│ ├── GeneBurdenApp.java # Main Swing window (5 tabs)
│ ├── ParameterPanel.java # Parameter form with file choosers
│ ├── ResultsPanel.java # Sortable/filterable JTable
│ ├── PlotPanel.java # Embedded PNG viewer
│ └── LogPanel.java # Real-time log console
├── src/cumulativegeneburden/
│ ├── CumulativeGeneBurdenDB.tsv # 71-column gene annotation database
│ ├── HPO_to_genes.txt # HPO term → gene mapping
│ └── MSC_v1.6_95.txt # Gene-level CADD MSC 95th percentile
├── lib/ # Dependency JARs (for IDE/NetBeans use)
│ ├── htsjdk-3.0.5.jar
│ ├── commons-math3-3.6.1.jar
│ ├── jfreechart-1.5.3.jar
│ └── (7 transitive dependencies)
├── dist/
│ └── GeneBurden.jar # Pre-built self-contained executable
├── pom.xml # Maven build (mvn package → fat JAR)
├── build.xml # NetBeans/Ant build (alternative)
├── CHANGELOG.md # Full version history
└── README.md
| Library | Version | Purpose |
|---|---|---|
| HTSJDK | 3.0.5 | VCF/BCF reading |
| Apache Commons Math 3 | 3.6.1 | Statistical distributions |
| JFreeChart | 1.5.3 | Chart rendering (QQ plot, bar chart) |
All dependencies are bundled in dist/GeneBurden.jar — no separate installation required.
This tool implements the core analysis strategy of CoCoRV while offering a self-contained executable with a GUI, richer pathogenicity filtering, and a built-in annotation database.
| Feature | Gene Burden Analyzer | CoCoRV |
|---|---|---|
| Dominant model | Yes | Yes |
| Recessive model | Yes | Yes |
| Expected rate formula | AF convolution | AF convolution |
| Statistical test | Binomial exact | Binomial exact |
| FDR correction | Benjamini–Hochberg | Benjamini–Hochberg |
| Callable-sites BED | Yes | Yes |
| Joint AF | Yes (default) | Yes |
| REVEL / CADD filter | Yes (both, OR logic) | REVEL only |
| Gene annotation DB | Built-in 71 columns | External |
| HPO prioritization | Yes | No |
| Visualizations | QQ plot + bar chart | QQ plot |
| GUI | Full Swing GUI | No |
| Platform | Any JVM | Python/conda |
For academic and research use within the Polychronakos Laboratory.
Please contact the authors for other usage.
If you use this tool in your research, please cite:
Eidi M. Gene Burden Analyzer. Polychronakos Laboratory, McGill University, 2025.
https://github.com/MiladEidi/GeneBurdenAnalyzer