Setting Up a Genome Browser for Your Own DNA

By the end of this you will have a browser running on your own machine that renders your alignments, your variant calls, and your RNA-seq coverage on GRCh38, alongside gene annotation, with everything indexed so that jumping to a locus takes milliseconds instead of minutes. You will also have a scripted path to the same data, because most of the questions you will want to ask are better answered by a one-line bcftools or tabix query than by squinting at a track. You need the following.

• a CRAM or BAM of your WGS
• a VCF (ideally a gVCF so you can distinguish “reference” from “not sequenced”)
• a BAM or quantification from RNA-seq
• the exact reference FASTA your data was aligned to
• roughly 150 GB of free disk if you are keeping a CRAM plus indices plus derived tracks

The tools are listed here.

• samtools
• bcftools
• htslib (bgzip, tabix)
• bedtools
• deeptools
• IGV
• optionally the JBrowse 2 CLI

1. Pin down your reference build before anything else

Almost every failure in this pipeline traces back to a reference mismatch. Find out which FASTA your CRAM was aligned to and keep that exact file. Check the header:

samtools view -H sample.cram | grep '^@SQ' | head -3
# @SQ SN:chr1 LN:248956422 M5:2648ae1bacce4ec4b6cf337dcae37816 UR:...

Three things matter. The build (chr1 length 248956422 is GRCh38, and 249250621 is GRCh37). The contig naming (chr1 in UCSC-style files, 1 in Ensembl/GRC-style files). The M5 checksum, which identifies the exact FASTA, including whether it is the analysis set with decoy and HLA contigs.

Our recommendation: use GRCh38, specifically the “no-alt analysis set” plus decoys, with chr-prefixed names. It is what most callers and most public track files assume, and keeping ALT contigs out of the primary alignment target avoids a class of multi-mapping artifacts that show up as phantom heterozygous calls in segmental duplications.

If your data is on GRCh37, do not mix builds in one browser session. Convert with CrossMap or bcftools +liftover and accept that a small fraction of positions will not map cleanly. Assembly-to-assembly comparison of the human genome identified thousands of structural differences between independently produced assemblies, which is the reason coordinate conversion is lossy rather than mechanical.¹

Index the reference once:

samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa

2. Index every file you plan to browse

A genome browser is a random-access client. Without indices it either refuses to load or reads the whole file.

# alignments
samtools index -@ 8 sample.cram            # produces sample.cram.crai
samtools index -@ 8 rna.bam                # rna.bam.bai

# variants: must be bgzip-compressed, coordinate-sorted
bcftools sort -Oz -o sample.vcf.gz sample.vcf
bcftools index -t sample.vcf.gz            # .tbi

Two gotchas. gzip is not bgzip, and a file compressed with plain gzip will fail tabix with not a BGZF file. And for files with contigs longer than 512 Mb, you need CSI indices: bcftools index -c, samtools index -c. That is relevant if you ever browse plant or amphibian genomes. Human contigs do not reach that length.

Verify the index works before opening a browser:

tabix sample.vcf.gz chr7:117,480,000-117,490,000 | head
samtools view -c sample.cram chr7:117,480,000-117,490,000

If samtools view on a CRAM errors with a reference lookup failure, pass the FASTA explicitly: samtools view -T GRCh38...fa -c sample.cram chr7:.... CRAM stores sequence as a diff against the reference, so the reference is not optional. Set REF_PATH or REF_CACHE if you want to avoid the flag every time.

3. Turn RNA-seq into a coverage track

Loading a 40 GB RNA BAM into a browser works, but rendering per-base coverage genome-wide is slow. Precompute a bigWig:

bamCoverage -b rna.bam -o rna.bw \
  --binSize 10 \
  --normalizeUsing BPM \
  --numberOfProcessors 8 \
  --ignoreDuplicates

BPM (bins per million) makes tracks from different timepoints comparable on the same y-axis. Keep the BAM around too: when a gene looks interesting you will want to see splice junctions and read-level evidence, and the bigWig throws that away. In IGV, load the bigWig as the always-on track and the BAM only when you zoom into a gene.

For quantification, a per-transcript TSV from salmon or kallisto is more useful than anything the browser shows. Browse to confirm what the numbers imply (is the signal at the annotated exons? is one isoform dominant?), and compute in the table.

4. Understand BED, because everything you build will be BED

A BED file is tab-separated plain text. Minimum three columns:

chr7	117479963	117668665
chr17	43044294	43170327

The full 12-column form:

chr7	117479963	117668665	CFTR	0	+	117480095	117667106	0,0,255	27	185,111,109,...	0,24672,25049,...

The columns come in this order.

• chrom
• chromStart
• chromEnd
• name
• score (0–1000)
• strand
• thickStart
• thickEnd
• itemRgb
• blockCount
• blockSizes
• blockStarts

Blocks are exons, and block starts are relative to chromStart.

The single most important fact is that BED is 0-based and half-open. The interval chr7 100 105 covers bases 101–105 in 1-based terms, five bases. VCF, GFF/GTF, SAM, and the coordinate box in every genome browser are 1-based, fully closed. If your annotations are consistently off by one at the left edge, this is why. Convert a VCF SNV position p to BED with start = p-1, end = p.

