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

N50 and Related Metrics

N50: Length of the shortest contig at which 50% of total assembly is contained in contigs of that length or longer

Related Metrics:

• L50: Number of contigs needed to reach N50
• N90: More stringent than N50 (90% coverage)
• NG50: N50 relative to genome size (better for comparisons)

Coverage and Depth

Coverage: Percentage of reference bases covered by at least one read Depth: Average number of reads covering each base

Recommended Depths:

• Genome assembly (HiFi): 30-50x
• Variant calling: 30x minimum
• RNA-seq: 20-40 million reads
• Hi-C scaffolding: 50-100x genomic coverage

See reference.md for complete coverage calculations, BUSCO interpretation, QV scores, and assembly quality metrics.

File Format Quick Reference

FASTA

>sequence_id description
ATCGATCGATCG

• Header line starts with >
• No quality scores

FASTQ

@read_id
ATCGATCGATCG
+
IIIIIIIIIIII

• Four lines per read
• Quality scores (Phred+33 encoding typical)

BED

chr1    1000    2000    feature_name    score    +

• 0-based coordinates
• Half-open interval [start, end)

AGP

chr1    1    5000    1    W    contig_1    1    5000    +
chr1    5001 5100    2    U    100    scaffold    yes    proximity_ligation

• Tab-delimited genome assembly format
• 1-based closed coordinates [start, end]
• Object and component lengths must match: obj_end - obj_beg + 1 == comp_end - comp_beg + 1

See reference.md for complete AGP specification, coordinate systems, validation rules, and processing patterns.

See common-issues.md for AGP coordinate debugging and unloc processing issues.

Coordinate Systems

1-based (SAM, VCF, GFF, AGP): First base is position 1. Interval [2,5] includes positions 2,3,4,5. 0-based (BED, BAM binary): First base is position 0. Interval [2,5) includes positions 2,3,4 (excludes 5).

Conversion: BED_start = SAM_start - 1; BED_end = SAM_end.

Best Practices

General

1. Always check data type: Paired-end vs single-end determines filtering strategy 2. Understand your sequencing technology: Hi-C behaves differently than HiFi 3. Filter in the right order: Proper pairs BEFORE region filtering 4. Validate outputs: Check file sizes, read counts, flagstat 5. Use appropriate MAPQ thresholds: Too stringent = lost data, too permissive = noise

For Hi-C Data

1. Expect distant read pairs: Don't be surprised by different scaffolds 2. Preserve proper pairs: Critical for downstream scaffolding 3. Use paired-aware tools: Standard filters may break pairs 4. Don't over-filter on MAPQ: Hi-C often has lower MAPQ than DNA-seq

For HiFi Data

1. Single-end processing: No pair concerns 2. High quality expected: Can use strict filters 3. Use appropriate presets: minimap2 map-hifi or map-pb 4. Consider read length distribution: HiFi reads vary in length

For Tool Testing

1. Create self-contained datasets: Both mates in selected region 2. Maintain proper pairs: Essential for realistic testing 3. Use representative data: Subsample proportionally, not randomly 4. Verify file sizes: Too small = overly filtered

Related Skills

• vgp-pipeline - VGP workflows process Hi-C and HiFi data
• galaxy-tool-wrapping - Galaxy tools work with SAM/BAM and sequencing data formats
• galaxy-workflow-development - Workflows process sequencing data

Supporting Documentation

• reference.md: Detailed format specifications (SAM/BAM complete reference, CIGAR operations, AGP format, FASTQ encoding, tool command reference, sequencing technology specs, assembly quality metrics, coverage calculations, coordinate systems)
• common-issues.md: Comprehensive troubleshooting guide (empty outputs, paired-end issues, quality/mapping problems, format conversion, Hi-C and HiFi specific issues, AGP processing errors, diagnostic commands)
• genomeark-data-access.md: GenomeArk AWS S3 data access patterns (directory structure evolution, QC data locations for GenomeScope/BUSCO/Merqury, fetching strategies, path normalization, assembly date extraction)
• genomic-analysis-patterns.md: Domain-specific analysis patterns (karyotype data curation, haploid vs diploid chromosome counts, phylogenetic tree species mapping, BED/telomere analysis, NCBI data integration strategies)

Version History

• v1.2.0: Split into SKILL.md + supporting files for maintainability; moved GenomeArk access, karyotype/chromosome analysis, phylogenetic mapping, AGP details, and telomere/NCBI patterns to supporting files
• v1.1.1: Added BED file processing patterns for telomere analysis and NCBI data integration strategies
• v1.1.0: Added comprehensive AGP format documentation including coordinate validation, unloc processing, and common error patterns
• v1.0.0: Initial release with SAM/BAM, Hi-C, HiFi, common filtering patterns

Common Bioinformatics Issues and Solutions

Troubleshooting guide for frequent problems in bioinformatics workflows, organized by symptom and data type.

---

Table of Contents

1. Empty or Missing Output Files 2. Paired-End Data Issues 3. Quality and Mapping Issues 4. Format and Conversion Issues 5. Hi-C Specific Issues 6. HiFi/Long-Read Specific Issues

---

Empty or Missing Output Files

Issue: Empty BAM After Region Filtering (Paired-End Data)

Symptom:

• Input BAM has reads
• After filtering by regions, output is empty
• Happens with paired-end data (especially Hi-C)

Root Cause: Region filtering (samtools view -L) operates on each read independently: 1. R1 maps to selected region → kept 2. R2 maps outside selected region → discarded 3. R1 loses "proper pair" flag (0x2) because mate is missing 4. Subsequent filtering for proper pairs removes all remaining reads 5. Result: Empty file

Solution:

Option A: Filter for proper pairs BEFORE region filtering

# Correct order: proper pair flag first, then region filtering
samtools view -b -f 2 -L regions.bed input.bam > output.bam

This ensures flag filtering happens before region filtering can break pairs.

Option B: Accept that one mate may be outside regions

# Keep pairs where at least one mate is in regions
# (mate may be anywhere in genome)
samtools view -b -L regions.bed input.bam > output.bam
# Don't filter for proper pairs afterward

Option C: Keep only pairs where BOTH mates are in regions (custom solution)

# Extract pairs in regions
samtools view -b -L regions.bed input.bam | \
# Filter for properly paired
bamtools filter -isPaired -isProperPair -out temp.bam

# Check if both mates are in selected regions (requires custom script)

When This Happens:

• Hi-C data (read pairs map to different scaffolds)
• Large-scale structural variation data
• Subsampling by genomic region
• Creating test datasets from larger assemblies

Prevention: Always think about filter order for paired-end data: 1. Pair-based filters first (proper pair, both mapped) 2. Region-based filters second 3. Quality filters can go either before or after

---

Issue: Empty FASTQ After BAM Conversion

Symptom:

• BAM file has reads (verified with samtools flagstat)
• samtools fastx produces empty FASTQ files
• Or R1 and R2 have different numbers of reads

Possible Causes:

Cause 1: Filtering Too Strict

# This might filter out ALL reads if none meet criteria
samtools fastx -1 R1.fq -2 R2.fq \
  --i1-flags 2,64 \  # Requires BOTH proper pair AND read1
  --i2-flags 2,128   # Requires BOTH proper pair AND read2
  input.bam

Solution: Relax filters or check that reads actually have these flags

# Check what flags are present
samtools flagstat input.bam

# Extract without strict filtering
samtools fastx -1 R1.fq -2 R2.fq input.bam

Cause 2: Single-End Data with Paired-End Extraction

# Wrong: trying to extract R1/R2 from single-end data
samtools fastx -1 R1.fq -2 R2.fq single_end.bam  # Empty output

Solution: Use correct extraction mode

# For single-end data
samtools fastx -0 output.fq single_end.bam

# Check if data is paired first
samtools flagstat input.bam | grep "paired in sequencing"

Cause 3: All Reads Are Unmapped

# Flags require mapped reads, but all are unmapped
samtools fastx --i1-flags 2 unmapped.bam  # Empty because flag 2 requires mapping

Solution:

# Check mapping rate
samtools flagstat input.bam

# Extract unmapped reads if needed
samtools fastx -f 4 -0 unmapped.fq input.bam

---

Paired-End Data Issues

Issue: R1 and R2 Files Have Different Read Counts

Symptom:

wc -l R1.fastq  # 1000000 lines (250000 reads)
wc -l R2.fastq  # 800000 lines (200000 reads)

Root Cause: Reads were filtered in a way that kept one mate but not the other.

Diagnosis:

# Check BAM for orphaned reads (mate unmapped)
samtools view -c -f 8 input.bam  # Count reads with unmapped mate

# Check for proper pairs
samtools view -c -f 2 input.bam  # Count proper pairs

Solution:

Option A: Filter BAM for proper pairs before extraction

# Only keep properly paired reads
samtools view -b -f 2 -F 12 input.bam | \
samtools fastx -1 R1.fq -2 R2.fq

Option B: Use paired-aware extraction

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

Option C: Repair pairs using external tools

# Use BBTools repair.sh
repair.sh in1=R1.fq in2=R2.fq out1=R1_fixed.fq out2=R2_fixed.fq outs=singletons.fq

---

Issue: Proper Pair Flag Missing After Alignment

Symptom:

samtools flagstat aligned.bam
# Shows:
# 1000000 paired in sequencing
# 0 properly paired (0.00%)  # <- Problem!

