How to Set Up a Genome View of Your Own Sequencing Data

By the end of this you will have a local genome browser pointed at your own data. It will show aligned reads from whole-genome sequencing, your variant calls layered on top, and RNA-seq junctions in a second track. Public annotation sits underneath, namely RefSeq genes, ClinVar, gnomAD allele frequencies and ENCODE candidate regulatory elements. You will also have a scripted path that produces PNG screenshots of any region you name, so you can put a picture of a locus into a notebook without clicking. You need the following.

• your aligned reads (.bam or .cram)
• the matching index
• a VCF
• the exact reference FASTA your files were aligned to
• roughly 100 GB of free disk if you keep a CRAM plus a decompressed reference
• 16 GB of RAM (8 GB will work if you are careful with track count)
• samtools, bcftools, and htslib from bioconda

A genome view is a coordinate system with tracks stacked on it. Almost every problem you will hit is a coordinate system mismatch.

1. Inventory what you have and confirm the reference build

Before installing anything, find out which genome your files live on. A BAM header carries the contig names and lengths used at alignment time:

samtools view -H sample.bam | grep '^@SQ' | head -3
samtools view -H sample.bam | grep '^@PG'

Look at two things. First, contig naming: SN:chr1 (UCSC style) or SN:1 (Ensembl/GRC style). Second, the length of chromosome 1. LN:248956422 means GRCh38. LN:249250621 means GRCh37/hg19. If your provider handed you a VCF from one build and a report with coordinates from another, you will silently look at the wrong base. GRCh37 and GRCh38 coordinates for the same variant can differ by thousands to millions of bases depending on the chromosome.

If you have a CRAM instead of a BAM, the header also records the reference checksums (M5: on each @SQ line). CRAM is reference-compressed: the file does not store bases that match the reference, so without the exact FASTA it cannot be decoded at all. Keep that FASTA. A 30x WGS CRAM is roughly 15-20 GB against 60-90 GB for the equivalent BAM, which is why most providers ship it now, and why “I lost the reference” is a genuinely destroyed dataset.

Verify file integrity before anything else:

samtools quickcheck -v sample.bam || echo "TRUNCATED"
md5sum -c checksums.md5

quickcheck reads the EOF block. A truncated download looks completely normal in a file listing and then produces reads that stop at chr14.

2. Pin one reference FASTA and index it

Use the GRCh38 no-alt analysis set unless your data says otherwise. The alt contigs in the full GRCh38 primary assembly steal reads from their primary counterparts and drop mapping quality to zero in those regions, which shows up in IGV as a white (MAPQ 0) pileup in genes like the HLA cluster. The no-alt analysis set from the 1000 Genomes FTP mirror is the de facto standard for short-read human alignment.

# index for random access
samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa
# confirm it matches your BAM's contig names and lengths
cut -f1,2 GRCh38_full_analysis_set_plus_decoy_hla.fa.fai | head -3

If your BAM uses 1 and your annotation files use chr1, fix the annotation. Do not change the BAM. Rewriting a BAM header is possible with samtools reheader but you will regret having two BAMs that differ only in naming. For VCFs:

bcftools annotate --rename-chrs chr_map.txt in.vcf.gz -Oz -o out.vcf.gz

where chr_map.txt is two whitespace-separated columns, old name and new name.

3. Index everything, in the right order

IGV does random access by byte range. Every track needs an index, and an index older than its data file produces reads at the wrong offsets, which manifests as garbage or nothing.

samtools index -@ 4 sample.bam            # produces sample.bam.bai
samtools index -@ 4 sample.cram           # produces sample.cram.crai

bgzip -@ 4 variants.vcf                   # NOT gzip
tabix -p vcf variants.vcf.gz              # produces .tbi

# for genomes with contigs > 512 Mb, or if tabix complains:
bcftools index --csi variants.vcf.gz

gzip will not work. Tabix requires BGZF, a gzip-compatible container written in independently decompressible ~64 KB blocks. A plain gzip’d VCF loads in zcat and fails in every genome browser with a message about an invalid index.

For BED, GFF3, and interval files, sort by coordinate first or tabix will error on an out-of-order line:

sort -k1,1 -k2,2n annotations.bed | bgzip > annotations.bed.gz
tabix -p bed annotations.bed.gz

4. Install IGV desktop and raise the heap

Download IGV from the Broad’s site (the “with Java included” bundle avoids Java version arguments). Then edit the memory setting before you open a whole genome. On macOS, the setting is in igv.args inside the app bundle. On Linux, it is in the igv.sh launcher.

-Xmx8g
-Dapple.laf.useScreenMenuBar=true

