Skip to content

Latest commit

 

History

2 Commits

Folders and files

NameName
Last commit message
Last commit date
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

CBS-Frontier

Live app: https://cbs-frontier.vercel.app/

Table of Contents


Overview

CBS-Frontier consists of two main components:

  1. Preprocessing Pipeline (Python): Converts raw genomic data (BAM, GFF3, FAI) into files optimized for visualization in a browser
  2. Web Application (React): Interactive browser-based visualization of k-mer density across the genome

Key Features

  • Multi-resolution genome browsing (1kb, 10kb, 50kb)
  • MAPQ-based quality filtering
  • Primary vs secondary alignment filtering
  • Gene and repeat annotation overlays
  • Searchable, sortable feature lists
  • Side-by-side chromosome comparison
  • Changepoint detection (Binary Segmentation) to identify density-transition boundaries
  • PNG and CSV export

System Requirements

Preprocessing Pipeline

The pipeline shells out to bedtools and samtools; it needs no third-party Python packages.

Requirement Version
Python 3.8+ (standard library only)
bedtools 2.30+ (system install)
samtools 1.15+ (system install)
AGAT optional, only for --standardize-genes (conda install -c bioconda agat)

Web Application:

This application has only been tested in Google Chrome. Other browsers (Firefox, Safari, Edge) may work but are not officially supported.


Preparing Your Data

Before running the preprocessing pipeline, you need:

Required Files

File Description Format
Genome index FASTA index file .fai
K-mer alignments K-mers aligned to genome .bam (indexed)

Important: BAM files must be indexed before preprocessing. If you don't have a .bai index file.

Optional Files

File Description Format Flag
Gene annotations Gene features GFF3 (RAGNAROK output) --genes
Repeat annotations Repeat/TE features GFF3 (EDTA v2.2.2 output) --repeats
Custom track Any other feature set — RepeatMasker output, curated regions, QTL intervals BED6 --bed + --bed-name

Forming a BED6 file for --bed

Anything the native GFF3 paths can't ingest can be added as a custom BED6 track. It is tab-delimited with six columns and 0-based, half-open coordinates:

chrom    start    end    name    score    strand
Chr01    15283    16209  myFeature   .        +
Col Field Notes
1 chrom Must match your .fai contig names exactlychr1 vs 1 is the most common failure
2 start 0-based, integer, start < end
3 end half-open, integer
4 name the searchable label in the web app; blank names fall back to chrom:start-end
5 score unused — . is fine
6 strand +, -, or .

Note on GFF3 Compatibility: The current preprocessing pipeline is optimized for specific GFF3 formats:

  • Gene annotations: Optimized for RAGNAROK output format
  • Repeat annotations: The native repeat path currently supports EDTA v2.2.2 GFF3 output only.

For gene GFF3 files from other annotation tools, pass --standardize-genes to normalize the file with AGAT first — it rebuilds the gene→mRNA→exon/CDS hierarchy so non-standard files parse correctly. Repeat GFF3 files have no equivalent normalization step; convert non-EDTA repeat annotations to BED6 and load them via --bed instead (see below).

Using RepeatMasker .out files: RepeatMasker output is not accepted by the native repeat path. To visualize it, convert it to BED6 (e.g. with rmsk2bed from BEDOPS, then cut -f1-6) and load it as an optional BED track. You must choose what goes in column 4 — either the repeat name or the class/family — since a BED6 name column holds only one value (unlike the native EDTA path, which tracks both). Ensure the BED chromosome names match your reference .fai/BAM exactly.


Preprocessing Pipeline

The preprocessing script generates all files needed for the web application.

Basic Usage

cd preprocessing

python preprocess_genome.py \
    --fai /path/to/genome.fa.fai \
    --kmers /path/to/kmers.bam \
    --output-dir /path/to/output/

Full Usage with Annotations

python preprocess_genome.py \
    --fai /path/to/genome.fa.fai \
    --kmers /path/to/kmers.bam \
    --genes /path/to/genes.gff3 \
    --repeats /path/to/repeats.gff3 \
    --output-dir /path/to/output/ \
    --sparse \
    --mapq-bins "0-19,20-39,40-60"

Recommended: genes + repeats, lowest runtime

The best balance of coverage and speed for most projects — full gene and repeat context, without the per-MAPQ-score track sets:

python preprocess_genome.py \
    --fai genome.fa.fai \
    --kmers kmers.bam \
    --genes genes.gff3 \
    --repeats repeats.gff3 \
    --output-dir output/ \
    --skip-mapping-scores

