annotatERs: Connects ERs to genes using junction data, then classifies...
In eolagbaju/ODER: Optimising the Definition of Expressed Regions

Description Usage Arguments Value Examples

Finds the overlap between junctions and ERs, then adds gene info and junction info as metadata columns. Then, uses a gtf file or a Txdb passed in to generate a genomic state used to label each ER as to whether they are exonic, intronic, intergenic or none.

1	annotatERs(opt_ers, junc_data, genom_state, gtf, txdb)

`opt_ers`	optimally defined ERs (the product of the ODER function)
`junc_data`	junction data that should match the ERs passed into opt_ers
`genom_state`	a genomic state object
`gtf`	gtf in a GRanges object, pre-imported using `rtracklayer::import` . This is used to provide the gene information for annotation.
`txdb`	TxDb-class (txdb object) to create genomic state. This is used to annotate the expressed regions as exonic, intronic or intergenic.

annotated ERs

gtf_url <- paste0(
    "http://ftp.ensembl.org/pub/release-103/gtf/",
    "homo_sapiens/Homo_sapiens.GRCh38.103.chr.gtf.gz"
)
# file_cache is an internal function to download a bigwig file from a link
# if the file has been downloaded recently, it will be retrieved from a cache
gtf_path <- file_cache(gtf_url)

gtf_gr <- rtracklayer::import(gtf_path)

ex_opt_ers <- GenomicRanges::GRanges(
    seqnames = S4Vectors::Rle(c("chr21"), c(5)),
    ranges = IRanges::IRanges(
        start = c(5032176, 5033408, 5034717, 5035188, 5036577),
        end = c(5032217, 5033425, 5034756, 5035189, 5036581)
    )
)

junctions <- SummarizedExperiment::rowRanges(dasper::junctions_example)
chrs_to_keep <- c("21", "22")

#### preparing the txdb and genomstate object(s)

hg38_chrominfo <- GenomeInfoDb::getChromInfoFromUCSC("hg38")
new_info <- hg38_chrominfo$size[match(
    chrs_to_keep,
    GenomeInfoDb::mapSeqlevels(hg38_chrominfo$chrom, "Ensembl")
)]
names(new_info) <- chrs_to_keep
gtf_gr_tx <- GenomeInfoDb::keepSeqlevels(gtf_gr,
    chrs_to_keep,
    pruning.mode = "tidy"
)
GenomeInfoDb::seqlengths(gtf_gr_tx) <- new_info
GenomeInfoDb::seqlevelsStyle(gtf_gr_tx) <- "UCSC"
rtracklayer::genome(gtf_gr_tx) <- "hg38"

ucsc_txdb <- GenomicFeatures::makeTxDbFromGRanges(gtf_gr_tx)
genom_state <- derfinder::makeGenomicState(txdb = ucsc_txdb)
ens_txdb <- ucsc_txdb
GenomeInfoDb::seqlevelsStyle(ens_txdb) <- "Ensembl"

annot_ers1 <- annotatERs(
    opt_ers = ex_opt_ers, junc_data = junctions,
    gtf = gtf_gr, txdb = ens_txdb, genom_state = genom_state
)

annot_ers1