Nothing
R/ASEset-class.R
In AllelicImbalance: Investigates Allele Specific Expression
Defines functions
.calcFrequencyFromAlleleCounts
.returnMaternalPhaseFrequency
.unknownStrandCountsExists
.minusStrandCountsExists
.plusStrandCountsExists
#'@include initialize-methods.R
NULL
#' ASEset objects
#'
#' Object that holds allele counts, genomic positions and map-bias for a set of
#' SNPs
#'
#' An ASEset object differs from a regular RangedSummarizedExperiment object
#' in that the assays contains an array instead of matrix. This array has
#' ranges on the rows, sampleNames on the columns and variants in the third
#' dimension.
#'
#' It is possible to use the commands barplot and locationplot on an ASEset
#' object see more details in \code{\link{barplot}} and
#' \code{\link{locationplot}}.
#'
#' Three different alleleCount options are available. The simples one is the
#' * option, and is for experiments where the strand information is not
#' known e.g. non-stranded data. The unknown strand could also be for strand
#' specific data when the aligner could not find
#' any strand associated with the read, but this should normally not happen,
#' and if it does probably having an extremely low mapping quality.
#' Then there are an option too add plus
#' and minus stranded data. When using this, it is essential to make sure that
#' the RNA-seq experiment under analysis has in fact been created so that
#' correct strand information was obtained. The most functions will by default
#' have their strand argument set to '*'.
#'
#' The phase information is stored by the convention of
#' 'maternal chromosome|paternal chromosome', with 0 as reference allele and 1
#' as alternative allele. '|' when the phase is known and '/' when the phase is
#' unknown. Internally the information will be stored as an three dimensional
#' array, dim 1 for SNPs, dim 2 for Samples and dim 3 which is fixed and stores
#' maternal chromosome, paternal chromosome and phased (1 equals TRUE).
#'
#'
#' @name ASEset-class
#' @rdname ASEset-class
#' @aliases ASEset-class ASEset alleleCounts mapBias fraction arank
#' frequency genotype genotype<- phase phase<- alleleCounts,ASEset-method mapBias,ASEset-method
#' fraction,ASEset-method arank,ASEset-method
#' frequency,ASEset-method genotype,ASEset-method genotype<-,ASEset-method
#' alleleCounts<- alleleCounts<-,ASEset-method phase,ASEset-method phase<-,ASEset-method
#'
#' @docType class
#' @param x ASEset object
#' @param strand which strand of '+', '-' or '*'
#' @param verbose makes function more talkative
#' @param return.type return 'names', rank or 'counts'
#' @param return.class return 'list' or 'array'
#' @param top.fraction.criteria 'maxcount', 'ref' or 'phase'
#' @param value replacement variable
#' @param ... additional arguments
#' @return An object of class ASEset containing location information and allele
#' counts for a number of SNPs measured in a number of samples on various
#' strand, as well as mapBias information. All data is stored in a manner
#' similar to the SummarizedExperiment class.
#'
#' @section Table: table(...)
#'
#' \describe{
#' Arguments: \item{...}{An \code{ASEset object} that contains the
#' variants of interest}
#'
#' The generics for table does not easily allow more than one argument
#' so in respect to the different strand options, \code{table} will
#' return a SimpleList with length 3, one element for each strand.
#' }
#'
#' @section Frequency: frequency(x,
#' return.class = "list", strand = "*",
#' threshold.count.sample = 15)
#'
#' \describe{
#' Arguments: \item{x}{An \code{ASEset object} that contains the
#' variants of interest}
#'
#' \item{x}{threshold.count.samples} if sample has fewer counts the function
#' return NA.
#'
#' }
#' @section Constructor: ASEsetFromCountList(rowRanges, countListNonStranded =
#' NULL, countListPlus = NULL, countListMinus = NULL, countListUnknown = NULL,
#' colData = NULL, mapBiasExpMean = array(), verbose=FALSE, ...)
#'
#' \describe{
#'
#' Arguments: \item{rowRanges}{A \code{GenomicRanges object} that contains the
#' variants of interest} \item{countListNonStranded}{A \code{list} where each
#' entry is a matrix with allele counts as columns and sample counts as rows}
#' \item{countListPlus}{A \code{list} where each entry is a matrix with allele
#' counts as columns and sample counts as rows} \item{countListMinus}{A
#' \code{list} where each entry is a matrix with allele counts as columns and
#' sample counts as rows} \item{countListUnknown}{A \code{list} where each
#' entry is a matrix with allele counts as columns and sample counts as rows}
#' \item{colData}{A \code{DataFrame} object containing sample specific data}
#' \item{mapBiasExpMean}{A 3D \code{array} describing mapping bias. The SNPs
#' are in the 1st dimension, samples in the 2nd dimension and variants in the
#' 3rd dimension.} \item{verbose}{Makes function more talkative}
#' \item{...}{arguments passed on to SummarizedExperiment constructor} }
#'
#'
#' @author Jesper R. Gadin, Lasse Folkersen
#' @keywords class ASEset
#' @examples
#'
#'
#' #make example countList
#' set.seed(42)
#' countListPlus <- list()
#' snps <- c('snp1','snp2','snp3','snp4','snp5')
#' for(snp in snps){
#' count<-matrix(rep(0,16),ncol=4,dimnames=list(
#' c('sample1','sample2','sample3','sample4'),
#' c('A','T','G','C')))
#' #insert random counts in two of the alleles
#' for(allele in sample(c('A','T','G','C'),2)){
#' count[,allele]<-as.integer(rnorm(4,mean=50,sd=10))
#' }
#' countListPlus[[snp]] <- count
#' }
#'
#' #make example rowRanges
#' rowRanges <- GRanges(
#' seqnames = Rle(c('chr1', 'chr2', 'chr1', 'chr3', 'chr1')),
#' ranges = IRanges(1:5, width = 1, names = head(letters,5)),
#' snp = paste('snp',1:5,sep='')
#' )
#'
#' #make example colData
#' colData <- DataFrame(Treatment=c('ChIP', 'Input','Input','ChIP'),
#' row.names=c('ind1','ind2','ind3','ind4'))
#'
#' #make ASEset
#' a <- ASEsetFromCountList(rowRanges, countListPlus=countListPlus,
#' colData=colData)
#'
#'
#' #example phase matrix (simple form)
#' p1 <- matrix(sample(c(1,0),replace=TRUE, size=nrow(a)*ncol(a)),nrow=nrow(a), ncol(a))
#' p2 <- matrix(sample(c(1,0),replace=TRUE, size=nrow(a)*ncol(a)),nrow=nrow(a), ncol(a))
#' p <- matrix(paste(p1,sample(c("|","|","/"), size=nrow(a)*ncol(a), replace=TRUE), p2, sep=""),
#' nrow=nrow(a), ncol(a))
#'
#' phase(a) <- p
#'
#'
#' #generate ASEset from array
#' snps <- 999
#' samples <-5
#' ar <-array(rep(unlist(lapply(1:snps,
#' function(x){(sample(c(TRUE,FALSE,TRUE,FALSE), size = 4))})), samples),
#' dim=c(4,snps,samples))
#' ar2 <- array(sample(50:300, 4*snps*samples,replace=TRUE), dim=c(4,snps,samples))
#' ar2[ar] <- 0
#' ar2 <- aperm(ar2, c(2, 3, 1))
#' dimnames(ar2) <- list(paste("snp",1:snps,sep=""),paste("sample",1:samples,sep=""),
#' c("A","C","G","T"))
#' gr <- GRanges(seqnames=c("chr2"), ranges=IRanges(start=1:dim(ar2)[1], width=1), strand="*")
#' a <- ASEsetFromArrays(gr, countsUnknown=ar2)
#'
#'
#' @exportClass ASEset
#' @exportMethod alleleCounts alleleCounts<- mapBias fraction arank
#' frequency genotype genotype<- phase phase<-
#'
#' @export
setClass("ASEset", contains = "RangedSummarizedExperiment",
representation(variants = "vector"))
#' @rdname ASEset-class
#' @export
setGeneric("alleleCounts", function(x, strand = "*", return.class="list") {
standardGeneric("alleleCounts")
})
#' @rdname ASEset-class
#' @export
setMethod("alleleCounts", signature(x = "ASEset"), function(x, strand = "*",
return.class="list") {
#check that strand parameter is correct set
if (!sum(strand %in% c("+", "-", "*", "both")) > 0) {
stop("strand parameter has to be either '+', '-', '*' or 'both' ")
}
#check if the assays are present
if(strand=="+" ){
if(!"countsPlus" %in% assayNames(x)){
stop(paste("strand",strand,"is not present in ASEset object",sep=""))
}
}
if(strand=="-" ){
if(!"countsMinus" %in% assayNames(x)){
stop(paste("strand",strand,"is not present in ASEset object",sep=""))
}
}
if(strand=="*" ){
if(!("countsMinus" %in% assayNames(x) | "countsPlus" %in% assayNames(x))){
stop(paste("for strand '*' at least one of '+' or '-' is required to be",
" present in ASEset object. Old ASEsets set the '*' strand in",
" countsUnknown, but that has been deprecated. If so, move your counts",
" to the plus or minus slot instead ",
" assays(yourASEset, withDimnames=FALSE)[['countsPlus']] <- assays(yourASEset)[['countsUnknown']]",
sep=""))
}
}
if(strand=="both" ){
if(!"countsMinus" %in% assayNames(x) & "countsPlus" %in% assayNames(x)){
stop(paste("for strand 'both' both '+' and '-' is required to be",
" present in the ASEset object",sep=""))
}
}
#access the correct assays
if (strand == "+") {
ar <- assays(x)[["countsPlus"]]
} else if (strand == "-") {
ar <- assays(x)[["countsMinus"]]
} else if (strand == "*") {
if(!("countsMinus" %in% assayNames(x))){
ar <- assays(x)[["countsPlus"]]
}else if(!("countsPlus" %in% assayNames(x))){
ar <- assays(x)[["countsMinus"]]
}else{
ar <- assays(x)[["countsMinus"]] + assays(x)[["countsPlus"]]
}
} else if (strand == "both") {
ar <- array(c(assays(x)[["countsPlus"]], assays(x)[["countsMinus"]]), dim=c(dim(assays(x)[["countsPlus"]]), 2))
} else {
stop("not existing strand option")
}
if(return.class=="array"){
dimnames(ar)[[3]]<- x@variants
if (strand == "both") {
dimnames(ar)[[4]]<- c("+","-")
}
ar
}else if(return.class=="list"){
if(strand=="both"){
strands <- 2
alleleCountList2 <- list()
}else{
strands <- 1
}
for( j in 1:strands){
alleleCountList <- list()
for (i in 1:nrow(ar)) {
if(strand=="both"){
mat <- ar[i, , ,j]
}else{
mat <- ar[i, , ]
}
if (class(mat)[1] == "integer") {
mat <- t(as.matrix(mat))
}
if (class(mat)[1] == "numeric") {
mat <- t(mat)
}
colnames(mat) <- x@variants
rownames(mat) <- colnames(x)
alleleCountList[[i]] <- mat
}
# add snp id
names(alleleCountList) <- rownames(x)
if(strand=="both"){
alleleCountList2[[j]] <- alleleCountList
}
}
if(strand=="both"){
names(alleleCountList2) <- c("+","-")
alleleCountList2
}else{
alleleCountList
}
}else{
stop("return.class has to be 'list' or 'array'")
}
})
### -------------------------------------------------------------------------
### helpers for alleleCounts
###
.plusStrandCountsExists <- function(x){
"countsPlus" %in% assayNames(x)
}
.minusStrandCountsExists <- function(x){
"countsMinus" %in% assayNames(x)
}
.unknownStrandCountsExists <- function(x){
"countsUnknown" %in% assayNames(x)
}
#' @rdname ASEset-class
#' @export
setGeneric("alleleCounts<-", function(x, ..., value) {
standardGeneric("alleleCounts<-")
})
#' @rdname ASEset-class
#' @export
setReplaceMethod("alleleCounts", "ASEset",
function(x, strand="*", return.class="array", ..., value) {
#for dev
#message("valueX: ", paste0(value, collapse=", "))
if (!sum(strand %in% c("+", "-", "*")) > 0) {
stop("strand parameter has to be either '+', '-', '*' ")
}
if (strand == "+") {
el <- "countsPlus"
} else if (strand == "-") {
el <- "countsMinus"
} else if (strand == "*") {
el <- "countsUnknown"
} else {
stop("not existing strand option")
}
#check that value has the right dimensions (make validObject method)
#if(!all(dim(value)==dim(assays(x, withDimnames=FALSE)[["countsUnknown"]]))){
# stop("dimensions in replacement object is not correct")
#}
assays(x, withDimnames=FALSE)[[el]][] <- value
#return object
x
})
#' @rdname ASEset-class
#' @export
setGeneric("mapBias", function(x, ...) {
standardGeneric("mapBias")
})
#' @rdname ASEset-class
#' @export
setMethod("mapBias", signature(x = "ASEset"), function(x,
return.class="list") {
# assume alleleCount information is stored as element 1
if(return.class=="array"){
return(assays(x)[["mapBias"]])
}else if(return.class=="list"){
mapBiasList <- list()
for (i in 1:nrow(x)) {
mat <- assays(x)[["mapBias"]][i, , ]
if (class(mat)[1] == "numeric") {
dim(mat) <- c(1, 4)
rownames(mat) <- colnames(x)
}
colnames(mat) <- x@variants
mapBiasList[[i]] <- mat
}
# add snp id
names(mapBiasList) <- rownames(x)
return(mapBiasList)
}
})
#' @rdname ASEset-class
#' @export
setGeneric("fraction", function(x, ...) {
standardGeneric("fraction")
})
#' @rdname ASEset-class
#' @export
setMethod("fraction", signature(x = "ASEset"), function(x, strand = "*",
top.fraction.criteria="maxcount", verbose = FALSE, ...) {
if (!sum(strand %in% c("+", "-", "*")) > 0) {
stop("strand parameter has to be either '+', '-', '*' ")
}
#calculate frequency
fr <- frequency(x, strand=strand, return.class="array", ...)
#check and use top.fraction.criteria=="phase"
if(top.fraction.criteria=="phase"){
if(!"phase" %in% assayNames(x)){
stop("the phase slot has not been initialized")
}
if(is.null(assays(x)[["phase"]])){
stop("the phase slot cannot be empty if 'top.fraction.criteria=\"phase\"'")
}
#take out reference allele information
ref <- ref(x)
#select only ref rows (here is mcols() ref required)
ar <- .arrayFromAlleleVector(x@variants, ref,nc=ncol(x))
#use the array to subset the ref frequencies
ret <- .subsetArrayToMatrix(fr, ar)
#use maternal phase freq, by reverse values in mat that are not ref (0)
ret <- .returnMaternalPhaseFrequency(phase(x,return.class="array")[,,1], ret)
}else if(top.fraction.criteria=="ref"){
#take out reference allele information
ref <- ref(x)
#select only ref rows (here is mcols() ref required)
ar <- .arrayFromAlleleVector(x@variants, ref,nc=ncol(x))
#use the array to subset the ref frequencies
ret <- .subsetArrayToMatrix(fr, ar)
}else if(top.fraction.criteria=="maxcount"){
#use output from rank as maxcount, select only 1st rank
ar <- .arrayFromAlleleVector(x@variants,
arank(x, strand = strand, return.class="matrix")[,1], ncol(x))
#use the array to subset the ref frequencies
ret <- .subsetArrayToMatrix(fr, ar)
}
# return object matrix
colnames(ret) <- colnames(x)
rownames(ret) <- rownames(x)
ret
})
### -------------------------------------------------------------------------
### helpers for fraction
###
#ph is an array, rf is the fraction matrix for ref frequencies
.returnMaternalPhaseFrequency <- function(ph,rf){
rf[ph==1] <- 1 - rf[ph==1]
rf
}
#' @rdname ASEset-class
#' @export
setGeneric("arank", function(x, return.type = "names",
return.class = "list", strand = "*", ...) {
standardGeneric("arank")
})
setMethod("arank", signature(x = "ASEset"), function(x, return.type = "names",
return.class = "list", strand = "*", ...) {
if(return.class=="matrix"){
if (return.type == "names") {
ar <- t(apply(apply(alleleCounts(x,
strand=strand,return.class="array"),c(1,3),sum),
1, function(x){rank(x,ties.method="first")}))
ar2 <- t(apply(ar,1,function(x){
x <- sort(x,index.return=TRUE,decreasing=TRUE)$ix
x
}))
mat <- matrix(x@variants[ar2],ncol=4, nrow=nrow(x), byrow=FALSE,
dimnames=list(dimnames(ar)[[1]],c(1,2,3,4)))
return(mat)
}else if (return.type == "rank") {
return(t(apply(apply(alleleCounts(x,
strand=strand,return.class="array"),c(1,3),sum),
1, function(x){rank(x,ties.method="first")})))
}else if (return.type == "counts") {
return(apply(alleleCounts(x,strand=strand,return.class="array"),c(1,3),sum))
}else{stop("return.type is not valid")}
}else if(return.class=="list"){
acounts <- alleleCounts(x, strand = strand)
if (return.type == "names") {
lapply(acounts, function(x) {
names(sort(apply(x, 2, sum), decreasing = TRUE))
})
}else if (return.type == "counts") {
lapply(acounts, function(x) {
as.numeric(sort(apply(x, 2, sum), decreasing = TRUE))
})
}else if (return.type == "rank") {
stop("rank for return.class 'list' is not implemented yet")
}else{stop("return.type is not valid")}
}
})
#setGeneric("table")
#setGeneric("table", function(x, strand = "*", sortBy="none", ...) {
# standardGeneric("table")
#})
# @rdname ASEset-class
#setMethod("table", signature(... = "ASEset"), function(...) {
#
# args <- list(...)
# if (length(args) > 1)
# stop("Only one argument in '...' supported")
# x <- args[[1L]]
#
# #because the generis of table is rubbish we have to return a list for each strand
# #retList <- list()
#
# #for(strand in c("+","-","*")){
# # df <- data.frame(row.names=rownames(x))
# # df[,c("chromosome","position")] <- c(as.character(seqnames(x)),start(x))
# # df <- cbind(df,as.data.frame(arank(x, return.type="counts",
# # return.class="matrix",strand=strand)))
#
# # #if only one sample add fraction to table()
# # if(ncol(x)==1){
# # df[,"fraction"] <- as.vector(fraction(x,strand=strand))
# # }
#
# # retList[[strand]] <- df
#
# }
# #return(SimpleList(retList))
#
# df()
#
#})
#' @rdname ASEset-class
#' @export
setGeneric("frequency")
setMethod("frequency", signature(x = "ASEset"), function(x,
return.class = "list", strand = "*",
threshold.count.sample = 1) {
#threshold.count.sample cannot be zero (further down f=1/0 woulf fail)
if(threshold.count.sample<1){stop("threshold.count.sample needs to be >1")}
#calculate frequency from Allele counts
ar2 <- .calcFrequencyFromAlleleCounts(alleleCounts(x, strand=strand, return.class="array"), threshold.count.sample)
if(return.class=="array"){
return(ar2)
}else if(return.class=="list"){
return(.array2MatrixList(ar2))
}else{
stop("return.class has to be 'array' or 'list'")
}
})
### -------------------------------------------------------------------------
### helpers for frequemcy
###
.calcFrequencyFromAlleleCounts <- function(ar, th){
actot <- apply(ar, c(1,2), sum, na.rm=TRUE)
actot[actot < th] <- NaN
ar * array(as.vector(1/actot),dim=dim(ar))
}
#' @rdname ASEset-class
#' @export
setGeneric("genotype", function(x, ...){
standardGeneric("genotype")
})
#' @rdname ASEset-class
#' @export
setMethod("genotype", signature(x = "ASEset"), function(x,
return.class="matrix"){
if (!("phase" %in% assayNames(x))) {
stop(paste("phase is not present in assays in",
" ASEset object, see '?inferGenotypes' "))
}
if(!("ref" %in% names(mcols(x)))){
stop(paste("ref allele is not present in mcols in",
" ASEset object, see '?ASEset' "))
}
if(!("alt" %in% names(mcols(x)))){
stop(paste("alt allele is not present in mcols in",
" ASEset object, see '?inferAltAllele' "))
}
phase2genotype(phase(x, return.class="array"), ref(x), alt(x), return.class=return.class)
})
#' @rdname ASEset-class
#' @export
setGeneric("genotype<-", function(x, value){
standardGeneric("genotype<-")
})
#' @rdname ASEset-class
#' @export
setMethod("genotype<-", signature(x = "ASEset"), function(x, value){
#check dimensions
if(!nrow(x)==nrow(value)){
stop("nrow(x) is not equal to nrow(value)")
}
if(!ncol(x)==ncol(value)){
stop("ncol(x) is not equal to ncol(value)")
}
#check prexence of ref
if(!"ref" %in% names(mcols(x))){
stop("ref() is not set in ASEset")
}
#It would be better to not allow withDimnames=FALSE
#I suspect that we need to change something in the init object to fix that.
assays(x, withDimnames=FALSE)[["phase"]] <- genotype2phase(value, ref(x))
x
})
#' @rdname ASEset-class
#' @export
setGeneric("countsPerSnp", function(x, ...){
standardGeneric("countsPerSnp")
})
#' @rdname ASEset-class
#' @export
setMethod("countsPerSnp", signature(x = "ASEset"), function(x,
return.class = "matrix", return.type="mean", strand = "*") {
if(return.class=="matrix"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum))
}else if(return.class=="vector"){
if(return.type=="all"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), 1, sum))
}else if(return.type=="mean"){
return(apply(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum), 1, mean))
}
}else{
stop("return.class has to be 'vector' or 'matrix'")
}
})
#' @rdname ASEset-class
#' @export
setGeneric("countsPerSample", function(x, ...){
standardGeneric("countsPerSample")
})
#' @rdname ASEset-class
#' @export
setMethod("countsPerSample", signature(x = "ASEset"), function(x,
return.class = "matrix", return.type="mean", strand = "*") {
if(return.class=="matrix"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum))
}else if(return.class=="vector"){
if(return.type=="all"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), 2, sum))
}else if(return.type=="mean"){
return(apply(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum), 2, mean))
}
}else{
stop("return.class has to be 'vector' or 'matrix'")
}
})
#' @rdname ASEset-class
#' @export
setGeneric("phase", function(x, ...){
standardGeneric("phase")
})
#' @rdname ASEset-class
#' @export
setMethod("phase", signature(x = "ASEset"), function(x,
return.class = "matrix" ) {
if(return.class=="matrix"){
mat <- phaseArray2phaseMatrix(assays(x)[["phase"]])
colnames(mat) <- colnames(x)
rownames(mat) <- rownames(x)
mat
}else if(return.class=="array"){
assays(x)[["phase"]]
}
})
#' @rdname ASEset-class
#' @export
setGeneric("phase<-", function(x, value){
standardGeneric("phase<-")
})
#' @rdname ASEset-class
#' @export
setMethod("phase<-", signature(x = "ASEset"), function(x,value) {
if(class(value)[1]=="matrix") {
if(!identical(dim(x),dim(value))){
stop("dimension of value does not correspond to the values of object ASEset")
}
assays(x, withDimnames=FALSE)[["phase"]] <- phaseMatrix2Array(value,dimnames=dimnames(x))
}else if(class(value)[1]=="array"){
assays(x, withDimnames=FALSE)[["phase"]] <- value
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("mapBias<-", function(x, value){
standardGeneric("mapBias<-")
})
#' @rdname ASEset-class
#' @export
setMethod("mapBias<-", signature(x = "ASEset"), function(x,value) {
if(class(value)[1]=="array") {
assays(x, withDimnames=FALSE)[["mapBias"]] <- value
}else {
stop("class has to be array")
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("refExist", function(x){
standardGeneric("refExist")
})
#' @rdname ASEset-class
#' @export
setMethod("refExist", signature(x = "ASEset"), function(x) {
!is.null(ref(x))
})
#' @rdname ASEset-class
#' @export
setGeneric("ref")
#' @rdname ASEset-class
#' @export
setMethod("ref", signature(x = "ASEset"), function(x) {
mcols(x)[["ref"]]
})
#' @rdname ASEset-class
#' @export
setGeneric("ref<-")
#' @rdname ASEset-class
#' @export
setMethod("ref<-", signature(x = "ASEset"), function(x, value) {
if(class(value)[1]=="character") {
mcols(x)[["ref"]] <- value
}else{
stop("wrong class")
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("altExist", function(x){
standardGeneric("altExist")
})
#' @rdname ASEset-class
#' @export
setMethod("altExist", signature(x = "ASEset"), function(x) {
!is.null(ref(x))
})
#' @rdname ASEset-class
#' @export
setGeneric("alt")
#' @rdname ASEset-class
#' @export
setMethod("alt", signature(x = "ASEset"), function(x) {
mcols(x)[["alt"]]
})
#' @rdname ASEset-class
#' @export
setGeneric("alt<-")
#' @rdname ASEset-class
#' @export
setMethod("alt<-", signature(x = "ASEset"), function(x, value) {
if(class(value)[1]=="character") {
mcols(x)[["alt"]] <- value
}else{
stop("wrong class")
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("aquals", function(x, ...){
standardGeneric("aquals")
})
#' @rdname ASEset-class
#' @export
setMethod("aquals", signature(x = "ASEset"), function(x) {
assays(x)[["aquals"]]
})
#' @rdname ASEset-class
#' @export
setGeneric("aquals<-", function(x, value){
standardGeneric("aquals<-")
})
#' @rdname ASEset-class
#' @export
setMethod("aquals<-", signature(x = "ASEset"), function(x,value) {
if(class(value)[1]=="array") {
if(!identical(dim(x),dim(value)[1:2])){
stop("dimension of value does not correspond to the values of object ASEset")
}
assays(x, withDimnames=FALSE)[["aquals"]] <- value
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("maternalAllele", function(x, ...){
standardGeneric("maternalAllele")
})
#' @rdname ASEset-class
#' @export
setMethod("maternalAllele", signature(x = "ASEset"),
function(x) {
mat <- phase(x,return.class="array")[,,1]
ref <- mcols(x)[["ref"]]
alt <- mcols(x)[["alt"]]
apply(t(mat),1,function(y, ref, alt){
vec <- rep(NA,length(y))
if(any(y == 1)){
vec[y == 1] <- alt[y == 1]
}
if(any(y == 0)){
vec[y == 0] <- ref[y == 0]
}
vec
}, ref=ref, alt=alt)
})
#' @rdname ASEset-class
#' @export
setGeneric("paternalAllele", function(x, ...){
standardGeneric("paternalAllele")
})
#' @rdname ASEset-class
#' @export
setMethod("paternalAllele", signature(x = "ASEset"),
function(x) {
mat <- phase(x,return.class="array")[,,2]
ref <- mcols(x)[["ref"]]
alt <- mcols(x)[["alt"]]
apply(t(mat),1,function(y, ref, alt){
vec <- rep(NA,length(y))
if(any(y == 1)){
vec[y == 1] <- alt[y == 1]
}
if(any(y == 0)){
vec[y == 0] <- ref[y == 0]
}
vec
}, ref=ref, alt=alt)
})
##' @rdname ASEset-class
##' @export
#setGeneric("testet<-", signature="x", function(x, ..., value){
# standardGeneric("testet<-")
#})
##we need to set a normal generic to be able to use the replacement generic
#setGeneric("testet", signature="x", function(x, ...){
# standardGeneric("testet")
#})
#
##' @rdname ASEset-class
##' @export
#setReplaceMethod("testet", c("ASEset"),
# function(x, ... ,value){
# message("value0: ", paste0(value, collapse=", "))
# .set_ASEset_ref(x, value)
# }
#)
#
#setMethod("testet","ASEset",
# function(x) ref(x)
#)
#
#.set_ASEset_ref <- function(x, value)
#{
# message("value1: ", paste0(value, collapse=", "))
# ## Use 'x@width[]' instead of 'x@width' so the right value is recycled.
# ref(x)[] <- value
# x
#}
#
#ref(x)
#testet(x) <- c("O","P","Q")
#testet(x)
#testet(x)[1:2] <- c("S","I")
#
#
##another test
#
#setGeneric("getLetters", function(x)standardGeneric("getLetters"))
#setGeneric("getLetters<-", function(x, value)standardGeneric("getLetters<-"))
#
#setMethod("getLetters", "ASEset", function(x) {
# x@variants
#})
#
#setReplaceMethod("getLetters", c("ASEset", "character"), function(x, value) {
# message("value: ", paste0(value, collapse=", "))
# x@variants <- value
# x
#})
#
## [1] "<0x3716b40>"
#getLetters(x)
## [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J"
#getLetters(x)[1:2] <- letters[1:2]
#
#
#
Try the AllelicImbalance package in your browser
Any scripts or data that you put into this service are public.
AllelicImbalance documentation built on Nov. 8, 2020, 6:52 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .calcFrequencyFromAlleleCounts .returnMaternalPhaseFrequency .unknownStrandCountsExists .minusStrandCountsExists .plusStrandCountsExists
#'@include initialize-methods.R
NULL
#' ASEset objects
#'
#' Object that holds allele counts, genomic positions and map-bias for a set of
#' SNPs
#'
#' An ASEset object differs from a regular RangedSummarizedExperiment object
#' in that the assays contains an array instead of matrix. This array has
#' ranges on the rows, sampleNames on the columns and variants in the third
#' dimension.
#'
#' It is possible to use the commands barplot and locationplot on an ASEset
#' object see more details in \code{\link{barplot}} and
#' \code{\link{locationplot}}.
#'
#' Three different alleleCount options are available. The simples one is the
#' * option, and is for experiments where the strand information is not
#' known e.g. non-stranded data. The unknown strand could also be for strand
#' specific data when the aligner could not find
#' any strand associated with the read, but this should normally not happen,
#' and if it does probably having an extremely low mapping quality.
#' Then there are an option too add plus
#' and minus stranded data. When using this, it is essential to make sure that
#' the RNA-seq experiment under analysis has in fact been created so that
#' correct strand information was obtained. The most functions will by default
#' have their strand argument set to '*'.
#'
#' The phase information is stored by the convention of
#' 'maternal chromosome|paternal chromosome', with 0 as reference allele and 1
#' as alternative allele. '|' when the phase is known and '/' when the phase is
#' unknown. Internally the information will be stored as an three dimensional
#' array, dim 1 for SNPs, dim 2 for Samples and dim 3 which is fixed and stores
#' maternal chromosome, paternal chromosome and phased (1 equals TRUE).
#'
#'
#' @name ASEset-class
#' @rdname ASEset-class
#' @aliases ASEset-class ASEset alleleCounts mapBias fraction arank
#' frequency genotype genotype<- phase phase<- alleleCounts,ASEset-method mapBias,ASEset-method
#' fraction,ASEset-method arank,ASEset-method
#' frequency,ASEset-method genotype,ASEset-method genotype<-,ASEset-method
#' alleleCounts<- alleleCounts<-,ASEset-method phase,ASEset-method phase<-,ASEset-method
#'
#' @docType class
#' @param x ASEset object
#' @param strand which strand of '+', '-' or '*'
#' @param verbose makes function more talkative
#' @param return.type return 'names', rank or 'counts'
#' @param return.class return 'list' or 'array'
#' @param top.fraction.criteria 'maxcount', 'ref' or 'phase'
#' @param value replacement variable
#' @param ... additional arguments
#' @return An object of class ASEset containing location information and allele
#' counts for a number of SNPs measured in a number of samples on various
#' strand, as well as mapBias information. All data is stored in a manner
#' similar to the SummarizedExperiment class.
#'
#' @section Table: table(...)
#'
#' \describe{
#' Arguments: \item{...}{An \code{ASEset object} that contains the
#' variants of interest}
#'
#' The generics for table does not easily allow more than one argument
#' so in respect to the different strand options, \code{table} will
#' return a SimpleList with length 3, one element for each strand.
#' }
#'
#' @section Frequency: frequency(x,
#' return.class = "list", strand = "*",
#' threshold.count.sample = 15)
#'
#' \describe{
#' Arguments: \item{x}{An \code{ASEset object} that contains the
#' variants of interest}
#'
#' \item{x}{threshold.count.samples} if sample has fewer counts the function
#' return NA.
#'
#' }
#' @section Constructor: ASEsetFromCountList(rowRanges, countListNonStranded =
#' NULL, countListPlus = NULL, countListMinus = NULL, countListUnknown = NULL,
#' colData = NULL, mapBiasExpMean = array(), verbose=FALSE, ...)
#'
#' \describe{
#'
#' Arguments: \item{rowRanges}{A \code{GenomicRanges object} that contains the
#' variants of interest} \item{countListNonStranded}{A \code{list} where each
#' entry is a matrix with allele counts as columns and sample counts as rows}
#' \item{countListPlus}{A \code{list} where each entry is a matrix with allele
#' counts as columns and sample counts as rows} \item{countListMinus}{A
#' \code{list} where each entry is a matrix with allele counts as columns and
#' sample counts as rows} \item{countListUnknown}{A \code{list} where each
#' entry is a matrix with allele counts as columns and sample counts as rows}
#' \item{colData}{A \code{DataFrame} object containing sample specific data}
#' \item{mapBiasExpMean}{A 3D \code{array} describing mapping bias. The SNPs
#' are in the 1st dimension, samples in the 2nd dimension and variants in the
#' 3rd dimension.} \item{verbose}{Makes function more talkative}
#' \item{...}{arguments passed on to SummarizedExperiment constructor} }
#'
#'
#' @author Jesper R. Gadin, Lasse Folkersen
#' @keywords class ASEset
#' @examples
#'
#'
#' #make example countList
#' set.seed(42)
#' countListPlus <- list()
#' snps <- c('snp1','snp2','snp3','snp4','snp5')
#' for(snp in snps){
#' count<-matrix(rep(0,16),ncol=4,dimnames=list(
#' c('sample1','sample2','sample3','sample4'),
#' c('A','T','G','C')))
#' #insert random counts in two of the alleles
#' for(allele in sample(c('A','T','G','C'),2)){
#' count[,allele]<-as.integer(rnorm(4,mean=50,sd=10))
#' }
#' countListPlus[[snp]] <- count
#' }
#'
#' #make example rowRanges
#' rowRanges <- GRanges(
#' seqnames = Rle(c('chr1', 'chr2', 'chr1', 'chr3', 'chr1')),
#' ranges = IRanges(1:5, width = 1, names = head(letters,5)),
#' snp = paste('snp',1:5,sep='')
#' )
#'
#' #make example colData
#' colData <- DataFrame(Treatment=c('ChIP', 'Input','Input','ChIP'),
#' row.names=c('ind1','ind2','ind3','ind4'))
#'
#' #make ASEset
#' a <- ASEsetFromCountList(rowRanges, countListPlus=countListPlus,
#' colData=colData)
#'
#'
#' #example phase matrix (simple form)
#' p1 <- matrix(sample(c(1,0),replace=TRUE, size=nrow(a)*ncol(a)),nrow=nrow(a), ncol(a))
#' p2 <- matrix(sample(c(1,0),replace=TRUE, size=nrow(a)*ncol(a)),nrow=nrow(a), ncol(a))
#' p <- matrix(paste(p1,sample(c("|","|","/"), size=nrow(a)*ncol(a), replace=TRUE), p2, sep=""),
#' nrow=nrow(a), ncol(a))
#'
#' phase(a) <- p
#'
#'
#' #generate ASEset from array
#' snps <- 999
#' samples <-5
#' ar <-array(rep(unlist(lapply(1:snps,
#' function(x){(sample(c(TRUE,FALSE,TRUE,FALSE), size = 4))})), samples),
#' dim=c(4,snps,samples))
#' ar2 <- array(sample(50:300, 4*snps*samples,replace=TRUE), dim=c(4,snps,samples))
#' ar2[ar] <- 0
#' ar2 <- aperm(ar2, c(2, 3, 1))
#' dimnames(ar2) <- list(paste("snp",1:snps,sep=""),paste("sample",1:samples,sep=""),
#' c("A","C","G","T"))
#' gr <- GRanges(seqnames=c("chr2"), ranges=IRanges(start=1:dim(ar2)[1], width=1), strand="*")
#' a <- ASEsetFromArrays(gr, countsUnknown=ar2)
#'
#'
#' @exportClass ASEset
#' @exportMethod alleleCounts alleleCounts<- mapBias fraction arank
#' frequency genotype genotype<- phase phase<-
#'
#' @export
setClass("ASEset", contains = "RangedSummarizedExperiment",
representation(variants = "vector"))
#' @rdname ASEset-class
#' @export
setGeneric("alleleCounts", function(x, strand = "*", return.class="list") {
standardGeneric("alleleCounts")
})
#' @rdname ASEset-class
#' @export
setMethod("alleleCounts", signature(x = "ASEset"), function(x, strand = "*",
return.class="list") {
#check that strand parameter is correct set
if (!sum(strand %in% c("+", "-", "*", "both")) > 0) {
stop("strand parameter has to be either '+', '-', '*' or 'both' ")
}
#check if the assays are present
if(strand=="+" ){
if(!"countsPlus" %in% assayNames(x)){
stop(paste("strand",strand,"is not present in ASEset object",sep=""))
}
}
if(strand=="-" ){
if(!"countsMinus" %in% assayNames(x)){
stop(paste("strand",strand,"is not present in ASEset object",sep=""))
}
}
if(strand=="*" ){
if(!("countsMinus" %in% assayNames(x) | "countsPlus" %in% assayNames(x))){
stop(paste("for strand '*' at least one of '+' or '-' is required to be",
" present in ASEset object. Old ASEsets set the '*' strand in",
" countsUnknown, but that has been deprecated. If so, move your counts",
" to the plus or minus slot instead ",
" assays(yourASEset, withDimnames=FALSE)[['countsPlus']] <- assays(yourASEset)[['countsUnknown']]",
sep=""))
}
}
if(strand=="both" ){
if(!"countsMinus" %in% assayNames(x) & "countsPlus" %in% assayNames(x)){
stop(paste("for strand 'both' both '+' and '-' is required to be",
" present in the ASEset object",sep=""))
}
}
#access the correct assays
if (strand == "+") {
ar <- assays(x)[["countsPlus"]]
} else if (strand == "-") {
ar <- assays(x)[["countsMinus"]]
} else if (strand == "*") {
if(!("countsMinus" %in% assayNames(x))){
ar <- assays(x)[["countsPlus"]]
}else if(!("countsPlus" %in% assayNames(x))){
ar <- assays(x)[["countsMinus"]]
}else{
ar <- assays(x)[["countsMinus"]] + assays(x)[["countsPlus"]]
}
} else if (strand == "both") {
ar <- array(c(assays(x)[["countsPlus"]], assays(x)[["countsMinus"]]), dim=c(dim(assays(x)[["countsPlus"]]), 2))
} else {
stop("not existing strand option")
}
if(return.class=="array"){
dimnames(ar)[[3]]<- x@variants
if (strand == "both") {
dimnames(ar)[[4]]<- c("+","-")
}
ar
}else if(return.class=="list"){
if(strand=="both"){
strands <- 2
alleleCountList2 <- list()
}else{
strands <- 1
}
for( j in 1:strands){
alleleCountList <- list()
for (i in 1:nrow(ar)) {
if(strand=="both"){
mat <- ar[i, , ,j]
}else{
mat <- ar[i, , ]
}
if (class(mat)[1] == "integer") {
mat <- t(as.matrix(mat))
}
if (class(mat)[1] == "numeric") {
mat <- t(mat)
}
colnames(mat) <- x@variants
rownames(mat) <- colnames(x)
alleleCountList[[i]] <- mat
}
# add snp id
names(alleleCountList) <- rownames(x)
if(strand=="both"){
alleleCountList2[[j]] <- alleleCountList
}
}
if(strand=="both"){
names(alleleCountList2) <- c("+","-")
alleleCountList2
}else{
alleleCountList
}
}else{
stop("return.class has to be 'list' or 'array'")
}
})
### -------------------------------------------------------------------------
### helpers for alleleCounts
###
.plusStrandCountsExists <- function(x){
"countsPlus" %in% assayNames(x)
}
.minusStrandCountsExists <- function(x){
"countsMinus" %in% assayNames(x)
}
.unknownStrandCountsExists <- function(x){
"countsUnknown" %in% assayNames(x)
}
#' @rdname ASEset-class
#' @export
setGeneric("alleleCounts<-", function(x, ..., value) {
standardGeneric("alleleCounts<-")
})
#' @rdname ASEset-class
#' @export
setReplaceMethod("alleleCounts", "ASEset",
function(x, strand="*", return.class="array", ..., value) {
#for dev
#message("valueX: ", paste0(value, collapse=", "))
if (!sum(strand %in% c("+", "-", "*")) > 0) {
stop("strand parameter has to be either '+', '-', '*' ")
}
if (strand == "+") {
el <- "countsPlus"
} else if (strand == "-") {
el <- "countsMinus"
} else if (strand == "*") {
el <- "countsUnknown"
} else {
stop("not existing strand option")
}
#check that value has the right dimensions (make validObject method)
#if(!all(dim(value)==dim(assays(x, withDimnames=FALSE)[["countsUnknown"]]))){
# stop("dimensions in replacement object is not correct")
#}
assays(x, withDimnames=FALSE)[[el]][] <- value
#return object
x
})
#' @rdname ASEset-class
#' @export
setGeneric("mapBias", function(x, ...) {
standardGeneric("mapBias")
})
#' @rdname ASEset-class
#' @export
setMethod("mapBias", signature(x = "ASEset"), function(x,
return.class="list") {
# assume alleleCount information is stored as element 1
if(return.class=="array"){
return(assays(x)[["mapBias"]])
}else if(return.class=="list"){
mapBiasList <- list()
for (i in 1:nrow(x)) {
mat <- assays(x)[["mapBias"]][i, , ]
if (class(mat)[1] == "numeric") {
dim(mat) <- c(1, 4)
rownames(mat) <- colnames(x)
}
colnames(mat) <- x@variants
mapBiasList[[i]] <- mat
}
# add snp id
names(mapBiasList) <- rownames(x)
return(mapBiasList)
}
})
#' @rdname ASEset-class
#' @export
setGeneric("fraction", function(x, ...) {
standardGeneric("fraction")
})
#' @rdname ASEset-class
#' @export
setMethod("fraction", signature(x = "ASEset"), function(x, strand = "*",
top.fraction.criteria="maxcount", verbose = FALSE, ...) {
if (!sum(strand %in% c("+", "-", "*")) > 0) {
stop("strand parameter has to be either '+', '-', '*' ")
}
#calculate frequency
fr <- frequency(x, strand=strand, return.class="array", ...)
#check and use top.fraction.criteria=="phase"
if(top.fraction.criteria=="phase"){
if(!"phase" %in% assayNames(x)){
stop("the phase slot has not been initialized")
}
if(is.null(assays(x)[["phase"]])){
stop("the phase slot cannot be empty if 'top.fraction.criteria=\"phase\"'")
}
#take out reference allele information
ref <- ref(x)
#select only ref rows (here is mcols() ref required)
ar <- .arrayFromAlleleVector(x@variants, ref,nc=ncol(x))
#use the array to subset the ref frequencies
ret <- .subsetArrayToMatrix(fr, ar)
#use maternal phase freq, by reverse values in mat that are not ref (0)
ret <- .returnMaternalPhaseFrequency(phase(x,return.class="array")[,,1], ret)
}else if(top.fraction.criteria=="ref"){
#take out reference allele information
ref <- ref(x)
#select only ref rows (here is mcols() ref required)
ar <- .arrayFromAlleleVector(x@variants, ref,nc=ncol(x))
#use the array to subset the ref frequencies
ret <- .subsetArrayToMatrix(fr, ar)
}else if(top.fraction.criteria=="maxcount"){
#use output from rank as maxcount, select only 1st rank
ar <- .arrayFromAlleleVector(x@variants,
arank(x, strand = strand, return.class="matrix")[,1], ncol(x))
#use the array to subset the ref frequencies
ret <- .subsetArrayToMatrix(fr, ar)
}
# return object matrix
colnames(ret) <- colnames(x)
rownames(ret) <- rownames(x)
ret
})
### -------------------------------------------------------------------------
### helpers for fraction
###
#ph is an array, rf is the fraction matrix for ref frequencies
.returnMaternalPhaseFrequency <- function(ph,rf){
rf[ph==1] <- 1 - rf[ph==1]
rf
}
#' @rdname ASEset-class
#' @export
setGeneric("arank", function(x, return.type = "names",
return.class = "list", strand = "*", ...) {
standardGeneric("arank")
})
setMethod("arank", signature(x = "ASEset"), function(x, return.type = "names",
return.class = "list", strand = "*", ...) {
if(return.class=="matrix"){
if (return.type == "names") {
ar <- t(apply(apply(alleleCounts(x,
strand=strand,return.class="array"),c(1,3),sum),
1, function(x){rank(x,ties.method="first")}))
ar2 <- t(apply(ar,1,function(x){
x <- sort(x,index.return=TRUE,decreasing=TRUE)$ix
x
}))
mat <- matrix(x@variants[ar2],ncol=4, nrow=nrow(x), byrow=FALSE,
dimnames=list(dimnames(ar)[[1]],c(1,2,3,4)))
return(mat)
}else if (return.type == "rank") {
return(t(apply(apply(alleleCounts(x,
strand=strand,return.class="array"),c(1,3),sum),
1, function(x){rank(x,ties.method="first")})))
}else if (return.type == "counts") {
return(apply(alleleCounts(x,strand=strand,return.class="array"),c(1,3),sum))
}else{stop("return.type is not valid")}
}else if(return.class=="list"){
acounts <- alleleCounts(x, strand = strand)
if (return.type == "names") {
lapply(acounts, function(x) {
names(sort(apply(x, 2, sum), decreasing = TRUE))
})
}else if (return.type == "counts") {
lapply(acounts, function(x) {
as.numeric(sort(apply(x, 2, sum), decreasing = TRUE))
})
}else if (return.type == "rank") {
stop("rank for return.class 'list' is not implemented yet")
}else{stop("return.type is not valid")}
}
})
#setGeneric("table")
#setGeneric("table", function(x, strand = "*", sortBy="none", ...) {
# standardGeneric("table")
#})
# @rdname ASEset-class
#setMethod("table", signature(... = "ASEset"), function(...) {
#
# args <- list(...)
# if (length(args) > 1)
# stop("Only one argument in '...' supported")
# x <- args[[1L]]
#
# #because the generis of table is rubbish we have to return a list for each strand
# #retList <- list()
#
# #for(strand in c("+","-","*")){
# # df <- data.frame(row.names=rownames(x))
# # df[,c("chromosome","position")] <- c(as.character(seqnames(x)),start(x))
# # df <- cbind(df,as.data.frame(arank(x, return.type="counts",
# # return.class="matrix",strand=strand)))
#
# # #if only one sample add fraction to table()
# # if(ncol(x)==1){
# # df[,"fraction"] <- as.vector(fraction(x,strand=strand))
# # }
#
# # retList[[strand]] <- df
#
# }
# #return(SimpleList(retList))
#
# df()
#
#})
#' @rdname ASEset-class
#' @export
setGeneric("frequency")
setMethod("frequency", signature(x = "ASEset"), function(x,
return.class = "list", strand = "*",
threshold.count.sample = 1) {
#threshold.count.sample cannot be zero (further down f=1/0 woulf fail)
if(threshold.count.sample<1){stop("threshold.count.sample needs to be >1")}
#calculate frequency from Allele counts
ar2 <- .calcFrequencyFromAlleleCounts(alleleCounts(x, strand=strand, return.class="array"), threshold.count.sample)
if(return.class=="array"){
return(ar2)
}else if(return.class=="list"){
return(.array2MatrixList(ar2))
}else{
stop("return.class has to be 'array' or 'list'")
}
})
### -------------------------------------------------------------------------
### helpers for frequemcy
###
.calcFrequencyFromAlleleCounts <- function(ar, th){
actot <- apply(ar, c(1,2), sum, na.rm=TRUE)
actot[actot < th] <- NaN
ar * array(as.vector(1/actot),dim=dim(ar))
}
#' @rdname ASEset-class
#' @export
setGeneric("genotype", function(x, ...){
standardGeneric("genotype")
})
#' @rdname ASEset-class
#' @export
setMethod("genotype", signature(x = "ASEset"), function(x,
return.class="matrix"){
if (!("phase" %in% assayNames(x))) {
stop(paste("phase is not present in assays in",
" ASEset object, see '?inferGenotypes' "))
}
if(!("ref" %in% names(mcols(x)))){
stop(paste("ref allele is not present in mcols in",
" ASEset object, see '?ASEset' "))
}
if(!("alt" %in% names(mcols(x)))){
stop(paste("alt allele is not present in mcols in",
" ASEset object, see '?inferAltAllele' "))
}
phase2genotype(phase(x, return.class="array"), ref(x), alt(x), return.class=return.class)
})
#' @rdname ASEset-class
#' @export
setGeneric("genotype<-", function(x, value){
standardGeneric("genotype<-")
})
#' @rdname ASEset-class
#' @export
setMethod("genotype<-", signature(x = "ASEset"), function(x, value){
#check dimensions
if(!nrow(x)==nrow(value)){
stop("nrow(x) is not equal to nrow(value)")
}
if(!ncol(x)==ncol(value)){
stop("ncol(x) is not equal to ncol(value)")
}
#check prexence of ref
if(!"ref" %in% names(mcols(x))){
stop("ref() is not set in ASEset")
}
#It would be better to not allow withDimnames=FALSE
#I suspect that we need to change something in the init object to fix that.
assays(x, withDimnames=FALSE)[["phase"]] <- genotype2phase(value, ref(x))
x
})
#' @rdname ASEset-class
#' @export
setGeneric("countsPerSnp", function(x, ...){
standardGeneric("countsPerSnp")
})
#' @rdname ASEset-class
#' @export
setMethod("countsPerSnp", signature(x = "ASEset"), function(x,
return.class = "matrix", return.type="mean", strand = "*") {
if(return.class=="matrix"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum))
}else if(return.class=="vector"){
if(return.type=="all"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), 1, sum))
}else if(return.type=="mean"){
return(apply(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum), 1, mean))
}
}else{
stop("return.class has to be 'vector' or 'matrix'")
}
})
#' @rdname ASEset-class
#' @export
setGeneric("countsPerSample", function(x, ...){
standardGeneric("countsPerSample")
})
#' @rdname ASEset-class
#' @export
setMethod("countsPerSample", signature(x = "ASEset"), function(x,
return.class = "matrix", return.type="mean", strand = "*") {
if(return.class=="matrix"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum))
}else if(return.class=="vector"){
if(return.type=="all"){
return(apply(alleleCounts(x, strand=strand, return.class="array"), 2, sum))
}else if(return.type=="mean"){
return(apply(apply(alleleCounts(x, strand=strand, return.class="array"), c(1,2), sum), 2, mean))
}
}else{
stop("return.class has to be 'vector' or 'matrix'")
}
})
#' @rdname ASEset-class
#' @export
setGeneric("phase", function(x, ...){
standardGeneric("phase")
})
#' @rdname ASEset-class
#' @export
setMethod("phase", signature(x = "ASEset"), function(x,
return.class = "matrix" ) {
if(return.class=="matrix"){
mat <- phaseArray2phaseMatrix(assays(x)[["phase"]])
colnames(mat) <- colnames(x)
rownames(mat) <- rownames(x)
mat
}else if(return.class=="array"){
assays(x)[["phase"]]
}
})
#' @rdname ASEset-class
#' @export
setGeneric("phase<-", function(x, value){
standardGeneric("phase<-")
})
#' @rdname ASEset-class
#' @export
setMethod("phase<-", signature(x = "ASEset"), function(x,value) {
if(class(value)[1]=="matrix") {
if(!identical(dim(x),dim(value))){
stop("dimension of value does not correspond to the values of object ASEset")
}
assays(x, withDimnames=FALSE)[["phase"]] <- phaseMatrix2Array(value,dimnames=dimnames(x))
}else if(class(value)[1]=="array"){
assays(x, withDimnames=FALSE)[["phase"]] <- value
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("mapBias<-", function(x, value){
standardGeneric("mapBias<-")
})
#' @rdname ASEset-class
#' @export
setMethod("mapBias<-", signature(x = "ASEset"), function(x,value) {
if(class(value)[1]=="array") {
assays(x, withDimnames=FALSE)[["mapBias"]] <- value
}else {
stop("class has to be array")
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("refExist", function(x){
standardGeneric("refExist")
})
#' @rdname ASEset-class
#' @export
setMethod("refExist", signature(x = "ASEset"), function(x) {
!is.null(ref(x))
})
#' @rdname ASEset-class
#' @export
setGeneric("ref")
#' @rdname ASEset-class
#' @export
setMethod("ref", signature(x = "ASEset"), function(x) {
mcols(x)[["ref"]]
})
#' @rdname ASEset-class
#' @export
setGeneric("ref<-")
#' @rdname ASEset-class
#' @export
setMethod("ref<-", signature(x = "ASEset"), function(x, value) {
if(class(value)[1]=="character") {
mcols(x)[["ref"]] <- value
}else{
stop("wrong class")
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("altExist", function(x){
standardGeneric("altExist")
})
#' @rdname ASEset-class
#' @export
setMethod("altExist", signature(x = "ASEset"), function(x) {
!is.null(ref(x))
})
#' @rdname ASEset-class
#' @export
setGeneric("alt")
#' @rdname ASEset-class
#' @export
setMethod("alt", signature(x = "ASEset"), function(x) {
mcols(x)[["alt"]]
})
#' @rdname ASEset-class
#' @export
setGeneric("alt<-")
#' @rdname ASEset-class
#' @export
setMethod("alt<-", signature(x = "ASEset"), function(x, value) {
if(class(value)[1]=="character") {
mcols(x)[["alt"]] <- value
}else{
stop("wrong class")
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("aquals", function(x, ...){
standardGeneric("aquals")
})
#' @rdname ASEset-class
#' @export
setMethod("aquals", signature(x = "ASEset"), function(x) {
assays(x)[["aquals"]]
})
#' @rdname ASEset-class
#' @export
setGeneric("aquals<-", function(x, value){
standardGeneric("aquals<-")
})
#' @rdname ASEset-class
#' @export
setMethod("aquals<-", signature(x = "ASEset"), function(x,value) {
if(class(value)[1]=="array") {
if(!identical(dim(x),dim(value)[1:2])){
stop("dimension of value does not correspond to the values of object ASEset")
}
assays(x, withDimnames=FALSE)[["aquals"]] <- value
}
x
})
#' @rdname ASEset-class
#' @export
setGeneric("maternalAllele", function(x, ...){
standardGeneric("maternalAllele")
})
#' @rdname ASEset-class
#' @export
setMethod("maternalAllele", signature(x = "ASEset"),
function(x) {
mat <- phase(x,return.class="array")[,,1]
ref <- mcols(x)[["ref"]]
alt <- mcols(x)[["alt"]]
apply(t(mat),1,function(y, ref, alt){
vec <- rep(NA,length(y))
if(any(y == 1)){
vec[y == 1] <- alt[y == 1]
}
if(any(y == 0)){
vec[y == 0] <- ref[y == 0]
}
vec
}, ref=ref, alt=alt)
})
#' @rdname ASEset-class
#' @export
setGeneric("paternalAllele", function(x, ...){
standardGeneric("paternalAllele")
})
#' @rdname ASEset-class
#' @export
setMethod("paternalAllele", signature(x = "ASEset"),
function(x) {
mat <- phase(x,return.class="array")[,,2]
ref <- mcols(x)[["ref"]]
alt <- mcols(x)[["alt"]]
apply(t(mat),1,function(y, ref, alt){
vec <- rep(NA,length(y))
if(any(y == 1)){
vec[y == 1] <- alt[y == 1]
}
if(any(y == 0)){
vec[y == 0] <- ref[y == 0]
}
vec
}, ref=ref, alt=alt)
})
##' @rdname ASEset-class
##' @export
#setGeneric("testet<-", signature="x", function(x, ..., value){
# standardGeneric("testet<-")
#})
##we need to set a normal generic to be able to use the replacement generic
#setGeneric("testet", signature="x", function(x, ...){
# standardGeneric("testet")
#})
#
##' @rdname ASEset-class
##' @export
#setReplaceMethod("testet", c("ASEset"),
# function(x, ... ,value){
# message("value0: ", paste0(value, collapse=", "))
# .set_ASEset_ref(x, value)
# }
#)
#
#setMethod("testet","ASEset",
# function(x) ref(x)
#)
#
#.set_ASEset_ref <- function(x, value)
#{
# message("value1: ", paste0(value, collapse=", "))
# ## Use 'x@width[]' instead of 'x@width' so the right value is recycled.
# ref(x)[] <- value
# x
#}
#
#ref(x)
#testet(x) <- c("O","P","Q")
#testet(x)
#testet(x)[1:2] <- c("S","I")
#
#
##another test
#
#setGeneric("getLetters", function(x)standardGeneric("getLetters"))
#setGeneric("getLetters<-", function(x, value)standardGeneric("getLetters<-"))
#
#setMethod("getLetters", "ASEset", function(x) {
# x@variants
#})
#
#setReplaceMethod("getLetters", c("ASEset", "character"), function(x, value) {
# message("value: ", paste0(value, collapse=", "))
# x@variants <- value
# x
#})
#
## [1] "<0x3716b40>"
#getLetters(x)
## [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J"
#getLetters(x)[1:2] <- letters[1:2]
#
#
#
Try the AllelicImbalance package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.