How to Visualize RNA-Seq Data from BAM Files in IGV

By Dr. Zubair Khalid, DVM, MS, PhD ·

How to Visualize RNA-Seq Data from BAM Files in IGV

By the end of this walkthrough you will have a coordinate-sorted RNA-seq BAM open in IGV, a coverage track you can describe accurately in a figure legend, a splice junction track that shows exon skipping at a glance, and a saved session plus snapshot images you can hand to a collaborator. You will also know how to confirm that what IGV draws matches what samtools reports, which is the difference between a picture and evidence.

You need five things on hand before you start: an IGV installation (version 2.16 or later), a Java runtime (version 11 or later, which IGV requires), samtools (version 1.16 or later is a safe baseline), your BAM file, and its .bai index. If you are working with a non-model organism, or with an assembly your aligner used directly rather than a hosted genome, you also need the reference FASTA and its .fai index. Everything below runs on Linux, macOS, or Windows with small path differences.

What IGV Actually Draws for an RNA-Seq BAM

A BAM file is a compressed binary form of SAM, the format that stores each sequencing read, its alignment position against a reference, and the operations that describe how the read relates to that reference [1]. The critical field for RNA-seq is the CIGAR string. A CIGAR of 60M5000N60M means 60 bases match, 5,000 reference bases are skipped, and another 60 bases match. In a spliced alignment produced by STAR, HISAT2, or a similar aligner, that N operation is the intron.

IGV reads the CIGAR and draws three things. The alignment track shows each read as a row of boxes, with the skipped intron rendered as a thin connecting line. The coverage track shows per-base depth of the reads currently visible. The splice junction track draws arcs between the donor and acceptor positions of every junction that appears in at least the minimum number of reads.

That structure is why IGV is the standard first look at RNA-seq data. Alternative splicing work routinely combines RNA sequencing with targeted validation, and the junction arcs are the fastest way to spot a candidate event before committing to a wet-lab assay [2]. Understanding splicing patterns in RNA-seq is a common use of high-throughput sequencing, and visualization tools have become the default entry point for exploring coverage patterns across regions of interest [3].

IGV is a viewer, not a quantifier. It will never tell you that a gene is differentially expressed. It tells you whether the alignment looks like the biology you expected.

Files You Need and the Formats Behind Them

Bioinformatics file formats are containers with conventions, and mixing them up is the single most common reason a track fails to load. Here is what each file in this workflow does.

BAM stores aligned reads. It is binary, compressed, and requires random access, which is why it must be coordinate sorted. The companion BAI index is a lookup table that lets IGV jump to a genomic position without scanning the whole file. Without the index, IGV cannot display a BAM from a local or remote source.

FASTA stores sequences with no quality information. It is the format of a reference genome, a transcriptome, or a protein set. A FASTA index, the .fai file, records the byte offset and line length of every contig so tools can seek into the file. IGV can use a hosted genome such as GRCh38 or mm10, or you can load your own FASTA, which is the same approach used when setting up other genome browsers from a FASTA reference [4].

GTF or GFF3 stores gene annotation as features with genomic coordinates. Loading one gives IGV gene models to draw and makes the search box accept gene symbols. Annotation in GFF3 is the standard companion to a FASTA reference in browser setup instructions [4].

FASTQ vs FASTA

These two are frequently confused because both are plain text and both start with a punctuation character.

FASTA uses a > header line followed by sequence lines wrapped at a fixed width. A record is two or more lines and contains no per-base quality.

>chr1
NNNNNNNNNNNNNNNNNNNNACGTACGTACGTACGTACGT
ACGTACGTACGTACGTNNNNNNNNNNNNNNNNNNNNNNNN

FASTQ uses exactly four lines per read. The first line starts with @ and holds the read identifier. The second line is the sequence. The third line starts with +. The fourth line is one quality character per base, encoded as Phred+33, where a I character is a high-confidence base and a ! is essentially no confidence.

@A00123:45:H7L2MDSX2:1:1101:1000:1000 1:N:0:ACGTACGT
GATCGATCGTAGCTAGCTAGCTAACGTTAGCTAGCTAGCTA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

The practical consequence is simple. FASTQ is the raw output of the sequencer and the input to your aligner. FASTA is the reference you align against or the file you load into IGV to define coordinates. RNA-seq analysis flows from FASTQ to BAM, and IGV visualizes the BAM. There is no route in this workflow that loads a FASTQ into IGV, because IGV has no alignment to draw from unaligned reads.

Step-by-Step: Visualize RNA-Seq BAM Files in IGV