Possible Causes:

Cause 1: Insert Size Distribution Mismatch Aligner expected insert size 300-500 bp, but actual library is 100-200 bp.

Solution:

# Check actual insert size
samtools stats aligned.bam | grep "insert size average"

# Re-align with correct insert size expectations
bwa mem -I 200,50 ref.fa R1.fq R2.fq  # Mean 200, StdDev 50

Cause 2: Wrong Reference Genome Reads are from different assembly/species than reference.

Solution:

# Check mapping rate
samtools flagstat aligned.bam
# If very low (<50%), likely wrong reference

Cause 3: Reads Are Actually Properly Paired (Hi-C) For Hi-C data, distant pairs are EXPECTED and may not be marked as "proper pairs" by standard aligners.

Solution:

• This is normal for Hi-C
• Use Hi-C-specific pipelines (HiC-Pro, Juicer)
• Don't rely on "proper pair" flag for Hi-C
• Filter by MAPQ and mapping status instead

Cause 4: Incorrect Orientation Flags

# Check read orientation distribution
samtools view aligned.bam | awk '{print and($2,16), and($2,32)}' | sort | uniq -c

# Should see expected patterns for paired-end:
# Many reads: 0 32 (R1 forward, R2 reverse)
# Many reads: 16 0 (R1 reverse, R2 forward)

---

Issue: Lost Reads During Merging or Filtering

Symptom: Start with 10M reads, end with 100K after filtering pipeline.

Diagnosis Strategy:

Step 1: Count reads at each step

samtools view -c step1.bam  # 10000000
samtools view -c step2.bam  # 5000000  <- lost 50% here
samtools view -c step3.bam  # 100000   <- lost 98% here!

Step 2: Identify the problem step

# For step that lost reads, check why
samtools flagstat problem_step.bam
samtools stats problem_step.bam

Step 3: Check filter parameters

# Was MAPQ too stringent?
samtools view problem_step.bam | awk '{print $5}' | sort -n | uniq -c

# Were flag filters too strict?
samtools flagstat input_to_problem_step.bam

Common Culprits:

• MAPQ > 30 on Hi-C data (too strict)
• Proper pair requirement on scaffolded/fragmented assemblies
• Region filtering on paired-end data
• Duplicate marking removing too many reads

---

Quality and Mapping Issues

Issue: Low Mapping Rate

Symptom:

samtools flagstat aligned.bam
# 1000000 reads
# 100000 mapped (10.00%)  # <- Very low!

Possible Causes:

Cause 1: Wrong Reference Most common cause.

Solution:

# Verify reference
head -1 reference.fa  # Check sequence ID matches expectations
samtools view aligned.bam | head  # Check RNAME column matches reference

Cause 2: Adapter Contamination Reads still contain sequencing adapters.

Solution:

# Check for adapters
fastqc reads.fq  # Look for adapter content

# Trim adapters
cutadapt -a AGATCGGAAGAGC -o trimmed.fq reads.fq

# Re-align
bwa mem ref.fa trimmed.fq > aligned.sam

Cause 3: Low Quality Reads

# Check quality distribution
fastqc reads.fq
# Look for quality drop-off

# Filter low quality
fastp -i reads.fq -o filtered.fq -q 20

Cause 4: Contamination Reads from different organism.

Solution:

# Align to suspected contaminant reference
bwa mem contaminant.fa reads.fq > contam_check.sam
samtools flagstat contam_check.sam  # High mapping? Contamination!

Expected Mapping Rates:

• DNA-seq to same species: >95%
• DNA-seq to closely related species: 70-90%
• RNA-seq (with introns): 70-85%
• Hi-C: 60-80% (lower is normal due to chimeric reads)
• Metagenomics: Highly variable

---

Issue: Many Secondary/Supplementary Alignments

Symptom:

samtools flagstat aligned.bam
# 1000000 total reads
# 500000 secondary  # <- Very high!
# 300000 supplementary

Meaning:

• Secondary (flag 256): Alternative mapping locations (multi-mappers)
• Supplementary (flag 2048): Chimeric alignment (read maps to multiple locations)

When This Is Normal:

• Repetitive regions (expect high secondary)
• Structural variants (expect high supplementary)
• Long reads spanning multiple regions (supplementary normal for PacBio)

When This Is a Problem:

• Fragmented assembly (reads mapping everywhere)
• Wrong reference
• Poor quality reads

Solution:

For most analyses, exclude secondary and supplementary:

samtools view -b -F 2304 aligned.bam > primary_only.bam
# -F 2304 = exclude 256 (secondary) + 2048 (supplementary)

---

Format and Conversion Issues

Issue: BAM to FASTQ Loses Read Names

Symptom: Original FASTQ has read names like @INSTRUMENT:RUN:FLOWCELL... After BAM→FASTQ, names are truncated or modified.

Cause: SAM specification only stores first field of read name (up to first space).

Solution:

Preserve names during alignment:

# bwa mem preserves full names with -C flag
bwa mem -C ref.fa R1.fq R2.fq > aligned.sam

Or accept that descriptions are lost (usually fine for most analyses).

---

Issue: CRAM Conversion Fails

Symptom:

samtools view -C -T ref.fa input.bam > output.cram
# Error: [main_samview] failed to read the header

Cause: CRAM requires reference genome for compression.

Solutions:

Provide reference:

samtools view -C -T reference.fa -o output.cram input.bam

Check reference matches:

# Reference sequence names must match BAM @SQ headers
samtools view -H input.bam | grep "^@SQ"
grep "^>" reference.fa

Embed reference in CRAM:

samtools view -C --output-fmt-option embed_ref=1 -o output.cram input.bam

---

Hi-C Specific Issues

Issue: Very Low Valid Pairs Percentage

Symptom: Hi-C analysis shows <20% valid pairs (expected: 40-70%).

Diagnosis:

Check pair types:

• Self-circles: Ligation of same fragment
• Dangling ends: No ligation occurred
• Re-ligation: Fragments from same restriction site

Common Causes:

Cause 1: Restriction Enzyme Problem Wrong enzyme used in analysis vs library prep.

Solution:

# Verify restriction enzyme
# Common enzymes: MboI (GATC), DpnII (GATC), HindIII (AAGCTT)

# Check if cut sites are where expected
grep "GATC" assembly.fa | head  # For MboI/DpnII

Cause 2: Poor Library Quality Ligation step failed, over-digestion, etc.

Solution: This is a wet-lab issue - library needs to be re-made.

Cause 3: Mapping Issues Reference genome fragmented or incorrect.

Solution:

# Check reference contiguity
assembly-stats assembly.fa
# Many small contigs = will have low valid pairs

# Use more contiguous assembly if available

---

Issue: Hi-C Reads Not Properly Paired After Filtering

Symptom: After region filtering to create test data, Hi-C reads are empty or unpaired.

Root Cause: Hi-C pairs span different regions - region filtering breaks them. (See "Empty BAM After Region Filtering" above).

Solution: See detailed solution in main issues section above.

---

HiFi/Long-Read Specific Issues

Issue: Low HiFi Mapping Rate Despite High Quality

Symptom:

• HiFi reads have Q20+
• But <80% mapping rate to reference

Possible Causes:

Cause 1: Wrong Mapper Preset

# Wrong: using short-read aligner
bwa mem ref.fa hifi.fq  # Designed for short reads!

# Wrong: using CLR preset
minimap2 -ax map-pb ref.fa hifi.fq  # For older PacBio

# Correct: use HiFi preset
minimap2 -ax map-hifi ref.fa hifi.fq
# Or
minimap2 -ax map-pb -k 19 ref.fa hifi.fq  # Explicit HiFi params

Cause 2: Structural Differences Reference is from different haplotype/individual.

Solution: This may be expected - HiFi reveals real structural variation.

Cause 3: Low Complexity or Repetitive Reads

# Check for low-complexity sequences
seqtk comp hifi.fq | awk '$2 < 50 || $3/$2 > 0.4'  # Potential low-complexity

---

Issue: HiFi Reads Marked as Duplicates

Symptom: Many HiFi reads flagged as duplicates (flag 1024).

Cause: Duplicate marking tools designed for short reads don't work well with long reads.

Solution:

Don't mark duplicates on HiFi:

# Skip duplicate marking
# Or use long-read aware duplicate removal

# Remove duplicate flags if already marked
samtools view -b -F 1024 marked.bam > no_dups.bam

For long reads: True PCR duplicates are rare (each molecule is unique due to length).

---

Diagnostic Commands

Quick BAM Inspection

# Overall statistics
samtools flagstat file.bam

# Detailed statistics
samtools stats file.bam | less

# Check first few alignments
samtools view file.bam | head

# Check header
samtools view -H file.bam

# Count by flag
samtools view file.bam | awk '{print $2}' | sort | uniq -c | sort -rn

# MAPQ distribution
samtools view file.bam | awk '{print $5}' | sort -n | uniq -c

# Insert size distribution (paired-end)
samtools stats file.bam | grep "^IS" | head -20

# Coverage per chromosome
samtools idxstats file.bam

Quick FASTQ Inspection

# Count reads
echo $(cat file.fq | wc -l)/4 | bc

# For gzipped
echo $(zcat file.fq.gz | wc -l)/4 | bc

# Read length distribution
awk 'NR%4==2 {print length($0)}' file.fq | sort -n | uniq -c

