Large-Scale Genomic Data Processing

Step 1 — Install Homebrew (if not already installed)

# Paste this into Terminal — it installs Homebrew
/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

# Verify installation
brew --version

Step 2 — Install Core Bioinformatics Tools

# samtools — reads, indexes, and queries genome FASTA and BAM files
brew install samtools

# bcftools — reads, queries, and filters VCF variant call files
brew install bcftools

# htslib — provides bgzip (block gzip) and tabix (indexed compression)
brew install htslib

# bwa — aligns short DNA reads to a reference genome
brew install bwa

# Java — required by snpEff (annotation engine)
brew install openjdk@17
echo 'export PATH="/opt/homebrew/opt/openjdk@17/bin:$PATH"' >> ~/.zshrc
source ~/.zshrc

# Verify Java is working
java -version

Step 3 — Install snpEff (Variant Annotator)

snpEff is not on Homebrew — install it manually. It is a Java program distributed as a ZIP file.

# Create a working directory
mkdir -p ~/na12878_pipeline && cd ~/na12878_pipeline

# Download snpEff
curl -L https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip -o snpEff.zip

# Unzip it
unzip snpEff.zip

# Download the GRCh38 annotation database (one-time, ~3 GB)
java -jar snpEff/snpEff.jar download GRCh38.109

# Verify snpEff works
java -jar snpEff/snpEff.jar -version

Step 4 — Install Miniconda (for snpEff conda environment — optional)

# Download Miniconda installer for Apple Silicon
curl -O https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-arm64.sh

# Run installer
bash Miniconda3-latest-MacOSX-arm64.sh

# Restart terminal, then create a dedicated snpEff environment
conda create -n snpeff_env -c bioconda snpeff -y
conda activate snpeff_env

Step 5 — Install Python Visualization Tools

# Install pip packages for plotting
pip3 install matplotlib pandas numpy biopython

# Verify
python3 -c "import matplotlib; print('matplotlib', matplotlib.__version__)"

Full Tool Reference

Option A — Download GRCh38 from NCBI (Authoritative Source)

mkdir -p ~/na12878_pipeline && cd ~/na12878_pipeline

# Full primary assembly — all chromosomes — ~900 MB compressed (~3 GB unzipped)
curl -O https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/GCA_000001405.15_GRCh38_assembly_structure/Primary_Assembly/assembled_chromosomes/FASTA/chroms/GRCh38_major_release_seqs_for_alignment_pipelines.fa.gz

# OR — faster mirror from Ensembl
curl -O https://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# OR — UCSC hg38 (same genome, different chromosome naming: chr1 vs 1)
rsync -avzP rsync://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz .

Option B — Download NA12878 Real Variant Calls (GIAB Gold Standard)

cd ~/na12878_pipeline

# High-confidence variant VCF for NA12878 against GRCh38
curl -O https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/NA12878_HG001/NISTv3.3.2/GRCh38/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz

# Download the index file too (required for tabix random access)
curl -O https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/NA12878_HG001/NISTv3.3.2/GRCh38/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz.tbi

# Rename to something shorter for convenience
mv HG001_GRCh38_GIAB_highconf_*.vcf.gz na12878.vcf.gz
mv HG001_GRCh38_GIAB_highconf_*.vcf.gz.tbi na12878.vcf.gz.tbi

Index the Reference Genome for Fast Querying

Before you can query specific coordinates, the genome file must be indexed. This creates a lookup table so tools can jump directly to any position without reading the entire file.

# Step 1 — Block-gzip compress the reference (if not already bgzip'd)
bgzip -c Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz > ref_grch38.fna.bgz

# Step 2 — Build the FASTA index (.fai file)
samtools faidx ref_grch38.fna.bgz

# This creates ref_grch38.fna.bgz.fai — a lookup table of chromosome positions

# Step 3 — Index the VCF for random access
tabix -p vcf na12878.vcf.gz

# Step 4 — Index VCF with bcftools (alternative, creates .csi index)
bcftools index na12878.vcf.gz

File Sizes to Expect

1. Stream a Compressed File Without Extracting It

gzcat reads a .gz file and outputs its contents as a text stream — no decompression to disk. This is the foundation of every genomic query on macOS.

# Stream the entire genome into your terminal (DO NOT run alone — will print forever)
gzcat ref_grch38.fna.bgz

# Always pipe gzcat into another command that controls the output
gzcat ref_grch38.fna.bgz | head -n 5        # first 5 lines only
gzcat ref_grch38.fna.bgz | wc -l             # count total lines (slow — reads all 38M)
gzcat ref_grch38.fna.bgz | grep "^>"         # show only chromosome headers

