Use Case 1: Reference-Based Sequence Analysis

This type of analysis starts with DNA sequencing data in the form of "sequencing reads", typically stored in FASTQ files, which contain both nucleotide sequences and their per-base quality scores. These reads are then mapped (aligned) to a reference genome, a curated, representative genome sequence that serves as a coordinate system for the organism under study.

Workflow Overview

Read mapping (alignment) → Variant calling
FASTQ (+ reference FASTA) → BAM/SAM → VCF

Key Steps

1. Read Mapping (Alignment)

The mapping process produces an alignment file, commonly in SAM (or its binary equivalent, BAM) format, which records the individual reads and their positions on the reference genome. In addition to the aligned sequences, these files store two kinds of quality measurements:

The original base quality scores from the sequencing
A new mapping quality score, reflecting the confidence of a read being placed at the correct genomic location

2. Variant Identification

Using this information, differences between the sequenced sample and the reference genome can be identified. These differences may include:

Naturally occurring variation
Deliberately introduced genetic changes (e.g., targeted mutations or deletions generated by genome editing experiments)

Such changes in the genetic code include:

Single Nucleotide Variants (SNVs) or Single Nucleotide Polymorphisms (SNPs) when common in a population
Insertions or deletions (indels)

3. Variant Calling

Identifying these variants is a dedicated analysis step known as "variant calling", which evaluates the aligned reads in the BAM/SAM file. The result of variant calling is typically stored in a "Variant Call Format" (VCF) file, which summarizes:

The detected variations
Their genomic positions
The observed alleles
Quality and confidence metrics
How many sequencing reads support each variant

Applications

Reference-based variant analysis is widely used to:

Identify specific mutations
Determine the genotype of an organism associated with a particular phenotype
Investigate how genetic changes alter biological function and lead to observable traits, diseases, or experimental outcomes

Common applications include:

Microbial strain identification
Evolutionary genetics
Genetic variant discovery

Popular Tools

BLAST - Basic sequence alignment
BWA, BWA-MEM - Fast read alignment
STAR - RNA-seq aligner
Minimap2 - Long-read alignment
samtools - SAM/BAM file manipulation

Code Example: Standard Reference-Based Workflow

Click to expand Bash workflow

#!/bin/bash
# Reference-based variant calling workflow

# Step 1: Quality control of raw reads
fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o qc_results/

# Step 2: Align reads to reference genome using BWA-MEM
# Index reference genome (only needed once)
bwa index reference_genome.fasta

# Align paired-end reads
bwa mem -t 8 reference_genome.fasta \
    sample_R1.fastq.gz \
    sample_R2.fastq.gz \
    > aligned_reads.sam

# Step 3: Convert SAM to BAM and sort
samtools view -bS aligned_reads.sam | samtools sort -o aligned_reads.sorted.bam
samtools index aligned_reads.sorted.bam

# Step 4: Mark duplicate reads
picard MarkDuplicates \
    INPUT=aligned_reads.sorted.bam \
    OUTPUT=aligned_reads.dedup.bam \
    METRICS_FILE=metrics.txt

# Step 5: Call variants using bcftools
bcftools mpileup -f reference_genome.fasta aligned_reads.dedup.bam | \
    bcftools call -mv -Oz -o variants.vcf.gz

# Step 6: Filter variants
bcftools filter -s LowQual -e 'QUAL<20 || DP<10' variants.vcf.gz > variants.filtered.vcf

# Step 7: Summarize results
bcftools stats variants.filtered.vcf > variant_statistics.txt

echo "Variant calling complete! Results in variants.filtered.vcf"

Python Example: Parsing VCF Files

Click to expand Python script

#!/usr/bin/env python3
"""
Parse VCF file and extract variant information
"""

