Bioinformatics Fundamentals

Bioinformatics Fundamentals

Foundation knowledge for genomics and bioinformatics workflows. Provides essential understanding of file formats, sequencing technologies, and common data processing patterns.

When to Use This Skill

• Working with sequencing data (PacBio HiFi, Hi-C, Illumina)
• Debugging SAM/BAM alignment or filtering issues
• Processing AGP files for genome assembly curation
• Validating AGP coordinate systems and unloc assignments
• Understanding paired-end vs single-end data
• Interpreting quality metrics (MAPQ, PHRED scores)
• Troubleshooting empty outputs or broken read pairs
• Accessing GenomeArk QC data (GenomeScope, BUSCO, Merqury)
• Curating karyotype data or chromosome count analysis
• General bioinformatics data analysis

SAM/BAM Format Essentials

SAM Flags (Bitwise)

Flags are additive - a read can have multiple flags set simultaneously.

Common Flags:

• 0x0001 (1): Read is paired in sequencing
• 0x0002 (2): Each segment properly aligned (proper pair)
• 0x0004 (4): Read unmapped
• 0x0008 (8): Mate unmapped
• 0x0010 (16): Read mapped to reverse strand
• 0x0020 (32): Mate mapped to reverse strand
• 0x0040 (64): First in pair (R1/forward)
• 0x0080 (128): Second in pair (R2/reverse)
• 0x0100 (256): Secondary alignment
• 0x0400 (1024): PCR or optical duplicate
• 0x0800 (2048): Supplementary alignment

Flag Combinations:

• Properly paired R1: 99 (0x63 = 1 + 2 + 32 + 64)
• Properly paired R2: 147 (0x93 = 1 + 2 + 16 + 128)
• Unmapped read: 4
• Mate unmapped: 8

See reference.md for complete flag tables, CIGAR operations, optional tags, and SAM mandatory fields.

Proper Pair Flag (0x0002)

What "proper pair" means:

• Both R1 and R2 are mapped
• Mapping orientations are correct (typically R1 forward, R2 reverse)
• Insert size is reasonable for the library
• Pair conforms to aligner's expectations

Important: Different aligners have different criteria for proper pairs!

MAPQ (Mapping Quality)

Formula: MAPQ = -10 * log10(P(mapping is wrong))

Common Thresholds:

• MAPQ >= 60: High confidence (error probability < 0.0001%)
• MAPQ >= 30: Good quality (error probability < 0.1%)
• MAPQ >= 20: Acceptable (error probability < 1%)
• MAPQ >= 10: Low confidence (error probability < 10%)
• MAPQ = 0: Multi-mapper or unmapped

Note: MAPQ=0 can mean either unmapped OR equally good multiple mappings.

CIGAR String

Represents alignment between read and reference:

• M: Match or mismatch (alignment match)
• I: Insertion in read vs reference
• D: Deletion in read vs reference
• S: Soft clipping (bases in read not aligned)
• H: Hard clipping (bases not in read sequence)
• N: Skipped region (for RNA-seq splicing)

Example: 100M = perfect 100bp match Example: 50M5I45M = 50bp match, 5bp insertion, 45bp match

Sequencing Technologies

PacBio HiFi (High Fidelity)

Characteristics:

• Long reads: 10-25 kb typical
• High accuracy: >99.9% (Q20+)
• Circular Consensus Sequencing (CCS)
• Single-end data (though from circular molecules)
• Excellent for de novo assembly

Best Mappers:

• minimap2 presets: map-pb, map-hifi
• BWA-MEM2 can work but optimized for short reads

Typical Use Cases:

• De novo genome assembly
• Structural variant detection
• Isoform sequencing (Iso-Seq)
• Haplotype phasing

Hi-C (Chromatin Conformation Capture)

Characteristics:

• Paired-end short reads (typically 100-150 bp)
• Read pairs capture chromatin interactions
• R1 and R2 often map to different scaffolds/chromosomes
• Requires careful proper pair handling
• Used for scaffolding and 3D genome structure

Best Mappers:

• BWA-MEM2 (paired-end mode)
• BWA-MEM (paired-end mode)