2. Peek at the File Without Opening It

# Show first 10 lines
gzcat ref_grch38.fna.bgz | head -n 10

# Show last 10 lines (reads entire file — slow on 3 GB)
gzcat ref_grch38.fna.bgz | tail -n 10

# Show first 20 lines of a specific chromosome (chr22)
gzcat ref_grch38.fna.bgz | sed -n '/^>chr22/,/^>/p' | head -n 20

# Count how many chromosomes are in the file
gzcat ref_grch38.fna.bgz | grep -c "^>"

3. samtools faidx — Query Exact Coordinates (Fastest Method)

samtools faidx uses the index to jump directly to any genomic coordinate in milliseconds — no reading from the start. This is how you query specific genes or regions.

# Pull a single base at a specific position
samtools faidx ref_grch38.fna.bgz chr17:43094464-43094464

# Pull the entire BRCA1 gene region
samtools faidx ref_grch38.fna.bgz chr17:43044295-43125483

# Pull a region from chromosome 22
samtools faidx ref_grch38.fna.bgz chr22:1-50000

# Pull multiple regions at once
samtools faidx ref_grch38.fna.bgz chr17:43044295-43125483 chr22:1-100000

# List all chromosomes and their sizes (from the .fai index)
cut -f1,2 ref_grch38.fna.bgz.fai

# Output format: FASTA header + sequence
# >chr17:43044295-43125483
# ATCGATCGATCG...

4. bcftools — Query the Variant Database (VCF File)

The VCF file contains 4.2 million positions where NA12878's DNA differs from the reference. bcftools is the primary tool for querying it.

# ─── BASIC QUERIES ───────────────────────────────

# Show VCF header (metadata about the file)
bcftools view -h na12878.vcf.gz | less

# Show first 10 variant records
bcftools view na12878.vcf.gz | grep -v "^#" | head -10

# Count total number of variants
bcftools view na12878.vcf.gz | grep -vc "^#"

# ─── FILTER BY CHROMOSOME / REGION ──────────────

# All variants on chromosome 22
bcftools view -r chr22 na12878.vcf.gz | less

# Variants in a specific gene region (BRCA1)
bcftools view -r chr17:43044295-43125483 na12878.vcf.gz

# Save region to a new file
bcftools view -r chr17:43044295-43125483 na12878.vcf.gz -o brca1_variants.vcf

# Multiple regions at once
bcftools view -r chr17:43044295-43125483,chr22:1-5000000 na12878.vcf.gz

# ─── FILTER BY VARIANT TYPE ──────────────────────

# Show only SNPs (single nucleotide changes)
bcftools view -v snps na12878.vcf.gz | head -20

# Show only indels (insertions + deletions)
bcftools view -v indels na12878.vcf.gz | head -20

# Count SNPs
bcftools view -v snps na12878.vcf.gz | grep -vc "^#"

# Count indels
bcftools view -v indels na12878.vcf.gz | grep -vc "^#"

# ─── STATISTICS ──────────────────────────────────

# Full stats report (ts/tv ratio, indel sizes, quality distribution)
bcftools stats na12878.vcf.gz > na12878_stats.txt
less na12878_stats.txt

# Quick summary stats
bcftools stats na12878.vcf.gz | grep "^SN"

# ─── ADVANCED QUERIES ────────────────────────────

# Extract only the chromosome, position, ref, alt columns
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' na12878.vcf.gz | head -20

# Count insertions (ALT longer than REF)
bcftools view -v indels na12878.vcf.gz | awk '!/^#/ && length($5) > length($4)' | wc -l

# Count deletions (REF longer than ALT)
bcftools view -v indels na12878.vcf.gz | awk '!/^#/ && length($4) > length($5)' | wc -l

# Find the largest indels (most disruptive)
bcftools view -v indels na12878.vcf.gz | awk '!/^#/ {print length($4)-length($5), $0}' | sort -n | tail -10

# Check if a specific position has a variant
bcftools view -r chr17:43094464-43094464 na12878.vcf.gz

5. tabix — Indexed Random Access (Jump to Any Position Instantly)

tabix is like a GPS for a compressed genomic file. Once a file is indexed, it can seek directly to any chromosome:position without reading from the start.

# Build a tabix index on a VCF (creates .tbi file)
tabix -p vcf na12878.vcf.gz

# Query a region directly — millisecond access on a 450 MB file
tabix na12878.vcf.gz chr17:43044295-43125483

# Query chromosome 22 entirely
tabix na12878.vcf.gz chr22

# List all chromosomes in the file
tabix -l na12878.vcf.gz

# Build a tabix index on a FASTA file
tabix -p fasta genome.fasta.gz

