# How to Plot Variant Density in Python: Step-by-Step Guide By the end of this guide you will be able to take a VCF file, compute the density of variants across every chromosome in windows of your chosen size, and produce a clean, labeled, publication-ready density plot with matplotlib. You will also know how to normalize density by callable sites, how to lay out a whole-genome multi-panel figure, how to sanity-check your numbers, and how to fix the errors that trip up most people the first time they attempt this. You need Python 3.9 or newer, plus the packages `numpy`, `pandas`, `matplotlib`, and `pysam`. If you want to work with the reference genome directly you will also want `biopython`. Everything installs in one line with pip, and all four core packages are also available through Bioconda, the conda channel built for bioinformatics software that distributes many Python tools alongside command-line programs [1]. Python has become the default scripting language for this kind of work because it combines readable syntax with fast numeric libraries and mature plotting. Recent bioinformatics tooling papers keep reinforcing the pattern. A whole ecosystem of Python packages now exists for scalable analysis and visualization of large biological datasets, from phylogenomics platforms like ETE Toolkit and TreeProfiler [2] to sequence-design evaluation libraries like seqme that ship their own diagnostics and visualization functions [3]. The workflow below follows the same design principles those projects use. Parse once, compute once, plot last, and keep every step in a script you can rerun. ## What a Variant Density Plot Shows ![Worked calculation showing a 1 Mb bin with 340 SNVs giving a density of 340 SNV per Mb, and the general formula count divided by bin width in Mb.](/article-images/illustrations/variant-density-plot-python-density-formula.webp) A variant density plot answers one question: how densely are variants packed along a genome or chromosome? The input is usually a VCF file, which stands for Variant Call Format, the standard text format for storing genetic variants such as SNVs (single nucleotide variants, changes of one DNA letter) and indels (short insertions or deletions). The output is a line or bar track showing variants per megabyte (1 Mb equals 1,000,000 bases) across genomic coordinates. Three quantities matter: - **Variant count** per bin. A bin, also called a window, is a fixed coordinate interval, for example 1,000,000 bases. - **Bin width**. This controls smoothing. Wide bins give smooth curves that hide local structure. Narrow bins give noisy tracks that reveal hotspots. - **Density**, which is count divided by bin size, usually expressed as variants per Mb. If a 1 Mb bin contains 340 SNVs, the density is 340 SNV/Mb. Density plots are used to spot regions with unusually few variants (often low mappability or strong conservation), regions with excess variation (often repetitive sequence, paralogs, or contamination), and to compare variant distributions between callers, filters, or samples before deeper analysis. In genome-wide association work, the same windowing logic underlies Manhattan-style genome-wide visualizations, where Python pipelines integrate genome-wide plotting with downstream LD and fine-mapping analyses [4]. ## Setting Up Your Environment Use a dedicated environment so you do not break other projects: ```bash python -m venv venv-variantdensity source venv-variantdensity/bin/activate pip install numpy pandas matplotlib pysam biopython ``` Verify the install: ```bash python -c "import numpy, pandas, matplotlib, pysam, Bio; print('ok')" ``` Expected output: ``` ok ``` If you use conda instead, install from the bioconda channel, which is the standard distribution route for many bioinformatics Python tools [1]. For test data, you can download any published VCF, or simulate one so you can verify every step before touching real data. A quick simulation: ```python import numpy as np rng = np.random.default_rng(42) chrom_len = 50_000_000 # 50 Mb toy chromosome n_variants = 40_000 # Poisson-distributed counts per 10 kb window create realistic-looking clumping bin_size = 10_000 n_bins = chrom_len // bin_size counts = rng.poisson(lam=n_variants / n_bins, size=n_bins) positions = [] for i, c in enumerate(counts): start = i * bin_size positions.extend(start + rng.integers(1, bin_size, size=c)) positions = np.sort(np.array(positions)) np.save("sim_positions.npy", positions) print(len(positions), "variants simulated") ``` Expected output: ``` 40000 variants simulated ``` Keep this simulation script around. It is the single most useful tool for debugging a density plot pipeline, because you know the ground truth. ## Step-by-Step: From VCF to Density Plot ### Step 1: Inspect and prepare the VCF VCF coordinates are 1-based, meaning the first base of a chromosome is position 1. This matters later when you bin with numpy, whose arrays are 0-based. Most pipelines also want a compressed, indexed VCF so random access works: ```bash bgzip input.vcf bcftools index input.vcf.gz ``` A useful pre-filter keeps only confident biallelic SNVs, which avoids double counting and simplifies interpretation: ```bash bcftools view -v snps -m2 -M2 input.vcf.gz -Oz -o snps.biallelic.vcf.gz bcftools index snps.biallelic.vcf.gz ``` The flags mean: `-v snps` keeps SNVs only, `-m2 -M2` keeps sites with a minimum and maximum of two alleles (biallelic), and `-Oz` writes compressed output. Record how many variants you started with, because you will check against this number later. ### Step 2: Parse variant positions with pysam `pysam` is the Python interface to htslib, the same C library behind samtools and bcftools. It streams the file rather than loading it all into memory, which matters when your VCF has tens of millions of rows. ```python import pysam import numpy as np def parse_positions(vcf_path): """Return a dict mapping chromosome name to a sorted numpy array of 1-based POS values.""" vcf = pysam.VariantFile(vcf_path) per_chrom = {} for rec in vcf: per_chrom.setdefault(rec.chrom, []).append(rec.stop) return {c: np.sort(np.array(p, dtype=np.int64)) for c, p in per_chrom.items()} positions = parse_positions("snps.biallelic.vcf.gz") for chrom, pos in sorted(positions.items()): print(f"{chrom}\t{len(pos)} variants\tlast pos {pos[-1]}") ``` Expected output for a human callset looks like: ``` chr1 1985210 variants last pos 248946422 chr2 1291044 variants last pos 242193529 chrM 5123 variants last pos 16569 ... ``` A note on `rec.stop`: for a simple SNV, `rec.start` and `rec.stop` differ by one because pysam internally stores 0-based half-open intervals, while `rec.pos` is the 1-based position you see in the VCF. For density plotting either works as long as you are consistent. Using `rec.pos` is the least confusing choice for beginners. Pick one and document it in your script. ### Step 3: Get chromosome lengths You need true chromosome lengths so trailing empty bins are not dropped. Two good sources: From the VCF header itself: ```python import pysam vcf = pysam.VariantFile("snps.biallelic.vcf.gz") lengths = {r.name: r.length for r in vcf.header.records if r.type == "CONTIG"} print(lengths["chr1"]) ``` Or from a reference FASTA with Biopython, which is the classic biopython tutorial use case for sequence handling: ```python from Bio import SeqIO lengths = {} for rec in SeqIO.parse("genome.fasta", "fasta"): lengths[rec.id] = len(rec.seq) ``` Biopython shines at sequence and annotation file handling, and it pairs naturally with the numeric work here. If your VCF header contig lengths disagree with your reference FASTA, stop and resolve that before plotting, because it usually means a version mismatch between the callset and the reference build. ### Step 4: Bin the positions The core computation is a histogram over genomic coordinates: ```python import numpy as np def bin_density(pos, chrom_len, bin_size): """Return (counts per bin, bin right-edges) for 1-based positions.""" edges = np.arange(0, chrom_len + bin_size, bin_size, dtype=np.int64) counts, edges = np.histogram(pos, bins=edges) density = counts / (bin_size / 1_000_000) # variants per Mb return density, edges pos = positions["chr1"] density, edges = bin_density(pos, lengths["chr1"], bin_size=1_000_000) print(density[:10].round(1)) print("total check:", density.sum() * 1.0, "vs", len(pos), "variants") ``` Expected output on the simulated 50 Mb chromosome binned at 10 kb: ``` [3. 4.1 4.2 3.8 4.5 4.0 3.9 4.4 3.6 4.3] ``` The formula inside the function is the whole method: density = (variants in bin) / (bin width in Mb). With a 1 Mb bin the density equals the raw count, which is why many papers plot raw counts for 1 Mb windows and call it density. For any other bin width you must divide, otherwise changing the bin width silently changes the y-axis scale. ### Step 5: Plot the first figure ```python import matplotlib.pyplot as plt def plot_density(density, edges, chrom, bin_size): midpoints = (edges[:-1] + edges[1:]) / 2 fig, ax = plt.subplots(figsize=(12, 3)) ax.fill_between(midpoints, density, step="mid", color="#3b6ea5", alpha=0.85) ax.set_xlabel(f"{chrom} position (Mb)") ax.set_ylabel(f"Variants per {bin_size/1_000_000:g} Mb") ax.set_xlim(midpoints, midpoints[-1]) ax.margins(y=0.05) fig.tight_layout() fig.savefig(f"{chrom}_density.png", dpi=300) plt.show() plot_density(density, edges, "chr1", 1_000_000) ``` What you should see: a filled step-style track running left to right across the chromosome, with the x-axis in megabases and the y-axis in variants per megabase. Real human data typically shows a rolling baseline of a few hundred SNVs per Mb in accessible regions, deep dips at centromeres and other low-mappability zones, and occasional sharp spikes. If your simulated data looks like low-amplitude noise around the Poisson rate, the pipeline is correct. ### Step 6: Smooth with a rolling mean Raw windowed counts at fine bin sizes are spiky. A centered rolling mean is the standard fix: ```python import pandas as pd def rolling_density(density, edges, bin_size, smooth_mb=5): s = pd.Series(density) smoothed = s.rolling(window=smooth_mb, center=True, min_periods=1).mean() midpoints = (edges[:-1] + edges[1:]) / 2 / 1_000_000 return midpoints, smoothed.to_numpy() mid, smooth = rolling_density(density, edges, 1_000_000, smooth_mb=5) ``` A 5 Mb smoothing window on 1 Mb bins is a reasonable starting point for human-scale chromosomes. Always state the bin width and smoothing window in your figure legend, because the same underlying data produces visibly different curves under different choices. ### Step 7: Lay out every chromosome A whole-genome figure stacks one track per chromosome: ```python import matplotlib.pyplot as plt chroms = [c for c in positions if c.startswith("chr") and c != "chrM"] chroms = sorted(chroms, key=lambda c: int(c[3:]) if c[3:].isdigit() else 99) fig, axes = plt.subplots(len(chroms), 1, figsize=(12, 1.2 * len(chroms)), sharex=False) for ax, chrom in zip(axes, chroms): d, e = bin_density(positions[chrom], lengths[chrom], 1_000_000) mid = (e[:-1] + e[1:]) / 2 / 1_000_000 ax.fill_between(mid, d, step="mid", color="#3b6ea5", lw=0) ax.set_ylabel(chrom, rotation=0, ha="right", va="center") ax.set_yticks([]) ax.set_xlim(0, lengths[chrom] / 1_000_000) axes[-1].set_xlabel("Position (Mb)") fig.suptitle("SNV density, 1 Mb bins") fig.tight_layout() fig.savefig("genome_density.png", dpi=300) ``` Two details make or break this figure. Sort chromosomes in karyotype order with a explicit key, because default string sorting puts chr10 between chr1 and chr2. And drop unplaced contigs or move them to a separate panel, otherwise they compress the real chromosomes into a thin band. ### Step 8: Compare two callsets on one axis Overlaying a tumor and its matched normal, or a caller A versus caller B, is the most common follow-up: ```python d1, e1 = bin_density(pos_tumor, lengths["chr1"], 1_000_000) d2, e2 = bin_density(pos_normal, lengths["chr1"], 1_000_000) mid = (e1[:-1] + e1[1:]) / 2 / 1_000_000 plt.figure(figsize=(12, 3)) plt.fill_between(mid, d1, step="mid", alpha=0.4, label="tumor") plt.plot(mid, d2, color="k", lw=1, label="normal") plt.legend() plt.xlabel("chr1 position (Mb)") plt.ylabel("SNV per Mb") ``` This pattern generalizes to any number of tracks, and it is the same overlay logic used by comparative visualization tools in other imaging domains, where multi-specimen comparison on a shared scale is the central design goal [5]. ## Choosing Parameters: A Summary Table | Choice | Typical values | Effect | Practical rule | |--|--|--|--| | Bin width | 10 kb to 10 Mb | Narrow bins show hotspots, wide bins smooth | Start at 1 Mb for whole genomes, 10 kb for a single gene region | | [Density units](/blog/research-skills/density-units-explained-si-cgs-and-conversions) | count, per Mb, per callable Mb | Normalization changes comparability | Divide by callable bases whenever masks are available | | Smoothing | rolling mean 3 to 11 bins | Suppresses single-bin noise | State the window in the legend | | Fill vs line | fill_between, step plot | Fills emphasize area, lines ease overlay | Line for overlays, fill for single tracks | | Log y-axis | on or off | Compresses extreme spikes | Use only if you annotate that you did | | Per-sample layout | one panel per sample | Exposes sample-specific dropout | Keep identical y-limits across panels | For two-dimensional views, a different toolkit applies. Per-variant scatter plots colored by density (hexbin or 2D histogram) work well when you want position plus a quality metric, and heatmap layouts are standard when the unit of observation is a window crossed with a sample, the same layout sortscore uses for tiled variant-effect data in Python [6]. ## Verifying the Result Run these four checks before trusting any density plot: 1. **Count conservation.** `density.sum()` times the reciprocal of your per-Mb normalization must equal the number of parsed variants. If you normalized by callable sites, the sum should be slightly lower than the raw count. A mismatch almost always means variants fell outside your bin edges, usually positions beyond the chromosome length you supplied. 2. **Edge check.** Print `counts` from the first and last bins. If the last bin is empty but the VCF reported variants near the chromosome end, your `edges` array is too short. 3. **Cross-check with a second tool.** Count variants in one 10 Mb region with `bcftools view -H -r chr1:1-10000000 file.vcf.gz | wc -l` and compare to your histogram sum over those bins. 4. **Known-feature check.** Centromeres and assembly gaps should produce near-zero density. If a gap region shows dense variant calls, you are probably plotting unfiltered low-quality calls, or your coordinate system is off by the length of an offset contig. If all four pass, the plot is very likely correct. If one fails, fix it before interpreting anything biological. ## Common Mistakes and Limitations ![Six cards on density plot mistakes: bin size, coordinate systems, multiallelic double counting, gaps read as biology, KDE bleeding, and sort order.](/article-images/illustrations/variant-density-plot-python-common-mistakes.webp) These are the errors seen most often, with exact symptoms and fixes. **Plotting counts when bins differ in size.** If you trim the last bin to the actual chromosome end (a 250 Mb chromosome in 60 Mb bins leaves a 10 Mb remainder), the raw counts are not comparable. Symptom: an artificial dip or spike at the right edge of every chromosome. Fix: always divide by actual bin width, or make all bins equal and accept one partial bin. **Mixing 0-based and 1-based coordinates.** VCF is 1-based, BED is 0-based half-open, and pysam exposes both conventions through different attributes. Symptom: every peak shifted by exactly one bin, or off-by-one alignment against annotation tracks. Fix: standardize on `rec.pos` (1-based) inside the parsing function and convert only at the boundary with external files. **Counting multiallelic sites more than once.** A site with alleles A,T,C is one row but three potential variant relationships. If you iterate records and also iterate `rec.alts`, you can double or triple count. Symptom: total parsed count exceeds the line count of the VCF. Fix: filter to biallelic sites up front, or count `rec.pos` once per record, which is what the parsing function above does. **Interpreting density gaps as biology.** A zero-density region may mean strong conservation, but it more often means the region is not callable: centromeric repeats, segmental duplications, or low mappability where the caller simply has no power. Symptom: reporting "conserved" regions that are actually just unmappable. Fix: overlay a mappability or gap mask, and normalize density by callable bases per bin. Callable-normalized density is variants divided by the number of reference bases in the bin that survived masking, not by total bin width. **Letting seaborn KDE invent density where there is none.** Kernel density estimation smooths point data into a continuous curve, and near a chromosome edge or inside a large gap it will report nonzero density where zero variants exist. Symptom: smooth curves bleeding over centromeres. Fix: for positional data, stick to windowed histograms and rolling means. Reserve KDE for distributions where smoothing is genuinely appropriate. **Default chromosome sort order.** Python sorts strings lexicographically, so chr1, chr10, chr11, chr2. Symptom: a genome figure with scrambled panel order. Fix: the numeric sort key shown in Step 7. **Loading a whole large VCF into a Python list.** A callset with 50 million variants stored as Python integers can consume many gigabytes. Symptom: memory errors or extreme slowdown in `parse_positions`. Fix: stream per chromosome, store as `numpy.int64` arrays, or convert the VCF to a plain two-column position file once with bcftools and read it with `np.loadtxt` in chunks. **Forgetting the y-axis label.** A density plot without units, bin width, and smoothing window cannot be compared to anything, including your own rerun next month. Fix: put all three in the axis label or caption, every time. One further limitation to state plainly: density is a summary, and it compresses away allele frequencies, genotypes, and variant types. Two callsets can have identical density profiles and completely different biology. Density plots are for quality control and pattern discovery, never for the final claim. ## Variant Density Plot Python Tips for Speed and Scale A few habits keep this workflow fast as data grows. **Convert once, reuse forever.** Parsing VCF text is the slow step. Dump positions to a compressed `.npy` file per chromosome and reload those in every later session. Reloading a 50 Mb chromosome of int64 positions takes milliseconds. **Use `np.histogram`, not Python loops.** A loop appending bins one variant at a time is orders of magnitude slower than a single vectorized histogram call. This is the same principle behind the lightweight matrix-file decomposition used in Python visualization frameworks for large multi-channel datasets, where pre-splitting data into simple, widely supported formats keeps later steps fast [5]. **Downsample for exploration.** While iterating on figure style, plot every tenth bin. Switch to full resolution only for the final export. **Decouple computation from plotting.** Save density arrays to a small `.npz` file, then plot from that file. This separation makes figure revisions take seconds instead of minutes and keeps the compute step identical across figures. **Package the script for reuse.** Even a single `density.py` with a `main()` function and argparse options will save you time across projects. Reproducibility is the recurring theme in current bioinformatics tooling: recent Python frameworks emphasize scripted, interfaces-based access precisely so analyses can run on remote and high-performance systems without GUI steps [7], and end-to-end pipelines increasingly wrap such scripts in containerized workflows with companion interactive viewers for exploration [8]. ## Where This Fits in the Larger Python for Bioinformatics Workflow A density plot is rarely the end goal. It is usually the first genome-wide look that gates everything downstream. In variant calling QC, density across a tumor-normal pair flags contamination, sample swaps, and caller artifacts before you spend compute on annotation. In population genetics, per-window density alongside diversity statistics shows where recombination and selection shaped the genome. In association studies, the visualization layer feeds directly into locus selection, and published Python pipelines now integrate genome-wide visualization with LD estimation and fine-mapping in one scripted workflow [4]. In variant-effect screens, the windowed-track aesthetic carries over directly, and Python packages built for multiplexed assays ship heatmap visualizations of variant scores as a core feature [6]. The wider ecosystem is worth knowing as your needs grow. For annotation-heavy figures on phylogenetic or phylogenomic data, the ETE Toolkit ecosystem provides Python-native tree and dataset visualization at scale [2]. For sequence design projects, seqme provides property and embedding metrics with built-in diagnostic plots in the same scientific Python stack [3]. For single-cell work adjacent to variant analysis, distribution-shift testing tools are packaged for pip and conda installation in the standard scientific Python mold [1]. And structured, step-by-step processing chapters increasingly demonstrate both programmatic Python workflows and GUI alternatives side by side so users at different skill levels can follow the same protocol [9]. Learning the density plot workflow taught here transfers to all of these, because the pattern is identical: parse, bin, normalize, plot, verify. ## Frequently Asked Questions ### What is a variant density plot used for? It shows how variants are distributed along a genome, in counts or variants per megabase across fixed windows. Analysts use it to find low-variant regions (often low mappability or conservation), detect hotspots (often repeats or artifacts), compare callsets, and QC a variant caller before deeper analysis. ### What bin size should I use for a variant density plot? Start with 1 Mb bins for whole-genome views and 10 kb bins for regional views, then adjust. Narrower bins reveal local structure but get noisy, and wider bins smooth away hotspots. Whatever you choose, report the bin width in the figure legend. ### Should I plot raw counts or variants per Mb? Plot variants per Mb whenever your bins are not all exactly 1 Mb, or when you compare datasets with different bin widths. Raw counts are acceptable only when every bin is the same size and you say so. For cross-dataset comparisons, normalize by callable bases per bin. ### Can I make a variant density plot with Biopython alone? You can parse sequences and get chromosome lengths with Biopython, and older Biopython VCF support has been retired, so the common modern approach is pysam or cyvcf2 for VCF parsing combined with Biopython for reference and annotation handling. The binning and plotting then use numpy and matplotlib. ### Why does my density plot show spikes at centromeres or assembly gaps? Gaps and repetitive centromeric sequence usually produce near-zero density because nothing maps there, so a spike in those regions means unfiltered low-quality calls, coordinate offsets, or contigs from a mismatched reference build. Check the count conservation and coordinate-system checks before drawing biological conclusions. ### How do I plot density for multiple samples or VCFs? Parse each VCF with the same function, bin each with identical bin sizes and chromosome lengths, then either overlay tracks with matching y-limits on one axis or draw one subplot panel per sample. Keep bin width, normalization, and smoothing identical across all of them. ### How do I handle a VCF that is too large for memory? Stream the file with pysam rather than building a full list, store positions as numpy int64 arrays per chromosome, and cache them to compressed `.npy` files so you parse only once. Converting the VCF to a plain position list with bcftools and reading it in chunks also works well. ### How do I save a publication-quality figure? Use `fig.savefig("name.png", dpi=300)` for raster or `fig.savefig("name.pdf")` for vector output, call `fig.tight_layout()` before saving, and set figure size in inches to match your target column width. Vector formats keep text sharp and are preferred by most journals. ## Sources 1. [Quantifying distribution shifts in single-cell data with scXMatch.](https://pubmed.ncbi.nlm.nih.gov/42830441/) 2. [An Integrated Ecosystem for Analysis and Visualization of Large-Scale Phylogenomic Data Combining ETE Toolkit, TreeProfiler, and PhyloCloud for Scalable Phylogenomic Analysis.](https://pubmed.ncbi.nlm.nih.gov/42961754/) 3. [seqme: a Python library for evaluating biological sequence design from generative models.](https://pubmed.ncbi.nlm.nih.gov/42633207/) 4. [Large-Scale Data Analysis of Post-Traumatic Stress Disorder (PTSD) Through GWAS Fine-Mapping and Systems Biology.](https://pubmed.ncbi.nlm.nih.gov/42645920/) 5. [ScaleBarOn and Muad'Data: simple Python tools for elemental imaging data visualization and comparative scaling.](https://pubmed.ncbi.nlm.nih.gov/42657225/) 6. [sortscore: Sort-seq MAVE scoring and visualization using Python.](https://pubmed.ncbi.nlm.nih.gov/42684039/) 7. [Fast-nnt: fast, reproducible, and scalable neighbour network analysis in R, Python, and CLI.](https://pubmed.ncbi.nlm.nih.gov/42578716/) 8. [PerSeqPIPE: Nextflow workflow for per-sequence analysis of small RNA-sequencing data.](https://pubmed.ncbi.nlm.nih.gov/42893276/) 9. [Computational Tools for LC-IMS-MS Data Processing in Metabolomics.](https://pubmed.ncbi.nlm.nih.gov/42680970/)