The sequence below assumes you have an unsorted or unknown BAM named sample.bam and want a clean, indexed file you can trust.

Step 1. Verify the BAM is complete and coordinate sorted

Run a quick integrity check first. samtools quickcheck reads the file header and the end-of-file marker, which catches truncated downloads and half-written files.

samtools quickcheck -v sample.bam
echo $?

Expected output on a healthy file is nothing at all, followed by 0 from the exit status. Any printed line means the file failed the check. A common message is that the file was not identified as a BAM file, which usually means the download is incomplete or the file is still a gzip of FASTQ under a misleading name.

Next, inspect the header to confirm the sort order and the contig naming.

samtools view -H sample.bam | head -n 5
@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:248956422
@SQ	SN:chr2	LN:242193529

The field that matters is SO:coordinate. If you see SO:unsorted or the field is absent, sort the file. Note the contig names too, because chr1 and 1 are different strings to a computer.

samtools sort -@ 8 -m 2G -o sample.sorted.bam sample.bam
samtools index -@ 4 sample.sorted.bam
ls -l sample.sorted.bam sample.sorted.bam.bai

Expected output is a BAM and a .bai of similar age and reasonable size. The BAM will typically be a few gigabytes for a human RNA-seq library depending on depth.

Step 2. Prepare the reference and its index

If the genome you aligned to is hosted in IGV (GRCh38, hg19, mm10, and many others), skip this step. If you aligned to a custom assembly, create the FASTA index.

samtools faidx GRCh38.fa
ls -l GRCh38.fa.fai
cut -f1,2 GRCh38.fa.fai | head -n 3

The .fai file is a tab-separated table of contig name, length, offset, line bases, and line width. The first two columns are the ones most tools want.

Contig names in your FASTA, your GTF, and your BAM header must match exactly. A BAM aligned with 1, 2, 3 will show nothing when loaded against a FASTA that names contigs chr1, chr2, chr3.

Step 3. Launch IGV with enough memory

Start IGV from the command line on Linux or macOS, or from igv.bat on Windows. Alignment-dense RNA-seq loci can exhaust the default Java heap. If IGV stalls at a highly expressed gene, raise the -Xmx value in the launcher script or pass a larger memory setting at launch.

igv.sh

If IGV fails to start with a message about an unsupported class file version, your Java is too old. Install Java 11 or later and point JAVA_HOME at it.

Step 4. Load the genome

From the menu, choose Genomes > Load Genome from File for a custom FASTA, or Genomes > Load Genome from Server for a hosted assembly. The genome selector in the top-left toolbar shows the current assembly. Confirm it matches the assembly your aligner used, not just the species. hg19 and hg38 coordinates for the same gene differ by megabases.

If you have a GTF or GFF3, load it as a track with File > Load from File. RefSeq and Ensembl gene models are available from the hosted genome list and are usually sufficient.

Step 5. Load the BAM and its index

Use File > Load from File and select sample.sorted.bam. IGV looks for sample.sorted.bam.bai in the same directory and loads it automatically. If the index has a different name, IGV will prompt you for it.

For a BAM stored on a web server or in cloud object storage, use File > Load from URL and paste the BAM URL. IGV will ask for the index URL separately. Both files must be reachable and the server must allow range requests. Remote loading from a slow network share is the most common cause of a track that appears frozen.

Step 6. Navigate to a locus

Type a coordinate or a gene symbol in the search box at the top. Coordinates use commas or not, both work.

chr7:55,190,000-55,200,000

Pressing Enter re-centers the view. If the box refuses a gene symbol, your annotation track is not loaded or it uses a different gene naming scheme from the one you typed.

Step 7. Configure the coverage track

The coverage track is the histogram at the top of the BAM panel. Right-click it to open the settings menu. The two controls that change your interpretation most are the data range and the windowing function.

Autoscale rescales the y-axis to the visible window every time you move. It is convenient and misleading, because a peak that looks identical to a neighboring sample may be twice as deep in absolute terms. For comparisons, right-click the coverage track, choose Set Data Range, and enter a fixed maximum used across all samples.

Windowing function controls how depth is summarized when you zoom out. Mean smooths over coverage gaps. Max exaggerates single-base pileups from mismapping or PCR duplicates. For most RNA-seq questions, mean is the honest default.

Log scale compresses the dynamic range so you can see low-level expression next to a highly expressed gene. Use it when you care about whether a transcript is present at all, not about fold change.

Step 8. Configure the alignment track

Right-click the alignment track and work through these settings in order.