6. grep — Pattern Search Inside Any Genomic File

# Find all chromosome headers in FASTA
gzcat ref_grch38.fna.bgz | grep "^>"

# Find a specific chromosome header
gzcat ref_grch38.fna.bgz | grep "^>chr22"

# Show chr22 header + the next 5 lines of sequence
gzcat ref_grch38.fna.bgz | grep "^>chr22" -A 5

# Count how many lines contain a specific motif (e.g. BRCA1 gene ID)
bcftools view na12878.vcf.gz | grep -c "BRCA1"

# Search VCF for a specific rs-ID (dbSNP variant identifier)
bcftools view na12878.vcf.gz | grep "rs28897672"

# Count lines NOT starting with # (count actual variants, skip header)
bcftools view na12878.vcf.gz | grep -vc "^#"

7. sed — Stream Filter by Range

sed is the tool that isolates one chromosome from the 38-million-line genome stream.

# Isolate chromosome 22 from start to the next chromosome header
gzcat ref_grch38.fna.bgz | sed -n '/^>chr22/,/^>chr[^2]/p' | head -20

# How it works:
#   -n          suppress default printing
#   /^>chr22/   start printing when this pattern is matched
#   ,/^>chr/    stop printing when ANY new chromosome header is matched
#   p           print the matching range

# Isolate lines between two patterns and save to file
gzcat ref_grch38.fna.bgz | sed -n '/^>chr22/,/^>chrX/p' > chr22_only.fasta

8. awk — Column Math on Genomic Data

awk treats each line as columns separated by tabs. VCF files have 8+ tab-separated columns — awk lets you calculate on them.

# VCF column structure:
# $1=CHROM  $2=POS  $3=ID  $4=REF  $5=ALT  $6=QUAL  $7=FILTER  $8=INFO

# Print chromosome and position of every variant
bcftools view na12878.vcf.gz | awk '!/^#/ {print $1, $2}'

# Calculate indel size (positive = deletion, negative = insertion)
bcftools view -v indels na12878.vcf.gz | awk '!/^#/ {print length($4)-length($5), $1, $2}'

# Show only variants with quality score above 50
bcftools view na12878.vcf.gz | awk '!/^#/ && $6 > 50'

# Count variants per chromosome
bcftools view na12878.vcf.gz | awk '!/^#/ {count[$1]++} END {for(c in count) print c, count[c]}' | sort -V

9. snpEff — Annotate What Variants Actually Do

# Annotate chr22 only (memory-safe, ~3 GB RAM)
bcftools view -r chr22 na12878.vcf.gz | \
  java -Xmx8g -jar snpEff/snpEff.jar GRCh38.109 - > chr22_annotated.vcf

# Annotate full genome (requires 16 GB free RAM)
java -Xmx16g -jar snpEff/snpEff.jar GRCh38.109 na12878.vcf.gz > na12878_annotated.vcf 2>snpeff.log

# Monitor progress while running
tail -f snpeff.log

# After annotation — compress and index the output
bgzip na12878_annotated.vcf
tabix -p vcf na12878_annotated.vcf.gz

# Query annotated file for missense variants only
bcftools view na12878_annotated.vcf.gz | grep "missense_variant"

# Query annotated file for BRCA1 gene
grep BRCA1 snpEff_genes.txt

# Open the HTML summary report in browser
open snpEff_summary.html

10. Full Pipeline — From Download to Annotation in One Flow

# ─── COMPLETE PIPELINE ───────────────────────────────────────────
mkdir -p ~/na12878_pipeline && cd ~/na12878_pipeline

# 1. Download reference genome
curl -O https://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 2. Download NA12878 variant calls
curl -O https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/NA12878_HG001/NISTv3.3.2/GRCh38/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz
curl -O https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/NA12878_HG001/NISTv3.3.2/GRCh38/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz.tbi
mv HG001_GRCh38_GIAB_highconf_*.vcf.gz na12878.vcf.gz
mv HG001_GRCh38_GIAB_highconf_*.vcf.gz.tbi na12878.vcf.gz.tbi

# 3. Compress and index the reference
bgzip -c Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz > ref_grch38.fna.bgz
samtools faidx ref_grch38.fna.bgz

# 4. Index the VCF
bcftools index na12878.vcf.gz

# 5. Quick sanity check — list chromosomes, count variants
tabix -l na12878.vcf.gz
bcftools view na12878.vcf.gz | grep -vc "^#"

# 6. Annotate chr22 (memory-safe first run)
bcftools view -r chr22 na12878.vcf.gz | \
  java -Xmx8g -jar snpEff/snpEff.jar GRCh38.109 - > chr22_annotated.vcf

