Module 2.6 SnpEff, Combining VCFs, and SNP Ordination - mestato/EPP622_2026 GitHub Wiki

01. SnpEff: Structurally annotating your variants

  1. Variants are not created equal.
  • Variation in non-coding regions (e.g., introns and intergenic regions) is useful for population structure, kinship estimation, molecular dating, etc.
  • Variation in regulatory and coding regions (e.g., regulatory elements and exons) can help us predict impacts on function (sometimes).
  • Variation in transposable elements is often difficult and unreliable using short-read data (i.e., Illumina).
  1. Let's install SnpEff (short for: "SNP Effect") to annotate where our variants occur.
conda create -n snpeff bioconda::snpeff
  1. Next, we need to set up the directory structure expected by snpEff and copy over our reference files.
mkdir -p data data/Fraxinus_pennsylvanica_PE-57
cp /lustre/isaac24/proj/UTK0505/module2_data/reference_genome/Fraxinus_pennsylvanica_PE-57.fasta data/Fraxinus_pennsylvanica_PE-57/sequences.fa
cp /lustre/isaac24/proj/UTK0505/module2_data/reference_genome/Fraxinus_pennsylvanica_PE-57.gff3 data/Fraxinus_pennsylvanica_PE-57/genes.gff

# Copy VCF
cp /lustre/isaac24/proj/UTK0505/module2_data/group1_results/04_bcftools/group1.filteredMQ40DP10.vcf.gz
  1. To build the database, we need to create a snpEff.config file. We can accomplish this by opening our text editor and copying the following text into the file. It must be named snpEff.config.
# snpEff Config

data.dir = ./data/
Fraxinus_pennsylvanica_PE-57.genome : Fraxinus_pennsylvanica_PE-57
  1. Now, we can build our snpEff database.
# Run snpEff database build
# Note: -noCheckProtein and -noCheckCds should be used with caution in actual data analysis.
# It is preferable to include a CDS and protein file, so exon boundaries can be verified.
snpEff build \
        -c ${PWD}/snpEff.config \
        -gff3 \
        -v Fraxinus_pennsylvanica_PE-57 \
        -noCheckProtein \
        -noCheckCds
  1. Finally, we can annotate our VCF with our structural gene annotations and relevant amino acid changes induced by our SNPs.
# Now, annotate the VCF using your database.
snpEff \
        Fraxinus_pennsylvanica_PE-57 \
        <YOUR_GROUP>.filteredMQ40DP10.vcf.gz \
        > <YOUR_GROUP>.filteredMQ40DP10_annotated.vcf.gz
  1. Now, let's do a quick sanity check to look at a variant within a genic region.
zgrep -v 'intergenic' group1.filteredMQ40DP10_annotated.vcf.gz | head -n 80 | tail -n 1

Syntax explanation: grep - lets you search for strings in a file. zgrep - lets you search for strings in a gzipped file. -v - this flag searches for rows NOT matching the search pattern.

  1. snpEff also produces a nice HTML report that you can view here: snpEff_summary.html

02. Combine VCFs

Sometimes, you will sequence additional samples for a project. In that case, you may find you want to only call SNPs on the new individuals and merge those calls with the previous VCF. Typically, I would suggest recalling all of your SNPs together to reduce batch effects from false monomorphic loci--loci that vary in the population, but not in the samples you are calling variants on.

As a caveat, the 5 sample sets (one per group) have too few samples to robustly combine between runs. The more samples you have in a run, the more accurately you can detect variation across a population, and the lower the chance that you will have false negative SNP errors (site where variation occurs, but the variant caller did not capture it due to insufficient individuals, read coverage, or read depth).

Merging VCFs after filtering is also typically a very bad idea because you can introduce additional filtering batch effects. For some genotyping methods--particularly GATK which performs a Bayesian joint-genotyping procedure across all samples--it is advisable to recall SNPs on your entire cohort as you add new samples, not merge VCFs.

TL;DR: Merge VCFs with the utmost caution. Caveat emptor.

  1. First, make sure you're in your conda environment!
conda activate snps
  1. Let's grab the raw VCFs from all 5 groups.