8 GB of heap handles a 30x WGS BAM plus a dozen annotation tracks. IGV is not streaming everything into memory. The alignment track does cache the visible window plus a margin, and sashimi plots on deep RNA-seq are memory hungry.

Then set preferences that matter (View > Preferences > Alignments):

• Visibility range threshold: 30 kb. Above this, IGV will not render reads. Raising it to 300 kb on a 30x genome will hang the UI.
• Downsample reads: on, 100 reads per 50 base window. Leave it on for browsing. Turn it off when you are judging allele fraction at a specific site, because downsampling is not allele-aware and will distort the ratio you see.
• Show soft-clipped bases: on. Soft clips are the single most informative visual cue for structural variants and for mismapping.
• Mapping quality threshold: 0 for exploration. Set it to 20 when you want to know whether a variant survives mapping-quality filtering.

5. Load your data and look at one known locus

File > Load from File for the BAM and the VCF. IGV follows the .bai next to the .bam automatically. For CRAM, set Genomes > Load Genome from File to your exact FASTA first, otherwise decoding fails with a cryptic reference mismatch error.

Type a gene name in the locus box (MTHFR, HFE, CYP2C19) and hit enter. IGV resolves gene symbols against its own annotation for the loaded genome. Now set up the alignment track the way it is useful:

• Right-click the track > Color alignments by > insert size and pair orientation. Normal pairs are gray. Red means insert size too large, blue too small, and the pair-orientation colors (green, teal, dark blue) flag inversions and duplications.
• Group alignments by > phase (HP tag) if you have phased reads. Long-read or linked-read data will have HP:i:1 / HP:i:2 and the two haplotypes split into separate stacks, which makes compound heterozygosity visible at a glance.
• View as pairs, on, when you are looking at a candidate structural variant.

For a heterozygous SNV you should see roughly half the reads carrying the alternate base, colored by base, with reads on both strands and both mate orientations. What you are doing here is visual quality control on a caller’s output, and it remains the fastest way to spot a false positive. Large-scale sequencing studies have shown that a substantial fraction of the interpretive burden in human genome analysis sits in exactly this gap between a call being emitted and a call being trustworthy ¹.

6. Add annotation tracks

Gene models alone will not tell you much. Add, in this order:

1. RefSeq genes (loaded by default with the IGV hg38 genome).
2. ClinVar VCF from NCBI’s FTP, bgzip’d and tabix’d. Note that ClinVar’s VCF uses 1, not chr1, so rename it as in step 2 if your genome is UCSC-styled.
3. gnomAD allele frequencies. The full genomes VCF is enormous, so pull only the INFO fields you want. Strip everything except AF, AF_popmax and nhomalt with bcftools annotate -x, and the file drops by an order of magnitude.
4. ENCODE candidate cis-regulatory elements (cCREs) as a BED from the SCREEN portal. This is what turns a variant that sits in “nothing” into a variant that sits in a distal enhancer-like signature. The pilot phase of ENCODE established that the majority of the genome is transcribed or associated with a regulatory feature, and that functional elements are not confined to protein-coding exons ².

# trim gnomAD to just the fields you will read
bcftools annotate \
  -x '^INFO/AF,INFO/AF_popmax,INFO/nhomalt' \
  gnomad.genomes.v4.1.sites.chr1.vcf.bgz \
  -Oz -o gnomad.chr1.slim.vcf.gz
bcftools index --tbi gnomad.chr1.slim.vcf.gz

Save the whole arrangement: File > Save Session. The .xml session records track order, colors, thresholds, and file paths. Check it into a repo. It is the reproducible artifact. The screenshot does not play that role.

7. Layer RNA-seq on top of DNA

If you have RNA-seq aligned with STAR, load the coordinate-sorted BAM as a second track and the SJ.out.tab splice junctions as a third. In IGV, right-click the RNA track and choose Sashimi Plot. You get exon coverage with arcs between junctions, labeled with the number of reads supporting each junction.

This is where a genome view earns its keep. A synonymous or intronic variant near a splice site is uninterpretable from DNA alone. With RNA on the same coordinates you can see whether an alternate junction is being used and roughly how often. You can also see allele-specific expression: navigate to a heterozygous coding SNV, and if the DNA track shows a 50/50 split while the RNA track shows 90/10, one allele is being expressed far less than the other.

Two caveats. RNA-seq coverage is not uniform, so low coverage at a site means low expression in that tissue and sample. Low coverage does not mean absence. And a single blood RNA sample tells you about blood.

8. Script the view so you stop clicking

IGV desktop listens on TCP port 60151 when Enable port is set in preferences. Batch scripts drive it:

