How to Visualize Your Own Genome and Transcriptome Data

Oak
genomics visualization igv rna-seq bioinformatics
A laboratory instrument projects a glowing rotating ring of banded light above its glass deck, with luminous horizontal track ribbons below.

By the end of this guide you will have a small set of reproducible figures over your own data: a locus-level view in IGV that shows aligned reads and called variants together, a genome-wide copy-number profile, a coverage track you can scan for dropouts, an expression MA plot and volcano plot from RNA sequencing, and a single circular summary that puts variants, coverage, and expression on one coordinate system. You need your data in standard formats: a coordinate-sorted BAM or CRAM file with its index (.bai or .crai), a bgzip-compressed VCF with a .tbi index, a reference FASTA matching the build your files were aligned to (GRCh38 in almost all current pipelines), and for RNA a quantification matrix from salmon or featureCounts. Plan for roughly 100 GB of disk for a 30x whole genome in CRAM plus indexes and intermediates, 16 GB of RAM as a practical floor, and samtools, bcftools, bedtools, mosdepth, CNVkit, IGV, and an R or Python environment installed. Everything below runs on a laptop except the alignment itself, which we assume has already happened.

1. Confirm your coordinate system before you plot anything

Nearly every confusing genomic figure traces back to a mismatch between the coordinate system of the data and the coordinate system of the annotation. Genome browsers place everything on a shared linear axis of chromosome and position, and that axis is only meaningful if every file agrees on the reference build and the contig naming convention 1. GRCh37 and GRCh38 differ by megabases of coordinate shift in places, and a variant plotted on the wrong build lands in the wrong gene.

Start by reading the header of your alignment file and the contigs of your VCF:

samtools view -H sample.cram | grep '^@SQ' | head -5
# @SQ  SN:chr1  LN:248956422  M5:...
bcftools view -h sample.vcf.gz | grep '^##contig' | head -5

A chr1 of length 248,956,422 is GRCh38. A 1 of length 249,250,621 is GRCh37. If one file uses chr1 and the other uses 1, fix it once at the file level rather than patching it in each tool. bcftools annotate --rename-chrs takes a two-column map:

printf '1\tchr1\n2\tchr2\n' > chr_map.txt   # extend to all contigs
bcftools annotate --rename-chrs chr_map.txt sample.vcf.gz -Oz -o sample.chr.vcf.gz
bcftools index -t sample.chr.vcf.gz

Also record the reference FASTA’s MD5 or the M5 tags from the BAM header somewhere in your project notes. CRAM stores sequence by reference difference, so a CRAM file read against the wrong FASTA will silently produce wrong bases rather than failing loudly. This is the single most common way to generate a plot that looks fine and is wrong.

2. Build a coverage track and look at it genome-wide first

Before examining any specific gene, look at coverage across the whole genome. Coverage is the variable that limits every downstream interpretation: a variant call in a region with 4x depth means something very different from the same call at 35x, and a gene that appears deleted is often a region that simply failed to capture or sequence. mosdepth is the tool we use because it computes binned depth from a BAM or CRAM in a few minutes on a single genome and writes BED and bigWig-ready output directly.

mosdepth --by 1000 --fast-mode --no-per-base \
  -f GRCh38.fa -t 4 sample_1kb sample.cram

This writes sample_1kb.regions.bed.gz with mean depth per 1 kb bin and sample_1kb.mosdepth.summary.txt with per-chromosome means. Convert the BED to bigWig so browsers can render it at any zoom without re-reading the whole file:

zcat sample_1kb.regions.bed.gz | cut -f1-4 | sort -k1,1 -k2,2n > sample_1kb.bg
bedGraphToBigWig sample_1kb.bg GRCh38.chrom.sizes sample_1kb.bw

The bigWig format matters here for a structural reason. It stores a precomputed pyramid of summary zoom levels, so a browser drawing chromosome 1 at whole-chromosome scale reads a few kilobytes rather than 250,000 bins. Any format without that pyramid (plain bedGraph, wiggle) forces the browser to stream everything and will make an interactive session unusable at genome scale.

