---
title: "Coordinate-level pharmacological annotation with drugannotR"
author: "drugannotR authors"
package: drugannotR
output:
  BiocStyle::html_document:
    toc_float: true
vignette: >
  %\VignetteIndexEntry{Coordinate-level pharmacological annotation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, message=FALSE}
library(drugannotR)
library(GenomicRanges)
```

# Scope

`drugannotR` provides coordinate-level pharmacological annotation. It reports
that a region overlaps a gene feature or an explicitly mapped variant connected
to a drug relationship. It does not infer that the region is functional,
predict response, rank drugs, or make treatment recommendations.

The evidence chain is:

```text
region -> genomic feature -> gene -> drug-gene relationship -> source record
```

Every result retains the genomic coordinates, controlled relationship category,
coordinate evidence scope, and source provenance.

# Build a relationship table

The canonical relationship categories are `target`, `pharmacodynamics`,
`metabolism`, `transport`, `carrier`, `efficacy`, `toxicity`, `dosage`,
`Other`, `resistance`, and `exposure`.

This local example uses a small user-curated table:

```{r build}
records <- data.frame(
    drug_id = c("drug:clopidogrel", "drug:clopidogrel"),
    drug_name = c("clopidogrel", "clopidogrel"),
    gene_symbol = c("CYP2C19", "P2RY12"),
    relationship_category = c("metabolism", "target"),
    evidence_level = c("1A", NA),
    evidence_scope = c("exact_variant", "gene_association"),
    source = c("example PGx source", "example target source"),
    source_record_id = c("pgx:1", "target:1"),
    variant_id = c("rs-example", NA)
)

drug_sets <- buildDrugGeneSets(
    sources = "custom",
    source_data = list(custom = records)
)
drug_sets
```

ClinPGx/PharmGKB, DGIdb, and DrugCentral public tables are fetched and cached
when no local entry is supplied. DGIdb and DrugCentral claims remain in the
`target` category; target interaction data are not promoted to efficacy or
toxicity:

```{r}
drug_sets <- buildDrugGeneSets(
    sources = c("ClinPGx", "DGIdb", "DrugCentral")
)
```

A local `summaryAnnotations.zip` can be supplied through `source_data` for a
fixed snapshot. ClinPGx annotations with multiple supported phenotype
categories are expanded into distinct relationship rows. ClinPGx phenotype
categories that are not assigned to a supported pharmacological relationship
are retained as `Other`.

# Map relationships to genomic features

A gene model can be a `TxDb` or a `GRanges`. An untyped `GRanges` represents
gene bodies. Here a compact model keeps the example self-contained:

```{r map}
gene_model <- GRanges(
    seqnames = c("chr10", "chr3"),
    ranges = IRanges(
        start = c(94760000, 151337000),
        end = c(94860000, 151384000)
    ),
    strand = c("+", "-"),
    gene_symbol = c("CYP2C19", "P2RY12")
)

variant_ranges <- GRanges(
    seqnames = "chr10",
    ranges = IRanges(94781859, 94781859),
    variant_id = "rs-example"
)

drug_ranges <- mapDrugGenes(
    drug_sets,
    genome = "GRCh38",
    gene_model = gene_model,
    features = c("gene_body", "promoter"),
    variant_ranges = variant_ranges
)
drug_ranges
```

Promoter mappings receive `evidence_scope = "regulatory_overlap"`. Gene bodies,
coding exons, UTRs, and introns receive `gene_association`. Only the supplied
variant coordinate receives `exact_variant`; merely mentioning an rsID in a
source record is not sufficient.

For `TxDb` inputs, the available feature modes are:

```r
c("coding_exon", "utr", "intron", "gene_body", "promoter")
```

Promoters use transcript starts and default to 2,000 bases upstream and 500
bases downstream. Distant nearest-gene assignment is not performed.

# Annotate regions

```{r annotate}
regions <- GRanges(
    seqnames = c("chr10", "chr10", "chr3"),
    ranges = IRanges(
        start = c(94781000, 94781859, 151383500),
        end = c(94782000, 94781859, 151384500)
    ),
    region_id = c("region_1", "region_2", "region_3")
)

annotations <- annotateDrugRegions(regions, drug_ranges)
mcols(annotations)[, c(
    "region_id",
    "gene_symbol",
    "drug_name",
    "relationship_category",
    "feature_type",
    "evidence_scope",
    "overlap_bp",
    "source"
)]
```

BED file paths are imported with `rtracklayer`, so zero-based BED coordinates
are converted correctly. Data-frame inputs are explicitly interpreted as
one-based closed genomic coordinates.

Rows remain feature- and source-record-specific. Disjoint transcript pieces
within one logical feature mapping are collapsed, and their overlap is unioned
so duplicated isoforms do not inflate `overlap_bp`.

# Summarise without scoring drugs

```{r summarise}
summariseDrugAnnotations(
    annotations,
    by = c("region", "drug", "relationship")
)
```

Summaries include applicable counts, union overlap bases, represented feature
classes, evidence levels, and source provenance. Columns for dimensions already
used in `by` are omitted.

# Session information

```{r session-info}
sessionInfo()
```