metaplot

Visualize m6A (or any other RNA modification) peaks across 5'UTR, CDS and 3'UTR regions, and inspect where a motif of interest sits relative to the peak centre.

You can install the development version of metaplot like so:

# install.packages("devtools")
devtools::install_github("junjunlab/metaplot")

# or
remotes::install_github("junjunlab/metaplot")

Everything below is shipped in inst/extdata, so the examples in this README can be copy-pasted and run as-is right after installing the package.

Three calls are enough for the standard meta-gene plot:

library(metaplot)

gtf_file <- system.file("extdata", "hg38.ncbiRefSeq.gtf.gz", package = "metaplot")
m6a_bed  <- system.file("extdata", "m6a_peaks.bed", package = "metaplot")

# 1. annotation, computed once per genome
features_anno <- prepare_features(gtf_file, features_name = c("5UTR", "CDS", "3UTR"))

# 2. + 3. map peaks and plot
res <- metaPlot(
  bed_file      = list(m6a_bed),
  group_names   = "m6A",
  features_anno = features_anno,
  scale_region  = TRUE
)

res$plot

The rest of this document walks through every function and every argument. The examples reuse features_anno from above and these additional paths:

m1a_bed    <- system.file("extdata", "m1a.sorted2.bed", package = "metaplot")
helako_bed <- system.file("extdata", "helaA30ko-1.m1a.sorted2.bed", package = "metaplot")
m6a_fa     <- system.file("extdata", "m6a_peaks.fa", package = "metaplot")

prepare_features() reads a GTF, keeps one transcript per gene (the one with the largest summed exon length), and returns the 5'UTR / CDS / 3'UTR segments of those transcripts as a GRanges.

features_anno <- prepare_features(
  gtf_file      = gtf_file,
  features_name = c("5UTR", "CDS", "3UTR")   # UCSC/RefSeq labels
)

features_name must be a character vector of length 3, given in the order 5'UTR, CDS, 3'UTR. Different annotation sources use different labels:

What you get back (abbreviated):

features_anno
#> GRanges object with 233837 ranges and 7 metadata columns:
#>          tx_len        type     f_len      gene_id    gene_name  transcript_id
#> 1           956         CDS       956          ND1          ND1        rna-ND1
#> 2          1042         CDS      1042          ND2          ND2        rna-ND2
#> ...         ...         ...       ...          ...          ...            ...
#> 233837      421        UTR3       421 LOC112268260 LOC112268260 XM_047436352.1
#>         trans_len
#> 1             956
#> ...          ...
#> 233837      1819

Each row is one exon-level segment; f_len is the summed length of that feature type across the whole transcript, and tx_len is the cumulative length up to and including this segment (counted in transcript orientation, so the minus strand is handled correctly).

as.data.frame(features_anno) |> dplyr::count(type)
#>   type      n
#> 1  CDS 181031
#> 2 UTR3  19790
#> 3 UTR5  33016

19,539 transcripts remain after the longest-transcript filter. The CDS is the largest region (median 2.1 kb per transcript), while the median 5'UTR is only 0.38 kb — which is exactly why scale_region exists (see section 3.2).

features_anno only needs to be computed once per genome and can be reused across every metaPlot() call, so it is worth saving:

saveRDS(features_anno, "hg38_features_anno.rds")

metaPlot() is a thin plotting layer on top of this function, which is exported so you can inspect or reuse the numbers behind the figure.

rel_df <- calculate_relative_position(
  bed_file      = m6a_bed,
  features_anno = features_anno,
  scale_region  = TRUE
)

head(rel_df)
#>    type  rel_pos
#> 1   CDS 1.129853
#> 2  UTR3 2.011796
#> 3  UTR3 2.071849
#> 4   CDS 1.284665
#> ...

Every feature type is normalised to a length of 1 and stacked along the x-axis:

     5'UTR            CDS            3'UTR
  |-----------|----------------|-----------|
  0           1                2           3

• rel_pos in [0, 1] → 5'UTR
• rel_pos in [1, 2] → CDS
• rel_pos in [2, 3] → 3'UTR

A peak is counted only when its midpoint falls inside a feature, which is why 20,435 input peaks produce 11,833 rows here — the rest fall in introns or intergenic space.

With scale_region = TRUE the 5'UTR and 3'UTR are not drawn at full width but compressed/expanded by their median length relative to the CDS, so the panel widths reflect real transcript proportions.

res <- metaPlot(
  bed_file      = list(m6a_bed),
  group_names   = "m6A",
  features_anno = features_anno,
  scale_region  = TRUE
)

res$plot

bed_file accepts a list of BED paths, or a single BED path / data.frame / GRanges object without wrapping it in a list.

res_equal  <- metaPlot(list(m6a_bed), "m6A", features_anno, scale_region = FALSE)
res_scaled <- metaPlot(list(m6a_bed), "m6A", features_anno, scale_region = TRUE)

With scale_region = FALSE (A) all three regions occupy the same width, so a 5'UTR that is 0.38 kb and a CDS that is 2.1 kb look identical. With scale_region = TRUE (B) the region widths follow the median lengths observed in the annotation, and the tick labels are placed at the centre of each region.

The scaling factors are estimated once from the whole annotation, not from the subset of features that happens to be covered by your peaks, so every sample in a multi-sample plot shares exactly the same x-axis.

If you would rather fix the proportions yourself, pass a ratio of length 3 (5'UTR : CDS : 3'UTR) instead of relying on the estimated median lengths:

res_cut <- metaPlot(
  bed_file      = list(m6a_bed),
  group_names   = "m6A",
  features_anno = features_anno,
  scale_region  = TRUE,
  cut_ratio     = c(0.1, 0.8, 0.1)
)