Plot the distribution of bin depths in R or Python before looking at the track. A clean 30x genome gives a roughly unimodal distribution centered near 30 with a long right tail from segmental duplications and a spike at zero from centromeres, telomeres, and unmappable regions. A bimodal distribution with a second mode near half the primary mode is a signal worth investigating, as is a mean that drifts systematically with GC content.

3. Set up IGV for locus-level inspection

The locus view is where you go when you have a specific question about a specific position, and it is the one visualization that cannot be replaced by a summary statistic. IGV was designed around exactly this workload: streaming indexed alignment and variant data from disk or over HTTP so that you can move from whole-genome to single-base resolution without loading everything into memory 2. Load your reference genome first (Genomes → Human GRCh38), then drag in sample.cram, sample.vcf.gz, and sample_1kb.bw.

Configure the alignment track deliberately, because the defaults hide the things you most want to see. Right-click the track and set:

  • Group alignments by → none, Color alignments by → insert size and pair orientation. This is what makes structural variants visible: a deletion appears as pairs with abnormally large insert size drawn in red, an inversion as pairs with the wrong orientation drawn in blue or green.
  • Show soft-clipped bases → on. Soft clips clustered at one position are the signature of an indel or breakpoint that the aligner declined to represent in the CIGAR string.
  • Shade base by quality → on, and set the coverage allele-fraction threshold (Preferences → Alignments) to 0.15 rather than the default 0.2 if you are looking for low-fraction events.

For reproducible screenshots, drive IGV in batch mode instead of clicking. A batch script is a plain text file of commands:

new
genome hg38
load sample.cram
load sample.vcf.gz
snapshotDirectory ./figures
goto chr7:117,559,590-117,559,700
sort base
squish
snapshot CFTR_exon10.png
exit

Run it with igv.sh --batch script.txt. Committing these scripts alongside your analysis code is what makes a figure re-creatable six months later, and it matters more than it sounds: the manual zoom level and sort order you chose are part of the figure and are otherwise lost.

4. Produce a copy-number profile

Copy number is the layer people most often skip and most often need, because deletions and duplications that span whole exons are invisible in a variant call file built for single-nucleotide changes. The standard visualization is a scatter plot of log2 coverage ratio against genomic position, with segmented mean lines drawn over the points. CNVkit computes this from coverage in fixed bins, corrects for GC content and mappability bias, and emits both the plot and the underlying segment table 3.

For a whole genome, run in whole-genome mode against a flat reference if you do not have a matched panel of normals:

cnvkit.py batch sample.cram -m wgs \
  -f GRCh38.fa --annotate refFlat.txt \
  -n -p 8 --output-dir cnv/
cnvkit.py scatter cnv/sample.cnr -s cnv/sample.cns -o cnv/sample_scatter.pdf
cnvkit.py diagram cnv/sample.cnr -s cnv/sample.cns -o cnv/sample_diagram.pdf

Read the scatter plot with the bias corrections in mind. Without a panel of normals (-n with no argument builds a flat reference), systematic coverage biases are only partially removed, and the residual wave pattern along each chromosome is technical rather than biological. Segments of a few bins at log2 ratios between -0.3 and 0.3 are noise. A convincing heterozygous deletion sits near -1.0 in log2 space across many consecutive bins, and a convincing duplication near +0.58. Zoom into any candidate with cnvkit.py scatter -c chr17:0-10000000 and then look at the same interval in IGV, because a real event usually shows supporting discordant read pairs and clipped reads at the boundaries.

Interpretation of a copy-number finding in a clinically relevant gene is a conversation with a clinical geneticist, not something to settle from a scatter plot. The plot tells you where to look and how confident the measurement is.

5. Visualize the transcriptome: MA plot, then volcano, then the gene itself

RNA sequencing gives you tens of thousands of measurements per sample, and the useful visualizations compress that into shape rather than listing values. Assume you have transcript-level quantification from salmon and have imported it with tximport into a DESeq2 object comparing two conditions or two timepoints from your own longitudinal sampling.

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ timepoint)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]
dds <- DESeq(dds)
res <- results(dds, contrast = c("timepoint", "t2", "t1"))

