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 106 ranges and 5 metadata columns:
#> seqnames ranges strand | name source_annotation
#> <Rle> <IRanges> <Rle> | <character> <character>
#> [1] chr1 10001018-10001351 * | EH38E1317660 dELS: distal Enhance..
#> [2] chr1 10013307-10013640 * | EH38E1317663 dELS: distal Enhance..
#> [3] chr1 10014094-10014261 * | EH38E1317664 dELS: distal Enhance..
#> [4] chr1 10018931-10019279 * | EH38E1317668 dELS: distal Enhance..
#> [5] chr1 10019632-10019888 * | EH38E1317669 CA-CTCF: chromatin a..
#> ... ... ... ... . ... ...
#> [102] chr1 10045614-10045768 * | EH38E3954155 dELS: distal Enhance..
#> [103] chr1 10064932-10065090 * | EH38E3954156 dELS: distal Enhance..
#> [104] chr1 10080547-10080876 * | EH38E3954157 TF: TF binding
#> [105] chr1 10088561-10088910 * | EH38E3954158 TF: TF binding
#> [106] chr1 10099598-10099944 * | EH38E3954159 CA: chromatin access..
#> type source source_url
#> <character> <character> <character>
#> [1] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [2] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [3] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [4] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [5] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> ... ... ... ...
#> [102] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [103] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [104] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [105] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> [106] candidate cis regula.. ENCODE_SCREEN (ccREs) https://data.igvf.or..
#> -------
#> 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 6
#> 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]