def parse_vcf(vcf_file):
    """
    Parse VCF file and extract variant positions and types
    
    Args:
        vcf_file: Path to VCF file
        
    Returns:
        List of dictionaries containing variant information
    """
    variants = []
    
    with open(vcf_file, 'r') as f:
        for line in f:
            # Skip header lines
            if line.startswith('#'):
                continue
                
            # Parse variant line
            fields = line.strip().split('\t')
            chrom = fields[0]
            pos = int(fields[1])
            ref = fields[3]
            alt = fields[4]
            qual = float(fields[5]) if fields[5] != '.' else None
            
            # Determine variant type
            if len(ref) == 1 and len(alt) == 1:
                var_type = 'SNV'
            elif len(ref) < len(alt):
                var_type = 'Insertion'
            elif len(ref) > len(alt):
                var_type = 'Deletion'
            else:
                var_type = 'Complex'
            
            variants.append({
                'chromosome': chrom,
                'position': pos,
                'reference': ref,
                'alternate': alt,
                'quality': qual,
                'type': var_type
            })
    
    return variants

# Usage example
if __name__ == '__main__':
    variants = parse_vcf('variants.filtered.vcf')
    
    # Print summary
    print(f"Total variants found: {len(variants)}")
    
    # Count variant types
    from collections import Counter
    var_types = Counter(v['type'] for v in variants)
    print("\nVariant type distribution:")
    for var_type, count in var_types.items():
        print(f"  {var_type}: {count}")

Expected Outputs

After running the complete workflow, you should have:

Quality reports - FastQC HTML reports showing read quality metrics
Alignment file - Sorted, indexed BAM file with aligned reads
Variant file - VCF file containing called variants with quality scores
Statistics - Summary statistics on alignment rates and variant counts

Interpreting Results

Alignment rate - Should typically be >90% for good quality samples with correct reference
Variant count - Depends on organism and expected polymorphism (bacterial: 10-1000s, human: 3-5 million SNPs)
Quality scores - PHRED scores >20 indicate >99% base call accuracy
Coverage depth - Aim for ≥30x for reliable variant calling in diploid organisms

Computational Requirements

Typical Resources Needed

Step	CPU Cores	RAM	Storage	Time (100M paired-end reads)
FastQC	2	2 GB	Minimal	10-20 min
BWA alignment	8-16	8-16 GB	~50 GB	2-4 hours
Variant calling	4-8	16-32 GB	~20 GB	1-3 hours

Where to Run

Local workstation - Suitable for small genomes (bacteria) or pilot studies
LRZ Linux Cluster - Recommended for human genomes or multiple samples
See: Infrastructure Guide

Common Issues & Troubleshooting

Common Problems

Low alignment rate (<80%)

Check if you're using the correct reference genome
Verify sample contamination or quality issues
Consider adapter contamination

Few or no variants called

Insufficient coverage depth
Too stringent filtering parameters
Sample may be very similar to reference

High memory usage

Use BWA-MEM instead of older BWA algorithms
Process chromosomes separately
Increase swap space or use cluster resources

Slow performance

Increase thread count (-t parameter)
Use indexed reference genomes
Consider pre-filtering low-quality reads

Sample Data & Practice

To practice this workflow, you can use publicly available datasets:

E. coli test data - ENA Browser
Human chromosome 22 - 1000 Genomes Project
NCBI SRA - Search for specific organisms

After variant calling, you might want to:

Variant annotation - Use tools like SnpEff or VEP to predict variant effects
Comparative analysis - Compare variants across multiple samples
Population genetics - Analyze allele frequencies and genetic diversity
Functional validation - Confirm important variants experimentally

Related Use Cases:

Transcriptomic Analysis - Study gene expression changes
De Novo Assembly - When no reference is available

Key Considerations

Quality control is essential - Always check read quality before and after alignment
Reference genome selection - Ensure you're using the appropriate reference for your organism
Coverage depth - Higher coverage provides more confident variant calls (≥30x recommended)
Filtering criteria - Adjust quality and depth thresholds based on your specific needs
Biological validation - Confirm critical variants with Sanger sequencing or other methods

Workflow Overview​

Key Steps​

1. Read Mapping (Alignment)​

2. Variant Identification​

3. Variant Calling​

Applications​

Popular Tools​

Code Example: Standard Reference-Based Workflow​

Python Example: Parsing VCF Files​

Expected Outputs​

Interpreting Results​

Computational Requirements​

Typical Resources Needed​

Where to Run​

Common Issues & Troubleshooting​

Sample Data & Practice​

Next Steps & Related Analyses​

Key Considerations​