Critical Concept: Hi-C read pairs intentionally map to distant loci. Region filtering can easily break pairs!

Typical Use Cases:

• Genome scaffolding (connecting contigs)
• 3D chromatin structure analysis
• Haplotype phasing
• Assembly quality assessment

Illumina Short Reads

Characteristics:

• Short reads: 50-300 bp
• Paired-end or single-end
• High throughput
• Well-established quality scores

Best Mappers:

• BWA-MEM2, BWA-MEM (general purpose)
• Bowtie2 (fast, local alignment)
• STAR (RNA-seq spliced alignment)

See reference.md for detailed technical specs, error profiles, and quality metrics per technology.

Common Tools and Their Behaviors

samtools view

Purpose: Filter, convert, and view SAM/BAM files

Key Flags:

• -b: Output BAM format
• -h: Include header
• -f INT: Require flags (keep reads WITH these flags)
• -F INT: Filter flags (remove reads WITH these flags)
• -q INT: Minimum MAPQ threshold
• -L FILE: Keep reads overlapping regions in BED file

Important Behavior:

• -L (region filtering) checks each read individually, not pairs
• Can break read pairs if mates map to different regions
• Flag filters (-f, -F) are applied before region filters (-L)

Example - Proper pairs in regions (correct order):

samtools view -b -f 2 -L regions.bed input.bam > proper_pairs_in_regions.bam

bamtools filter

Purpose: Advanced filtering with complex criteria

Common Filters:

• isPaired: true - Read is from paired-end sequencing
• isProperPair: true - Read is part of proper pair
• isMapped: true - Read is mapped
• mapQuality: >=30 - Mapping quality threshold

Important Difference from samtools:

• isProperPair is more strict than samtools -f 2
• Checks pair validity more thoroughly

samtools fastx

Purpose: Convert SAM/BAM to FASTQ/FASTA

Critical: Use appropriate filters to ensure R1/R2 files match!

See reference.md for complete tool command reference with all options and examples.

Common Patterns and Best Practices

Pattern 1: Filtering Paired-End Data by Regions

WRONG WAY (breaks pairs):

# Region filter first -> breaks pairs when mates are in different regions
samtools view -b -L regions.bed input.bam | bamtools filter -isPaired -isProperPair
# Result: Empty output (all pairs broken)

RIGHT WAY (preserves pairs):

# Proper pair filter FIRST, then region filter
samtools view -b -f 2 -L regions.bed input.bam > output.bam

Pattern 2: Extracting FASTQ from Filtered BAM

For Paired-End:

samtools fastx -1 R1.fq.gz -2 R2.fq.gz \
  --i1-flags 2 \  # Require proper pair
  input.bam

For Single-End:

samtools fastx -0 output.fq.gz input.bam

Pattern 3: Quality Filtering

Conservative (high quality):

samtools view -b -q 30 -f 2 -F 256 -F 2048 input.bam
# MAPQ >= 30, proper pairs, no secondary/supplementary

Permissive (for low-coverage data):

samtools view -b -q 10 -F 4 input.bam
# MAPQ >= 10, mapped reads

Common Issues Summary

Issue 1: Empty Output After Region Filtering (Hi-C Data)

Region filter (samtools view -L) breaks read pairs. One mate in region, other outside. Proper pair flag lost. Apply proper pair filter BEFORE region filtering:

samtools view -b -f 2 -L regions.bed input.bam > output.bam

Issue 2: R1 and R2 Files Have Different Read Counts

Improper filtering broke some pairs. Require proper pairs during extraction:

samtools fastx -1 R1.fq -2 R2.fq --i1-flags 2 input.bam

Issue 3: Low Mapping Rate for Hi-C Data

This is normal for Hi-C due to chimeric reads. Use Hi-C-specific pipelines (HiC-Pro, Juicer). Don't filter too aggressively on MAPQ.

Issue 4: Proper Pairs Lost After Mapping

Check insert size distribution, reference mismatch, or incorrect orientation flags.

samtools stats input.bam | grep "insert size"
samtools flagstat input.bam

See common-issues.md for comprehensive troubleshooting with detailed solutions, including AGP processing issues, HiFi-specific problems, and diagnostic commands.

Quality Metrics