# 1. Dispersion plot: the QC check nobody skips twice
plotDispEsts(dds)

# 2. MA plot on shrunken effect sizes
library(apeglm)
resLFC <- lfcShrink(dds, coef = "timepoint_t2_vs_t1", type = "apeglm")
plotMA(resLFC, ylim = c(-4, 4))

Look at the dispersion plot first. You expect a cloud of gene-wise estimates scattered around a smooth fitted curve that decreases with mean expression, with shrunken final estimates pulled toward that curve. A cloud with no trend, or a large group of points far above the curve, means the model is not describing your data and every p-value downstream is suspect.

The MA plot (log2 fold change against mean normalized count) exists to show you that the fold changes are centered on zero and that the funnel of variance narrows as expression rises. Use shrunken fold changes for it. Unshrunken estimates at low counts produce enormous apparent fold changes from a handful of reads, and those points dominate the visual field without carrying information.

The volcano plot is then a straightforward scatter of the same fold changes against significance:

library(ggplot2)
d <- as.data.frame(resLFC)
d$sig <- !is.na(d$padj) & d$padj < 0.05 & abs(d$log2FoldChange) > 1
ggplot(d, aes(log2FoldChange, -log10(pvalue), color = sig)) +
  geom_point(size = 0.6, alpha = 0.5) +
  scale_color_manual(values = c("grey70", "firebrick")) +
  geom_vline(xintercept = c(-1, 1), linetype = 2) +
  theme_minimal()

Plot the raw p-value on the y-axis and use the adjusted value for coloring. Plotting adjusted p-values produces a horizontal shelf of tied values that misrepresents the underlying evidence. Then, for any gene you care about, go back to IGV with the RNA BAM loaded in Sashimi view (right-click → Sashimi Plot) to see whether the expression change is a whole-gene shift or a change in which exons are included. That distinction is invisible in gene-level counts and often the more interesting finding.

6. Assemble a genome-wide summary on one coordinate system

The final figure worth building puts several data types on a shared axis so you can see coincidences: a region of low coverage that also shows a copy-number loss and reduced expression of the genes inside it. Circular layouts are the conventional choice for this because a linear genome axis at readable resolution is about three meters long, and wrapping it into a circle lets multiple tracks share one axis in one screen 4.

In R, circlize is the shortest path:

library(circlize)
circos.initializeWithIdeogram(species = "hg38")
cov <- read.table("sample_1kb.bg", col.names = c("chr","start","end","depth"))
circos.genomicTrack(cov, numeric.column = 4, track.height = 0.12,
  panel.fun = function(region, value, ...) {
    circos.genomicLines(region, value, type = "l", col = "grey30")
  })
seg <- read.table("cnv/sample.cns", header = TRUE)
circos.genomicTrack(seg[, c("chromosome","start","end","log2")], numeric.column = 4,
  panel.fun = function(region, value, ...) {
    circos.genomicRect(region, value, ytop.column = 1, ybottom = 0,
      col = ifelse(value[[1]] > 0, "firebrick", "steelblue"), border = NA)
  })

Downsample the coverage track before plotting. Two and a half million 1 kb bins will render, eventually, into a figure where each bin occupies a fraction of a pixel. Bin to 100 kb with bedtools map and plot the median, which gives roughly 30,000 points and a readable line.

If you want the same view interactively in a browser, Gosling is the tool we would reach for. It is a declarative grammar in which you describe the data source, the mark type, and the channel mappings, and the runtime handles tiled loading and zooming across scales via GPU rendering 5. The cost is that your data must be in tiled formats (bigWig, or Higlass multivec for matrices), so you pay a preprocessing step to get fluid interaction.

Common problems

Coverage looks fine but variants in a gene are missing entirely. Check whether the gene has a high-identity paralog. Reads from segmental duplications get mapping quality 0, and most variant callers and browser tracks filter those by default. In IGV, turn off the mapping-quality threshold temporarily and see whether reads appear with white (MAPQ 0) shading. Multiple-alignment context is the right way to reason about these regions, and browser tracks built from whole-genome alignments across species exist precisely to expose such duplicated and conserved structure 6.

