---
title: "Assigning a direction to CpG sites"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Assigning a direction to CpG sites}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
# The demonstration chunks need the package's data layers. They are present
# in the GitHub (fat) build, in local-data mode, and once the cpgdirectionData
# resources are available from ExperimentHub; on a machine without them the
# code is shown but not evaluated.
has_data <- cpgdirection::cpgd_has_data("lookup_blood_hg19") &&
  cpgdirection::cpgd_has_data("smr_directions") &&
  cpgdirection::cpgd_has_data("epicv2_probe_gene_annotation")
knitr::opts_chunk$set(eval = has_data)
```

*The code below is evaluated only when the package's data layers are
available -- bundled, via `options(cpgdirection.data_dir = )`, or through
the cpgdirectionData package. The text and code are complete either way.*

## The problem

An epigenome-wide study gives you CpG sites. It rarely gives you matched RNA, so
for most sites you do not know whether higher methylation goes with **higher**
or **lower** expression of the nearby gene. That sign is not fixed: methylation
is not uniformly repressive, and pooled estimates across studies can average to
zero simply because opposing signs cancel.

This package assigns the sign, and tells you how much to trust it.

## One call

```{r one-call}
library(cpgdirection)
my_cpgs <- c("Zcg00335286_TC21_MC2R", "cg02079741_TC21_POMC", "cg00000029")
res <- cpg_expression_direction(my_cpgs, universal = FALSE, verbose = FALSE)
res[, c("cpg_id", "best_direction", "best_evidence"), with = FALSE]
```

`my_cpgs` can be a character vector, a file path or a data frame column.
Identifiers carrying array suffixes or analysis prefixes are reduced to their
canonical form automatically, so `Zcg00335286_TC21_MC2R` is understood, and the
gene embedded in the name is used as the target. (`universal = FALSE` skips the
distance-only fallback layer here, which needs the Bioconductor gene
annotation packages on first use.)

## The two columns that matter

| column | meaning |
|---|---|
| `best_direction` | `+1`, `-1`, or `NA` where the evidence does not support a call |
| `best_evidence` | which layer supplied it |

**Always read them together.** A direction worth 0.96 and a direction worth 0.62
are both a `+1`. `best_expected_accuracy` carries the validated figure.

| `best_evidence` | accuracy | what it means |
|---|---|---|
| `measured` | ~1.00 | a catalogue measured this exact pair |
| `smr_high` | 0.96 | Mendelian randomisation, several instruments agreeing |
| `catalogue_consensus` | 0.77–0.87 | genomic-context models, tissues agreeing |
| `smr_moderate` | 0.85 | Mendelian randomisation, one instrument |
| `catalogue_single` | 0.55–0.87 | one tissue only |
| `smr_weak` | 0.70 | instruments disagree; a warning, not evidence |
| `distance_only` | 0.60–0.65 | distance to the TSS, nothing else |
| `no_evidence` | — | nothing available; a correct answer, not a failure |

## Every target gene, one row per pair

A CpG is not obliged to regulate one gene. `cpg_gene_pairs()` discovers every
supported target automatically -- bare identifiers suffice -- and resolves
direction separately for every CpG x gene pair, excluding probes flagged as
unreliable (cross-hybridizing or degenerately mapped) by default:

```{r gene-pairs}
p <- cpg_gene_pairs(c("cg26261055", "cg06495038"), verbose = FALSE)
p[, c("cpg_id", "target_gene", "mapping_primary", "best_direction",
      "best_evidence", "usable"), with = FALSE]
```

`mapping_sources` records how each gene entered the candidate set (manifest
annotation, catalogue, measured eQTM, SMR); `best_evidence` records what the
direction rests on. They are different facts and are kept apart.

## Abstention is deliberate

Where the evidence is insufficient the package declines rather than guessing.
That is why accuracy is high on the sites it does call. If you need a complete
column, name the tissue your question is about:

```{r target-tissue, eval = FALSE}
res <- cpg_expression_direction(my_cpgs, target_tissue = "solid_tissue")
```

Declare that tissue on biological grounds **before** looking at the answers.
Choosing the tissue that gives the result you wanted is selecting on the
outcome.

## Checking a result that looks one-sided

```{r direction-balance}
cpgd_direction_balance(res)
```

Inverse associations are the majority everywhere, so a mostly-`-1` result is
expected. This tests whether yours is more one-sided than the base rate explains.

## The limitation that governs everything

Every accuracy figure is **conditional on the CpG being an eQTM at all**. The
catalogues behind the models contain only significant associations, so nothing
here tells you whether a CpG affects expression. Feed it a site with no real
effect and you will get a confident direction for something that is not there.

## Session info

```{r session-info, eval = TRUE}
sessionInfo()
```