---
title: "Case study: hare cytochrome b fragments"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Case study: hare cytochrome b fragments}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6)
library(genoaligner)
```

Public sequence databases give you *fragments*, not tidy columns of
equal-length strings. Records differ in length and in which part of the gene
they cover, and a few are not the gene you searched for at all. This case
study takes real cytochrome b (cytb) records for three hare species —
*Lepus granatensis*, *L. timidus* and *L. europaeus*, all of which the
[EcoSeek genominer](https://genominer.ecoseek.org) catalogue flags as having
population-genetic data — plus the rabbit *Oryctolagus cuniculus* as an
outgroup, and asks:

> Can bounded edit distance recover the species clusters, and what does the
> `resolved`/`smax` semantics buy us along the way?

The dataset ships with the package, exactly as downloaded from NCBI
(see `?lepus_cytb` for provenance):

```{r}
data(lepus_cytb)
table(lepus_cytb$species)
range(nchar(lepus_cytb$sequence))
```

Thirty-seven records, 400–1200 bp, mixed coverage — the raw material of every
barcode or marker survey.

## QC pass: which records actually are cytb?

A first `align_sw()` against the complete reference gene (accession
HQ596476.1, *L. granatensis* voucher EBD-L094) answers two questions at once:
does this record overlap the locus at all, and where?

```{r}
ref <- lepus_cytb$sequence[lepus_cytb$accession == "HQ596476.1"]
qc <- align_sw(ref, lepus_cytb$sequence, scoring = c(2, -3, 5, 2))
lepus_cytb$sw_score <- qc$score
lepus_cytb$span_start <- qc$start_i   # 0-based, on the reference
lepus_cytb$span_end <- qc$end_i

lepus_cytb[order(lepus_cytb$sw_score),
           c("accession", "species", "sw_score", "span_start", "span_end")][1:4, ]
```

The two lowest scores are not low — they are *nothing*: spans of ~10 bp with
score 22. XM_017337704.3 and XM_051846530.2 are predicted **nuclear** mRNAs
that a "cytochrome b" title search returns anyway. A near-zero local score is
the flag: whatever these records are, they are not this locus. Every pipeline
has this problem; here it is visible instead of silently folded into a
distance matrix.

```{r}
cytb <- lepus_cytb[lepus_cytb$sw_score > 500, ]
nrow(cytb)
```

## A common window, honestly chosen

Fragments that do survive QC still cover different parts of the gene:

```{r}
plot(cytb$span_start, cytb$span_end, pch = 19,
     xlab = "aligned start on reference (bp)", ylab = "aligned end (bp)")
abline(h = max(cytb$span_start), v = min(cytb$span_end), lty = 2)
```

Two records start around position 708 — much later than everything else —
so the intersection over *all* records is a useless ~24 bp. Dropping the
3'-end pair leaves a 617 bp window shared by the remaining 33 records:

```{r}
keep <- cytb[cytb$span_start < 400, ]
w <- c(max(keep$span_start), min(keep$span_end))
w
dropped <- setdiff(cytb$accession, keep$accession)
dropped
```

Cropping each record to that window is itself a local alignment: align the
window to the record and take the aligned span on the record side.

```{r}
window_seq <- substr(ref, w[1] + 1, w[2] + 1)   # 0-based coords -> 1-based substr
m <- align_sw(window_seq, keep$sequence, scoring = c(2, -3, 5, 2))
keep$frag <- mapply(substr, keep$sequence, m$start_j + 1, m$end_j + 1)
range(nchar(keep$frag))
```

## Edit distance recovers the species

Now the question of the study: pick one reference haplotype per species and
ask every fragment how far it sits from each.

```{r}
pick <- function(acc) keep$frag[keep$accession == acc]
refs <- c(granatensis = "JF299036.1", timidus = "LC132658.1",
          europaeus = "MN098896.1", cuniculus = "AY292717.1")

for (sp in names(refs))
    keep[[paste0("d_", sp)]] <-
        align_edit(keep$frag, pick(refs[sp]), smax = 511, with_cigar = FALSE)$score

d <- aggregate(cbind(d_granatensis, d_timidus, d_europaeus, d_cuniculus) ~ species,
               keep, median)
print(d, right = FALSE)
```

The diagonal is the within-species scale (0–8 edits); off-diagonal entries are
the between-species scale (~55–64 within *Lepus*, ~100+ to the rabbit). The
**barcode gap** — the empty interval between "same species" and "different
species" distances — is the whole reason cytb works as a marker, and here it
is computed with three calls.

```{r}
ids <- keep$accession
within <- unlist(lapply(seq_len(nrow(keep)), function(i)
    align_edit(keep$frag[i], keep$frag[keep$species == keep$species[i] & ids != ids[i]],
          smax = 511, with_cigar = FALSE)$score))
between <- unlist(lapply(seq_len(nrow(keep)), function(i)
    align_edit(keep$frag[i], keep$frag[keep$species != keep$species[i]],
          smax = 511, with_cigar = FALSE)$score))

br <- seq(0, max(c(within, between)) + 5, by = 5)
hw <- hist(within, breaks = br, plot = FALSE)
hb <- hist(between, breaks = br, plot = FALSE)
plot(hb, col = adjustcolor("tomato", 0.6), border = "white",
     ylim = c(0, max(hw$counts, hb$counts)),
     main = "edit distance over the 617 bp window", xlab = "edits")
plot(hw, col = adjustcolor("steelblue", 0.6), border = "white", add = TRUE)
legend("topleft", c("between species", "within species"),
       fill = c(adjustcolor("tomato", 0.6), adjustcolor("steelblue", 0.6)))
```

Two honest caveats the numbers surface for free: the *Oryctolagus* within
group splits into two mtDNA lineages (~0 and ~55 edits — a documented
polymorphism in the rabbit), and *timidus* isolate ALT2 sits ~48 edits from
the European *timidus* references, a divergent Asian lineage. Edit distance
does not paper over structure that is really there.

## `smax` is a claim about relatedness

At the default `smax = 64` the outgroup comparisons simply do not resolve:

```{r}
align_edit(keep$frag, pick(refs["cuniculus"]), smax = 40, with_cigar = FALSE)$resolved |> table()
```

Twenty-eight unresolved pairs is not a failure — it is the library saying
"these are more than 40 edits apart" instead of returning 28 wrong numbers.
Raise the bound and they resolve, with the true scores:

```{r}
align_edit(keep$frag, pick(refs["cuniculus"]), smax = 200, with_cigar = FALSE)$score |> range()
```

The bound is a hypothesis: *if two sequences are related at the scale you care
about, they resolve; if not, you hear about it.* Choosing `smax` is choosing
where the barcode gap ends — cheap to compute, impossible to get silently
wrong.

## Takeaway

Starting from a naive NCBI dump — mixed lengths, mixed coverage, two
mislabelled records — three `align_sw`/`align` calls produced: a locus-QC
filter, a defensible common window, and a distance table in which the species
clusters are unambiguous. The flags `resolved`, `rescore_ok` and
`wellformed_ok` travel with every row, so downstream code never has to guess
whether a number is real.