How genevintage works

Detecting species

detect_species() first works out what fraction of the input each candidate species’ identifier prefix explains:

ids <- c(
  "ENSMUSG00000051951", "ENSMUSG00000089699", "ENSMUSG00000051951",
  "TP53"
)
cand <- detect_species(ids)
head(cand[, c("species", "prefix", "frac")])
#>                  species  prefix frac
#> 1           mus_musculus ENSMUSG 0.75
#> 2      mus_musculus_cbaj ENSMUSG 0.75
#> 3     mus_musculus_fvbnj ENSMUSG 0.75
#> 4    mus_musculus_c3hhej ENSMUSG 0.75
#> 5 mus_musculus_nodshiltj ENSMUSG 0.75
#> 6     mus_musculus_dba2j ENSMUSG 0.75

Three of the four inputs (one a repeat) carry the ENSMUSG prefix, so every mouse strain that shares it ties at frac = 0.75. Ties are broken in order by:

1. Prefix length. A longer matching prefix ranks higher.
2. reach. Could the species have issued the input’s highest identifier number? The largest number in the input is compared with the largest qmax (the highest identifier number in a gene set) ever indexed for the species, within QMAX_TOLERANCE.
3. fit. How close is the input’s gene count to a real gene set? It is the smallest log-scale gap between the number of unique IDs and any indexed release’s gene count. It is only informative for a full gene set.

Any remaining tie is broken by species name, for a stable order; such species are indistinguishable from the identifiers alone. reach and fit are used for ranking only and are not returned.

Determining the annotation release

detect_release() compares the input against the fingerprint index and rules out releases that could not have produced it, using two facts:

1. Identifier numbers (usually) only go up, so the highest number in the input is a floor, i.e. no earlier release would have it
2. A gene’s version only goes up too, so the highest version in the input is also a floor.

QMAX_TOLERANCE lets a release’s qmax fall slightly short of the input’s highest number. Ensembl sometimes retires its top-numbered genes, which lowers qmax without changing the assembly. qmax_violations() lists every case in the shipped index:

qmax_violations()
#> [1] "camarhynchus_parvulus 11: qmax falls sharply at release 108 -- check the assembly labels" 
#> [2] "corvus_moneduloides 1: qmax falls sharply at release 108 -- check the assembly labels"    
#> [3] "melopsittacus_undulatus 1: qmax falls sharply at release 108 -- check the assembly labels"
#> [4] "ovis_aries 31: qmax falls sharply at release 114 -- check the assembly labels"            
#> [5] "takifugu_rubripes 4: qmax falls sharply at release 49 -- check the assembly labels"

Where a violation exists, detect_release() skips the floor for that species, since the check would fail against the species’ own history.

Among the releases left feasible, distance ranks the rest:

ids2 <- readLines(system.file("extdata", "example_ids.txt.gz",
  package = "genevintage"
))
head(detect_release(ids2), 3)
#>        species  source release assembly feasible   dist n_index
#> 1 homo_sapiens gencode      44       38     TRUE 1.1551   62700
#> 2 homo_sapiens gencode      46       38     TRUE 1.1582   63086
#> 3 homo_sapiens gencode      45       38     TRUE 1.1599   63187

dist is the absolute difference in log1p(gene count) between the input and each candidate release, so a thousand genes matter less in human (tens of thousands of genes) than in a species with a few hundred. For species whose identifier numbers are ordered, it adds the same comparison on the 90th percentile of identifier numbers, which separates releases that share a gene count but use different numbers.

Most non-Ensembl schemes carry no order in their identifier numbers. Such schemes skip the numeric floor and the quantile term entirely:

head(detect_release(c("YAL001C", "YAL002W", "YAL003W"), species = "yeast"), 5)
#> sgd identifiers are not numerically ordered; ranking by gene count alone.
#>                    species  source release assembly feasible   dist n_index
#> 1 saccharomyces_cerevisiae ensembl      94        4       NA 7.4110    6615
#> 2 saccharomyces_cerevisiae ensembl      95        4       NA 7.4110    6615
#> 3 saccharomyces_cerevisiae ensembl      96        4       NA 7.4110    6615
#> 4 saccharomyces_cerevisiae ensembl     100        4       NA 7.4246    6706
#> 5 saccharomyces_cerevisiae ensembl     101        4       NA 7.4246    6706

Handling GENCODE and Ensembl

Human and mouse are indexed against GENCODE as well as Ensembl, since their RNA-seq is often quantified against GENCODE, whose gene set differs in size from Ensembl’s at the same release. detect_release() ranks across both and reports which one each row came from in source.

Caution around new releases

The index covers a fixed range of releases. If your best match is the oldest or newest release in that range, treat it with caution. Your IDs may come from a release just outside the range. detect_release() warns when this happens.

head(detect_release(c("ENSG00000141510", "ENSG00000012048"), species = "human"), 3)
#>        species  source release assembly feasible   dist n_index
#> 1 homo_sapiens ensembl      48       36     TRUE 9.7936   32561
#> 2 homo_sapiens ensembl      49       36     TRUE 9.7936   32561
#> 3 homo_sapiens ensembl      47       36     TRUE 9.8311   33761