# 7. View results
open snpEff_summary.html

The CEPH Collection (1980s)

NA12878's story begins with the Centre d'Etude du Polymorphisme Humain (CEPH), a French genetics institute founded in Paris in 1983. CEPH recruited large three-generation families and collected DNA samples from grandparents, parents, and children. These families, called the CEPH pedigrees, gave scientists the ability to track how genetic variants are inherited across generations — a process called linkage analysis.

NA12878 belongs to CEPH family 1463. The "NA" prefix stands for North American. The "12878" is her permanent, anonymized lab identifier. She is one member of a multi-generation family that donated samples together, allowing scientists to verify which variants she inherited from which parent.

Genome in a Bottle (GIAB) — The Gold Standard

In 2011, the National Institute of Standards and Technology (NIST) launched the Genome in a Bottle Consortium. Their mission: create a human genome reference so deeply verified that it could be trusted as a universal benchmark. NA12878 was selected as the first subject.

Her genome was sequenced using every major sequencing technology available — Illumina short reads, Pacific Biosciences long reads, Oxford Nanopore, Complete Genomics, Ion Torrent, and more — and the results were cross-validated against each other to produce a high-confidence variant call set: a list of every position in her genome where she differs from the reference, with the highest confidence rating in the field.

Her Scientific Legacy

NA12878's genome has been downloaded millions of times. It is used to:

• Validate every new DNA sequencing machine before it goes to market
• Test and benchmark bioinformatics software (aligners, variant callers, annotation tools)
• Train machine learning models that detect cancer mutations and rare diseases
• Calibrate clinical laboratory pipelines before they process real patient samples
• Serve as a teaching resource — exactly as in this project

Every time a hospital sequences a cancer patient's tumor to find the right treatment, the pipeline that produced that result was likely validated against NA12878 first. Her contribution is immeasurable.

The Three Build Components

What the Central Database Contains

The variant calls in na12878.vcf.gz are cross-referenced against major public databases:

Autosomal Ancestry Distribution (Assigned Segments Only)

NA12878 vs. Anonymized Profile — Key Differences

How the Comparison Pipeline Works

When comparing this profile's chromosome segments to the central database, the process follows these steps:

1. Phased segment extraction: 23andMe's internal algorithm phases the ~650,000 SNP positions into maternal and paternal haplotypes (Copy 1 / Copy 2)
2. Reference panel matching: Each phased segment is compared against 31 reference populations using a Hidden Markov Model (HMM) — the same statistical method used in scientific population genetics
3. Boundary detection: Recombination breakpoints are identified where population signatures change — these are the segment boundaries in the coordinate ledger
4. Confidence thresholding: At 0.5 threshold, segments are only assigned when the top reference panel match has at least 50% posterior probability — tighter than the default 0.3 setting
5. Gap classification: Regions where no panel exceeds the threshold are marked Unassigned — these are the 8 gap windows in the dataset

Chromosome 22 Variant Classification — From Least to Most Significant

snpEff assigns each variant a functional consequence based on where it falls in the genome. Severity increases from top to bottom.

Visual Outputs Generated

Critical Gene Comparison Results

The pipeline also extracted and compared variants in four clinically significant genes against the reference genome:

Commands: View Results & Scale to Full Genome

# Open all generated PNG graphs
open ~/na12878_pipeline/*.png

# ─────────────────────────────────────────────────────
# To run snpEff on the FULL 3.6M variant dataset:
# Restart Mac first to clear RAM, then:
# ─────────────────────────────────────────────────────
cd ~/na12878_pipeline

/opt/homebrew/Caskroom/miniconda/base/bin/conda run -n snpeff_env \
  snpEff -Xmx16g GRCh38.99 na12878.vcf.gz > na12878_annotated.vcf 2>snpeff.log

# -Xmx16g = allow snpEff to use up to 16 GB RAM
# Requires at least 16 GB free before running
# Monitor progress: tail -f snpeff.log

The zcat Problem on macOS

# Error you will see on macOS with the wrong command:
# zcat: can't stat: genome.fasta.gz (genome.fasta.gz.Z): No such file or directory

# Correct macOS commands:
gzcat genome.fasta.gz | head
gunzip -c genome.fasta.gz | head

Convert This Report to PDF

# Install pandoc and a minimal LaTeX engine
brew install pandoc basictex

# Convert the Markdown report to a PDF
pandoc GENOMIC_PIPELINE_REPORT.md -o Genomic_Operations_Report.pdf

Or open this HTML file in any browser and use Cmd + P → Save as PDF for a direct PDF export of this full reference page.