# Quality score distribution
awk 'NR%4==0' file.fq | perl -ne 'chomp; @q=split(//); foreach(@q){print ord($_)-33, "\n";}' | sort -n | uniq -c

# Check for adapters (simple check)
grep "AGATCGGAAGAGC" file.fq | wc -l

Comparing BAM Files

# Compare read counts
samtools view -c file1.bam
samtools view -c file2.bam

# Compare flagstats
diff <(samtools flagstat file1.bam) <(samtools flagstat file2.bam)

# Check if same reads
samtools view file1.bam | cut -f1 | sort > reads1.txt
samtools view file2.bam | cut -f1 | sort > reads2.txt
diff reads1.txt reads2.txt

---

Prevention Checklist

Before running a complex filtering pipeline:

• [ ] Understand input data type (paired-end vs single-end)
• [ ] Know sequencing technology (Hi-C, HiFi, Illumina)
• [ ] Plan filter order (pairs first, then regions)
• [ ] Set reasonable MAPQ thresholds for data type
• [ ] Validate each step with flagstat/stats
• [ ] Check output file sizes (empty = problem)
• [ ] Verify read counts match expectations
• [ ] Test on small subset first

---

AGP Processing Issues

Issue: Incorrect Object Coordinates When Creating Unlocs

Symptom:

ERROR: object coordinates (1, 19328398) and component coordinates (19274039, 19328398)
do not have the same length

Cause: When converting a region of a scaffold into an unlocalized sequence (unloc), the object coordinates must represent the length of the extracted region, not the original component end coordinate.

Wrong Approach:

# Setting object end to component end coordinate
agp_df.loc[index, 'chr_end'] = agp_df.loc[index, 'scaff_end']  # WRONG

Correct Approach:

# Calculate actual length from component coordinates
agp_df.loc[index, 'chr_end'] = int(agp_df.loc[index, 'scaff_end']) - int(agp_df.loc[index, 'scaff_start']) + 1  # CORRECT

---

Issue: Component Numbering Not Reset for New Objects

Symptom: Unloc scaffolds have component numbers > 1 when they should start at 1.

Cause: When creating a new object (unloc scaffold), component numbering wasn't reset.

Solution:

# When creating unlocs, reset component number
agp_df.loc[index, '#_scaffs'] = 1  # Column 4: component number

---

Issue: AGPcorrect Accumulating Coordinates

Symptom: Unloc sequences inherit cumulative coordinates from parent scaffolds.

Cause: AGPcorrect adjusts coordinates based on sequence length corrections. When scaffolds are later split into unlocs, the accumulated corrections need to be recalculated based on actual component spans.

Solution: Always recalculate object coordinates from component spans when creating new objects (unlocs).

---

AGP Coordinate Debugging Pattern

When encountering coordinate errors:

# For each AGP line, verify:
obj_length = int(obj_end) - int(obj_beg) + 1
comp_length = int(comp_end) - int(comp_beg) + 1

assert obj_length == comp_length, f"Length mismatch: obj={obj_length}, comp={comp_length}"

# For sequential component numbers:
assert comp_num == expected_num, f"Component number gap: got {comp_num}, expected {expected_num}"

AGP Processing Best Practices

Creating Unlocalized Sequences (Unlocs)

# When extracting a region to create an unloc:
# 1. Calculate the actual length of the region
length = int(comp_end) - int(comp_start) + 1

# 2. Set object coordinates for the new unloc
obj_start = 1  # Always starts at 1
obj_end = length  # Equals the length

# 3. Reset component number
component_num = 1  # New object, new numbering

# 4. Rename the object
new_object_name = f"{parent_scaffold}_unloc_{unloc_number}"

Validating AGP Files

Use NCBI's AGP validator:

agp_validate assembly.agp

Common validation checks:

• Object/component length match
• Sequential component numbering
• No coordinate overlaps
• Gap specifications valid
• Orientation values (+, -, ?, 0, na)

Handling Haplotype-Split Assemblies

When splitting diploid assemblies into haplotypes: 1. Identify haplotype markers in sequence names (H1/hap1, H2/hap2) 2. Maintain proper pairing information 3. Process unlocs separately per haplotype 4. Remove haplotig duplications 5. Track gaps appropriately (especially proximity ligation gaps)

AGP File Structure by Assembly Stage

1. Raw Assembly AGP:

• Direct representation from assembler
• May have incorrect sequence lengths
• Needs coordinate correction (AGPcorrect)

2. Corrected AGP:

• Sequence lengths match actual FASTA
• Coordinates adjusted for length discrepancies
• Ready for haplotype splitting

3. Haplotype-Split AGP:

• Separate files per haplotype
• Unlocs identified but not separated
• Haplotigs marked but not removed

4. Final Curated AGP:

• Unlocs separated into individual objects
• Haplotigs removed to separate file
• Proximity ligation gaps cleaned
• Ready for database submission

---

Getting Help

When asking for help with bioinformatics issues, include:

1. Data type: Paired-end? Single-end? Hi-C? HiFi? 2. Commands used: Exact command lines 3. Input stats: samtools flagstat input.bam 4. Output stats: samtools flagstat output.bam 5. Expected vs actual: What did you expect? What happened? 6. File sizes: Are outputs unexpectedly small/large?

This helps diagnose issues quickly!

GenomeArk AWS S3 Data Access

Patterns and strategies for accessing VGP genome assembly data and QC metrics from GenomeArk public S3 buckets.

---

Overview