With a custom BED6 track

Adds one custom annotation track alongside the native gene and repeat paths.

python preprocess_genome.py \
    --fai genome.fa.fai \
    --kmers kmers.bam \
    --genes genes.gff3 \
    --repeats repeats.gff3 \
    --bed repeatmasker.bed \
    --bed-name "RepeatMasker" \
    --output-dir output/ \
    --skip-mapping-scores

Command Line Options

Option Required Default Description
--fai Yes - Genome FASTA index file (.fai)
--kmers Yes - K-mer alignments in BAM format (must be indexed)
--output-dir Yes - Directory for output files
--genes No - Gene annotations in GFF3 format. GFF3 is strongly preferred over GTF — GTF from some tools (e.g. BRAKER) yields duplicate gene models
--repeats No - Repeat annotations in GFF3 format (EDTA v2.2.2)
--standardize-genes No false Normalize --genes with AGAT (agat_convert_sp_gxf2gxf.pl) before processing, rebuilding the gene→mRNA→exon/CDS hierarchy so non-standard GFF3 parses correctly. Requires AGAT on PATH
--agat-path No agat_convert_sp_gxf2gxf.pl Path or command for the AGAT conversion script, if it isn't on PATH
--bed No - Optional custom BED6 track (see Optional BED6 Track). Requires --bed-name
--bed-name No - Display name for the custom BED6 track (required with --bed)
--size-threshold No 10,000,000 Min chromosome size for genome-wide view
--min-chrom-size No 5,000,000 Min chromosome size to include in analysis
--min-kmer-count No 50 Min k-mer count for features to appear in the web app's searchable lists (lower includes more low-signal features; 0 includes all). Does not affect all_kmer_counts/, which is always complete
--sparse true|false No true Compressed sparse gzipped-JSON output (default; ~84% smaller, much faster to load. Set --sparse false for plain .bedgraph (interoperable with IGV/bedtools, but far larger)
--low-memory No false Process chromosomes one at a time (slower, uses less RAM)
--mapq-bins No - Group MAPQ scores into bins (e.g., "0-19,20-39,40-60")
--skip-mapping-scores No false Skip per-MAPQ-score coverage files

How do I reduce the preprocessing runtime?

  • --skip-mapping-scores. The per-score tracks are the single biggest runtime cost. The cost scales with the number of distinct scores, so on a typical BWA BAM (dozens of values) this dominates everything else. Skip it unless you specifically need MAPQ-stratified tracks; if you do need them, use --mapq-bins "0-19,20-39,40-60" rather than the default per-score behaviour, which collapses dozens of file sets into three.
  • Leave --low-memory off unless you actually run out of RAM. It re-runs samtools index and samtools view once per chromosome per track, so on a 16-chromosome assembly it multiplies the samtools work roughly 16-fold.
  • Keep sparse output on (the default).

Optional BED6 Track

In addition to the native gene (GFF3) and repeat (EDTA GFF3) paths, you can supply one custom track as a plain BED6 file with --bed. This is the way to visualize annotation sets the native paths can't ingest — most notably RepeatMasker / RepeatModeler output, or any curated / subset region list.

BED6 specification — tab-delimited, exactly (at least) 6 columns, coordinates 0-based, half-open:

Col Field Notes
1 chrom Must match the genome .fai contig names exactly (e.g. chr1 vs 1) — the #1 failure mode
2 start 0-based, integer, start < end
3 end half-open, integer
4 name feature label (used as the searchable index key; empty names fall back to chrom:start-end)
5 score any value / . (unused)
6 strand +, -, or .

This produces custom_kmer_density_{1kb,10kb,50kb}, custom_density_{1kb,10kb,50kb}, custom_index.json, and the custom_track.json manifest.

Complete k-mer count tables (all_kmer_counts/)

The gene_index.json and repeat_index.json files the web app loads are filtered by --min-kmer-count (default 50), so they are not a complete record of all overlaps. Every run therefore also writes an all_kmer_counts/ subdirectory containing the same information unfiltered, as plain TSVs for downstream analysis in R, pandas, or Excel: gene_kmer_counts.tsv (one row per gene), repeat_family_kmer_counts.tsv and repeat_locus_kmer_counts.tsv (one row per TE family and per individual genomic copy, since the index only records family-wide totals), and custom_kmer_counts.tsv when --bed is given. Each row carries the feature's chromosome, coordinates, strand, k-mer count, length, and density, sorted with the highest counts first. Counts are distinct k-mers rather than overlap events — a k-mer spanning two records of the same feature is counted once, which is also why a repeat family's total can be less than the sum of its loci.


Web Application

Loading Data

  1. Click the file upload area
  2. Select all files from your preprocessing output directory
  3. Wait for files to load (progress bar shows status)

Required Files for Web App

File Description
metadata.json Genome information
kmer_density_1kb.bedgraph* 1kb resolution k-mer coverage
kmer_density_10kb.bedgraph* 10kb resolution k-mer coverage
kmer_density_50kb.bedgraph* 50kb resolution k-mer coverage

Optional Files

Gene group — required together if gene_index.json is present:

File Description
gene_index.json Gene feature index
gene_kmer_density_{1kb,10kb,50kb}.bedgraph* K-mers overlapping genes
gene_density_{1kb,10kb,50kb}.bedgraph* Coverage of the gene annotations themselves

Repeat group — required together if repeat_index.json is present:

File Description
repeat_index.json Repeat feature index
repeat_classifications.json Repeat superfamily breakdown
repeat_kmer_density_{1kb,10kb,50kb}.bedgraph* K-mers overlapping repeats
repeat_density_{1kb,10kb,50kb}.bedgraph* Coverage of the repeat annotations themselves

Custom BED6 group — gated on the manifest:

File Description
custom_track.json Custom BED6 track manifest (name + gate — see note below)
custom_index.json Custom BED6 feature index
custom_kmer_density_*.bedgraph* K-mers overlapping the custom intervals
custom_density_*.bedgraph* Coverage of the custom intervals themselves

Independently optional — each can be uploaded or omitted on its own:

File Description
mapping_scores.json Available MAPQ scores
alignment_flags.json Alignment flag groups
*_score*.bedgraph* Per-MAPQ coverage files
*_primary.bedgraph* Primary alignment coverage
*_secondary.bedgraph* Secondary alignment coverage

Navigation

  • Genome View: Click on a chromosome to zoom in
  • Chromosome View: Click and drag to select a region
  • Region View: Use zoom/pan buttons or keyboard shortcuts
    • Arrow keys: Pan left/right
    • +/-: Zoom in/out

Filtering

  • K-mer Filter: Show all k-mers, gene-associated, repeat-associated, or specific feature
  • Mapping Score: Select a MAPQ score (or bin, if --mapq-bins was used)
  • Alignment Flag: Show primary or secondary alignments

Changepoint Analysis

Detect statistically distinct regions of k-mer density along a chromosome using Binary Segmentation. This is useful for locating boundaries such as sex-determining regions or other density transitions.

Open the Changepoint Analysis panel in the view controls (available at 1 kb resolution in chromosome or region view) and set:

Parameter Description
Breakpoints Number of change points to detect (nBkps)
Min segment size Minimum segment length in bins (minSize)
Cost model L2 (detects shifts in mean) or Normal (detects shifts in mean and variance)

Click Run Changepoint to compute results (Clear removes them). The results panel reports:

  • Segments: start/end coordinates, size, mean, median, max density, and bin count
  • Pairwise significance: Mann-Whitney U test between adjacent segments to validate each boundary

Exporting

Available from the Export menu in the navigation bar:

Export Output
PNG Current view as a high-resolution image (10x scale)
JPG Same capture at 10x scale, JPEG quality 0.95
PDF Current view scaled to fit a single A4 page (4x capture)
Genes CSV Genes visible in the current view
Repeats CSV Repeats visible in the current view
Custom CSV Custom BED6 features visible in the current view (only with a custom track loaded)

CSV exports include feature name, coordinates, length, classification, k-mer count and k-mer density, filtered to whatever the current view shows.

Attribution & Development

CBS-Frontier was developed by Laramie Akozbek.

This project was developed with the assistance of Anthropic's Claude models (Opus and Sonnet). The models assisted with the design, implementation, and refinement of both the Python preprocessing pipeline and the React web application. All code was directed, reviewed, and integrated by the author.

Funding

This work was funded by National Science Foundation IOS-PGRP CAREER #2239530 (Alex Harkess) and by the National Science Foundation Graduate Research Fellowship Program (Laramie Akozbek).


Citation

If you use CBS-Frontier in your work, please cite the software via CITATION.cff — GitHub's "Cite this repository" button will generate BibTeX or APA from it.


License

Released under the Apache License 2.0.

Copyright © 2025-2026 HudsonAlpha Institute for Biotechnology.

About

Genome-wide k-mer density browser: Python preprocessing pipeline + React visualization web app

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Contributors

Languages