ln -s /lustre/isaac24/proj/UTK0505/module2_data/group?_results/04_bcftools/*.raw.vcf.gz* .
  1. Now, we need to index them with bcftools, which we can do with the following in-line loop. Note that there is NO semicolon after do! By echoing the file name in our loop, we can track the loop's progress.
for VCF in *vcf.gz; do echo $VCF; bcftools index $VCF; done
  1. To merge these VCFs, we'll also use bcftools.
bcftools merge -m both -Oz -o merged.vcf.gz *.raw.vcf.gz
  • -m both makes sure to merge any variants (SNPs or indels) into a single locus at each position.
  • -Oz ensures the output is gzipped.
  • -o specifies your outfile name.
  1. Next, we would filter all of the samples in the merged VCF.
bcftools filter -i 'type="snp" && QUAL>=40 && FORMAT/DP>=10 && F_MISSING>=0.90' merged.vcf.gz > merged.QUAL40-DP10-FM90.vcf.gz
  1. To finish up, you should always index your VCFs!
bcftools index merged.vcf.gz
  1. Now, let's look at our VCF to see the perils of merging. We'll use less, so we can interactively scroll through without loading the entire file into memory and not have to gunzip the VCF.
less merged.QUAL40-DP10-FM90.vcf.gz
  • Image Example: Notice how the presence/absence of called alleles conspicuously tend to come in presence-absence blocks of 10 individuals (the number of individuals we had in each original VCF?)
image

03. SNP Ordination

Slides: https://docs.google.com/presentation/d/10CeNgG8AQG81MJdBoaLDNawafleOTdqpWTqCmpP0FDI/edit?slide=id.p#slide=id.p

  1. Time to make a conda environment! This will create the environment and install all dependencies in one step. This might take a few minutes...
conda create -n pca scikit-allel scikit-learn matplotlib numpy pandas plotly
  1. Link over your VCF from your group's directories. Replace <YOUR_GROUP> with your group.
ln -s /lustre/isaac24/proj/UTK0505/module2_data/<YOUR_GROUP>_results/04_bcftools/*.vcf.gz* .
  1. Now, we can run a PCA. The simplest way to accomplish this is using python. I've prepared a simplified form of my usual plotting technique with Claude for teaching purposes.

Open a text-editor on the server via nano and copy this script into a file called pcas.py.

import argparse
import sys
import numpy as np
import pandas as pd

import allel  # scikit-allel
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (needed for 3D projection)

# ---------------------------------------------------------------------------
# Parse arguments
# ---------------------------------------------------------------------------
parser = argparse.ArgumentParser(description="Basic 2D/3D PCA on a VCF file.")
parser.add_argument("vcf", help="Path to VCF or VCF.gz file")
parser.add_argument("--out", default="pca", help="Output file prefix (default: pca)")
parser.add_argument("--n-components", type=int, default=3,
                     help="Number of PCs to compute (default: 3, need >=3 for 3D plot)")
parser.add_argument("--keep-multiallelic", action="store_true",
                     help="Keep multiallelic/non-SNP sites instead of filtering them out")
parser.add_argument("--labels", default=None,
                     help="Optional CSV/TSV with sample,group columns for coloring points")
args = parser.parse_args()

# ---------------------------------------------------------------------------
# Load genotype matrix from VCF
# ---------------------------------------------------------------------------
print(f"Reading VCF: {args.vcf}")
callset = allel.read_vcf(
    args.vcf,
    fields=['samples', 'calldata/GT', 'variants/ALT', 'variants/is_snp']
)

if callset is None:
    sys.exit("No data read from VCF. Check the file path/format.")

samples = callset['samples']
gt = allel.GenotypeArray(callset['calldata/GT'])

# Optionally restrict to biallelic SNPs (ALT has exactly 1 allele)
if not args.keep_multiallelic:
    n_alt = (callset['variants/ALT'] != '').sum(axis=1)
    is_biallelic_snp = callset['variants/is_snp'] & (n_alt == 1)
    gt = gt[is_biallelic_snp]
    print(f"Kept {is_biallelic_snp.sum()} biallelic SNPs out of {len(is_biallelic_snp)} variants")

# Convert genotypes to alternate-allele counts per sample (0, 1, 2), missing -> -1
ac_matrix = gt.to_n_alt(fill=-1).astype(float)  # shape: (n_variants, n_samples)

# ---------------------------------------------------------------------------
# Clean and impute missing genotypes
# ---------------------------------------------------------------------------
ac_matrix[ac_matrix == -1] = np.nan

# Drop variants with no variance (uninformative for PCA) or all-missing
variant_std = np.nanstd(ac_matrix, axis=1)
valid = (variant_std > 0) & (~np.isnan(variant_std))
ac_matrix = ac_matrix[valid]
print(f"Retained {ac_matrix.shape[0]} variable variants for PCA")

# Mean-impute remaining missing values per variant
row_means = np.nanmean(ac_matrix, axis=1, keepdims=True)
nan_mask = np.isnan(ac_matrix)
ac_matrix[nan_mask] = np.take(row_means, np.where(nan_mask)[0])

# ---------------------------------------------------------------------------
# Run PCA
# ---------------------------------------------------------------------------
# Ensure enough PCs for the 3D plot, but never request more components than samples/features allow.
max_possible = min(ac_matrix.shape[1], ac_matrix.shape[0])  # samples x variants bound
n_components = max(args.n_components, 3, 10)
n_components = min(n_components, max_possible)

X = ac_matrix.T  # samples x variants

# Standardize each variant (zero mean, unit variance) - common for genotype PCA
X = (X - X.mean(axis=0)) / (X.std(axis=0) + 1e-9)

pca = PCA(n_components=n_components)
coords = pca.fit_transform(X)
explained = pca.explained_variance_ratio_ * 100

# ---------------------------------------------------------------------------
# Optional labels for coloring
# ---------------------------------------------------------------------------
labels = None
if args.labels:
    df = pd.read_csv(args.labels, sep=None, engine='python')
    df = df.set_index(df.columns[0])
    group_col = df.columns[0]
    labels = df.reindex(samples)[group_col].fillna('unknown').values

# ---------------------------------------------------------------------------
# 2D plot (PC1 vs PC2)
# ---------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(7, 6))

if labels is not None:
    for grp in pd.unique(labels):
        mask = labels == grp
        ax.scatter(coords[mask, 0], coords[mask, 1], label=str(grp), s=40, alpha=0.8)
    ax.legend(title="Group", fontsize=8)
else:
    ax.scatter(coords[:, 0], coords[:, 1], s=40, alpha=0.8)

ax.set_xlabel(f"PC1 ({explained[0]:.1f}% var)")
ax.set_ylabel(f"PC2 ({explained[1]:.1f}% var)")
ax.set_title("2D PCA of VCF Genotypes")
fig.tight_layout()
fig.savefig(f"{args.out}_2d.png", dpi=200)
print(f"Saved: {args.out}_2d.png")
plt.close(fig)

# ---------------------------------------------------------------------------
# 3D plot (PC1 vs PC2 vs PC3)
# ---------------------------------------------------------------------------
fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(111, projection='3d')

if labels is not None:
    for grp in pd.unique(labels):
        mask = labels == grp
        ax.scatter(coords[mask, 0], coords[mask, 1], coords[mask, 2],
                   label=str(grp), s=40, alpha=0.8)
    ax.legend(title="Group", fontsize=8)
else:
    ax.scatter(coords[:, 0], coords[:, 1], coords[:, 2], s=40, alpha=0.8)

ax.set_xlabel(f"PC1 ({explained[0]:.1f}%)")
ax.set_ylabel(f"PC2 ({explained[1]:.1f}%)")
ax.set_zlabel(f"PC3 ({explained[2]:.1f}%)")
ax.set_title("3D PCA of VCF Genotypes")
fig.tight_layout()
fig.savefig(f"{args.out}_3d.png", dpi=200)
print(f"Saved: {args.out}_3d.png")
plt.close(fig)


# ---------------------------------------------------------------------------
# Print explained variance for the first 10 PCs
# ---------------------------------------------------------------------------
n_report = min(10, len(explained))
print(f"\nExplained variance ratio, first {n_report} PCs (%):")
for i in range(n_report):
    print(f"  PC{i+1}: {explained[i]:.2f}%")

if len(explained) < 10:
    print(f"\nNote: only {len(explained)} PCs were available (limited by data dimensions).")
  1. Below, you'll find an explanation of how this script works. Python programming is outside the scope of this course, but feel free to ask questions if you're interested in learning more.
Basic 2D and 3D PCA on a VCF file.

Dependencies:
    pip install scikit-allel scikit-learn matplotlib numpy pandas

Usage:
    python pcas.py input.vcf.gz --out pca_output --n-components 3

Notes:
    - Missing genotypes are mean-imputed per variant.
    - Only biallelic SNPs are used by default (set --keep-multiallelic to skip filtering).
    - If sample population/group labels exist in a separate file (sample, group columns),
      pass it via --labels to color the plots.
    - The final report always shows explained variance for the first 10 PCs
      (regardless of --n-components used for plotting). This is important for knowing whether your
      PCA is capturing a sufficient amount of the variation in your data.
  1. Now, let's run our PCAs! We're going to run a 2D PCA and a 3D PCA. In the following images, you'll notice 2 dimensions doesn't do an excellent job of describing our data, but 3 dimensions does a better job.

I most commonly see 2D PCAs, but I generally see this as a mistake, unless it truly captures the dimensions of your data. In our case below, the second and third PC loadings explain nearly equivalent proportions of the variance of our data. This script will report the values of the first 10 PC loadings, so you can judge whether or not it's appropriate.

python3 pcas.py <your_group>.filteredMQ40DP10.vcf.gz
  • 2D PCA:
pca_2d
  • 3D PCA:
pca_3d
  • Variance Explained:
Explained variance ratio, first 10 PCs (%):
  PC1: 23.24%
  PC2: 17.26%
  PC3: 14.08%
  PC4: 9.60%
  PC5: 8.87%
  PC6: 8.13%
  PC7: 7.23%
  PC8: 6.20%
  PC9: 5.37%
  PC10: 0.00%

Bonus: Interactive PCAs with Plotly

  1. This is a bonus sections showing you how I generally prefer to visualize my SNP data. Here's the script to generate an interactive, HTML pca. The syntax for this script is the same, except it uses plotly instead of matplotlib.

I saved the following script to a file called pcas_plotly.py.

import argparse
import sys
import numpy as np
import pandas as pd

import allel  # scikit-allel
from sklearn.decomposition import PCA
import plotly.express as px
import plotly.graph_objects as go

# ---------------------------------------------------------------------------
# Parse arguments
# ---------------------------------------------------------------------------
parser = argparse.ArgumentParser(description="Basic 2D/3D PCA on a VCF file.")
parser.add_argument("vcf", help="Path to VCF or VCF.gz file")
parser.add_argument("--out", default="pca", help="Output file prefix (default: pca)")
parser.add_argument("--n-components", type=int, default=3,
                     help="Number of PCs to compute (default: 3, need >=3 for 3D plot)")
parser.add_argument("--keep-multiallelic", action="store_true",
                     help="Keep multiallelic/non-SNP sites instead of filtering them out")
parser.add_argument("--labels", default=None,
                     help="Optional CSV/TSV with sample,group columns for coloring points")
args = parser.parse_args()

# ---------------------------------------------------------------------------
# Load genotype matrix from VCF
# ---------------------------------------------------------------------------
print(f"Reading VCF: {args.vcf}")
callset = allel.read_vcf(
    args.vcf,
    fields=['samples', 'calldata/GT', 'variants/ALT', 'variants/is_snp']
)

if callset is None:
    sys.exit("No data read from VCF. Check the file path/format.")

samples = callset['samples']
gt = allel.GenotypeArray(callset['calldata/GT'])

# Optionally restrict to biallelic SNPs (ALT has exactly 1 allele)
if not args.keep_multiallelic:
    n_alt = (callset['variants/ALT'] != '').sum(axis=1)
    is_biallelic_snp = callset['variants/is_snp'] & (n_alt == 1)
    gt = gt[is_biallelic_snp]
    print(f"Kept {is_biallelic_snp.sum()} biallelic SNPs out of {len(is_biallelic_snp)} variants")

# Convert genotypes to alternate-allele counts per sample (0, 1, 2), missing -> -1
ac_matrix = gt.to_n_alt(fill=-1).astype(float)  # shape: (n_variants, n_samples)

# ---------------------------------------------------------------------------
# Clean and impute missing genotypes
# ---------------------------------------------------------------------------
ac_matrix[ac_matrix == -1] = np.nan

# Drop variants with no variance (uninformative for PCA) or all-missing
variant_std = np.nanstd(ac_matrix, axis=1)
valid = (variant_std > 0) & (~np.isnan(variant_std))
ac_matrix = ac_matrix[valid]
print(f"Retained {ac_matrix.shape[0]} variable variants for PCA")

# Mean-impute remaining missing values per variant
row_means = np.nanmean(ac_matrix, axis=1, keepdims=True)
nan_mask = np.isnan(ac_matrix)
ac_matrix[nan_mask] = np.take(row_means, np.where(nan_mask)[0])

# ---------------------------------------------------------------------------
# Run PCA
# ---------------------------------------------------------------------------
# Ensure enough PCs for the 3D plot AND for the final 10-PC report,
# but never request more components than samples/features allow.
max_possible = min(ac_matrix.shape[1], ac_matrix.shape[0])  # samples x variants bound
n_components = max(args.n_components, 3, 10)
n_components = min(n_components, max_possible)

X = ac_matrix.T  # samples x variants

# Standardize each variant (zero mean, unit variance) - common for genotype PCA
X = (X - X.mean(axis=0)) / (X.std(axis=0) + 1e-9)

pca = PCA(n_components=n_components)
coords = pca.fit_transform(X)
explained = pca.explained_variance_ratio_ * 100

# ---------------------------------------------------------------------------
# Optional labels for coloring
# ---------------------------------------------------------------------------
labels = None
if args.labels:
    df = pd.read_csv(args.labels, sep=None, engine='python')
    df = df.set_index(df.columns[0])
    group_col = df.columns[0]
    labels = df.reindex(samples)[group_col].fillna('unknown').values

# ---------------------------------------------------------------------------
# Save coordinates to CSV
# ---------------------------------------------------------------------------
pc_cols = [f"PC{i+1}" for i in range(coords.shape[1])]
out_df = pd.DataFrame(coords, columns=pc_cols)
out_df.insert(0, "sample", samples)
if labels is not None:
    out_df.insert(1, "group", labels)
out_df.to_csv(f"{args.out}_coordinates.csv", index=False)
print(f"Saved: {args.out}_coordinates.csv")

# ---------------------------------------------------------------------------
# 2D interactive plot (PC1 vs PC2)
# ---------------------------------------------------------------------------
plot_df = out_df.copy()
color_col = "group" if labels is not None else None

fig_2d = px.scatter(
    plot_df, x="PC1", y="PC2",
    color=color_col,
    hover_name="sample",
    labels={
        "PC1": f"PC1 ({explained[0]:.1f}% var)",
        "PC2": f"PC2 ({explained[1]:.1f}% var)",
    },
    title="2D PCA of VCF Genotypes",
)
fig_2d.update_traces(marker=dict(size=9, opacity=0.8))
fig_2d.write_html(f"{args.out}_2d.html")
print(f"Saved: {args.out}_2d.html")

# ---------------------------------------------------------------------------
# 3D interactive plot (PC1 vs PC2 vs PC3)
# ---------------------------------------------------------------------------
fig_3d = px.scatter_3d(
    plot_df, x="PC1", y="PC2", z="PC3",
    color=color_col,
    hover_name="sample",
    labels={
        "PC1": f"PC1 ({explained[0]:.1f}%)",
        "PC2": f"PC2 ({explained[1]:.1f}%)",
        "PC3": f"PC3 ({explained[2]:.1f}%)",
    },
    title="3D PCA of VCF Genotypes",
)
fig_3d.update_traces(marker=dict(size=5, opacity=0.8))
fig_3d.write_html(f"{args.out}_3d.html")
print(f"Saved: {args.out}_3d.html")

# ---------------------------------------------------------------------------
# Print explained variance for the first 10 PCs
# ---------------------------------------------------------------------------
n_report = min(10, len(explained))
print(f"\nExplained variance ratio, first {n_report} PCs (%):")
for i in range(n_report):
    print(f"  PC{i+1}: {explained[i]:.2f}%")

if len(explained) < 10:
    print(f"\nNote: only {len(explained)} PCs were available (limited by data dimensions).")
  1. Now, we can run this script with the same syntax.
python3 pcas_plotly.py <your_group>.filteredMQ40DP10.vcf.gz
  1. And voila, here are the output Files for your convenience:
⚠️ **GitHub.com Fallback** ⚠️