Nothing
.makeObj.feat <- function(isPointer) {
feat <- list()
if (isPointer) {
attr(feat, "class") <- c("feat", "list")
} else attr(feat, "class") <- c("feat", "data.frame")
feat
}
##' Creates a copy of a features object
##'
##' If x is stored in R (as it is by default), then this is no different
##' than x2 <- x. But if it is stored as a pointer to a structure in C,
##' then this is the only way to make an explicity copy of the features.
##' @title Features copy
##' @param x an object of type \code{feat}
##' @return a features object which can be modified independently from the
##' original object
##' @author Melissa J. Hubisz and Adam Siepel
##' @export
copy.feat <- function(x) {
if (is.null(x$externalPtr)) return(x)
result <- .makeObj.feat(TRUE)
result$externalPtr <- .Call.rphast("rph_gff_copy", x$externalPtr)
result
}
##' Read a features object from a file
##'
##' The function will guess the format of the input file automatically.
##'
##' @title Read a Feature File (GFF, BED, or GenePred)
##' @param filename the name of the file (can be GFF, BED, GenePred, or wig:
##' rphast will auto-detect)
##' @param pointer.only Whether to store object by reference instead of a
##' data.frame
##' @return If \code{pointer.only==FALSE}, a data.frame with columns corresponding
##' to the GFF specification. Otherwise, an object which is a pointer to
##' an object stored in C.
##' @seealso \code{\link{feat}} for more description of features objects.
##'
##' \code{\link{msa}} for more explanation of the pointer.only option.
##'
##' \url{http://www.sanger.ac.uk/resources/software/gff/spec.html}
##' for a detailed description of GFF file format. The columns in features
##' objects mirror the GFF column definitions.
##'
##' \url{http://genome.ucsc.edu/FAQ/FAQformat} for descriptions
##' of BED and GenePred formats.
##' @keywords GFF
##' @keywords Genepred
##' @keywords BED
##' @keywords features
##' @author Melissa J. Hubisz and Adam Siepel
##' @example inst/examples/read-feat.R
##' @export
read.feat <- function(filename, pointer.only=FALSE) {
filename <- check.arg(filename, "filename", "character", null.OK=FALSE,
min.length=1L, max.length=1L)
feat <- .makeObj.feat(TRUE)
feat$externalPtr <- .Call.rphast("rph_gff_read", filename)
if (!pointer.only) {
feat <- from.pointer.feat(feat)
}
feat
}
##' Create a new features object
##'
##' See \url{http://www.sanger.ac.uk/resources/software/gff/spec.html}
##' for more detailed description of each parameter.
##'
##' All arguments which are provided should be vectors of equal length.
##'
##' If pointer.only==FALSE, the new object is a data frame, with
##' columns mirroring the GFF Specification
##' Otherwise, it is a list containing a single element, which is
##' a pointer to an object stored in C.
##' @title Features Objects
##' @param seqname a character vector containing the name of the sequence. If
##' the features correspond to regions of an alignment, then seqname
##' should be the name of the sequence in the alignment that is used
##' as the frame of reference in the features. To use the entire alignment
##' as a frame of reference, set seqname to "MSA".
##' @param src The source of the feature
##' @param feature The feature type name
##' @param start The start of the feature. Sequence numbering begins at 1.
##' @param end The end of the feature. This is the last coordinate included
##' in the feature.
##' @param score The feature score, or NA if there is no score.
##' @param strand A character string which is either "+", "-", or "." (if
##' strand is not available or relevant).
##' @param frame A 0, 1, or 2, which specifies whether the feature is in frame.
##' @param attribute A feature attribute (character string).
##' @param pointer.only Whether to store object as a pointer to an object
##' in C, rather than as a data.frame in R.
##' @return If pointer.only==FALSE, returns a data.frame whose format
##' mirrors the GFF specification. Otherwise, returns a list with a single
##' object, which is a external pointer to a C structure representing a
##' features object.
##' @seealso \code{\link{read.feat}}
##'
##' \code{\link{msa}} for more details on the pointer.only option.
##' @keywords features
##' @author Melissa J. Hubisz and Adam Siepel
##' @example inst/examples/feat.R
##' @export
feat <- function(seqname="default", src=".", feature=".",
start, end, score=NULL,
strand=NULL, frame=NULL, attribute=NULL,
pointer.only=FALSE) {
check.arg(start, "start", "integer", null.OK=FALSE, min.length=NULL,
max.length=NULL)
len <- length(start)
seqname <- rep(seqname, length.out = len)
check.arg(seqname, "seqname", "character", null.OK=FALSE, min.length=len,
max.length=len)
src <- rep(src, length.out = len)
check.arg(src, "src", "character", null.OK=FALSE, min.length=len,
max.length=len)
feature <- rep(feature, length.out = len)
check.arg(feature, "feature", "character", null.OK=FALSE, min.length=len,
max.length=len)
check.arg(end, "end", "integer", null.OK=FALSE, min.length=len,
max.length=len)
if (!is.null(score)) score <- rep(score, length.out=len)
check.arg(score, "score", "numeric", null.OK=TRUE, min.length=len,
max.length=len)
if (!is.null(strand)) strand <- rep(strand, length.out=len)
check.arg(strand, "strand", "character", null.OK=TRUE, min.length=len,
max.length=len)
if (!is.null(frame)) frame <- rep(frame, length.out = len)
check.arg(frame, "frame", "integer", null.OK=TRUE, min.length=len,
max.length=len)
if (!is.null(attribute)) attribute <- rep(attribute, length.out = len)
check.arg(attribute, "attribute", "character", null.OK=TRUE,
min.length=len, max.length=len)
check.arg(pointer.only, "pointer.only", "logical", null.OK=FALSE)
# NOTE: I'd like to impose this, but in some genepred files elements have length -1
# if (sum(end - start < 0) > 0)
# stop("Error: end should be >= start for all features")
if (pointer.only) {
if (!is.null(score)) score <- as.numeric(score)
if (!is.null(strand)) strand <- as.character(strand)
if (!is.null(frame)) frame <- as.integer(frame)
if (!is.null(attribute)) attribute <- as.character(attribute)
ptr <- .Call.rphast("rph_gff_new", as.character(seqname),
as.character(src), as.character(feature),
as.integer(start), as.integer(end),
score, strand, frame, attribute)
x <- .makeObj.feat(TRUE)
x$externalPtr <- ptr
} else {
x <- data.frame(seqname=seqname, src=src, feature=feature,
start=start, end=end)
if (!is.null(score)) x <- cbind(x, score)
if (!is.null(strand)) x <- cbind(x, strand)
if (!is.null(frame)) x <- cbind(x, frame)
if (!is.null(attribute)) x <- cbind(x, attribute)
attr(x, "class") <- c("feat", "data.frame")
}
x
}
##' Take a set of features stored in R and return one stored by reference
##' @title Features To Pointer
##' @param x an object of type \code{feat} stored as a data frame in R
##' @return an object of type \code{feat} stored by reference as a pointer to an
##' object created in C.
##' @seealso \code{\link{feat}} for more details on features storage
##' options.
##' @author Melissa J. Hubisz
##' @example inst/examples/as-pointer-feat.R
##' @export
as.pointer.feat <- function(x) {
if (!is.null(x$externalPtr))
return(x)
feat(x$seqname, x$src, x$feature,
x$start, x$end, x$score,
x$strand, x$frame, x$attribute,
pointer.only=TRUE)
}
##' Prints a features object.
##' @title Printing a features Object
##' @param x an object of type \code{feat}
##' @param ... further arguments to be passed to or from other methods
##' @keywords features
##' @seealso \code{\link{write.feat}}
##' @author Melissa J. Hubisz and Adam Siepel
##' @method print feat
##' @export
##' @export print.feat
print.feat <- function(x, ...) {
cat(paste("Features object\n"))
write.feat(x, NULL)
}
##' Write a features object to a file in GFF format.
##' @title Writing a features Object
##' @param x an object of type \code{feat}
##' @param file The name of the file to write to (will be overwritten)
##' @keywords features GFF
##' @author Melissa J. Hubisz and Adam Siepel
##' @example inst/examples/write-feat.R
##' @export
write.feat <- function(x, file) {
check.arg(file, "file", "character", null.OK=TRUE)
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
invisible(.Call.rphast("rph_gff_print", file, x$externalPtr))
}
##' Get the number of rows in a features object
##' @title Number of Features
##' @param x An object of type \code{feat}
##' @return An integer containing the number of rows in each features object
##' @export
##' @export nrow.feat
##' @method nrow feat
##' @keywords features
##' @example inst/examples/nrow-feat.R
##' @author Melissa J. Hubisz and Adam Siepel
nrow.feat <- function(x) {
if (is.data.frame(x)) {
attr(x, "class") <- "data.frame"
nrow(x)
# attr(x, "class") <- "data.frame"
# NextMethod("nrow", x)
# NextMethod()
# UseMethod("nrow", x)
} else {
if (is.null(x$externalPtr)) stop("invalid object sent to nrow.feat")
.Call.rphast("rph_gff_numrow", x$externalPtr)
}
}
##' Get the number of columns in a features object
##' @title Number of Columns in Features
##' @param x An object of type \code{feat}
##' @return An integer containing the number of columns in the features object
##' @note If the features object is stored as a pointer in C, the number
##' of columns is always 9.
##' @export
##' @export ncol.feat
##' @method ncol feat
##' @keywords features
##' @example inst/examples/ncol-feat.R
##' @author Melissa J. Hubisz
ncol.feat <- function(x) {
if (is.data.frame(x)) {
attr(x, "class") <- "data.frame"
# NextMethod("nrow")
ncol(x)
} else return(9) # feat objects stored in C always have 9 columns
}
##' Prints a brief summary of a features object.
##' @title Features Summary
##' @param object an object of type \code{feat}
##' @param ... further arguments to be passed to or from other methods
##' @export
##' @export summary.feat
##' @method summary feat
##' @keywords features
##' @example inst/examples/summary-feat.R
##' @author Melissa J. Hubisz
summary.feat <- function(object, ...) {
if (is.null(object$externalPtr)) {
as <- "stored as data frame"
object <- as.pointer.feat(object)
} else as <- "stored as a pointer to a C structure"
nrow <- nrow.feat(object)
cat(paste("features object with", nrow, "rows", as))
cat("\n")
}
##' Convert a features object to a data frame
##' @title Features to Data Frame
##' @param x an object of type \code{feat}
##' @param row.names optional names for each feature
##' @param optional logical, if \code{TRUE}, setting row names and
##' converting column names (to syntactic names: see
##' \code{\link{make.names}} is optional.
##' @param ... additional arguments to be passed to other methods
##' @return a data frame containing features data
##' @seealso \code{\link{feat}} for a description of features data frames,
##' and \code{\link{as.pointer.feat}} for conversion in the other
##' direction.
##' @export
##' @export as.data.frame.feat
##' @method as.data.frame feat
##' @example inst/examples/as-data-frame-feat.R
##' @author Melissa J. Hubisz and Adam Siepel
as.data.frame.feat <- function(x, row.names=NULL, optional=FALSE, ...) {
if (is.data.frame(x)) return(x)
if (!is.null(x$externalPtr)) {
x <- .Call.rphast("rph_gff_dataframe", x$externalPtr)
}
attr(x, "class") <- "list"
as.data.frame(x, row.names, optional, ...)
}
##' Convert a features object from C memory (external pointer) to R memory
##' @param x A features object stored as a pointer to C memory
##' @return A features object (inheriting from the data.frame class) stored
##' in R memory
##' @export
##' @author Melissa J. Hubisz
from.pointer.feat <- function(x) {
x <- as.data.frame.feat(x)
attr(x, "class") <- c("feat", "data.frame")
x
}
##' Get the dimensions of a features object
##' @title Feature dimensions
##' @param x an object of type \code{feat}
##' @return An integer vector of length two containing the number of
##' rows and number of columns in the features object.
##' @export
##' @export dim.feat
##' @method dim feat
##' @keywords features
##' @author Melissa J. Hubisz
##' @example inst/examples/dim-feat.R
dim.feat <- function(x) {
c(nrow.feat(x), ncol.feat(x))
}
##' Get the range of a features object
##' @title Features range
##' @param ... Objects of type \code{feat}
##' @param na.rm Whether to remove values of NA before calculating range.
##' @return A vector of size 2 indicating minimum and maximum coord in
##' the features object
##' @method range feat
##' @export range.feat
##' @export
##' @keywords features
##' @author Melissa J. Hubisz
range.feat <- function(..., na.rm=FALSE) {
feats <- list(...)
mins <- numeric(length(feats))
maxs <- numeric(length(feats))
for (i in 1:length(feats)) {
x <- feats[[i]]
if (is.data.frame(x)) {
r <- range(c(x$start, x$end), na.rm=na.rm)
} else {
r <- c(.Call.rphast("rph_gff_minCoord", x$externalPtr),
.Call.rphast("rph_gff_maxCoord", x$externalPtr))
}
mins[i] <- r[1]
maxs[i] <- r[2]
}
c(min(mins, na.rm=na.rm),
max(maxs, na.rm=na.rm))
}
##' plot features
##' @title Features plot
##' @param x an object of type \code{feat}
##' @param y the location of the plot on the y axis
##' @param height the height of the boxes
##' @param plottype either "r" for rectangles or "a" for arrows, "b"
##' for arrows within rectangles, or "l" for line segments only.
##' @param arrow.density If plottype=="a" or "b", then this gives the density
##' of arrows in arrows per inch. Otherwise it gives the density of
##' shading lines in the rectangles, and a value of \code{NULL} implies
##' no shading lines.
##' @param angle angle (in degrees) of the shading lines or arrows.
##' @param col color to draw the boxes/lines/arrows with.
##' @param fill.col Color to fill the rectangles with. If \code{NULL}
##' then do not fill.
##' @param lty line type for lines, arrows, borders, and shading
##' @param lwd line width for lines, arrows, borders and shading
##' @param add if \code{TRUE}, add to existing plot
## @param labels whether to label the boxes. If a vector of strings gives
## the lables for each element. If \code{TRUE}, use x$feature for labels.
##' @param xlim A numerical vector of length 2 giving the range for the x-axis.
##' @param ylim A numerical vector of length 2 giving the range for the y-axis.
##' @param ... graphical parameters to be passed to \code{plot}.
##' @method plot feat
##' @export plot.feat
##' @export
##' @keywords features plot
##' @example inst/examples/plot-feat.R
##' @author Melissa J. Hubisz
plot.feat <- function(x, y=0, height=1, plottype="r",
arrow.density=5,
angle=30,
col="black",
fill.col=if (plottype == "r") col else NULL,
lty=par("lty"),
lwd=par("lwd"),
add=FALSE,
# labels=FALSE,
xlim=range.feat(x),
ylim=c(y-height*3/4, y+height*3/4), ...) {
if (!is.null(x$externalPtr))
x <- from.pointer.feat(x)
if (is.null(x$start) || is.null(x$end))
stop("invalid features object")
if (plottype=="a" && is.null(x$strand))
stop("cannot plot feature arrows without strand data")
if (!add)
plot(c(0), c(0), type="n", xlim=xlim, ylim=ylim, xlab="Coordinate", ylab="", ...)
f <- x$end >= xlim[1] & x$start <= xlim[2]
if (sum(f)==0) return(invisible(NULL))
x <- x[f,]
if (plottype=="r" || plottype=="b") {
rect(x$start-0.5, y-height/2, x$end+0.5, y+height/2,
col=fill.col, border=col, lty=lty,
lwd=lwd)
}
if (plottype == "l")
segments(x$start-0.5, y, x$end+0.5, col=col, lty=lty, lwd=lwd)
if (plottype == "a" || plottype == "b") {
usr <- par("usr")
pin <- par("pin")
upi <- c(usr[2L] - usr[1L], usr[4L] - usr[3L])/pin
arrow.density <- upi[1L]/arrow.density #convert from lines per inch
# to coordinates per line
f <- x$strand == "+"
angle <- angle*pi/180 #convert to radians
xlen <- height/2/tan(angle)/upi[2]*upi[1]
if (sum(f) > 0) {
xplus <- x[f,]
for (i in 1:nrow.feat(xplus)) {
start <- xplus[i,]$start-0.5
end <- xplus[i,]$end+0.5
xlen <- height/2/tan(angle)/upi[2]*upi[1]
x0 <- seq(from = end + arrow.density*floor(xlen/arrow.density),
to = start-xlen, by=-arrow.density)
# x0 <- seq(from=end+xlen, to=start-xlen, by=-arrow.density)
x1 <- x0 - xlen
y0top <- rep(y, length(x0))
y0bottom <- rep(y, length(x0))
y1top <- rep(y + height/2, length(x0))
y1bottom <- rep(y - height/2, length(x0))
f <- x0 > end
if (sum(f) > 0L) {
y0top[f] <- y + height/2*(x0[f] - end)/xlen
y0bottom[f] <- y - height/2*(x0[f]-end)/xlen
x0[f] <- end
}
f <- x1 < start
if (sum(f) > 0L) {
y1top[f] <- y + height/2*(1.0-(start - x1[f])/xlen)
y1bottom[f] <- y - height/2*(1.0-(start - x1[f])/xlen)
x1[f] <- start
}
f <- x1 < x0
if (sum(f) > 0L) {
segments(x0[f], y0top[f], x1[f], y1top[f], col=col, lty=lty, lwd=lwd)
segments(x0[f], y0bottom[f], x1[f], y1bottom[f], col=col, lty=lty, lwd=lwd)
}
}
}
f <- x$strand == "-"
if (sum(f) > 0) {
xneg <- x[f,]
for (i in 1:nrow.feat(xneg)) {
start <- xneg[i,]$start-0.5
end <- xneg[i,]$end+0.5
x0 <- seq(from=start - arrow.density*floor(xlen/arrow.density),
to = end+xlen, by=arrow.density)
# x0 <- seq(from=start-xlen, to=end+xlen, by=arrow.density)
x1 <- x0 + xlen
y0top <- rep(y, length(x0))
y0bottom <- rep(y, length(x0))
y1top <- rep(y + height/2, length(x0))
y1bottom <- rep(y - height/2, length(x0))
f <- x1 > end
if (sum(f) > 0L) {
y1top[f] <- y + height/2*(1.0-(x1[f] - end)/xlen)
y1bottom[f] <- y - height/2*(1.0-(x1[f]-end)/xlen)
x1[f] <- end
}
f <- x0 < start
if (sum(f) > 0L) {
y0top[f] <- y + height/2*(start - x0[f])/xlen
y0bottom[f] <- y - height/2*(start - x0[f])/xlen
x0[f] <- start
}
f <- x0 < x1
if (sum(f) > 0L) {
segments(x0[f], y0top[f], x1[f], y1top[f], col=col, lty=lty, lwd=lwd)
segments(x0[f], y0bottom[f], x1[f], y1bottom[f], col=col, lty=lty, lwd=lwd)
}
}
}
}
invisible(NULL)
}
##' make gene plot
##' @title Gene plot
##' @param x An object of type \code{feat}
##' @param y the location of the plot on the y axis
##' @param height the height of the boxes
##' @param arrow.density The density of the arrows in arrows per inch
##' @param angle angle (in degrees) of the arrow heads
##' @param col color to use for plotting
##' @param lty line type for arrows, borders, and shading
##' @param lwd line width for arrows, borders and shading
##' @param add if \code{TRUE}, add to existing plot
## @param labels whether to label the boxes. If a vector of strings gives
## the lables for each element. If \code{TRUE}, use x$feature for labels.
##' @param xlim A numerical vector of length 2 giving the range for the x-axis.
##' @param ylim A numerical vector of length 2 giving the range for the y-axis.
##' @param ... graphical parameters to be passed to \code{plot}.
##' @method plot gene
##' @export plot.gene
##' @export
##' @keywords features plot
##' @example inst/examples/plot-gene.R
##' @author Melissa J. Hubisz
plot.gene <- function(x, y=0, height=1,
arrow.density=5,
angle = 30,
col="black",
lty=par("lty"),
lwd=par("lwd"),
add=FALSE,
# labels=FALSE,
xlim=range.feat(x),
ylim=c(y-height*3/4, y+height*3/4), ...) {
if (!is.null(x$externalPtr))
x <- from.pointer.feat(x)
if (is.null(x$start) || is.null(x$end))
stop("invalid features object")
seqname <- unique(x$seqname)
if (length(seqname) > 1L)
stop("feature has multiple sequences, can only plot 1")
if (!add)
plot(c(0), c(0), type="n", xlim=xlim, ylim=ylim, ...)
f <- x$end >= xlim[1] & x$start <= xlim[2]
if (sum(f) == 0) return(invisible(NULL))
x <- x[f,]
f <- x$feature == "intron"
if (sum(f) > 0L)
plot.feat(x[f,], plottype="l", y=y, lty=1, col=col, add=TRUE)
fexon <- x$feature == "exon"
fcds <- x$feature == "CDS"
if (sum(fexon) > 0L) {
if (sum(fcds) > 0L) {
region.bounds <- feat(seqname=seqname, start=min(x$start), end=max(x$end))
exon <- coverage.feat(x[fexon,], x[fcds,], region.bounds,
get.feats=TRUE,
not=c(FALSE, TRUE, FALSE))
} else exon <- x[fexon,]
if (nrow.feat(exon) > 0L)
plot.feat(exon, plottype="a", y=y, height=height/5, add=TRUE, col=col,
angle=angle, arrow.density=arrow.density)
}
if (sum(fcds) > 0L)
plot.feat(x[fcds,], plottype="b", y=y, height=height, add=TRUE, col=col,
angle=angle, arrow.density=arrow.density)
}
##' Features kernel density
##' @param x An object of type \code{feat}
##' @param type a character string, denoting the value to compute
##' the density for. Currently the only valid types are "length"
##' and "score"
##' @param ... additional arguments to be passed to \code{density}
##' @return A kernel density object as defined by \code{\link{density}}
##' @export
##' @export density.feat
##' @method density feat
##' @keywords features
##' @importFrom stats density
##' @author Melissa J. Hubisz
density.feat <- function(x, type="length", ...) {
if (type == "length") {
if (!is.null(x$externalPtr)) {
vals <- .Call.rphast("rphast_gff_lengths", x$externalPtr)
} else {
vals <- x$end - x$start
}
} else if (type == "score") {
if (!is.null(x$externalPtr)) {
vals <- .Call.rphast("rphast_gff_getScores", x$externalPtr)
} else {
vals <- x$score
}
} else stop("unknown type (should be \"length\" or \"score\")")
density(vals, ...)
}
##' plot histogram of feature lengths
##' @param x an object of type \code{feat}
##' @param type a character string, denoting the value to make the histogram with.
##' Currently the only valid types are "length" or "score"
##' @param ... additional arguments to be passed to \code{hist}
##' @method hist feat
##' @export hist.feat
##' @export
##' @keywords features
##' @author Melissa J. Hubisz
hist.feat <- function(x, type="length", ...) {
if (type == "length") {
if (!is.null(x$externalPtr)) {
starts <- .Call.rphast("rph_gff_starts", x$externalPtr)
ends <- .Call.rphast("rph_gff_ends", x$externalPtr)
vals <- ends - starts + 1
} else {
vals <- x$end - x$start + 1
}
} else if (type == "score") {
if (!is.null(x$externalPtr)) {
vals <- .Call.rphast("rph_gff_scores", x$externalPtr)
} else {
vals <- x$score
}
} else stop("unknown type (should be \"length\" or \"score\")")
hist(vals, ...)
}
##' Feature overlap
##'
##' Creates a features object containing all the features from one set which overlap
##' features from another.
##' @param x An object of type \code{feat} containing features to select
##' @param filter An object of type \code{feat} which determines which elements of x to select
##' @param numbase The number of bases of overlap between x and filter required to choose
##' a record. Use NULL to ignore (but then min.percent must be defined)
##' @param min.percent The minimum percent that a record must overlap with the combined records in filter
##' in order to be chosen
##' @param overlapping If \code{FALSE}, choose records with less than numbase overlapping bases,
##' and less than min.percent fraction overlap if min.percent is not \code{NULL}
##' @param get.fragments If \code{FALSE}, entire records are selected from x based on whether
##' they meet selection criteria.
##' If \code{TRUE}, return only the fragments of x that overlap
##' with filter. In this case, the same fragments may be output multiple times, if they are
##' selected by multiple entries in filter. numbase and min.percent apply in either case.
##' When this option is used, the return value is a list with two gffs. The first (named frags)
##' contains the overlapping fragments, and the second (filter.frags) contain the fragments
##' from filter which selected the overlapping fragments.
##' @param pointer.only If \code{TRUE}, the return object will only be a
##' pointer to an object stored in C (useful for very large features; advanced use only).
##' @return an object of type \code{feat} containing the selected entries from x (unless
##' get.fragments==TRUE, then it returns a list with two \code{feat} objects; see
##' get.fragments).
##' @note If either x or filter are feature objects stored as a pointer to C memory,
##' then this function may reorder the elements in these objects, but leave them
##' otherwise unchanged.
##' @export
##' @keywords features
##' @example inst/examples/overlap-feat.R
##' @author Melissa J. Hubisz
overlap.feat <- function(x, filter, numbase=1, min.percent=NULL,
overlapping=TRUE, get.fragments=FALSE, pointer.only=FALSE) {
check.arg(numbase, "numbase", "integer", null.OK=TRUE)
check.arg(min.percent, "min.percent", "numeric", null.OK=TRUE)
if (is.null(numbase) && is.null(min.percent)) stop("one of numbase or min.percent should not be NULL")
if (!is.null(numbase) && numbase < 0) stop("numbase should be at least 1 (if it isn't NULL)")
if (!is.null(min.percent) && (min.percent < 0 || min.percent > 1))
stop("min.percent should be NULL or in the range (0,1)")
check.arg(overlapping, "overlapping", "logical", null.OK=FALSE)
check.arg(get.fragments, "get.fragments", "logical", null.OK=FALSE)
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
if (is.null(filter$externalPtr))
filter <- as.pointer.feat(filter)
rv <- list()
rv[[1]] <- .Call.rphast("rph_gff_overlapSelect",
x$externalPtr, filter$externalPtr,
numbase, min.percent, !overlapping,
get.fragments)
if (get.fragments) {
rv <- rphast.simplify.list(rv)
if (!pointer.only) {
rv$frags <- from.pointer.feat(rv$frags)
rv$filter.frags <- from.pointer.feat(rv$filter.frags)
}
return(rv)
}
feat.rv <- .makeObj.feat(TRUE)
feat.rv$externalPtr <- rv[[1]]
if (!is.null(rv) && !pointer.only) {
feat.rv <- from.pointer.feat(feat.rv)
}
if (nrow.feat(feat.rv) == 0) return(NULL)
feat.rv
}
##' Get inverse features
##' @param x An object of type \code{feat}
##' @param region.bounds An object of type \code{feat} which defines
##' the boundaries of all relevant chromosomes in the first argument
##' @param pointer.only If \code{TRUE}, return a pointer to a structure stored
##' in C (advanced use only).
##' @return An object of type \code{feat} which contains all regions in
##' region.bounds that are not in the first argument
##' @note If x is stored as a pointer to C memory, then its elements will
##' be sorted by this function. region.bounds will not be changed.
##' @export
##' @keywords features
##' @author Melissa J. Hubisz
inverse.feat <- function(x, region.bounds, pointer.only=FALSE) {
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
if (is.null(region.bounds$externalPtr))
region.bounds <- as.pointer.feat(region.bounds)
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_inverse",
x$externalPtr,
region.bounds$externalPtr)
if (!pointer.only)
rv <- from.pointer.feat(rv)
rv
}
##' Features coverage
##' @param ... objects of type \code{feat}
##' @param or if \code{TRUE}, get the coverage of union of feat arguments.
##' or is \code{FALSE} by default, which takes the intersection.
##' @param not If not \code{NULL}, a vector of logicals the same length
##' as the number of features
##' provided (or will be recycled to this length).
##' For each value which is \code{TRUE}, then any base *not* included in this
##' feature will be counted. (The negation is done before any other operation).
##' If \code{NULL}, do not negate any features. There must be at least one
##' feature which is not negated (so that boundaries can be established).
##' @param get.feats if \code{TRUE}, return an object of type \code{feat}
##' representing the intersection (or union of \code{or==TRUE}) of the
##' features
##' @param pointer.only (Only used if \code{get.feats==TRUE}). If \code{TRUE},
##' the features object returned will be stored as a pointer to an object
##' in C.
##' @return The number of bases covered by the feat arguments, or the
##' combined feat object if \code{get.feats==TRUE}.
##' @note Any features object passed into this function which is stored as a
##' pointer to an object stored in C may be reordered (sorted) by this function.
##' @export
##' @keywords features
##' @example inst/examples/coverage-feat.R
##' @author Melissa J. Hubisz
coverage.feat <- function(..., or=FALSE, not=FALSE, get.feats=FALSE,
pointer.only=FALSE) {
check.arg(or, "or", "logical", null.OK=FALSE)
check.arg(get.feats, "get.feats", "logical", null.OK=FALSE)
if (get.feats)
check.arg(pointer.only, "pointer.only", "logical", null.OK=FALSE)
featlist <- list(...)
check.arg(not, "not", null.OK=TRUE, min.length=1L,
max.length=length(featlist))
if (is.null(not)) {
not <- rep(FALSE, length(featlist))
} else {
not <- rep(not, length.out=length(featlist))
}
if (sum(not == FALSE) == 0L)
stop("at least one feature must have not==FALSE")
region.bounds <- featlist[[which(not==FALSE)[1]]]
for (i in 1:length(featlist)) {
x <- featlist[[i]]
if (not[i]) x <- inverse.feat(x, region.bounds)
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
featlist[[i]] <- x$externalPtr
}
if (get.feats) {
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_featureBits", featlist,
or, get.feats)
if (pointer.only) return(rv)
return(from.pointer.feat(rv))
}
.Call.rphast("rph_gff_featureBits", featlist, or, get.feats)
}
##' Add UTRs to features
##' @param x An object of type \code{feat}. CDS regions must be present with type "CDS", and
##' the transcript_id must be indicated in the attribute field.
##' @return An object of type \code{feat}, with all the entries of the original object, but
##' also with UTR annotations.
##' @note If x is stored as a pointer to an object stored in C, then UTRs will be added
##' to x.
##' @export
##' @keywords features
##' @example inst/examples/add-UTRs-feat.R
##' @author Melissa J. Hubisz and Adam Siepel
add.UTRs.feat <- function(x) {
if (is.null(x$externalPtr)) {
x <- as.pointer.feat(x)
getDataFrame <- TRUE
} else getDataFrame <- FALSE
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_add_UTRs", x$externalPtr)
if (getDataFrame) return(from.pointer.feat(rv))
rv
}
##' Add introns to features
##' @param x An object of type \code{feat}. CDS regions must be present with type "CDS", and
##' the transcript_id must be indicated in the attribute field.
##' @return An object of type \code{feat}, with all the entries of the original object, but
##' also with intron annotations.
##' @note If x is stored as a pointer to an object stored in C, introns will be
##' added to x.
##' @export
##' @keywords features
##' @example inst/examples/add-introns-feat.R
##' @author Melissa J. Hubisz and Adam Siepel
add.introns.feat <- function(x) {
if (is.null(x$externalPtr)) {
x <- as.pointer.feat(x)
getDataFrame <- TRUE
} else getDataFrame <- FALSE
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_add_introns", x$externalPtr)
if (getDataFrame) return(from.pointer.feat(rv))
rv
}
##' Add start/stop codon, 3'/5' splice signals to features
##' @param x An object of type \code{feat}. CDS regions must be present with type "CDS", and
##' the transcript_id must be indicated in the attribute field.
##' @return An object of type \code{feat}, with all the entries of the original object, but
##' also with stop codons, start, codons, 3' splice, and 5' splice sites annotated.
##' @note \itemize{
##' \item{If x is stored as a pointer to an object stored in C, signals will be
##' added to x.}
##' \item{Does not correctly handle case of splice site in middle of start
##' or stop codon.}}
##' @export
##' @keywords features
##' @example inst/examples/add-signals-feat.R
##' @author Melissa J. Hubisz and Adam Siepel
add.signals.feat <- function(x) {
if (is.null(x$externalPtr)) {
x <- as.pointer.feat(x)
getDataFrame <- TRUE
} else getDataFrame <- FALSE
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_add_signals", x$externalPtr)
if (getDataFrame) return(from.pointer.feat(rv))
rv
}
##' Fix start and stop signals
##' @param x An object of type \code{feat}. CDS regions must be present with type
##' "CDS", and the transcript_id must be indicated in the attribute field.
##' Start and stop codons should have feature type "start_codon" and "stop_codon"
##' (as produced by addSignals.feat).
##' @return An object of type \code{feat}, in which CDS regions are ensured to
##' include start codons and exclude stop codons, as required by the GTF2 standard.
##' @note \itemize{
##' \item{If x is stored as a pointer to an object stored in C, signals will be
##' added to x.}
##' \item{Assumes at most one start_codon and at most one stop_codon per transcript.}
##' }
##' @export
##' @keywords features
##' @example inst/examples/fix-start-stop-feat.R
##' @author Melissa J. Hubisz and Adam Siepel
fix.start.stop.feat <- function(x) {
if (is.null(x$externalPtr)) {
x <- as.pointer.feat(x)
getDataFrame <- TRUE
} else getDataFrame <- FALSE
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_fix_start_stop", x$externalPtr)
if (getDataFrame) return(from.pointer.feat(rv))
rv
}
##' concatenate feature objects
##' @param ... objects of type \code{feat} to be combined into a single object
##' @return An object of type \code{feat} containing entries from all
##' given features
##' @method rbind feat
##' @export rbind.feat
##' @export
##' @keywords features
##' @author Melissa J. Hubisz and Adam Siepel
rbind.feat <- function(...) {
feat <- .makeObj.feat(TRUE)
featlist <- list(...)
idx <- 1
p.featlist <- list()
for (i in 1:length(featlist)) {
currfeat <- featlist[[i]]
if (!is.null(currfeat) && nrow.feat(currfeat) > 0L) {
if (is.null(currfeat$externalPtr))
currfeat <- as.pointer.feat(currfeat)
p.featlist[[idx]] <- currfeat$externalPtr
idx <- idx+1
}
}
if (idx == 1) return(NULL)
feat$externalPtr <- .Call.rphast("rph_gff_append", p.featlist)
from.pointer.feat(feat)
}
##' Split features by length
##' @param x An object of type \code{feat}
##' @param f The maximum length of features in new object. Can be a
##' vector giving a different length for each row, or a single numeric
##' value. Values will be recycled to the same length
##' as \code{nrow.feat(x)}.
##' @param drop A logical value saying whether to drop "left-over" elements
##' which do not have exactly length f.
##' @param start.from A character string, current valid values are "left"
##' (start split at smallest coordinate for each feature), or "right"
##' (start splitting at the last coordinate and work down). Values will
##' be recycled to the length of \code{nrow.feat(x)}
##' @param pointer.only If \code{TRUE}, return an object which is a pointer
##' to a features object stored in C (advanced use only).
##' @param ... Currently not used (for S3 compatibility).
##' @return An object of type \code{feat} with the same features as x but
##' with all features of length > max.length broken into segments
##' (starting from the first position in feature). The last piece
##' of each split segment may be smaller than max.length
##' @method split feat
##' @export split.feat
##' @export
##' @keywords features
##' @author Melissa J. Hubisz
split.feat <- function(x, f, drop=FALSE, start.from="left",
pointer.only=FALSE, ...) {
featSize <- nrow.feat(x)
check.arg(f, "f", "integer", null.OK=FALSE,
min.length=1, max.length=NULL)
max.length <- f
check.arg(drop, "drop", "logical", null.OK=FALSE)
check.arg(start.from, "start.from", "character", null.OK=FALSE,
min.length=1L, max.length=nrow.feat(x))
if (sum(start.from!="left" & start.from!="right") > 0L)
stop("start.from invalid value (should be \"left\" or \"right\")")
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
splitFeat <- .makeObj.feat(TRUE)
splitFeat$externalPtr <- .Call.rphast("rph_gff_split", x$externalPtr,
max.length, drop,
ifelse(start.from=="left", 0, 1))
if (!pointer.only)
splitFeat <- from.pointer.feat(splitFeat)
splitFeat
}
##' Sort a GFF
##' @param x An object of type \code{feat}
##' @param decreasing Set to TRUE to sort from highest to lowest coordinates
##' @param ... Currently not used
##' @return An object of type \code{feat} sorted primarily by
##' seqname, then by start position, then by end position.
##' @method sort feat
##' @export sort.feat
##' @export
##' @note If x is stored as a pointer to an object in C, the object will
##' be modified to the return value.
##' @keywords features
##' @author Melissa J. Hubisz and Adam Siepel
sort.feat <- function(x, decreasing = FALSE, ...) {
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_sort", x$externalPtr)
rv <- from.pointer.feat(rv)
if (decreasing)
rv <- rv[dim(rv)[1]:1,]
rv
}
##' Composition of features with respect to annotations
##' @param x An object of type \code{feat}.
##' @param annotations An object of type \code{feat} containing
##' some annotations.
##' @return A data frame with two columns and a row for
##' each type of element in the annotations. The second
##' column gives the fraction of
##' x which fall in the corresponding annotation type.
##' Given non-overlapping annotations which cover the
##' entire range of interest, the second column should sum
##' to 1 (otherwise not).
##' @export
##' @note If x or annotations are passed to this function as pointers to objects
##' stored in C, they will be sorted after the function call.
##' @keywords features
##' @author Melissa J. Hubisz
composition.feat <- function(x, annotations) {
if (!is.null(annotations$externalPtr))
annotations <- from.pointer.feat(annotations)
if (!is.null(x$externalPtr))
x <- from.pointer.feat(x)
annTypes <- unique(annotations$feature)
rv <- list()
for (anntype in annTypes) {
annfeat <- annotations[annotations$feature == anntype,]
rv[[anntype]] <- coverage.feat(x, annfeat)/coverage.feat(x)
}
data.frame(type=names(rv), composition=as.numeric(rv), stringsAsFactors=TRUE)
}
##' Enrichment of features with respect to annotation types
##' @param x An object of type \code{feat}
##' @param annotations An object of type \code{feat} containing
##' some annotations.
##' @param region.bounds An object of type \code{feat} representing
##' the boundary coordinates of the regions of interest (such as chromosome
##' boundaries). All elements from the first two arguments should fall
##' entirely within region.bounds.
##' @return A data frame with two columns and a row for each type of
##' element in the annotations. The second column gives the
##' fold-enrichment
##' of x across the corresponding annotation type,
##' which is equal to the
##' fraction of x which fall within the annotation type,
##' divided by the fraction of the entire region covered by the
##' annotation type.
##' @export
##' @note If any of the arguments to this function are passed as pointers to
##' objects stored in C, they will be sorted after this function call.
##' @keywords features
##' @author Melissa J. Hubisz
enrichment.feat <- function(x, annotations, region.bounds) {
if (!is.null(annotations$externalPtr))
annotations <- from.pointer.feat(annotations)
if (!is.null(x$externalPtr))
x <- from.pointer.feat(x)
annTypes <- unique(annotations$feature)
rv <- list()
totalNumBase <- coverage.feat(region.bounds)
for (anntype in annTypes) {
annfeat <- annotations[annotations$feature == anntype,]
rv[[anntype]] <- as.numeric(coverage.feat(x, annfeat, region.bounds))*as.numeric(totalNumBase)/(as.numeric(coverage.feat(annfeat, region.bounds))*as.numeric(coverage.feat(x, region.bounds)))
}
data.frame(type=names(rv), enrichment=as.numeric(rv), stringsAsFactors=TRUE)
}
##' Remove overlapping genes
##' @param x An object of type \code{feat}, usually read from a genepred
##' file. Should have attributes labelled "transcript_id" which identify
##' features belonging to the same gene.
##' @param incomparables Not currently used (present for S3 compatibility).
##' @param ... Not currently used.
##' @return An object of type \code{feat} with overlapping genes removed.
##' If the features are scored, then the feature with the highest score
##' is kept; otherwise the feature with the longest length. If
##' x is a pointer to an object stored in C, the return value will also
##' be a pointer (and x will be altered to the return value).
##' @note \itemize{
##' \item{Long UTRs can have undesirable effects; may want to filter these
##' out first.}
##' \item{If x is a pointer to an object in C, it will be modified (to
##' the return value).}}
##' @method unique feat
##' @export unique.feat
##' @export
##' @keywords features
##' @note VERY IMPORTANT: this function is not currently implemented to look
##' at chromosomes (ie the seqname field of the feature). Therefore any
##' genes which have overlapping coordinates REGARDLESS OF THE CHROMOSOME
##' will be pruned to a single "non-overlapping" gene. To get around this,
##' first subset the features by chromosome and call uniq.feat on each
##' subset.
##' @note Also, this algorithm considers genes to be overlapping even if
##' they are on different strands. If this is undesirable, then subset
##' the features by strand as well as chromosome.
##' @author Melissa J. Hubisz and Adam Siepel
unique.feat <- function(x, incomparables=FALSE, ...) {
if (is.null(x$externalPtr)) {
x <- as.pointer.feat(x)
getDataFrame <- TRUE
} else getDataFrame <- FALSE
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_nonOverlapping_genes", x$externalPtr)
if (getDataFrame)
rv <- from.pointer.feat(rv)
rv
}
##' Extract value from tag-value formatted attribute in features object
##' @param x A features object of type \code{feat}. The attribute field
##' should be in tag-value format (either GFF 2 standard; ie,
##' "tag1 val1a val1b; tag2 val2 ; ...", or, GFF 3 standard; ie,
##' "tag1=val1a,val1b;tag2=val2; ...".
##' where vals are in quotes if they are strings.
##' @param tag The tag whose values are to be extracted.
##' @return If there is at most one relevant value for each feature,
##' a character vector of the same length as x will be returned, containing
##' the value for each feature, or NA where the tag does not exist for that
##' feature. If some elements have multiple values, then the return value
##' will be a list with the same length as x, each element being a character
##' vector containing the values for the corresponding element of x (or
##' NA for no value).
##' @export
##' @keywords features GFF
##' @example inst/examples/tagval-feat.R
##' @author Melissa J. Hubisz
tagval.feat <- function(x, tag) {
check.arg(tag, "tag", "character", null.OK=FALSE)
if (is.null(x$externalPtr)) {
if (is.null(x$attribute)) return (rep(NA, nrow(x)))
x <- as.pointer.feat(x)
}
rv <- rphast.simplify.list(.Call.rphast("rph_gff_one_attribute",
x$externalPtr, tag))
maxlen <- max(sapply(rv, length))
if (maxlen == 1L) rv <- as.character(rv)
f <- rv==""
if (sum(f) > 0L)
rv[f] <- NA
rv
}
##' Extract value from tag-value formatted attributes
##' @param x A vector of character strings in tag-val format (as
##' described in the GFF 2 standard; ie, "tag1 val1a val1b; tag2 val2 ; ...",
##' or in the GFF 3 format "tag1=val1a,val1b; tag2=val2; ..."),
##' where vals are in quotes if they are strings.
##' @param tag The tag whose values are to be extracted.
##' @return If there is at most one value per tag for each element of x,
##' a character vector of the same length as x will be returned, containing
##' the value for each element, or NA if the tag does not exist for that
##' element. If some elements have multiple values, then the return value
##' will be a list with the same length as x, each element being a character
##' vector containing the values for the corresponding element of x (or
##' NA for no value).
##' @export
##' @keywords GFF
##' @example inst/examples/tagval.R
##' @author Melissa J. Hubisz
tagval <- function(x, tag) {
check.arg(x, "x", "character", null.OK=FALSE, min.length=1L, max.length=NULL)
tagval.feat(x=feat(start=rep(1, length(x)),
end=rep(2, length(x)), attribute=x),
tag=tag)
}
##' Combine adjacent features with the same "feature" field
##' @param x An object of type feat
##' @param weightedAverageScore If TRUE, scores of merged features
##' are given by the average of those features, weighted by their
##' original lengths. Otherwise, scores of merged features are
##' simply summed.
##' @param minScore If TRUE, scores of merged features are given
##' by the minimum score of the combined features.
##' @return A features object in which adjacent features are
##' combined into one longer feature.
##' @note If x is stored as a pointer to a C structure, then
##' x will be modified to the return value.
##' @author Melissa J. Hubisz and Adam Siepel
##' @export
flatten.feat <- function(x, weightedAverageScore=FALSE, minScore=FALSE) {
if (minScore && weightedAverageScore) {
stop("flatten.feat: only one of minScore, weightedAverageScore can be TRUE")
}
isPointer <- !is.null(x$externalPtr)
if (!isPointer) x <- as.pointer.feat(x)
if (!is.logical(weightedAverageScore))
stop("Expected weightedAverageScore to be a logical")
sort.feat(x)
.Call.rphast("rph_gff_flatten", x$externalPtr, weightedAverageScore,
minScore)
if (!isPointer) return(from.pointer.feat(x))
x
}
##' Read a wig file
##'
##' Reads fixed or variable step wig files. Stores them as a
##' features object.
##' @param file The file to read
##' @param pointer.only If \code{TRUE}, store as a pointer to a
##' C structure
##' @return A GFF object representing data in wig file
##' @export
read.wig <- function(file, pointer.only=FALSE) {
read.feat(file, pointer.only=pointer.only)
}
##' Write a fixedStep wig file
##' @title Writing a wig file
##' @param chrom A character vector giving chromosome name for each point.
##' Will be recycled to length(start)
##' @param start An integer vector giving start coordinate for each point.
##' @param score A numeric vector giving score at each point
##' Will be recycled to length(start)
##' @param span An integer giving span (ie, length) of each element (all
##' elements must have the same length, so only a single value is allowed).
##' @param file The name of the file to write to (will be overwritten).
##' A value of NULL implies write to console.
##' @param append Whether to append to the file. If FALSE, file will be
##' overwritten.
##' @keywords wig fixedStep
##' @author Melissa J. Hubisz
##' @export
##' @example inst/examples/write-wig.R
write.wig <- function(chrom, start, score, span=1, file=NULL, append=FALSE) {
chrom <- check.arg(chrom, "chrom", "character", null.OK=FALSE, max.length=NULL)
start <- check.arg(start, "start", "integer", null.OK=FALSE, max.length=NULL)
score <- check.arg(score, "score", "numeric", null.OK=FALSE, max.length=NULL)
span <- check.arg(span, "span", "integer", null.OK=FALSE)
file <- check.arg(file, "file", "character", null.OK=TRUE)
append <- check.arg(append, "append", "logical", null.OK=FALSE)
if (span <= 0) stop("span should be integer >= 1")
if (length(chrom) > length(start))
stop("length(chrom) > length(start)")
if (length(score) > length(start))
stop("length(score) > length(start)")
invisible(.Call.rphast("rph_wig_print",
feat(seqname=chrom, start=start, end=start+span-1,
score=score, pointer.only=TRUE)$externalPtr,
file,
append))
}
##' Write a features object in fixedStep wig format
##' @param x An object of type feat
##' @param file The name of the file to write to. A value of NULL
##' implies write to console.
##' @param append If TRUE, append to the file. Otherwise overwrite.
##' @keywords wig fixedStep GFF features
##' @author Melissa J. Hubisz
##' @note Wig format only contains chromosome, coordinates, and score. Any
##' other data will be lost.
##' @note This function will quit with an error if the elements of x are
##' not all the same length (as required by fixedStep wig format).
##' @note If x is stored as a pointer to a C structure, the elements
##' will be sorted by this function.
##' @example inst/examples/write-wig-feat.R
##' @export
write.wig.feat <- function(x, file=NULL, append=FALSE) {
append <- check.arg(append, "append", "logical", null.OK=FALSE)
file <- check.arg(file, "file", "character", null.OK=TRUE)
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
invisible(.Call.rphast("rph_wig_print", x$externalPtr, file, append))
}
##' Convert coordinates from one frame of reference to another
##'
##' Converts coordinates of features in a GFF according to a multiple
##' alignment. Will map from the coordinate system of any sequence
##' to any other sequence; can also map to or from the coordinate system
##' of the entire alignment.
##' @param x A features object; if the from parameter is -1, the first
##' column should indicate which frame of reference is used for that row
##' (ie, the species name).
##' For coordinates in the frame of reference of the entire alignment, set
##' the first column to "MSA".
##' @param align An msa object containing the alignment
##' @param from A single character string or integer, used to indicate the
##' current frame of reference of the features. If the rows of the features are
##' not all in the same frame of reference, set this to -1 and indicate the
##' frame of reference in the 1st column of the features. Otherwise, it can
##' be specified here as a single integer (from 0 to nrow.msa(align)), with
##' 0 indicating the frame of reference of the entire alignment, and 1
##' indicating the 1st species, 2 the second, etc. Or it can be a single
##' character string giving the name of the species, with "MSA" indicating
##' the entire alignment.
##' @param to A single character string or integer, used to indicate the frame
##' of reference to convert to. This is specified in the same way as the
##' "from" argument, above, except that -1 is not an option.
##' @return A features object with elements in the frame of refernece
##' indicated by the "to" argument.
##' @note Ignores any offset in MSA. All coordinates should start with the
##' first position in the alignment as 1.
##' @note If the endpoints of an element have gaps in the "to" species,
##' the elements will be truncated
##' @example inst/examples/convert-coords-feat.R
##' @author Melissa J. Hubisz and Adam Siepel
##' @export
convert.coords.feat <- function(x, align, from=-1, to=1) {
x <- copy.feat(x)
if (is.null(x$externalPtr))
x <- as.pointer.feat(x)
if (!is.msa(align)) stop("align should be of type msa")
if (length(from) != 1L) {
cat("from=", from, "\n")
cat("length(from)=", length(from), "\n")
stop("from should be character or integer of length 1")
}
is.wholenumber <-
function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
if (is.character(from)) {
if (from=="MSA") from <- 0
else {
if (!is.element(from, names.msa(align))) {
stop("argument from does not match any names in the alignment")
}
from <- which(names.msa(align)==from)
}
} else if (is.numeric(from)) {
if (from < -1 || from > nrow.msa(align) || !is.wholenumber(from))
stop("argument from should be in the range -1, .., ", nrow.msa(align))
} else
stop("from should be a character or integer of length 1")
# now from is an integer of length 1
# do same for to
if (length(to) != 1L)
stop("to should be character or integer of length 1")
if (is.character(to)) {
if (to=="MSA") to <- 0
else {
if (!is.element(to, names.msa(align)))
stop("argument to does not match any names in the alignment")
to <- which(names.msa(align)==to)
}
} else if (is.numeric(to)) {
if (to < 0 || to > nrow.msa(align) || !is.wholenumber(to))
stop("argument to should be in the range 0, .., ", nrow.msa(align))
} else
stop("to should be a character or integer of length 1")
# now to is an integer of length 1
rv <- .makeObj.feat(TRUE)
rv$externalPtr <- .Call.rphast("rph_gff_convert_coords",
x$externalPtr,
as.pointer.msa(align)$externalPtr,
to)
from.pointer.feat(rv)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.