Sort, compress, and index BED for browsing:

sort -k1,1 -k2,2n regions.bed > regions.sorted.bed
bgzip regions.sorted.bed
tabix -p bed regions.sorted.bed.gz

For tracks with many thousands of features that you want to pan across smoothly, convert to bigBed:

cut -f1,2 GRCh38...fa.fai > hg38.chrom.sizes
bedToBigBed regions.sorted.bed hg38.chrom.sizes regions.bb

And the reason BED matters day to day: it is the interchange format for set operations on your own genome.

# which of my rare variants fall in exons of MANE Select transcripts?
bedtools intersect -a rare.bed -b mane_exons.bed -wa -wb -sorted -g hg38.chrom.sizes

-sorted with -g streams instead of loading both files into memory, which matters once you are intersecting 4 million variants against an annotation set.

5. Run IGV locally and drive it with batch scripts

IGV desktop is the right tool for read-level inspection. Install it, then edit the launcher to give the JVM room: -Xmx8g is a reasonable floor for CRAM plus a couple of tracks.

Preferences worth changing on day one:

• Alignments → visibility range threshold: raise from 30 kb to 100 kb if you have the RAM, otherwise leave it.
• Alignments → downsample reads: keep on (default 100 reads per 50-base window). Turn it off only when you are counting reads by eye at one locus, and remember to turn it back on.
• Alignments → show soft-clipped bases: turn on. Soft clips are the visual signature of structural variant breakpoints and of adapter contamination.

Then stop clicking. IGV takes batch scripts:

new
genome /data/ref/GRCh38_full_analysis_set_plus_decoy_hla.fa
load /data/sample.cram
load /data/sample.vcf.gz
load /data/rna.bw
snapshotDirectory /data/snapshots
goto chr7:117,559,590-117,559,610
sort base
snapshot CFTR_F508.png
goto chr16:88,713,000-88,717,000
snapshot region2.png
exit

igv.sh --batch inspect.txt

This is how you review 200 candidate loci in one pass rather than 200 manual navigations, and it gives you a reproducible artifact.

For an embedded, shareable version, igv.js takes the same index files over HTTP range requests. Your web server must support Range headers and CORS, or tracks silently fail to render.

6. Use JBrowse 2 when you want a persistent instance

If you want a browser you can leave running, with your tracks always loaded and reachable from a laptop, JBrowse 2 is the better choice. It is a static site plus a JSON config.

npm install -g @jbrowse/cli
jbrowse create ~/jb2
cd ~/jb2

jbrowse add-assembly /data/ref/GRCh38...fa --load symlink --name hg38

jbrowse add-track /data/sample.cram --load symlink \
  --name "WGS alignments" --assemblyNames hg38

jbrowse add-track /data/sample.vcf.gz --load symlink \
  --name "Germline variants" --assemblyNames hg38

jbrowse add-track /data/rna.bw --load symlink \
  --name "RNA coverage (BPM)" --assemblyNames hg38

# gene annotation: sorted, bgzipped, tabixed GFF3
jbrowse add-track /data/anno/gencode.v44.sorted.gff3.gz --load symlink \
  --name "GENCODE v44" --assemblyNames hg38

jbrowse text-index --tracks GENCODE_v44 --assemblies hg38

npx serve -l 8080 .

jbrowse text-index is what makes typing CFTR in the location box work. Without it you can only navigate by coordinates.

Keep it bound to localhost or behind auth. A VCF of your own genome is identifying, and once it is on a public URL you cannot take it back.

7. Query annotation programmatically instead of browsing for it

Most “look this up in a browser” tasks are faster as an API call or a local join.

Ensembl is free and open access, and the underlying code and database schema are published, which is why third-party tools can talk to it directly.² The core software distributes DNA sequence and genome annotation through a stable schema with Perl and REST access layers, and the same release cycle feeds the web browser and the API.² The 2025 release covers a large set of vertebrate genomes with regularly refreshed gene sets and variation data.³

A single-variant consequence lookup:

curl -s 'https://rest.ensembl.org/vep/human/region/7:117559590-117559592:1/A?content-type=application/json' \
  | jq '.[0].transcript_consequences[] | select(.canonical==1) | {gene_symbol, consequence_terms, hgvsc}'

Ensembl’s variation resources aggregate variant calls, frequencies, phenotype associations, and predicted consequences, which is what the VEP endpoint is reading from.⁴ For anything larger than a few hundred queries, stop hitting the REST API and use offline VEP with a cache, or the public MySQL mirror:

mysql -h ensembldb.ensembl.org -u anonymous -P 3306 -e \
  "SELECT stable_id, seq_region_start, seq_region_end FROM gene LIMIT 5" \
  homo_sapiens_core_112_38

Annotating a whole-genome VCF against the web API takes days and will get you rate-limited. Offline VEP with --cache --offline --fork 8 --assembly GRCh38 finishes a 4-million-variant VCF in under an hour on a laptop.

