///|
/// ChIPseeker peak annotation demo.
///
/// Loads sample ChIP-seq peaks and gene annotations, annotates each peak
/// with its nearest gene and genomic feature, prints the per-peak
/// annotations, summarizes the feature distribution, and demonstrates
/// filtering by feature and by distance to TSS.
fn main {
println("=== ChIPseeker Peak Annotation Demo ===")
println("Inspired by the Bioconductor ChIPseeker package")
println("")
// 1. Load sample data.
println("1. Loading sample data...")
let peaks = @bio.chipseeker_sample_peaks()
let genes = @bio.chipseeker_sample_genes()
println(" Peaks loaded: \{peaks.length()}")
println(" Genes loaded: \{genes.length()}")
println("")
// 2. Annotate peaks. Use the default promoter window of 3000bp.
println("2. Annotating peaks (promoter window = 3000bp)...")
let annotations = @bio.annotate_peaks(peaks, genes, 3000)
println(" Annotated \{annotations.length()} peaks")
println("")
// 3. Display the first few per-peak annotations.
println("3. Per-peak annotations (first 6):")
println(" -----------------------------------------------")
let mut i = 0
while i < annotations.length() && i < 6 {
let ann = annotations[i]
let p = ann.get_peak()
println(
" Peak \{i + 1}: \{p.chr}:\{p.start}-\{p.end} (summit=\{p.summit})",
)
println(" Gene : \{ann.gene_name} (\{ann.gene_id})")
println(" TSS : \{ann.tss}")
println(" Distance : \{ann.distance_to_tss} (signed, + = downstream)")
println(" Feature : \{ann.feature}")
println(" Annotation: \{ann.annotation}")
println("")
i = i + 1
}
// 4. Summary statistics.
println("4. Annotation summary:")
println(" --------------------")
let summary = @bio.annotation_summary(annotations)
println(@bio.chipseeker_summary_string(summary))
// 5. Filter by feature: show promoter peaks.
println("5. Filtering peaks by feature: 'Promoter'")
let promoters = @bio.filter_by_feature(annotations, "Promoter")
println(" Promoter peaks: \{promoters.length()}")
for ann in promoters {
let p = ann.get_peak()
println(
" - \{p.chr}:\{p.summit} -> \{ann.gene_name} (dist=\{ann.distance_to_tss})",
)
}
println("")
// 6. Filter by feature: show exon peaks.
println("6. Filtering peaks by feature: 'Exon'")
let exons = @bio.filter_by_feature(annotations, "Exon")
println(" Exon peaks: \{exons.length()}")
for ann in exons {
let p = ann.get_peak()
println(
" - \{p.chr}:\{p.summit} -> \{ann.gene_name} (dist=\{ann.distance_to_tss})",
)
}
println("")
// 7. Filter by distance: peaks within 3000bp of a TSS.
println("7. Filtering peaks by distance: |distance| <= 3000")
let close = @bio.filter_by_distance(annotations, 3000)
println(" Close peaks: \{close.length()}")
for ann in close {
let p = ann.get_peak()
println(
" - \{p.chr}:\{p.summit} -> \{ann.gene_name} (dist=\{ann.distance_to_tss}, \{ann.feature})",
)
}
println("")
// 8. Demonstrate promoter window customization. With a smaller window
// (100bp), the 5'UTR of GeneA becomes accessible because the 5'UTR
// no longer sits entirely inside the promoter window.
println("8. Promoter window customization (window=100bp):")
println(" Annotating a peak inside the 5'UTR of GeneA (summit=10300)...")
let peak_in_5utr = @bio.Peak::new(
chr="chr1",
start=10200,
end=10400,
summit=10300,
)
let ann_small_window = @bio.annotate_peak(peak_in_5utr, genes, 100)
println(" Feature with window=100 : \{ann_small_window.feature}")
let ann_large_window = @bio.annotate_peak(peak_in_5utr, genes, 3000)
println(" Feature with window=3000: \{ann_large_window.feature}")
println(" (smaller window exposes the 5'UTR assignment)")
println("")
// 9. Edge case: peak on a chromosome with no genes.
println("9. Edge case: peak on a chromosome with no genes (chrZ)")
let orphan_peak = @bio.Peak::new(
chr="chrZ",
start=1000,
end=2000,
summit=1500,
)
let orphan_ann = @bio.annotate_peak(orphan_peak, genes, 3000)
println(" Gene id : '\{orphan_ann.gene_id}' (empty -> no nearby gene)")
println(" Feature : \{orphan_ann.feature}")
println(" Annotation: \{orphan_ann.annotation}")
println("")
println("=== ChIPseeker demo complete ===")
}