.Rhistory

head(df)
txLoc
# Sanity checks
CheckClass(txLoc, "txLoc")
regions <- c("5'UTR", "CDS", "3'UTR")
objName <- deparse(substitute(txLoc))
if (!identical(GetRegions(txLoc), regions)) stop(sprintf(
"%s must contain the regions %s!",
objName,
paste(regions, collapse = ", ")))
# Get transcript region start coordinates along transcript and
# coordinates of start/stop codons
lst <- GetTxRegionCoordinates(GetRef(txLoc))
# Get transcript region coordinates of loci
df.loci <- do.call(
rbind, c(Map(
function(locus, reg)
setNames(cbind(
as.data.frame(locus$locus_in_tx_region)[, 1:3],
reg),
c("seqnames", "start", "end", "region")),
GetLoci(txLoc),
regions),
make.row.names = FALSE)
)
# Merge transcript region coordinates and widths
df.loci <- merge(
df.loci,
setNames(
cbind(
lst[["start_region"]][, "seqnames"],
stack(lst[["start_region"]], select = -seqnames)),
c("seqnames", "start_region", "region")),
by = c("seqnames", "region"))
# Convert transcript region to transcript coordinates
df.loci[, "start"] <- df.loci[, "start"] + df.loci[, "start_region"]
df.loci[, "end"] <- df.loci[, "end"] + df.loci[, "start_region"]
df.loci
# Store loci as `GRanges`
gr.loci <- GRanges(
seqnames = df.loci[, "seqnames"],
IRanges(df.loci[, "start"], df.loci[, "end"]),
strand = "*")
# Calculate distances to nearest start/stop codons
GetRelDistNearest(
list(start = gr.loci, stop = gr.loci),
list(start = lst[["gr.start"]], stop = lst[["gr.stop"]]))
source('~/.active-rstudio-document', echo=TRUE)
Sys.Date()
format(Sys.Date(), "%d-%m-%Y)"
format(Sys.Date(), "%d-%m-%Y")
format(Sys.Date(), "%d-%m-%Y")
format(Sys.Date(), "%d %M %Y")
?strptime
?strptime
format(Sys.Date(), "%d %b %Y")
m1A
GetLoci(m1A)
gr1 <- GetLoci(m1A)$locus_in_tx_region
gr1
GetLoci(m1A)
gr1 <- GetLoci(m1A)[[1]]$locus_in_tx_region
gr1
countOverlaps(gr1, gr1)
countOverlaps(gr1, sample(gr1, 10))
countOverlaps(gr1, sample(gr1, 100))
countOverlaps(gr1, sample(gr1, 1000))
PAS <- read.table("~/Desktop/Ulrike_m5C_PAS/gencode.v29.polyAs.gtf", skip = 2) %>%
filter(V3 == "polyA_site") %>%
transmute(
chr = V1,
start = V4 + 1,
end = V5,
strand = V7,
score = 0,
id = V3) %>%
makeGRangesFromDataFrame(keep.extra.columns = T) %>%
SmartMap(id = "PAS", refGenome = "hg38", showPb = FALSE)
PAS
motif <- GetMotifLoc(refGenome = "hg38", id = "PAS", showPb = FALSE)
motif
library(gplots)
txLoc1 <- PAS
txLoc2 <- motif
# Sanity check
CheckClassTxLocConsistency(txLoc1, txLoc2)
gr1 <- TxLoc2GRangesList(txLoc1)
#' @param method A character string; specifies whether coordinates
#' are genome (\code{method = "genome"}) or transcriptome coordinates
#' (\code{method = "tx_region"}).
#'
#' @return A \code{list} of \code{GRanges} objects. See 'Details'.
#'
#' @author Maurits Evers, \email{maurits.evers@@anu.edu.au}
#' @keywords internal
#'
#' @importFrom S4Vectors DataFrame
TxLoc2GRangesList <- function(txLoc, method = c("tx_region", "genome")) {
# Sanity checks
CheckClass(txLoc, "txLoc")
method <- match.arg(method)
# Return a `list` of `GRanges`
lapply(GetLoci(txLoc), function(locus) {
gr <- switch(
method,
"genome" = locus$locus_in_genome,
"tx_region" = locus$locus_in_tx_region)
mcols(gr) <- DataFrame(
source = GetId(txLoc),
tx_refseq = locus$tx_refseq,
tx_region = locus$tx_region,
id = locus$id)
gr
})
}
# Determine figure panel layout
if (length(gr1) < 4) {
par(mfrow = c(1, length(gr1)))
} else {
par(mfrow = c(ceiling(length(gr1) / 2), 2))
}
# Determine figure panel layout
if (length(GetRegions(txLoc1)) < 4) {
par(mfrow = c(1, length(GetRegions(txLoc1))))
} else {
par(mfrow = c(ceiling(length(GetRegions(txLoc1)) / 2), 2))
}
invisible(Map(
function(loci1, loci2, region) {
# Get transcriptome coordinates
gr1 <- loci1$locus_in_tx_region
gr2 <- loci2$locus_in_tx_region
# Make sure that seqlevels match
lvls <- intersect(seqlevels(gr1), seqlevels(gr2))
seqlevels(gr1, pruning.mode = "coarse") <- lvls
seqlevels(gr2, pruning.mode = "coarse") <- lvls
# Count overlap
m <- countOverlaps(gr1, gr2)
overlap <- sum(m > 0)
# Plot
grps <- list(
seq(1, length(gr1)),
seq(length(gr1) - overlap + 1, length.out = length(gr2)))
names(grps) <- c(
sprintf(
"%s (%3.2f%%)",
GetId(txLoc1), overlap / length(gr1) * 100),
sprintf(
"%s (%3.2f%%)",
GetId(txLoc2), overlap / length(gr2) * 100))
venn(grps)
mtext(names(gr1))
},
GetLoci(loci1),
GetLoci(loci2),
GetRegions(loci1)))
invisible(Map(
function(loci1, loci2, region) {
# Get transcriptome coordinates
gr1 <- loci1$locus_in_tx_region
gr2 <- loci2$locus_in_tx_region
# Make sure that seqlevels match
lvls <- intersect(seqlevels(gr1), seqlevels(gr2))
seqlevels(gr1, pruning.mode = "coarse") <- lvls
seqlevels(gr2, pruning.mode = "coarse") <- lvls
# Count overlap
m <- countOverlaps(gr1, gr2)
overlap <- sum(m > 0)
# Plot
grps <- list(
seq(1, length(gr1)),
seq(length(gr1) - overlap + 1, length.out = length(gr2)))
names(grps) <- c(
sprintf(
"%s (%3.2f%%)",
GetId(txLoc1), overlap / length(gr1) * 100),
sprintf(
"%s (%3.2f%%)",
GetId(txLoc2), overlap / length(gr2) * 100))
venn(grps)
mtext(names(gr1))
},
GetLoci(txLoc1),
GetLoci(txLoc2),
GetRegions(txLoc1)))
txLoc1
txLoc2
?venn
invisible(Map(
function(loci1, loci2, region) {
# Get transcriptome coordinates
gr1 <- loci1$locus_in_tx_region
gr2 <- loci2$locus_in_tx_region
# Make sure that seqlevels match
lvls <- intersect(seqlevels(gr1), seqlevels(gr2))
seqlevels(gr1, pruning.mode = "coarse") <- lvls
seqlevels(gr2, pruning.mode = "coarse") <- lvls
# Count overlap
m <- countOverlaps(gr1, gr2)
overlap <- sum(m > 0)
# Plot
grps <- list(
seq(1, length(gr1)),
seq(length(gr1) - overlap + 1, length.out = length(gr2)))
names(grps) <- c(
sprintf(
"%s (%3.2f%%)",
GetId(txLoc1), overlap / length(gr1) * 100),
sprintf(
"%s (%3.2f%%)",
GetId(txLoc2), overlap / length(gr2) * 100))
venn(grps, main = region)
mtext(names(gr1))
},
GetLoci(txLoc1),
GetLoci(txLoc2),
GetRegions(txLoc1)))
invisible(Map(
function(loci1, loci2, region) {
# Get transcriptome coordinates
gr1 <- loci1$locus_in_tx_region
gr2 <- loci2$locus_in_tx_region
# Make sure that seqlevels match
lvls <- intersect(seqlevels(gr1), seqlevels(gr2))
seqlevels(gr1, pruning.mode = "coarse") <- lvls
seqlevels(gr2, pruning.mode = "coarse") <- lvls
# Count overlap
m <- countOverlaps(gr1, gr2)
overlap <- sum(m > 0)
# Plot
grps <- list(
seq(1, length(gr1)),
seq(length(gr1) - overlap + 1, length.out = length(gr2)))
names(grps) <- c(
sprintf(
"%s (%3.2f%%)",
GetId(txLoc1), overlap / length(gr1) * 100),
sprintf(
"%s (%3.2f%%)",
GetId(txLoc2), overlap / length(gr2) * 100))
venn(grps)
title(region)
mtext(names(gr1))
},
GetLoci(txLoc1),
GetLoci(txLoc2),
GetRegions(txLoc1)))
##'
##' @param txLoc1 A \code{txLoc} object.
##' @param txLoc2 A \code{txLoc} object.
##'
##' @return \code{NULL}.
#'
##' @author Maurits Evers, \email{maurits.evers@@anu.edu.au}
##'
##' @import GenomicRanges IRanges
##' @importFrom gplots venn
PlotOverlap <- function(txLoc1, txLoc2) {
# Sanity check
CheckClassTxLocConsistency(txLoc1, txLoc2)
# Determine figure panel layout
if (length(GetRegions(txLoc1)) < 4) {
par(mfrow = c(1, length(GetRegions(txLoc1))))
} else {
par(mfrow = c(ceiling(length(GetRegions(txLoc1)) / 2), 2))
}
invisible(Map(
function(loci1, loci2, region) {
# Get transcriptome coordinates
gr1 <- loci1$locus_in_genome
gr2 <- loci2$locus_in_genome
# Make sure that seqlevels match
lvls <- intersect(seqlevels(gr1), seqlevels(gr2))
seqlevels(gr1, pruning.mode = "coarse") <- lvls
seqlevels(gr2, pruning.mode = "coarse") <- lvls
# Count overlap
m <- countOverlaps(gr1, gr2)
overlap <- sum(m > 0)
# Plot
grps <- list(
seq(1, length(gr1)),
seq(length(gr1) - overlap + 1, length.out = length(gr2)))
names(grps) <- c(
sprintf(
"%s (%3.2f%%)",
GetId(txLoc1), overlap / length(gr1) * 100),
sprintf(
"%s (%3.2f%%)",
GetId(txLoc2), overlap / length(gr2) * 100))
venn(grps)
title(region)
mtext(names(gr1))
},
GetLoci(txLoc1),
GetLoci(txLoc2),
GetRegions(txLoc1)))
}
PlotOverlap(PAS, motif)
##'
##' @param txLoc1 A \code{txLoc} object.
##' @param txLoc2 A \code{txLoc} object.
##'
##' @return \code{NULL}.
#'
##' @author Maurits Evers, \email{maurits.evers@@anu.edu.au}
##'
##' @import GenomicRanges IRanges
##' @importFrom gplots venn
PlotOverlap <- function(txLoc1, txLoc2) {
# Sanity check
CheckClassTxLocConsistency(txLoc1, txLoc2)
# Determine figure panel layout
if (length(GetRegions(txLoc1)) < 4) {
par(mfrow = c(1, length(GetRegions(txLoc1))))
} else {
par(mfrow = c(ceiling(length(GetRegions(txLoc1)) / 2), 2))
}
invisible(Map(
function(loci1, loci2, region) {
# Get transcriptome coordinates
gr1 <- loci1$locus_in_tx_region
gr2 <- loci2$locus_in_tx_region
# Make sure that seqlevels match
lvls <- intersect(seqlevels(gr1), seqlevels(gr2))
seqlevels(gr1, pruning.mode = "coarse") <- lvls
seqlevels(gr2, pruning.mode = "coarse") <- lvls
# Count overlap
m <- countOverlaps(gr1, gr2)
overlap <- sum(m > 0)
# Plot
grps <- list(
seq(1, length(gr1)),
seq(length(gr1) - overlap + 1, length.out = length(gr2)))
names(grps) <- c(
sprintf(
"%s (%3.2f%%)",
GetId(txLoc1), overlap / length(gr1) * 100),
sprintf(
"%s (%3.2f%%)",
GetId(txLoc2), overlap / length(gr2) * 100))
venn(grps)
title(region)
mtext(names(gr1))
},
GetLoci(txLoc1),
GetLoci(txLoc2),
GetRegions(txLoc1)))
}
PlotOverlap(m1A, motif)
PlotOverlap(PAS, motif)
1638 / (223821 + 1638)
1638 / (223821 + 1638) *100
WriteBED(PAS, file = "PAS_GENCODE.bed")
WriteBED(motif, file = "PAS_motif.bed")
PlotOverlap(PAS, motif)
?venn
library(RNAModR)
library(magrittr)
m5C <- ReadBED("~/Desktop/Ulrike_m5C_PAS/Ulrike_m5C_sites_hg38.bed") %>%
SmartMap(id = "m5C", refGenome = "hg38") %>%
FilterTxLoc(c("5'UTR", "CDS", "3'UTR"))
null <- GenerateNull(m5C, nt = "C")
ago <- ReadBED(system.file(
"extdata",
"PARCLIP_AGO1234_Hafner2010_hg38.bed",
package = "RNAModR")) %>%
SmartMap(id = "AGO1-4", refGenome = "hg38") %>%
FilterTxLoc(c("5'UTR", "CDS", "3'UTR"))
PlotRelDistEnrichment(m5C, null, ago)
txLoc <- m5C
txLoc1 <- m5C
txLoc2 <- null
txLocRef <- ago
# Sanity checks
CheckClass(txLoc1, "txLoc")
CheckClass(txLoc2, "txLoc")
CheckClass(txLocRef, "txLoc")
CheckClassTxLocConsistency(txLoc1, txLocRef)
CheckClassTxLocConsistency(txLoc2, txLocRef)
# Get ids
id1 <- GetId(txLoc1)
id2 <- GetId(txLoc2)
idRef <- GetId(txLocRef)
id1
id2
idREF
idRef
# Convert sites to `list` of `GRanges`
lst1 <- TxLoc2GRangesList(txLoc1, method = "tx_region")
#' @param method A character string; specifies whether coordinates
#' are genome (\code{method = "genome"}) or transcriptome coordinates
#' (\code{method = "tx_region"}).
#'
#' @return A \code{list} of \code{GRanges} objects. See 'Details'.
#'
#' @author Maurits Evers, \email{maurits.evers@@anu.edu.au}
#' @keywords internal
#'
#' @importFrom S4Vectors DataFrame
TxLoc2GRangesList <- function(txLoc, method = c("tx_region", "genome")) {
# Sanity checks
CheckClass(txLoc, "txLoc")
method <- match.arg(method)
# Return a `list` of `GRanges`
lapply(GetLoci(txLoc), function(locus) {
gr <- switch(
method,
"genome" = locus$locus_in_genome,
"tx_region" = locus$locus_in_tx_region)
mcols(gr) <- DataFrame(
source = GetId(txLoc),
tx_refseq = locus$tx_refseq,
tx_region = locus$tx_region,
id = locus$id)
gr
})
}
# Convert sites to `list` of `GRanges`
lst1 <- TxLoc2GRangesList(txLoc1, method = "tx_region")
lst2 <- TxLoc2GRangesList(txLoc2, method = "tx_region")
lstRef <- TxLoc2GRangesList(txLocRef, method = "tx_region")
lst1
# Calculate distances
lstDist1 <- GetRelDistNearest(lst1, lstRef)
#' @param lst1 A \code{list} of \code{GRanges} object.
#' @param lst2 A \code{list} of \code{GRanges} object.
#'
#' @return A list of \code{integer} vectors. See 'Details'.
#'
#' @author Maurits Evers, \email{maurits.evers@@anu.edu.au}
#' @keywords internal
#'
#' @import GenomicRanges
#' @importFrom S4Vectors queryHits subjectHits
GetRelDistNearest <- function(lst1, lst2) {
# Sanity checks
all(sapply(lst1, CheckClass, "GRanges"))
all(sapply(lst2, CheckClass, "GRanges"))
stopifnot(identical(names(lst1), names(lst2)))
# Calculate nearest distances from start
lst <- Map(
function(gr1, gr2) {
# Make sure that seqlevels match
lvls <- intersect(seqlevels(gr1), seqlevels(gr2))
seqlevels(gr1, pruning.mode = "coarse") <- lvls
seqlevels(gr2, pruning.mode = "coarse") <- lvls
# Collapse range of gr1 and gr2 to the start coordinate
end(gr1) <- start(gr1)
end(gr2) <- start(gr2)
# Calculate distance to nearest
hits <- distanceToNearest(gr1, gr2, ignore.strand = TRUE)
# Define and return distance d as
#   d > 0 : if pos(gr1) > pos(gr2)
#   d < 0 : if pos(gr1) < pos(gr2)
# In words: Negative distances => gr1 is upstream of gr2
#           Positive distances => gr1 is downstream of gr2
ifelse(
end(gr1)[queryHits(hits)] > start(gr2)[subjectHits(hits)],
mcols(hits)$distance,
-mcols(hits)$distance)
},
lst1, lst2)
# Return lst
lst
}
# Calculate distances
lstDist1 <- GetRelDistNearest(lst1, lstRef)
lstDist2 <- GetRelDistNearest(lst2, lstRef)
lstDist1
lstDist2
# Set breaks & binwidth
breaks <- seq(-flank, flank, by = binWidth)
flank <- 1000
binWidt <- 20
binWidth <- 20
# Set breaks & binwidth
breaks <- seq(-flank, flank, by = binWidth)
bwString <- sprintf("bw = %3.2f", binWidth)
invisible(Map(
function(dist1, dist2, region) {
# Filter distances that are within window [-flank, flank]
if (flank > 0) {
dist1 <- dist1[abs(dist1) <= flank]
dist2 <- dist2[abs(dist2) <= flank]
}
# Bin distances and count matrix
cts1 <- table(cut(dist1, breaks = breaks))
cts2 <- table(cut(dist2, breaks = breaks))
ctsMat <- as.matrix(rbind(cts1, cts2))
rownames(ctsMat) <- c("pos", "neg")
# Plot
title <- sprintf(
"%s\nN(d(%s,%s)) = %i, N(d(%s,%s)) = %i\nbw = %i nt, flank = %i nt",
region,
id1, idRef, sum(cts1),
id2, idRef, sum(cts2),
binWidth,
flank)
PlotEnrichment.Generic(
ctsMat,
title = title,
xlab = sprintf("Distance relative to %s [nt]", idRef),
x.las = 1, x.cex = 0.8, x.padj = 0,
xAxisLblFmt = 3)
},
lstDist1, lstDist2, names(lstDist1)))
dist1 <- lstDist1[[1]]
dist1
length(lstDist1)
names(lstDist1)
lstDist1
lstDist2