GenomeArk (s3://genomeark/) is a public AWS S3 bucket containing VGP genome assemblies and QC data. Access requires no credentials using --no-sign-request.

Critical Discovery: GenomeArk S3 structure has evolved over time (2022 -> 2024). Always try multiple path patterns for reliability.

Directory Structure Evolution

Base structure:

s3://genomeark/species/{Species_name}/{ToLID}/assembly_vgp_{type}_2.0/evaluation/

Key variations:

1. Case sensitivity:

• Table may store: assembly_vgp_hic_2.0
• S3 requires: assembly_vgp_HiC_2.0 (case-sensitive!)
• Always normalize: Replace hic -> HiC before fetching

2. Subspecies handling (CRITICAL, discovered Feb 2026):

• Dataset may contain trinomial names: "Elephas maximus indicus"
• GenomeArk uses only binomial: Elephas_maximus
• Always use only first two words (Genus + species) for S3 paths
• Example: "Elephas maximus indicus" -> s3://genomeark/species/Elephas_maximus/...

3. Assembly directory patterns (multiple generations exist):

Try in this order for maximum coverage (Feb 2026 - 8 patterns):

   path_patterns = [
       # Version 2.0 paths (original, most common)
       f"assembly_vgp_HiC_2.0/",
       f"assembly_vgp_standard_2.0/",
       f"assembly_vgp_hic_2.0/",

       # Without version suffix (newer assemblies)
       f"assembly_vgp_HiC/",
       f"assembly_vgp_standard/",

       # Alternative assembly types
       f"assembly_curated/",
       f"assembly_cambridge/",
       f"assembly_rockefeller/",
   ]

Coverage impact:

• 3 patterns (original): ~24% assembly coverage
• 8 patterns (expanded): ~35-45% assembly coverage
• Always try multiple patterns - don't assume structure

4. Species name construction:

   # Handle subspecies correctly
   species_parts = scientific_name.strip().split()
   if len(species_parts) >= 2:
       species_name = f"{species_parts[0]}_{species_parts[1]}"
   else:
       species_name = scientific_name.strip().replace(' ', '_')

   # Full path example
   s3_path = f"s3://genomeark/species/{species_name}/{tolid}/{assembly_type}/"

QC Data Locations and Formats

1. GenomeScope (Genome Size, Heterozygosity, Repeat Content)

Path: {assembly}/evaluation/genomescope/

Filename patterns (try in order): 1. {ToLID}_genomescope__Summary.txt (Pattern A: double underscore - most common) 2. {ToLID}_genomescope_Summary.txt (Pattern C: single underscore - EASILY MISSED) 3. {ToLID}_Summary.txt (Pattern B: no prefix - older assemblies)

CRITICAL: ALL THREE patterns must be checked! Pattern C (single underscore) was discovered in Feb 2026 during debugging - checking only patterns A and B causes ~30-40% of data to be missed!

Example of Pattern C:

• Missing: rPlaMeg1_genomescope__Summary.txt (not found)
• Found: rPlaMeg1_genomescope_Summary.txt (exists)

CRITICAL: Validate Data Quality

Failed GenomeScope runs show unrealistic ranges:

Heterozygous (ab)    0%    100%     <- FAILED RUN - DO NOT USE

Good runs show narrow ranges:

Heterozygous (ab)    0.49%    0.54%  <- VALID - use max value

Validation logic:

# Extract min and max percentages
percentages = [0.49, 0.54]  # Example from parsing
min_val, max_val = percentages[0], percentages[-1]
range_width = max_val - min_val

# Validate before using
if range_width <= 50.0 and max_val <= 95.0:
    heterozygosity = max_val  # ACCEPT
else:
    heterozygosity = None  # REJECT - failed run

Skip values if:

• Range width > 50% (indicates model failure)
• Max value > 95% (unrealistic for most genomes)
• Range is exactly 0%-100% (complete failure)

Summary.txt format:

GenomeScope version 2.0
...
property                      min               max
Genome Haploid Length         4,077,481,159 bp  4,095,803,536 bp
Heterozygous (ab)             1.43264%          1.47696%
Genome Repeat Length          2,528,408,288 bp  2,539,769,824 bp

Parsing:

• Genome size: Take max value (second number), remove commas
• Heterozygosity: Take max percentage (validate range first!)
• Repeat content: Calculate (repeat_length / genome_size) * 100

2. BUSCO (Assembly Completeness)

Path: {assembly}/evaluation/busco/{subdir}/

Subdirectories vary:

• c/, c1/ - primary results
• p/, p1/ - alternate results
• Search dynamically, don't hardcode

Files: *short_summary*.txt (case-insensitive search)

Filename patterns:

• HiC assemblies: {ToLID}_HiC__busco_hap1_busco_short_summary.txt
• Standard assemblies: {ToLID}_busco_short_summary.txt

Format:

# BUSCO version is: 5.2.2
# The lineage dataset is: vertebrata_odb10
...
	C:94.0%[S:92.4%,D:1.6%],F:2.7%,M:3.3%,n:3354

Parse line starting with `C:`: Extract 94.0 from C:94.0%

Expected coverage: ~20-30% of VGP assemblies have BUSCO data

3. Merqury (Assembly QV Scores)

TWO PATH PATTERNS (structure changed 2022 -> 2024):

Pattern A (Newer - Direct, 2024+):

{assembly}/evaluation/merqury/{ToLID}_qv/output_merqury.tabular

Pattern B (Older - Nested, 2022):

{assembly}/evaluation/merqury/{c,p}/{ToLID}_qv/output_merqury.tabular

Strategy: Try direct path first, then search for nested subdirectories

File format (tab-separated, may have header):

assembly	unique k-mers	common k-mers	QV	error rate
assembly_01	20197	2133011206	63.4592	4.50896e-07
assembly_02	19654	2304717679	63.9138	4.06084e-07
Both	39851	4437728885	63.6894	4.27623e-07

Parsing:

• Skip header line if starts with assembly\t
• QV is always column 4 (index 3)
• Take first data line (usually assembly_01 or Both)

Complete Fetching Strategy

def normalize_s3_path(s3_path):
    """Normalize path for GenomeArk (case sensitivity!)"""
    if not s3_path:
        return None
    # Critical: HiC capitalization
    s3_path = s3_path.replace('/assembly_vgp_hic_2.0/', '/assembly_vgp_HiC_2.0/')
    if not s3_path.endswith('/'):
        s3_path += '/'
    return s3_path

def fetch_genomescope_data(s3_path):
    """Fetch with validation"""
    s3_path = normalize_s3_path(s3_path)
    tolid = s3_path.rstrip('/').split('/')[-2]

    # Try ALL THREE filename patterns
    for filename in [
        f'{tolid}_genomescope__Summary.txt',   # Pattern A: double underscore
        f'{tolid}_genomescope_Summary.txt',    # Pattern C: single underscore
        f'{tolid}_Summary.txt'                  # Pattern B: no prefix
    ]:
        file_path = f"{s3_path}evaluation/genomescope/{filename}"
        result = subprocess.run(['aws', 's3', 'cp', file_path, '-', '--no-sign-request'],
                               capture_output=True, text=True, timeout=30)

        if result.returncode == 0 and result.stdout:
            # Parse and validate
            data = parse_genomescope(result.stdout)

            # Validate heterozygosity range
            if 'heterozygosity' in data:
                # Check if range is reasonable
                if heterozygosity_range > 50.0 or max_het > 95.0:
                    del data['heterozygosity']  # Skip invalid value

            if data:
                return data
    return None

def fetch_merqury_data(s3_path):
    """Fetch from direct or nested paths"""
    s3_path = normalize_s3_path(s3_path)
    tolid = s3_path.rstrip('/').split('/')[-2]

    # Try direct path first (newer structure)
    direct_path = f"{s3_path}evaluation/merqury/{tolid}_qv/output_merqury.tabular"
    result = subprocess.run(['aws', 's3', 'cp', direct_path, '-', '--no-sign-request'],
                           capture_output=True, text=True, timeout=30)

    if result.returncode == 0 and result.stdout:
        # Parse QV from column 4
        for line in result.stdout.split('\n'):
            if line.strip() and not line.startswith('assembly\t'):
                parts = line.split('\t')
                if len(parts) >= 4:
                    return {'qv': float(parts[3]), 'path_type': 'direct'}

    # Fallback: search nested subdirectories (older structure)
    # List subdirectories, try c/, p/, etc.
    ...

def fetch_busco_data(s3_path):
    """Search dynamic subdirectories"""
    s3_path = normalize_s3_path(s3_path)

    # List busco/ subdirectories
    list_result = subprocess.run(['aws', 's3', 'ls', f"{s3_path}evaluation/busco/", '--no-sign-request'],
                                capture_output=True, text=True, timeout=10)

    # Find subdirectories (lines with 'PRE')
    subdirs = [line.split('PRE')[1].strip().rstrip('/')
               for line in list_result.stdout.split('\n') if 'PRE' in line]

    # Try each subdirectory for short_summary files
    ...

Extracting Assembly Completion Dates

Problem: NCBI release_date doesn't reflect actual assembly completion due to curation/submission delays.

Solution: GenomeArk filenames contain YYYYMMDD timestamps showing actual assembly completion dates.

Filename pattern:

mLoxAfr1.HiC.hap1.20221209.fasta.gz
                  ^^^^^^^^
                  YYYYMMDD = Dec 9, 2022

Strategy:

import re
from datetime import datetime

def extract_assembly_year(tolid, scientific_name):
    """Extract assembly year from GenomeArk filenames"""

    # 1. Construct S3 path (try multiple patterns)
    species_name = '_'.join(scientific_name.split()[:2])  # Handle subspecies!

    for assembly_type in ['assembly_vgp_HiC_2.0', 'assembly_vgp_standard_2.0',
                          'assembly_vgp_HiC', 'assembly_curated', ...]:
        s3_path = f"s3://genomeark/species/{species_name}/{tolid}/{assembly_type}/"

        # 2. List files recursively (exclude _curated subdirs)
        cmd = ['aws', 's3', 'ls', s3_path, '--recursive', '--no-sign-request']
        result = subprocess.run(cmd, capture_output=True, text=True)

        if result.returncode != 0:
            continue  # Try next pattern

        # 3. Extract YYYYMMDD dates from filenames
        date_pattern = r'(\d{8})'
        dates = []

        for line in result.stdout.split('\n'):
            if '_curated' in line:  # Exclude curated versions (later dates)
                continue

            matches = re.findall(date_pattern, line)
            for date_str in matches:
                year = int(date_str[:4])
                month = int(date_str[4:6])
                day = int(date_str[6:8])

                # 4. Validate date (2000-2030, valid month/day)
                if 2000 <= year <= 2030 and 1 <= month <= 12 and 1 <= day <= 31:
                    dates.append(date_str)

        # 5. Return most recent year found
        if dates:
            most_recent = max(dates)
            return int(most_recent[:4])

    return None  # No dates found in any path

Coverage expectations:

• Newer assemblies (2020+): 80-90%
• Older assemblies (pre-2020): 50-60%
• Overall: 60-80%

Real examples of delays:

• mLoxAfr1: Assembly 2022, NCBI release 2023 (1-year delay)
• mLemCat1: Assembly 2021, NCBI release 2021 (same year)
• mEleMax1: Assembly 2022, NCBI release 2022 (same year)

Use case: Temporal trend analysis where accurate assembly dates are critical for identifying methodology vs. technology effects.

S3 Path Normalization

Always normalize paths:

def normalize_s3_path(s3_path):
    s3_path = s3_path.strip()
    s3_path = s3_path.replace('/assembly_vgp_hic_2.0/', '/assembly_vgp_HiC_2.0/')
    if not s3_path.endswith('/'):
        s3_path += '/'
    return s3_path

AWS CLI Usage

Public access (no credentials):

aws s3 ls s3://genomeark/... --no-sign-request
aws s3 cp s3://genomeark/.../file.txt - --no-sign-request

Prefer subprocess + aws s3 CLI with --no-sign-request over boto3 (which requires credential config even for public access).

# Simple and works
cmd = ['aws', 's3', 'cp', s3_path, '-', '--no-sign-request']
result = subprocess.run(cmd, capture_output=True, text=True, timeout=30)

Timeouts: Use 10-30s timeouts for robustness

Expected Performance

• S3 path inference: ~5-10 seconds per ToLID
• QC data fetching: ~1-2 minutes per assembly
• Full dataset (700+ assemblies): 2-3 hours total

Common Pitfalls

1. Case sensitivity: assembly_vgp_hic_2.0 in table -> assembly_vgp_HiC_2.0 in S3 2. Subspecies names: "Elephas maximus indicus" -> use only "Elephas_maximus" (CRITICAL) 3. Limited path patterns: Only trying _2.0 paths misses ~40% of assemblies - always try 8+ patterns 4. Directory evolution: Merqury moved from nested to direct structure 5. Failed QC runs: Always validate genomescope ranges before use 6. Subdirectory variations: BUSCO/Merqury use different subdir names (c vs c1 vs p) 7. File format variations: Merqury may/may not have header line 8. Haplotype-specific files: HiC assemblies have separate hap1/hap2 BUSCO results 9. Excluding `_curated` directories: These have later curation dates, not original assembly dates

Best Practices

1. Path normalization: Always fix case sensitivity 2. Handle subspecies: Extract only Genus + species (first 2 words) 3. Try multiple patterns: 8 assembly type patterns for maximum coverage 4. Validate data: Check ranges, detect failed analyses 5. Exclude curated versions: When extracting dates, skip _curated subdirectories 6. Dynamic discovery: List subdirectories, don't hardcode 7. Error handling: Continue on failures, report what succeeded 8. Timeouts: 10-30s per fetch, don't hang indefinitely 9. Rate limiting: 0.2s delay between fetches (respectful to AWS)

Testing Examples

Confirmed working paths:

# GenomeScope - Pattern A (double underscore)
aws s3 cp s3://genomeark/species/Gastrophryne_carolinensis/aGasCar1/assembly_vgp_HiC_2.0/evaluation/genomescope/aGasCar1_genomescope__Summary.txt - --no-sign-request

# GenomeScope - Pattern C (single underscore)
aws s3 cp s3://genomeark/species/Platysternon_megacephalum/rPlaMeg1/assembly_vgp_HiC_2.0/evaluation/genomescope/rPlaMeg1_genomescope_Summary.txt - --no-sign-request

# GenomeScope - Pattern B (no prefix - older)
aws s3 cp s3://genomeark/species/Spea_bombifrons/aSpeBom1/assembly_vgp_standard_2.0/evaluation/genomescope/aSpeBom1_Summary.txt - --no-sign-request

# BUSCO
aws s3 cp s3://genomeark/species/Gastrophryne_carolinensis/aGasCar1/assembly_vgp_HiC_2.0/evaluation/busco/c/aGasCar1_HiC__busco_hap1_busco_short_summary.txt - --no-sign-request

# Merqury - Direct path (2024+)
aws s3 cp s3://genomeark/species/Ia_io/mIaxIox2/assembly_vgp_HiC_2.0/evaluation/merqury/mIaxIox2_qv/output_merqury.tabular - --no-sign-request

# Merqury - Nested path (2022)
aws s3 cp s3://genomeark/species/Gastrophryne_carolinensis/aGasCar1/assembly_vgp_HiC_2.0/evaluation/merqury/aGasCar1_qv/output_merqury.tabular - --no-sign-request

Genomic Analysis Patterns

Domain-specific patterns for karyotype data curation, chromosome count analysis, phylogenetic tree mapping, telomere classification, and NCBI data integration.

---

Karyotype Data Curation and Literature Search

Overview

Karyotype data (diploid 2n and haploid n chromosome numbers) is critical for genome assembly validation but rarely available via APIs. Manual literature curation is required.

Search Strategy

Effective Search Terms

"{species_name} karyotype chromosome 2n"
"{species_name} diploid number karyotype"
"{genus} karyotype evolution"
"cytogenetic analysis {family_name}"
"{species_name} chromosome number diploid"

Best Reference Sources

1. PubMed/PMC: Primary cytogenetic studies 2. ResearchGate: Karyotype descriptions and figures 3. Specialized databases:

• Bird Chromosome Database: https://sites.unipampa.edu.br/birdchromosomedatabase/
• Animal Genome Size Database: http://www.genomesize.com/

4. Genome assembly papers: Often mention expected karyotype 5. Comparative cytogenetic studies: Family-level analyses

Search Time Estimates

• Model organisms, domestic species: 2-3 minutes
• Well-studied taxonomic groups: 5-10 minutes
• Rare/uncommon species: 10-20 minutes or not found

Taxonomic Conservation Patterns

Mammals

• Cetaceans: Highly conserved 2n = 44, n = 22 (exceptions: pygmy sperm whale, right whale, beaked whales = 2n = 42)
• Felidae: Conserved 2n = 38, n = 19
• Canidae: Conserved 2n = 78, n = 39
• Primates: Variable (great apes 2n = 48, macaques 2n = 42, marmosets 2n = 46)

Birds

• Anatidae (waterfowl): Highly conserved 2n = 80, n = 40 across ducks, geese, swans
• Galliformes (game birds): Typically 2n = 78, n = 39 (chicken, quail, grouse)
• Passerines: Variable 2n = 78-82, most common 2n = 80
• Ancestral avian karyotype: Putative 2n = 80
• General pattern: 50.7% of birds have 2n = 78-82; 21.7% have exactly 2n = 80

Reptiles

• Lacertidae (wall lizards): Often 2n = 38, n = 19

Genome Assembly Interpretation

Warning: Chromosome-level assemblies often report fewer chromosomes than actual diploid number.

Why: Assemblies typically capture only:

• Macrochromosomes (large chromosomes)
• Larger microchromosomes
• Small microchromosomes remain unassembled

Example: Waterfowl with 2n = 80 often have genome assemblies with 34-42 "chromosomes"

• True karyotype: 10 macro pairs + 30 micro pairs = 80
• Assembly: ~34-42 scaffolds (only macro + larger micros)

Using Conservation for Inference

When specific karyotype data is unavailable but genus/family patterns are strong:

1. High confidence inference (acceptable for publication):

• Multiple congeneric species confirmed
• Family-level conservation documented
• No known exceptions in genus

2. Document inference clearly:

   accession,taxid,species,2n,n,notes,reference
   GCA_XXX,123,Species name,80,40,Inferred from Anatidae conservation,https://family-level-study.url

3. Priority for direct confirmation:

• Species with conservation exceptions
• Type specimens or reference species
• Phylogenetically divergent lineages

VGP-Specific: Sex Chromosome Adjustment

When both sex chromosomes are in main haplotype (common in VGP assemblies):

• Expected scaffolds = n + 1 (not n)
• Reason: X+Y or Z+W = two distinct chromosomes
• Check: VGP metadata column "Sex chromosomes main haplotype"
• Patterns: "Has X and Y", "Has Z and W", "Has X1, X2, and Y"

Data Recording Format

CSV Structure:

accession,taxid,species_name,diploid_2n,haploid_n,notes,reference
GCA_XXXXXX,12345,Species name,80,40,Brief description,https://doi.org/...

Notes field examples:

• "Standard {family} karyotype"
• "Conserved {genus} karyotype"
• "Inferred from {family} conservation"
• "Unusual karyotype for family"
• "Geographic variation reported"

Prioritization for Literature Searches

TIER 1 (>90% success rate):

• Model organisms (zebrafish, mouse, medaka)
• Domestic species (chicken, goat, sheep)
• Game animals (waterfowl, deer)
• Laboratory species (fruit fly, nematode)

TIER 2 (70-90% success rate):

• Well-studied taxonomic groups (Podarcis lizards, corvids)
• Conservation focus species (raptors, large mammals)
• Commercial species (salmonids, oysters)

TIER 3 (50-70% success rate):

• Common but not economically important
• Widespread distribution
• Recent phylogenetic interest

Low priority (<50% success rate):

• Deep-sea species
• Rare/endangered without conservation genetics
• Recently described species
• Cryptic species complexes

---

Haploid vs Diploid Chromosome Counts in Assembly Analysis

The Critical Distinction

Genome assembly metadata typically includes both haploid and diploid chromosome counts:

• Haploid count (n): Number of chromosomes in a single genome copy
• Example: Human n=23 (22 autosomes + X or Y)
• Represents unique chromosome types
• Diploid count (2n): Number of chromosomes in diploid organism
• Example: Human 2n=46 (23 pairs)
• Represents total chromosomes in a diploid cell

Common Dataset Column Names

# Typical column names (exact names vary by dataset):
df['num_chromosomes']               # Often diploid (2n)
df['total_number_of_chromosomes']   # Often haploid (n)
df['karyotype']                     # Usually haploid (n)
df['num_chromosomes_haploid_adjusted']  # Haploid with sex chr adjustment

WARNING: Column names are NOT standardized across datasets - always verify which is which!

Which Count to Use When

Use HAPLOID (n) for:

• Per-assembly comparisons (scaffolds per assembly)
• Chromosome assignment ratios
• Expected vs observed chromosome counts
• Telomere counts (2 per chromosome x n chromosomes)
• Scaffold-to-chromosome mapping

Use DIPLOID (2n) for:

• Cell-level comparisons
• Comparing to diploid karyotypes
• Ploidy analyses
• Cytogenetic studies

Real-World Example: VGP Assembly Analysis

Problem: Used num_chromosomes (diploid) for per-assembly comparison

Result: All assemblies appeared to have 2x expected chromosomes

Fix: Changed to total_number_of_chromosomes (haploid)

Validation: Ratio now ~1.0 instead of ~2.0

# WRONG - uses diploid count
fig, ax = plt.subplots()
ax.scatter(df['num_chromosomes'], df['num_scaffolds_assigned'])
# Result: Everything appears at 2x diagonal

# CORRECT - uses haploid count
fig, ax = plt.subplots()
ax.scatter(df['total_number_of_chromosomes'], df['num_scaffolds_assigned'])
# Result: Expected 1:1 diagonal relationship

Sex Chromosome Adjustments

Some species have different haploid counts by sex:

• Male XY systems: n = autosomes + 2 (X and Y count separately)
• Female XX systems: n = autosomes + 1 (both X chromosomes count as one type)
• For telomere counts: Male XY may need +1 adjustment (X and Y both have telomeres)

Check for adjusted counts:

# Some datasets provide sex-adjusted haploid counts
# Example: Human male
# Karyotype n = 23 (22 autosomes + X or Y)
# But for telomere counting: 24 (22 autosomes + X + Y both have telomeres)

df['num_chromosomes_haploid_adjusted']  # May add +1 for male XY

Validation Checks

# Check if counts are haploid or diploid by testing known species
human_samples = df[df['species'] == 'Homo sapiens']
median_count = human_samples['column_name'].median()

if median_count > 40:
    print("Likely diploid (2n) - expect ~46 for humans")
elif median_count > 20:
    print("Likely haploid (n) - expect ~23 for humans")
else:
    print("Check data - values unexpectedly low")

# Verify ratios make biological sense
df['ratio'] = df['scaffolds_assigned'] / df['haploid_count']
assert 0.5 < df['ratio'].median() < 2.0, "Ratio should be near 1.0 for good assemblies"

# Check for systematic doubling
if df['ratio'].median() > 1.8:
    print("WARNING: May be using diploid count - ratios systematically doubled")

Common Pitfalls

1. Assuming column names are accurate

• num_chromosomes could be either n or 2n
• Always validate with known species

2. Not accounting for sex chromosomes

• Male XY vs Female XX can have different expected counts
• Telomere analyses need special handling

3. Mixing haploid and diploid across analyses

• Be consistent within each analysis
• Document which count you're using

4. Forgetting about polyploids

• Some species are naturally 3n, 4n, 6n, 8n
• Check literature for ploidy level

Key Takeaways

1. Always verify which count (n or 2n) a column contains 2. Don't trust column names - validate with known species 3. Use haploid (n) for per-assembly metrics 4. Add validation checks to catch errors early 5. Document which count you're using in code comments 6. Account for sex chromosomes when relevant

---

Phylogenetic Tree Species Mapping

Time Tree Species Replacement

Time Tree databases sometimes use proxy/replacement species when they don't have phylogenetic data for the exact species needed. This creates a mismatch between tree species names and dataset species names.

Pattern:

• Tree contains: Anniella_pulchra (proxy species with available data)
• Dataset contains: Anniella_stebbinsi (actual species being studied)
• Time Tree selected Anniella_pulchra as closest relative with data

Solution Workflow:

1. Document replacements in species_replacements.json:

{
  "actual_species_name": "tree_proxy_name",
  "Anniella_stebbinsi": "Anniella_pulchra",
  "Pelomedusa_somalica": "Pelomedusa_subrufa"
}

2. Update tree file to use actual dataset names:

• Read Newick tree file
• Replace proxy names with actual species names
• Ensures tree matches dataset exactly

3. Synchronize all config files using actual names:

• iTOL colorstrip configs
• Label configs
• Any taxonomic annotation files

4. Recover missing data if needed:

• Check deprecated datasets for actual species
• Proxy species indicates actual species likely exists in data
• Add to current dataset after recovery

Why This Matters:

• Prevents "missing species" that actually exist in dataset
• Ensures tree and dataset species names match exactly
• Required for iTOL visualization configs to work correctly
• Improves tree coverage metrics (e.g., 506->508 species)

Common Files Needing Synchronization:

• Tree_final.nwk - Main phylogenetic tree
• itol_taxonomic_colorstrip_final.txt - Taxonomic annotations
• species_*_methods.csv - Species classification configs
• All iTOL visualization config files

Tree Coverage Analysis Pattern

When reconciling phylogenetic trees with species datasets:

Coverage Metric:

Coverage = (Species in both tree AND dataset) / (Total species in tree) x 100%

Identifying Missing Species:

1. Extract species from tree (Newick format):

with open('Tree_final.nwk', 'r') as f:
    tree_content = f.read()
# Extract species names (underscored format)
tree_species = set(re.findall(r'([A-Z][a-z]+_[a-z]+)', tree_content))

2. Extract species from dataset:

df = pd.read_csv('species_methods.csv')
dataset_species = set(df['Species'].str.replace(' ', '_'))

3. Find missing species:

missing = tree_species - dataset_species

4. Categorize missing species:

• Recoverable: Time Tree replacements or in deprecated datasets
• Phylogenetic context: Tree-only species for evolutionary context
• Unknown curation: In dataset but cannot classify

Recovery Workflow:

# Check if missing species are Time Tree replacements
replacements = json.load(open('species_replacements.json'))
for species in missing:
    tree_name = species.replace('_', ' ')
    if tree_name in replacements.values():
        actual_name = [k for k,v in replacements.items() if v==tree_name][0]
        # Search deprecated datasets for actual_name
        # Recover and add to current dataset

Acceptable Coverage Levels:

• 100%: Ideal, all tree species have data
• 99%+: Excellent, few phylogenetic context species
• 95-99%: Good, some context species expected
• <95%: Investigate missing species for recovery opportunities

Example Results:

• Initial: 506/511 species (99.0%)
• After Time Tree mapping: 508/511 (99.4%)
• Remaining 3: Phylogenetic context only (acceptable)

---

BED File Processing and Telomere Analysis

Pattern: Classifying Scaffolds by Telomere Types

When analyzing telomere data from BED files to classify scaffolds:

File Structure:

• Terminal telomeres BED: columns include scaffold, start, end, orientation (p/q), accession
• Interstitial telomeres BED: similar structure with position markers (p/q/u for internal)

Best Practice - Use Python CSV Module:

import csv
from collections import defaultdict

# Use defaultdict for automatic initialization
telomere_counts = defaultdict(lambda: {'terminal': 0, 'interstitial': 0})

# Process with csv.reader (more portable than pandas)
with open('telomeres.bed', 'r') as f:
    reader = csv.reader(f, delimiter='\t')
    for row in reader:
        scaffold = row[0]
        accession = row[10]  # GCA accession
        key = (accession, scaffold)
        telomere_counts[key]['terminal'] += 1

Why CSV over pandas:

• No external dependencies (pandas may not be installed)
• Faster for simple tabular operations
• Lower memory footprint for large files
• Better portability across environments

Classification Categories: 1. Category 1: 2 terminal telomeres, 0 interstitial (complete chromosomes) 2. Category 2: 1 terminal telomere, 0 interstitial (partial) 3. Category 3: Has interstitial telomeres (likely assembly issues)

---

NCBI Data Integration Strategies

Check Existing Data Sources Before API Calls

Problem: Need chromosome counts for 400+ assemblies from NCBI.

Anti-pattern: Query NCBI datasets API for each accession

# DON'T: Query 400+ times
for accession in missing_data:
    result = subprocess.run(['datasets', 'summary', 'genome', 'accession', accession])
    # Takes 10+ minutes, hits API rate limits

Better Pattern: Check if data already exists in compiled tables

# DO: Look for existing compiled data first
# VGP table has multiple chromosome count columns:
# - num_chromosomes (column 54)
# - total_number_of_chromosomes (column 106)
# - num_chromosomes_haploid (column 122)

# Read from existing comprehensive table
with open('VGP-table.csv') as f:
    reader = csv.reader(f)
    header = next(reader)
    for row in reader:
        num_chr = row[53] if row[53] else row[105]  # Fallback strategy

Results: Filled 392/417 missing values instantly vs 10+ minutes of API calls.

Fallback Strategy for Multiple Columns:

# Try multiple sources in order of preference
num_chromosomes = row[53] if (len(row) > 53 and row[53]) else ''
if not num_chromosomes and len(row) > 105:
    num_chromosomes = row[105]  # Alternative column

When to use NCBI API:

• Data not in existing tables
• Need real-time/latest data
• Fetching assembly reports or sequence data
• Small number of queries (<20)

API Best Practices (when necessary):

• Use full path to datasets command (may be aliased)
• Add delays between calls (time.sleep(0.5))
• Set reasonable timeouts
• Handle errors gracefully

Bioinformatics Fundamentals - Reference Documentation

Detailed technical specifications, complete tables, and reference material for bioinformatics file formats and tools.

---

SAM/BAM Format Complete Reference

Complete SAM Flag Table

Common Flag Combinations

SAM Mandatory Fields

CIGAR Operations Complete Table

Optional Tags (Common)

---

FASTQ Format Specification

Format Structure

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

Line 1: Sequence identifier (starts with @) Line 2: Raw sequence Line 3: Separator (starts with +, optionally repeats identifier) Line 4: Quality scores (same length as sequence)

Quality Score Encodings

Phred+33 (Sanger, Illumina 1.8+)

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJ
|                                        |
0                                       41

Formula: Q = ASCII - 33 Range: 0-41 (sometimes extends to ~60)

Phred+64 (Illumina 1.3-1.7)

@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefgh
|                                        |
0                                       41

Formula: Q = ASCII - 64 Range: 0-41

Quality Score Interpretation

---

Tool Command Reference

samtools view

Complete Syntax:

samtools view [options] <in.bam>|<in.sam>|<in.cram> [region...]

Common Options:

-b          # Output BAM
-C          # Output CRAM
-h          # Include header
-H          # Header only
-c          # Count only
-o FILE     # Output file
-U FILE     # Output unselected reads
-f INT      # Required flags (include)
-F INT      # Filter flags (exclude)
-q INT      # Min MAPQ
-L FILE     # Regions from BED file
-r STR      # Read group
-R FILE     # Read group from file
-d TAG:VAL  # Tag filtering
-D TAG:FILE # Tag from file
-s FLOAT    # Subsample fraction
--threads N # Number of threads

Examples:

# Extract properly paired reads
samtools view -b -f 2 input.bam > proper_pairs.bam

# Exclude unmapped and secondary alignments
samtools view -b -F 4 -F 256 input.bam > mapped_primary.bam

# High quality mappings only
samtools view -b -q 30 input.bam > high_qual.bam

# Reads in specific region
samtools view -b input.bam chr1:1000-2000 > region.bam

# Subsample 10% of reads
samtools view -b -s 0.1 input.bam > subsample.bam

# Count mapped reads
samtools view -c -F 4 input.bam

samtools flagstat

Purpose: Count reads by flag status

Output Format:

12345 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
456 + 0 duplicates
11890 + 0 mapped (96.31% : N/A)
12345 + 0 paired in sequencing
6789 + 0 read1
5556 + 0 read2
11234 + 0 properly paired (90.99% : N/A)
11400 + 0 with itself and mate mapped
490 + 0 singletons (3.97% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools stats

Purpose: Comprehensive BAM statistics

Usage:

samtools stats [options] <in.bam>

# Common options
-r REF.fa   # Reference file
-c INT      # Coverage cap
-d INT      # Maximum coverage depth
--threads N # Threads

Key Output Sections:

• Summary numbers (SN)
• First fragment qualities (FFQ)
• Last fragment qualities (LFQ)
• GC content (GCC)
• Insert sizes (IS)
• Read lengths (RL)
• Indel distribution (ID)
• Coverage distribution (COV)

bamtools filter

Purpose: Advanced filtering

Usage:

bamtools filter -in <input.bam> [filter options]

# Filter options
-mapQuality ">30"
-isPaired true
-isProperPair true
-isMapped true
-isDuplicate false
-isReverseStrand true
-tag "RG:sample1"
-insertSize ">=100"

Filter File (JSON):

{
  "filters": [
    {
      "mapQuality": ">=30",
      "isPaired": true,
      "isProperPair": true
    }
  ]
}

bamtools filter -in input.bam -script filter.json -out output.bam

---

Sequencing Technology Details

PacBio HiFi Technical Specs

Platform: PacBio Sequel II, Sequel IIe, Revio

Chemistry:

• SMRTbell template preparation
• Circular Consensus Sequencing (CCS)
• Multiple passes over same molecule

Read Characteristics:

• Length: 10-25 kb (mode ~15 kb)
• Accuracy: >99.9% (Q20+), often Q30+
• Error mode: Random (not systematic)
• No GC bias
• Can sequence through modifications

Recommended Coverage:

• De novo assembly: 30-50x
• Variant calling: 20-30x
• Isoform sequencing: Depends on expression

Quality Metrics:

• Accuracy (CCS passes): More passes = higher quality
• Predicted accuracy in Phred scale (rq tag)
• Read length distribution

Hi-C Technical Specs

Protocol Steps: 1. Crosslink chromatin with formaldehyde 2. Digest with restriction enzyme 3. Fill in and label ends with biotin 4. Ligate (creates chimeric molecules) 5. Shear DNA 6. Pull down biotinylated junctions 7. Paired-end sequencing

Read Characteristics:

• Paired-end: 100-150 bp each end
• R1 and R2 from same ligation product
• Can be on different chromosomes/scaffolds
• Many "invalid pairs" (self-ligations, etc.)

Quality Metrics:

• Valid pairs percentage
• Cis vs trans ratio (intra vs inter-chromosomal)
• Contact distance distribution
• Coverage uniformity

Expected Pair Types:

• Valid pairs (useful): ~40-70%
• Self-circles: ~10-20%
• Dangling ends: ~10-20%
• Other invalid: ~10-30%

Illumina Technical Specs

Platforms: NovaSeq, NextSeq, HiSeq, MiSeq

Chemistry:

• Sequencing by synthesis (SBS)
• Clonal amplification (bridge PCR)
• Four-color imaging

Read Characteristics:

• Length: 50-300 bp (platform dependent)
• Paired-end or single-end
• Quality decreases toward 3' end
• Systematic errors possible (GGC motif)

Quality Metrics:

• Cluster density
• %PF (passing filter)
• Q30 percentage (% bases >Q30)
• Index balance (for multiplexing)

---

Assembly Quality Metrics

Contiguity Metrics

N50:

• Sort contigs by length (largest first)
• Sum lengths until reaching 50% of total assembly
• N50 = length of contig at 50% mark

L50:

• Number of contigs needed to reach N50

N90:

• Same as N50 but using 90% threshold
• More stringent

NG50:

• N50 relative to expected genome size
• Better for comparing assemblies

auN:

• Area under Nx curve
• Less sensitive to individual long contigs
• Better for comparing fragmented assemblies

Completeness Metrics

BUSCO (Benchmarking Universal Single-Copy Orthologs):

Complete: 95.2% (C:95.2%[S:94.1%,D:1.1%],F:2.3%,M:2.5%,n:3950)
- Complete and single-copy (S): 94.1%
- Complete and duplicated (D): 1.1%
- Fragmented (F): 2.3%
- Missing (M): 2.5%

Interpretation:

• >95% complete: Excellent
• 90-95% complete: Good
• <90% complete: May have issues

QV (Consensus Quality Value):

• Phred-scaled accuracy of consensus sequence
• QV30 = 99.9% accurate (1 error per 1000 bp)
• QV40 = 99.99% accurate (1 error per 10,000 bp)
• QV50 = 99.999% accurate (1 error per 100,000 bp)

Structural Metrics

LAI (LTR Assembly Index):

• For plant genomes with LTR retrotransposons
• Scale 0-100
• >20 = excellent continuity

BUSCO Structural:

• Checks for fragmentation of conserved genes
• High duplication may indicate haplotigs

---

Coverage and Depth Calculations

Theoretical Coverage

Formula:

Coverage (X) = (Number of reads × Read length) / Genome size

Example:

• 100 million reads
• 150 bp read length
• 3 Gb genome

Coverage = (100M × 150) / 3G = 15,000M / 3,000M = 5X

Effective Coverage

Account for duplicates, unmapped, low quality:

Effective coverage = Theoretical coverage × (1 - duplicate rate) × mapping rate

Recommended Coverage Levels

---

Coordinate Systems

0-based vs 1-based

1-based (SAM, VCF, GFF):

Sequence: A T C G A T C G
Position: 1 2 3 4 5 6 7 8

• First base is position 1
• Interval [2,5] includes bases at positions 2,3,4,5

0-based (BED, BAM binary):

Sequence: A T C G A T C G
Position: 0 1 2 3 4 5 6 7

• First base is position 0
• Interval [2,5) includes bases at positions 2,3,4 (excludes 5)

0-based half-open [start, end):

• BED format
• start included, end excluded
• Length = end - start

1-based closed [start, end]:

• SAM format
• Both start and end included
• Length = end - start + 1

Conversion

BED to SAM:

SAM_start = BED_start + 1
SAM_end = BED_end

SAM to BED:

BED_start = SAM_start - 1
BED_end = SAM_end

---

AGP Format Complete Reference

Overview

AGP (A Golden Path) is a tab-delimited text format that describes the assembly of larger sequence objects (chromosomes, scaffolds) from smaller components (contigs, scaffolds) and gaps.

Official Specification: https://www.ncbi.nlm.nih.gov/genbank/genome_agp_specification/

AGP Line Types

AGP files contain two types of lines: 1. Sequence lines (component_type = 'W'): Describe actual sequence components 2. Gap lines (component_type = 'N' or 'U'): Describe gaps between components

Sequence Line Format (9 columns)

Orientation Values:

• +: Component in same orientation as object
• -: Component reverse complemented relative to object
• ?: Unknown orientation
• 0: Unspecified (deprecated)
• na: Not applicable (for single-stranded sequences)

Gap Line Format (9+ columns)

Gap Types:

• scaffold: Gap within scaffold
• contig: Gap within contig (rare)
• centromere: Centromeric gap
• short_arm: Short arm of acrocentric chromosome
• heterochromatin: Heterochromatic gap
• telomere: Telomeric gap
• repeat: Gap due to repeat

Linkage Evidence Types:

• paired-ends: Paired read evidence
• align_genus: Alignment to related species
• align_xgenus: Alignment to different genus
• align_trnscpt: Alignment to transcript
• within_clone: Same clone
• clone_contig: Clone and contig evidence
• map: Genetic/physical map
• strobe: Strobe sequencing
• proximity_ligation: Hi-C or similar

Complete AGP Example

##agp-version 2.1
# ORGANISM: Genus species
# DESCRIPTION: Curated assembly
chr1    1       5000    1       W       contig_1        1       5000    +
chr1    5001    5100    2       U       100     scaffold        yes     proximity_ligation
chr1    5101    15000   3       W       contig_2        1       9900    -
chr1    15001   15100   4       N       100     scaffold        yes     paired-ends
chr1    15101   25000   5       W       contig_3        1       9900    +
chr1_unloc_1    1       3000    1       W       contig_4        1       3000    +
chr2    1       8000    1       W       contig_5        1       8000    -

Critical Validation Rules

Rule 1: Length Consistency

For sequence lines (type W):

object_end - object_beg + 1 == component_end - component_beg + 1

Valid:

chr1    1000    2999    1       W       ctg1    1       2000    +
# Object length: 2999 - 1000 + 1 = 2000 ✓
# Component length: 2000 - 1 + 1 = 2000 ✓

Invalid:

chr1    1000    5000    1       W       ctg1    1       2000    +
# Object length: 5000 - 1000 + 1 = 4001 ✗
# Component length: 2000 - 1 + 1 = 2000 ✗
# ERROR: Lengths don't match!

Rule 2: Sequential Part Numbers

Part numbers (column 4) must:

• Start at 1 for each new object
• Increment by 1 for each subsequent line of same object
• No gaps or duplicates

Valid:

chr1    1       1000    1       W       ctg1    1       1000    +
chr1    1001    1100    2       U       100     scaffold        yes     paired-ends
chr1    1101    2000    3       W       ctg2    1       900     -

Invalid:

chr1    1       1000    1       W       ctg1    1       1000    +
chr1    1001    1100    3       U       100     scaffold        yes     paired-ends  # ✗ Skipped 2
chr1    1101    2000    3       W       ctg2    1       900     -  # ✗ Duplicate 3

Rule 3: Coordinate Continuity

Object coordinates must be continuous with no gaps or overlaps:

• Next object_beg = Previous object_end + 1

Valid:

chr1    1       1000    1       W       ctg1    1       1000    +
chr1    1001    1100    2       U       100     scaffold        yes     paired-ends
chr1    1101    2000    3       W       ctg2    1       900     -

Invalid:

chr1    1       1000    1       W       ctg1    1       1000    +
chr1    1002    1100    2       U       100     scaffold        yes     paired-ends  # ✗ Gap at 1001

Rule 4: Component Usage

• Each component region (component_id:component_beg-component_end) should appear only once in the assembly
• Exception: Tandem repeats or duplications may appear multiple times if biologically accurate

AGP Processing Patterns

Pattern 1: Extracting Component Length

def get_component_length(agp_line):
    """Calculate length from AGP sequence line."""
    obj_beg, obj_end = int(agp_line[1]), int(agp_line[2])
    comp_beg, comp_end = int(agp_line[6]), int(agp_line[7])

    obj_length = obj_end - obj_beg + 1
    comp_length = comp_end - comp_beg + 1

    assert obj_length == comp_length, "Length mismatch!"
    return obj_length

Pattern 2: Creating New AGP Object

def create_unloc_line(parent_line, unloc_name):
    """Create AGP line for unlocalized scaffold."""
    # Extract component coordinates from parent
    comp_id = parent_line[5]
    comp_beg = int(parent_line[6])
    comp_end = int(parent_line[7])
    orientation = parent_line[8]

    # Calculate length
    length = comp_end - comp_beg + 1

    # Create new AGP line
    return [
        unloc_name,          # object
        1,                   # object_beg (always 1)
        length,              # object_end (equals length)
        1,                   # part_number (reset to 1)
        'W',                 # component_type
        comp_id,             # component_id
        comp_beg,            # component_beg (preserve)
        comp_end,            # component_end (preserve)
        orientation          # orientation (preserve)
    ]

Pattern 3: Validating AGP Coordinates

def validate_agp_line(line):
    """Validate AGP line coordinates."""
    if line[4] == 'W':  # Sequence line
        obj_beg, obj_end = int(line[1]), int(line[2])
        comp_beg, comp_end = int(line[6]), int(line[7])

        obj_length = obj_end - obj_beg + 1
        comp_length = comp_end - comp_beg + 1

        if obj_length != comp_length:
            raise ValueError(
                f"Length mismatch: object {obj_length} bp, "
                f"component {comp_length} bp"
            )

        if obj_beg < 1 or comp_beg < 1:
            raise ValueError("Coordinates must be >= 1 (1-based)")

    elif line[4] in ['N', 'U']:  # Gap line
        obj_beg, obj_end = int(line[1]), int(line[2])
        gap_length = int(line[5])

        obj_length = obj_end - obj_beg + 1

        if obj_length != gap_length:
            raise ValueError(
                f"Gap length mismatch: object span {obj_length} bp, "
                f"specified gap {gap_length} bp"
            )

Pattern 4: Splitting AGP by Object

def split_agp_by_object(agp_file):
    """Split AGP into separate files per object."""
    objects = {}

    with open(agp_file) as f:
        for line in f:
            if line.startswith('#'):
                continue

            parts = line.strip().split('\t')
            obj_name = parts[0]

            if obj_name not in objects:
                objects[obj_name] = []

            objects[obj_name].append(parts)

    return objects

Common AGP Errors and Fixes

Error: "object and component coordinates do not have the same length"

Cause: Object span ≠ component span

Fix:

# WRONG: Using component end coordinate as object end
obj_end = comp_end

# CORRECT: Calculate length and use that
length = comp_end - comp_beg + 1
obj_end = obj_beg + length - 1

Error: "part number is not sequential"

Cause: Part numbers have gaps or aren't incrementing

Fix:

# Track and reset part numbers per object
current_object = None
part_num = 1

for line in agp_lines:
    obj_name = line[0]

    if obj_name != current_object:
        current_object = obj_name
        part_num = 1  # Reset for new object

    line[3] = part_num
    part_num += 1

Error: "coordinates are not continuous"

Cause: Gap or overlap between adjacent lines

Fix:

# Ensure continuity
prev_end = 0

for line in agp_lines:
    if line[0] == current_object:  # Same object
        expected_beg = prev_end + 1
        line[1] = expected_beg
        line[2] = expected_beg + length - 1
        prev_end = line[2]
    else:  # New object
        line[1] = 1
        prev_end = line[2]

AGP Validation Tools

NCBI AGP Validator

# Download validator
wget https://ftp.ncbi.nlm.nih.gov/toolbox/ncbi_tools/converters/by_program/agp_validate/linux64.agp_validate.gz
gunzip linux64.agp_validate.gz
chmod +x linux64.agp_validate

# Validate AGP
./linux64.agp_validate -assembly assembly.agp -fasta assembly.fasta

# Common options
-o output.txt          # Write report to file
-euk                   # Eukaryotic assembly (default)
-prok                  # Prokaryotic assembly
-chr2scaf chr2scaf.txt # Chromosome to scaffold mapping

Quick Python Validation

def quick_validate_agp(agp_file):
    """Quick validation checks."""
    with open(agp_file) as f:
        prev_obj = None
        prev_end = 0
        part_num = 0

        for line_num, line in enumerate(f, 1):
            if line.startswith('#'):
                continue

            parts = line.strip().split('\t')
            obj, obj_beg, obj_end, part = parts[0:4]
            obj_beg, obj_end, part = int(obj_beg), int(obj_end), int(part)

            # Check object continuity
            if obj == prev_obj:
                if obj_beg != prev_end + 1:
                    print(f"Line {line_num}: Coordinate gap or overlap")
                if part != part_num + 1:
                    print(f"Line {line_num}: Part number not sequential")
            else:
                if obj_beg != 1:
                    print(f"Line {line_num}: Object doesn't start at 1")
                if part != 1:
                    print(f"Line {line_num}: Part number doesn't start at 1")

            # Check length consistency for sequence lines
            if parts[4] == 'W':
                comp_beg, comp_end = int(parts[6]), int(parts[7])
                obj_len = obj_end - obj_beg + 1
                comp_len = comp_end - comp_beg + 1
                if obj_len != comp_len:
                    print(f"Line {line_num}: Length mismatch ({obj_len} vs {comp_len})")

            prev_obj = obj
            prev_end = obj_end
            part_num = part

AGP Coordinate System Summary

• 1-based: Both object and component coordinates start at 1
• Inclusive: Both start and end positions are included
• Closed interval: [start, end] notation
• Length formula: end - start + 1

Related Formats

Comparison with other coordinate formats:

---

Error Rates and Sequencing Technology

Error Profiles

Base Quality Distributions

Illumina:

• High quality at 5' end (Q35-40)
• Gradual decline toward 3' end
• Last 10-20 bases often Q20-30

PacBio HiFi:

• Consistent across read
• Q30+ typical
• Length-independent quality

ONT:

• Variable by position
• Homopolymers lower quality
• Improving with newer chemistry

---

File Size Estimates

Compression Ratios

Example File Sizes (30x WGS, 3Gb genome)

---

Useful Regular Expressions

FASTQ Header Parsing

Illumina:

@INSTRUMENT:RUN:FLOWCELL:LANE:TILE:X:Y READ:FILTERED:CONTROL:BARCODE

Regex:

^@([^:]+):(\d+):([^:]+):(\d+):(\d+):(\d+):(\d+) ([12]):([YN]):(\d+):([ACTGN+]+)$

FASTA Header Parsing

^>(\S+)\s*(.*)$
# Group 1: Sequence ID
# Group 2: Description

CIGAR String Parsing

(\d+)([MIDNSHP=X])
# Group 1: Length
# Group 2: Operation

---

Additional Resources

• SAM Specification: https://samtools.github.io/hts-specs/SAMv1.pdf
• SAM Tags Specification: https://samtools.github.io/hts-specs/SAMtags.pdf
• VCF Specification: https://samtools.github.io/hts-specs/VCFv4.3.pdf
• FASTQ Format: https://en.wikipedia.org/wiki/FASTQ_format
• Phred Quality Scores: https://www.drive5.com/usearch/manual/quality_score.html
• BUSCO: https://busco.ezlab.org/
• samtools documentation: http://www.htslib.org/doc/