cut_ratio only takes effect together with scale_region = TRUE; if you set it while scale_region = FALSE, metaPlot() warns and ignores it.

res <- metaPlot(
  bed_file      = list(m6a_bed, m1a_bed, helako_bed),
  group_names   = c("m6A", "m1A", "m1A_helaA30KO"),
  features_anno = features_anno,
  scale_region  = TRUE
)

res$plot

The legend title comes from the internal column name; rename it on the returned ggplot object if you prefer something cleaner:

res$plot + ggplot2::labs(fill = "Group", colour = "Group")

Overlaid densities get hard to read beyond two or three groups; faceting gives each sample its own panel.

res_facet <- metaPlot(
  bed_file      = list(m6a_bed, m1a_bed, helako_bed),
  group_names   = c("m6A", "m1A", "m1A_helaA30KO"),
  features_anno = features_anno,
  scale_region  = TRUE,
  facet         = TRUE
)

res_facet$plot

Facets are ordered alphabetically by group_names. Convert the column to a factor first if you need the input order instead.

# thinner fill, thicker outline
metaPlot(list(m6a_bed, m1a_bed), c("m6A", "m1A"), features_anno,
         scale_region = TRUE, col_alpha = 0.15, line_width = 1.6)

# bandwidth controls how much detail survives
metaPlot(list(m6a_bed, m1a_bed), c("m6A", "m1A"), features_anno,
         scale_region = TRUE, geom_density_bw = 0.15)

For a single sample, line_color replaces the default "#CC0033":

metaPlot(list(m6a_bed), "m6A", features_anno,
         scale_region = TRUE, line_color = "#2C7FB8", col_alpha = 0.3)

metaPlot() returns a list with two elements, and the plot is a plain ggplot object, so anything ggplot2 offers keeps working:

res$data   # data frame with type / rel_pos / group_names
res$plot   # ggplot object

res$plot +
  ggplot2::labs(title = "m6A metaplot", y = "Peak density") +
  ggplot2::theme(legend.position = "bottom")

This part answers a different question: where inside the peak does the motif sit? Four helpers wrap small Python scripts through reticulate.

Requires a genome FASTA (not bundled — point it at your own reference):

getSeqFromBed(
  peak      = m6a_bed,
  genomefie = "hg38.fa",     # note the historical spelling of this argument
  outfasta  = "peaks.fa",
  type      = "bed",         # or "narrowpeak"
  a         = 25,            # extend 25 bp upstream of the peak
  b         = 25,            # extend 25 bp downstream of the peak
  minlen    = 5              # drop sequences shorter than this
)

The return value is the number of sequences written to outfasta, and every record is named peak_name|chrom|start|end|strand:

>1|chr19|58357367|58357517|-
...extracted sequence...

a and b are interpreted in transcript orientation: on the minus strand "upstream" means towards higher coordinates, and the extracted sequence is reverse-complemented. Both "bed" and "narrowpeak" are 0-based half-open, so neither needs a coordinate adjustment.

motif_hits <- findMotif(peaksfa = m6a_fa)   # m6a_fa ships with the package
length(motif_hits)
#> [1] 12220
head(sort(table(toupper(motif_hits)), decreasing = TRUE), 8)
#>
#> GGACT GAACT GGACA AGACT AAACT GAACA AAACA AGACA
#>  1114   959   853   850   822   695   680   646

The default motif is "[GA][GA]AC[ACT]|[TGA]GT[CT][CT]", i.e. "RRACH" plus its reverse complement "DGGYY" (both strands of the same consensus), and matching is case-insensitive. Pass any regular expression to search for something else:

findMotif(peaksfa = m6a_fa, motif = "GGAC")

rel_dist <- CalculatePeaksRelativeDist(peaksfa = m6a_fa)
length(rel_dist)
#> [1] 12220
range(rel_dist)
#> [1] -242  244

Each value is the distance, in bp, between a motif hit and the centre of the sequence it came from. mid = 3 (the default) shifts the hit from the start of the match to the base of interest inside the motif.

plotRelMotif(relpos = list(rel_dist), exp = "m6A")

The RRACH motif concentrates sharply at the peak centre — the expected result for m6A. relpos accepts a list, so several samples can be compared on one axis:

# rel_dist_m6a / rel_dist_m1a: one vector per sample, as returned above
plotRelMotif(
  relpos    = list(rel_dist_m6a, rel_dist_m1a),
  exp       = c("m6A", "m1A"),
  linewidth = 0.8
)

The function returns a ggplot object, so the legend can be dropped for a single sample with + ggplot2::theme(legend.position = "none").

getSeqFromBed() needs pyfaidx; findMotif() and CalculatePeaksRelativeDist() only use the standard library. If reticulate cannot find an interpreter, point it at one explicitly:

reticulate::install_miniconda()          # or use an existing Python
getSeqFromBed(..., pythonPath = "/path/to/python")

No peak falls within the annotated features — the peaks and the annotation do not agree on chromosome naming. chr1 vs 1 is the usual culprit; convert one side with GenomeInfoDb::seqlevelsStyle() before plotting.

No feature matched features_name — the GTF uses different type labels. Inspect them with table(rtracklayer::import(gtf_file)$type) and pass the right ones to features_name.

features_anno contains no CDS feature — scale_region = TRUE needs a CDS to scale against. Either pass cut_ratio explicitly or set scale_region = FALSE.

group_names has length ... but bed_file has ... sample(s) — the two must have the same length; group_names is matched to bed_file by position.

Please install python first! — reticulate found no interpreter. See Python setup.

Jun Zhang (2024). Visualizing m6A Distribution on RNA. https://github.com/junjunlab/metaplot

写个包 metaplot 绘制 m6A peak 分布

All figures in this README are regenerated by data-raw/make-readme-figures.R.