Color alignments by > read strand. This colors reads by the strand they aligned to. It is the fastest visual check for library strandedness. For a non-stranded library, both colors appear in roughly equal proportion across an exon. For a stranded library, one color dominates.

Sort alignments by > base. Sorting by the base at the current position turns the alignment panel into a visual pileup and makes allele-specific expression and RNA editing candidates easy to spot.

Group by > sample. When multiple BAMs are loaded, grouping keeps samples contiguous and prevents read rows from interleaving across tracks.

View as pairs connects paired reads with a line. For RNA-seq it is mostly a distraction, since pairs often span an intron, but it is useful when checking library prep problems.

Downsampling. IGV draws only a subset of reads when coverage is high, so that the interface stays responsive. The cap is set in View > Preferences > Alignments. If your conclusion depends on seeing every read, raise the cap or load a smaller region. Otherwise state in your figure legend that reads were downsampled.

Alignment filters. The filter panel lets you hide reads by mapping quality, by flag, or by insert size. Hiding mapping quality below 20 removes most ambiguous multi-mapping reads and changes the coverage profile you see. That is fine as long as you say so.

Step 9. Read the splice junction track

The splice junction track sits below the coverage track. Each arc spans one junction, and the number on the arc is the read count supporting it. Right-click the track to set the minimum junction coverage. The default hides low-support junctions, which is usually desirable and always worth knowing, because a junction supported by two reads and a junction supported by 200 reads look the same once they both clear the threshold.

Junctions whose coordinates match annotated exon boundaries support the annotated transcript model. Junctions that land in the middle of an annotated exon are candidates for unannotated splicing, mapping artifacts, or retained introns. Zoom in until the read boxes and the arc endpoints align, then confirm that the skipped region in the CIGAR matches the gap you see on screen.

Step 10. Save the session and export snapshots

Before you change anything else, save the session. File > Save Session writes an XML file that records the genome, the tracks, the loaded loci, and the colors. A reviewer or collaborator can reopen it and see exactly what you saw.

For figures, use File > Save Image on the IGV window. Set the panel height explicitly rather than relying on the window size, and include the locus in the title bar so the figure carries its own coordinates.

Step 11. Automate repetitive snapshots with a batch script

If you need the same view for 20 genes across 6 samples, drive IGV from a batch script. Write a plain text file with one command per line.

new
genome hg38
load sample1.sorted.bam
load sample2.sorted.bam
load sample3.sorted.bam
snapshotDirectory /home/lab/figures/snapshots
goto chr7:55,190,000-55,200,000
sort base
snapshot chr7_region_base.png
goto chr7:55,208,000-55,215,000
snapshot chr7_region_zoom.png
exit

Run it with the batch flag.

igv.sh -b snapshots.txt

Expected output is a set of PNG files in the snapshot directory, one per snapshot command, plus a log line per command in the terminal. The new command clears any previously loaded session, which prevents state from a prior run leaking into your figures.

Key IGV Settings at a Glance

SettingWhere to find itWhy it matters for RNA-seqStarting point
Reference genomeGenomes menuCoordinates must match the aligner's assemblyHosted assembly or your FASTA plus .fai
BAM indexAutomatic or prompted on loadWithout it, no track renders.bai beside the BAM, same basename
Visibility windowTrack settings, right-clickReads disappear when zoomed out past the windowKeep under a few tens of kb, or accept the warning
Downsampling capView > Preferences > AlignmentsLimits how many reads are drawnRaise only when the conclusion needs every read
Coverage data rangeRight-click coverage trackAutoscale hides cross-sample differencesFixed maximum shared by all samples
Windowing functionRight-click coverage trackMean vs max changes apparent peak heightMean
Color alignments byRight-click alignment trackReveals library strandednessRead strand
Sort alignments byRight-click alignment trackTurns reads into a pileup for allele checksBase
Minimum junction coverageRight-click junction trackControls which splicing events appearRaise for confident events, lower to explore
Mapping quality filterRight-click alignment trackRemoves ambiguous multi-mappers20 if you will report the filter
Session fileFile > Save SessionReproducibility for collaborators and reviewersSave before exporting images

Common Errors and Exact Fixes

The track fails to load and IGV reports that it cannot find an index file. The .bai is missing, misnamed, or older than the BAM. Regenerate it with samtools index sample.sorted.bam. An index older than its BAM is stale and IGV may refuse it or display wrong positions. Compare modification times with ls -l.

The track loads but is completely empty. This is almost always a contig naming mismatch. Your BAM header says 1 and the loaded genome says chr1. Either load the genome that matches the aligner, or rename contigs consistently and reindex. There is usually no error message, which is what makes this one expensive.