Why the Ensembl web browser feels slow: each region view assembles gene models, transcript support levels, regulatory features, comparative alignments, and variation tracks from a set of MySQL databases and renders them server-side. Ask for a 5 Mb window with variation tracks on and you are asking for millions of rows. The fixes are listed here in order of effectiveness. Narrow the window. Turn off the variation and comparative genomics tracks in the track configuration panel. Use the REST endpoint for anything you can express as a query, or download the GTF and query it locally. The project has been structured around downloadable data and programmatic access since its early releases for exactly this reason.⁵

8. Ask questions of your own data

Once the browser and the indices are in place, the workflow runs as follows. Form a hypothesis in the shell, confirm it visually, and go back to the shell.

# depth and quality at a specific position
samtools mpileup -f ref.fa -r chr7:117559590-117559590 sample.cram

# all variants in one gene, with genotype and depth
bcftools query -r chr7:117479963-117668665 \
  -f '%CHROM\t%POS\t%REF\t%ALT\t[%GT\t%DP\t%GQ]\n' sample.vcf.gz

# coverage summary over a region
samtools coverage -r chr17:43044294-43170327 sample.cram

Use a gVCF for the second query. A plain VCF that returns nothing for a region is ambiguous between “homozygous reference” and “zero reads covered this”. A gVCF distinguishes them, and the difference matters when you are checking whether an absent call is real.

Common problems

Contig naming mismatch. Your CRAM says chr1, your annotation BED says 1, and every intersect returns empty. Fix the annotation, not the alignment: bcftools annotate --rename-chrs map.txt, or awk 'BEGIN{OFS="\t"}{$1="chr"$1; print}' for BED. Check with cut -f1 file.bed | sort -u | head.

Off-by-one at interval starts. BED is 0-based half-open, GFF/VCF are 1-based closed. When a variant at a splice acceptor looks like it falls one base outside the exon, verify the conversion before you believe the biology.

CRAM fails to open. Missing or wrong reference. Pass -T, or set REF_CACHE=/path/cache/%2s/%2s/%s. If the M5 checksums in the header do not match your FASTA, you have the wrong FASTA and decoding will produce wrong bases. Decoding will not raise an error.

Apparent heteroplasmy in mitochondrial calls. Segments of mtDNA are present in the nuclear genome as nuclear mitochondrial sequences (NUMTs), and nuclear reads mismapping to chrM produce low-frequency variant calls that look like heteroplasmy. Annotated NUMT tracks exist for the UCSC browser and are worth loading before you interpret anything on chrM.⁶ Filter chrM calls by mapping quality and check whether the region overlaps a known NUMT.

IGV runs out of memory or hangs. You loaded a whole-genome BAM and navigated to a whole chromosome. Raise -Xmx, keep downsampling on, and use bigWig for anything you view above 100 kb.

A .dna file will not open. That extension is SnapGene, a plasmid or construct map, unrelated to human genome data. Open it with the free SnapGene Viewer, or parse it in Python: from Bio import SeqIO; rec = next(SeqIO.parse("construct.dna", "snapgene")), then write FASTA or GenBank. It will not contain your genome.

“Humans are 99% identical” does not mean your browser will look boring. That figure comes from single-nucleotide differences and understates the rest. Independent assemblies of the human genome differ by thousands of insertions, deletions, and inversions, many of them kilobases long, in regions that short-read alignment handles poorly.¹ The structural variants and copy-number differences are where a browser earns its keep, because they are visible as coverage steps, soft-clip clusters, and discordant pairs long before any caller labels them.

Interpretation. A browser shows you evidence. The tracks do not show you meaning. A variant that looks striking in IGV is a hypothesis, and the step from “this call is well supported” to “this matters for my health” is a clinical one. Take it to a physician or a certified genetic counselor, with confirmatory testing in a clinical lab, before acting on anything.

Footnotes

1. Razi Khaja, Junjun Zhang, Jeffrey R. MacDonald, et al. Genome assembly comparison identifies structural variants in the human genome. Nature Genetics, 2006. https://doi.org/10.1038/ng1921 ↩ ↩²

2. Magali Ruffier, Andreas Kähäri, Monika Komorowska, et al. Ensembl core software resources: storage and programmatic access for DNA sequence and genome annotation. Database, 2017. https://doi.org/10.1093/database/bax020 ↩ ↩²

3. Sarah Dyer, Olanrewaju Austine-Orimoloye, Andrey G Azov, et al. Ensembl 2025. Nucleic Acids Research, 2024. https://doi.org/10.1093/nar/gkae1071 ↩

4. Sarah Hunt, William McLaren, Laurent Gil, et al. Ensembl variation resources. Database, 2018. https://doi.org/10.1093/database/bay119 ↩

5. Tim Hubbard. The Ensembl genome database project. Nucleic Acids Research, 2002. https://doi.org/10.1093/nar/30.1.38 ↩

6. Francesco Maria Calabrese, Domenico Simone, Marcella Attimonelli. Primates and mouse NumtS in the UCSC Genome Browser. BMC Bioinformatics, 2012. https://doi.org/10.1186/1471-2105-13-s4-s15 ↩