Use with Bioconductor
Compute overlap with plyranges
We can query genomic elements using elements() and then compute overlaps with the plyranges package as below. See the plyranges package documentation for more examples.
library(plyranges)
rng <- data.frame(seqnames="chr1", start=10e6+1, end=10.1e6) |>
as_granges()
e <- elements(rng, limit=200L) |>
plyranges::filter(source != "ENCODE_EpiRaction")
e
#> GRanges object with 200 ranges and 5 metadata columns:
#> seqnames ranges strand | name
#> <Rle> <IRanges> <Rle> | <character>
#> [1] chr1 10001018-10001351 * | EH38E1317660
#> [2] chr1 10001823-10002161 * | EH38E3954141
#> [3] chr1 10002812-10003022 * | EH38E3954142
#> [4] chr1 10006462-10006688 * | EH38E3954143
#> [5] chr1 10007231-10007550 * | EH38E2785201
#> ... ... ... ... . ...
#> [196] chr1 10000983-10001482 * | genic_chr1_10000982_..
#> [197] chr1 10000983-10001482 * | genic_chr1_10000982_..
#> [198] chr1 10000983-10001482 * | genic_chr1_10000982_..
#> [199] chr1 10000983-10001482 * | genic_chr1_10000982_..
#> [200] chr1 10000983-10001482 * | genic_chr1_10000982_..
#> source_annotation type source
#> <character> <character> <character>
#> [1] dELS: distal Enhance.. candidate cis regula.. ENCODE
#> [2] TF: TF binding candidate cis regula.. ENCODE
#> [3] dELS: distal Enhance.. candidate cis regula.. ENCODE
#> [4] CA: chromatin access.. candidate cis regula.. ENCODE
#> [5] dELS: distal Enhance.. candidate cis regula.. ENCODE
#> ... ... ... ...
#> [196] genic accessible dna eleme.. ENCODE
#> [197] genic accessible dna eleme.. ENCODE
#> [198] genic accessible dna eleme.. ENCODE
#> [199] genic accessible dna eleme.. ENCODE
#> [200] genic accessible dna eleme.. ENCODE
#> source_url
#> <character>
#> [1] https://www.encodepr..
#> [2] https://www.encodepr..
#> [3] https://www.encodepr..
#> [4] https://www.encodepr..
#> [5] https://www.encodepr..
#> ... ...
#> [196] https://www.encodepr..
#> [197] https://www.encodepr..
#> [198] https://www.encodepr..
#> [199] https://www.encodepr..
#> [200] https://www.encodepr..
#> -------
#> seqinfo: 1 sequence from hg38 genome; no seqlengthstiles <- tile_ranges(rng, width=10e3) %>%
plyranges::select(-partition) %>%
plyranges::mutate(id = letters[seq_along(.)])
# count overlaps of central basepair of elements in tiles
e |>
plyranges::anchor_center() |>
plyranges::mutate(width = 1) |>
plyranges::join_overlap_left(tiles) |>
tibble::as_tibble() |>
dplyr::count(id)
#> # A tibble: 10 × 2
#> id n
#> <chr> <int>
#> 1 a 100
#> 2 b 13
#> 3 c 10
#> 4 d 15
#> 5 e 11
#> 6 f 15
#> 7 g 12
#> 8 h 12
#> 9 i 7
#> 10 j 5Plot variants with plotgardener
Below we show a simple example of plotting variants and elements in a gene context, using the plotgardener package.
First selecting some variants around the gene GCK (reminder that we only obtain the first limit number of variants, see ?catalog_queries for more details).
# up to 200 variants:
v <- gene_variants(gene_name = "GCK", limit=200L, verbose=TRUE) |>
dplyr::select(-c(gene, chr, source, source_url)) |>
tidyr::unnest_wider(sequence_variant) |>
dplyr::rename(seqnames = chr) |>
dplyr::mutate(start = pos + 1, end = pos + 1) |>
as_granges()We then load plotgardener and define some default parameters.
library(plotgardener)
par <- pgParams(
chrom = "chr7",
chromstart = 44.1e6,
chromend = 44.25e6,
assembly = "hg38",
just = c("left", "bottom")
)Renaming some columns:
To match with the correct gene annotation, we could explicitly define the transcript database using plotgardener::assembly(), where we would provide a database built by running GenomicFeatures::makeTxDbFromGFF() on a GENCODE GTF file. Here we use one of the standard hg38 TxDb with UCSC-style chromosome names.
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
#> Loading required package: GenomicFeatures
#> Loading required package: AnnotationDbi
#> Loading required package: Biobase
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#>
#> Attaching package: 'AnnotationDbi'
#> The following object is masked from 'package:plyranges':
#>
#> select
library(org.Hs.eg.db)
#> Below we build the page and populate it with annotation from the TxDb used by plotgardener and from IGVF Catalog.
pageCreate(width = 5, height = 4, showGuides = FALSE)
plotGenes(
params = par, x = 0.5, y = 3.5, width = 4, height = 1
)
#> genes[genes1]
plotGenomeLabel(
params = par,
x = 0.5, y = 3.5, length = 4,
just = c("left", "top")
)
#> genomeLabel[genomeLabel1]
mplot <- plotManhattan(
params = par, x = 0.5, y = 2.5, width = 4, height = 2,
v_for_plot, trans = "",
sigVal = -log10(5e-8), sigLine = TRUE, col = "grey", lty = 2
)
#> manhattan[manhattan1]
annoYaxis(
plot = mplot, at=0:4 * 4, axisLine = TRUE, fontsize = 8
)
#> yaxis[yaxis1]
annoXaxis(
plot = mplot, axisLine = TRUE, label = FALSE
)
#> xaxis[xaxis1]
plotText(
params = par,
label = "-log10(p-value)", x = 0.2, y = 2, rot = 90,
fontsize = 8, fontface = "bold",
default.units = "inches"
)
#> text[text1]