A variant appears in the VCF track but the reads underneath do not support it. Look at the soft-clip display and the base qualities. The usual causes are a nearby indel that the aligner represented differently in different reads, a strand bias where all supporting reads point the same direction, or a homopolymer run where the sequencer’s error mode produces consistent-looking noise. If the supporting reads all start and end at the same coordinates, you are probably looking at PCR duplicates that escaped marking.

The copy-number scatter shows a wave along every chromosome. This is GC and replication-timing bias that a flat reference cannot remove. Either build a reference from several samples processed identically, or restrict interpretation to segments whose log2 ratio exceeds the wave amplitude, which you can estimate from the standard deviation of bin-level log2 values in known-diploid regions.

RNA expression and DNA copy number disagree for a gene. This is often real. Transcript abundance is regulated by much more than gene dosage, and large integrative studies of matched DNA and RNA from the same samples show only partial correlation between copy number and expression across genes 7. Treat the disagreement as information about regulation rather than as a pipeline error, after you have checked that both measurements used the same gene model.

A heatmap of a contact matrix or any large matrix renders as a solid block. Matrix visualizations need explicit balancing and a clipped color scale. Raw contact or correlation matrices have dynamic ranges spanning several orders of magnitude, and the convention in Hi-C analysis is to normalize the matrix and then map color to a percentile-clipped range so that structure at moderate values is visible rather than swamped by the diagonal 8.

Interactive browsers become unusable at whole-genome zoom. The cause is always a track in a non-indexed, non-tiled format. Convert bedGraph to bigWig, BED to bigBed, plain VCF to bgzip plus tabix, and SAM to indexed CRAM. The historical lesson here is old and well documented: presenting genome-scale annotation interactively requires precomputed hierarchical summaries, not on-the-fly aggregation of raw records 9.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Cydney B Nielsen, Michael Cantor, Inna Dubchak, et al. Visualizing genomes: techniques and challenges. Nature Methods, 2010. https://doi.org/10.1038/nmeth.1422 ↩

  2. H. Thorvaldsdottir, J. T. Robinson, J. P. Mesirov. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Briefings in Bioinformatics, 2012. https://doi.org/10.1093/bib/bbs017 ↩

  3. Eric Talevich, A. Hunter Shain, Thomas Botton, et al. CNVkit: Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing. PLOS Computational Biology, 2016. https://doi.org/10.1371/journal.pcbi.1004873 ↩

  4. Zhonglin Qu, Chng Wei Lau, Quang Vinh Nguyen, et al. Visual Analytics of Genomic and Cancer Data: A Systematic Review. Cancer Informatics, 2019. https://doi.org/10.1177/1176935119835546 ↩

  5. Sehi LYi, Qianwen Wang, Fritz Lekschas, et al. Gosling: A Grammar-based Toolkit for Scalable and Interactive Genomics Data Visualization. IEEE Transactions on Visualization and Computer Graphics, 2022. https://doi.org/10.1109/tvcg.2021.3114876 ↩

  6. Mathieu Blanchette, W. James Kent, Cathy Riemer, et al. Aligning Multiple Genomic Sequences With the Threaded Blockset Aligner. Genome Research, 2004. https://doi.org/10.1101/gr.1933104 ↩

  7. The Cancer Genome Atlas Network. Comprehensive molecular portraits of human breast tumours. Nature, 2012. https://doi.org/10.1038/nature11412 ↩

  8. Suhas S.P. Rao, Miriam H. Huntley, Neva C. Durand, et al. A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping. Cell, 2015. https://doi.org/10.1016/j.cell.2015.07.024 ↩

  9. Ann E Loraine, Gregg A Helt. Visualizing the genome: techniques for presenting human genome data and annotations. BMC Bioinformatics, 2002. https://doi.org/10.1186/1471-2105-3-19 ↩