samtools index fails with a message about unsorted positions. The file is not coordinate sorted, or a previous sort was interrupted. Re-run samtools sort and index the output, not the input.

IGV shows a message asking you to zoom in before alignments appear. You are outside the track's visibility window. Zoom in, or raise the visibility window in the track settings if you deliberately want a wide view. The coverage track still renders at wide zoom, which is why people mistake a missing alignment panel for missing data.

IGV crashes with java.lang.OutOfMemoryError: Java heap space. The locus has more reads than the heap can hold. Raise the heap in the launcher, or load a narrower region by extracting it first.

samtools view -b sample.sorted.bam chr7:55190000-55200000 > chr7_window.bam
samtools index chr7_window.bam

Colors look inverted and nearly every read appears antisense. Your library is strand-specific but you are reading the colors as if it were not. dUTP-based protocols reverse the relationship between transcript strand and read strand. Confirm the protocol, then interpret the colors accordingly or flip your mental map.

Splice junctions are missing entirely. The junction track may be hidden, or the minimum coverage threshold may exceed the support in your data. Lower the threshold and confirm that junction arcs appear.

Loading from a network mount is unusably slow. Copy the BAM and BAI to local disk. Indexed random access over a high-latency share is slow by design, because every pan triggers new range requests.

How to Check That What You See Is Correct

A picture is only evidence if you can reproduce its numbers. Three checks catch nearly every configuration error.

Check the read count in a window. Count reads with the same filters IGV applies by default, which excludes unmapped, secondary, and supplementary alignments. The flag mask 0x904 covers those three categories.

samtools view -c -F 0x904 sample.sorted.bam chr7:55190000-55200000

The number printed should be in the same range as the read count IGV reports for the same window. A large discrepancy means the filters, the downsampling, or the reference are not what you think.

Check the mean depth. Compute it directly and compare it with the coverage track's mean windowing.

samtools depth -a -r chr7:55190000-55200000 sample.sorted.bam \
  | awk '{s+=$3} END {printf "%.1f\n", s/NR}'

If IGV's coverage peak disagrees with this value by more than the effect of downsampling, check the data range setting. A fixed data range that is too low clips the top of the histogram and makes a deep region look flat.

Check the junctions against the annotation. Pick a gene with a well-known structure, confirm that the arcs land on annotated exon boundaries, and confirm that the skipped distance matches the annotated intron length. If the arcs are offset by a constant amount, you are almost certainly viewing a liftover or a mismatched annotation build.

Run one positive control gene and one negative control region. A housekeeping gene should show strong, uniform exonic coverage. An intergenic desert should show near-zero coverage with only scattered mismapped reads. If the intergenic region looks busy, your mapping quality filter is too permissive or the library has genomic DNA contamination.

Common Mistakes and Limitations

Treating the coverage track as a quantification. Coverage height is a function of library size, gene length, mappability, and the filters you applied. Two samples with equal peaks can differ twofold in normalized expression. Use a dedicated quantifier such as featureCounts, HTSeq, or RSEM, then test with DESeq2 or edgeR. IGV answers "does this look right", not "how much".

Leaving autoscale on across a multi-sample figure. Each track rescales independently, so a weak sample and a strong sample both fill the panel. Fix the data range before you export.

Forgetting that the coverage track sums both strands. The coverage histogram does not separate sense from antisense signal. If your question involves antisense transcription, overlapping genes, or strand-specific library QC, you must inspect the alignment track colored by read strand, and even then you are looking at reads rather than a computed stranded coverage profile.

Ignoring multi-mapping reads. Reads from paralogs and pseudogenes pile up on the wrong locus. The NH tag records how many alignments a read had. Filtering on mapping quality removes most of the problem and also removes real signal from duplicated genes, so state which choice you made.

Leaving PCR duplicates in the picture. Duplicates inflate coverage at highly expressed transcripts. Marking or removing them changes the shape of your peaks, and the choice belongs in the methods, not in a caption footnote.

Assuming empty means unexpressed. A gene with zero visible reads in one sample may be unexpressed, or may be below the downsampling threshold at that zoom level, or may be missing because your annotation uses a different transcript. Check the raw count table before drawing a biological conclusion.

Reading alternative splicing off the junction track alone. A junction with low support can be real or can be noise, and a junction below the display threshold is invisible. Alternative splicing conclusions from RNA-seq are conventionally validated by an orthogonal assay, which is precisely why published protocols pair sequencing with targeted amplification [2].