new
genome hg38
load /data/sample.bam
load /data/variants.vcf.gz
load /data/gnomad.chr1.slim.vcf.gz
snapshotDirectory /data/snapshots
maxPanelHeight 1200
goto chr1:11,796,321-11,796,421
sort base
squish
snapshot MTHFR_c665.png
exit

Run headless or against a running instance:

igv.sh --batch loci.txt
# or, against a live IGV:
echo -e "goto chr7:117559590\nsnapshot cftr.png" | nc 127.0.0.1 60151

Generate the batch file from a variant list with awk and you have a hundred locus screenshots in a few minutes. Pad each region by 100 bp on each side, because a screenshot with the variant at the exact center edge is useless for judging context.

9. Put it on a web page with igv.js

For anything you want to share, or embed in a notebook, use igv.js. It is the same rendering engine in JavaScript, reading BAM/CRAM/VCF over HTTP range requests. No server-side code is required beyond a static file host that honors Range headers (S3 and most CDNs do).

<div id="igv"></div>
<script src="https://cdn.jsdelivr.net/npm/igv@3/dist/igv.min.js"></script>
<script>
igv.createBrowser(document.getElementById("igv"), {
  genome: "hg38",
  locus: "chr1:11,796,321-11,796,421",
  tracks: [
    { type: "alignment",
      format: "cram",
      url: "https://data.example.org/sample.cram",
      indexURL: "https://data.example.org/sample.cram.crai",
      name: "WGS 30x",
      colorBy: "unexpectedPair",
      visibilityWindow: 30000 },
    { type: "variant",
      format: "vcf",
      url: "https://data.example.org/variants.vcf.gz",
      indexURL: "https://data.example.org/variants.vcf.gz.tbi",
      name: "Calls" }
  ]
});
</script>

Two operational notes. CORS: the host must return Access-Control-Allow-Origin and Access-Control-Expose-Headers: Content-Length, Content-Range, or you get an empty track with no error in the UI. And if the data is yours, do not put it on a public bucket. Use presigned URLs with a short expiry, or run igv.js against localhost with a two-line Python static server. Your genome is the most re-identifiable file you own.

Common problems

“No reads shown” but the file is fine. There are three usual causes. You may be zoomed out past the visibility range threshold, where IGV prints a small message that is easy to miss. The contig naming may not match the loaded genome (chr1 vs 1). Or every read in the window may have MAPQ below your filter. Zoom to under 30 kb and drop the MAPQ threshold to 0 to distinguish them.

Index is stale. IGV does not check timestamps. If you regenerated a BAM and kept the old .bai, you get reads from wrong offsets or a BufferUnderflowException. Re-index whenever the data file changes.

tabix: the file is not BGZF. You used gzip. Do gunzip file.vcf.gz && bgzip file.vcf.

CRAM fails with a reference mismatch. The M5 checksums in the header do not match the FASTA you loaded. Compare samtools view -H x.cram | grep M5 against samtools dict ref.fa. There is no workaround other than finding the right FASTA.

Allele fraction looks wrong. Downsampling is on. Turn it off for the site you are judging, then turn it back on.

A “variant” sitting in a homopolymer or short tandem repeat. Indels at the edge of a run of the same base are the most common short-read artifact. Look for the indel appearing at slightly different positions across reads, soft clips nearby, and support that is strongly strand-biased. Illumina-style sequencing-by-synthesis chemistry has known error characteristics in these contexts ³, and the read-level view is where they become obvious.

Coverage cliffs in GC-extreme regions. First exons of many genes are GC-rich and routinely under-covered at 30x WGS. A gap in the pileup means you have no information there, which is different from a reference match. Depth and uniformity, not just total yield, determine which parts of your genome you can call on ⁴.

Coordinates from an old report land in the wrong gene. The report is GRCh37. Lift over with CrossMap or bcftools +liftover plus the UCSC chain file, and re-check by confirming the reference allele at the new position matches.

Finally: a genome view shows you evidence. The browser does not hand you a conclusion. Read-level inspection tells you whether a call is real. Whether a real call means anything for your health is a separate question, and one that belongs with a clinical geneticist or a genetic counselor working from a CLIA-validated result. Your research BAM is not the basis for that answer.

Footnotes

1. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩

2. Ewan Birney, Paul Flicek, Damian Keefe, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature, 2007. https://doi.org/10.1038/nature05874 ↩

3. Simon T. Bennett, C.L. Barnes, Anthony J. Cox, et al. Toward the $1000 Human Genome. Pharmacogenomics, 2005. https://doi.org/10.1517/14622416.6.4.373 ↩

4. Sang Tae Park, Jayoung Kim. Trends in Next-Generation Sequencing and a New Era for Whole Genome Sequencing. International Neurourology Journal, 2016. https://doi.org/10.5213/inj.1632742.371 ↩