Using IGV where a genome-wide summary is needed. IGV shows one locus at a time. When you need average signal across thousands of regions, profile plots and clustered heatmaps are the appropriate visualization, and tools built for that purpose accept the major genomic file formats directly [3]. Circular genome views are similarly a different tool for a different question [5].

Assuming every assay uses the same container. BAM is a genome-based representation that has been adapted well beyond DNA and RNA alignments, including proteomics data visualized in genome browsers [6]. Formats diverge quickly across assays, so confirm what your file actually contains before loading it.

Skipping the strandedness check. A stranded library analyzed as unstranded will produce plausible-looking but wrong antisense signal. Run a strandedness check on a known expressed gene before you interpret anything.

Biological interpretation still needs a second line of evidence. A junction arc or a coverage peak is a hypothesis. Confirm it with targeted amplification, with an independent dataset, or with a published annotation before it goes in a manuscript.

Sensible Next Steps

Once the BAM looks right in IGV, move the analysis forward rather than staring at more loci.

Quantify first. Generate a counts matrix with featureCounts using your GTF, then run a differential expression model. IGV becomes the tool you return to when a specific gene in the results table needs a visual check.

Build normalized signal tracks for comparison. Convert each BAM to a bigWig so coverage is comparable across samples and loads quickly in any browser. bigWig is a common format supported by many genome browsers, which makes it a good interchange format for collaborators [7].

bedtools genomecov -bg -split -ibam sample.sorted.bam > sample.bedgraph
cut -f1,2 GRCh38.fa.fai > GRCh38.chrom.sizes
bedGraphToBigWig sample.bedgraph GRCh38.chrom.sizes sample.bw

Use -split so spliced reads contribute coverage to their aligned blocks rather than to the intron they skip.

Move to a genome-wide view when the question is about patterns rather than one gene. Profile plots and heatmaps over sets of regions answer questions that a single-locus browser cannot [3].

Consider a different browser when you need to publish a shareable, web-accessible instance rather than a desktop session. JBrowse 2 is a modular browser that loads a FASTA reference and a GFF3 annotation track and runs on a standard Linux web server [4].

Check your upstream QC before you trust any of it. Adapter content, rRNA contamination, and read duplication all show up as visual artifacts, and they are cheaper to fix at the FASTQ stage than to explain in review.

Frequently Asked Questions

Can I open a BAM in IGV without an index?

No. IGV needs the .bai index to jump to a locus, so a BAM without one will not render as a track. Create it with samtools index after confirming the file is coordinate sorted.

Why is my BAM track empty when reads clearly aligned?

The usual cause is a contig naming mismatch between the BAM header and the loaded genome, such as 1 versus chr1. Confirm the header with samtools view -H and load the genome that matches your aligner.

Does IGV display spliced reads correctly?

Yes. IGV reads the N operation in the CIGAR string and draws the skipped intron as a gap, with a connecting arc in the splice junction track. That is the main reason it is useful for RNA-seq.

How do I compare two RNA-seq samples in IGV?

Load both BAMs and their indexes, then set a fixed data range on each coverage track instead of leaving autoscale on. Autoscale rescales per track and hides real differences in depth.

What is the difference between the coverage track and the alignment track?

The coverage track is a per-base depth histogram computed from the visible reads. The alignment track shows individual reads, their CIGAR operations, and their strand. Coverage summarizes, alignments explain.

Can I use CRAM instead of BAM?

IGV supports CRAM, but CRAM stores reads relative to a reference and therefore needs that reference available to decode. If you share CRAM files, share the matching reference or use a reference server.

How do I make a publication-quality figure from IGV?

Set the panel height explicitly, fix the coverage data range, load only the tracks you will discuss, and use File > Save Image. Batch scripts with the snapshot command generate the same figure across many loci without manual clicking.

Why does IGV show fewer reads than samtools counts?

IGV downsamples high-coverage regions by default and applies any alignment filters you have set. Check View > Preferences > Alignments and the filter panel, then compare with a count that uses the same flag mask.

Sources

  1. Online resources for genomic analysis using high-throughput sequencing.
  2. A Method to Quantitatively Examine Heat Stress-Induced Alternative Splicing in Plants by RNA-Seq and RT-PCR.
  3. SeqPlots - Interactive software for exploratory data analyses, pattern discovery and visualization in genomics.
  4. Setting Up the JBrowse 2 Genome Browser.
  5. Circleator: flexible circular visualization of genome-associated data with BioPerl and SVG.
  6. proBAMsuite, a Bioinformatics Framework for Genome-Based Representation and Analysis of Proteomics Data.
  7. Prioritising positively selected variants in whole-genome sequencing data using FineMAV.