Nothing
R/importData.R
In methimpute: Imputation-guided re-construction of complete methylomes from WGBS data
Defines functions
importRene
importBismark
importBSSeeker
importMethylpy
importBSMAP
Documented in
importBismark
importBSMAP
importBSSeeker
importMethylpy
importRene
#' Methimpute data import
#'
#' This page provides an overview of all \pkg{\link{methimpute}} data import functions.
#' @param file The file to import.
#' @param chrom.lengths A data.frame with chromosome names in the first, and chromosome lengths in the second column. Only chromosomes named in here will be returned. Alternatively a tab-separated file with such a data.frame (with headers).
#' @param skip The number of lines to skip. Usually 1 if the file contains a header and 0 otherwise.
#' @param contexts A character vector of the contexts that are to be assigned. Since some programs report 5-letter contexts, this parameter can be used to obtain a reduced number of contexts. Will yield contexts CG, CHG, CHH by default. Set \code{contexts=NULL} to obtain all available contexts.
#' @return A \code{\link{methimputeData}} object.
#' @name import
#' @examples
#'## Get an example file in BSSeeker format
#'file <- system.file("extdata","arabidopsis_bsseeker.txt.gz", package="methimpute")
#'data(arabidopsis_chromosomes)
#'bsseeker.data <- importBSSeeker(file, chrom.lengths=arabidopsis_chromosomes)
#'
#'## Get an example file in Bismark format
#'file <- system.file("extdata","arabidopsis_bismark.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'arabidopsis_chromosomes$chromosome <- sub('chr', '', arabidopsis_chromosomes$chromosome)
#'bismark.data <- importBismark(file, chrom.lengths=arabidopsis_chromosomes)
#'
#'## Get an example file in BSMAP format
#'file <- system.file("extdata","arabidopsis_BSMAP.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'bsmap.data <- importBSMAP(file, chrom.lengths=arabidopsis_chromosomes)
#'
#'## Get an example file in Methylpy format
#'file <- system.file("extdata","arabidopsis_methylpy.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'arabidopsis_chromosomes$chromosome <- sub('chr', '', arabidopsis_chromosomes$chromosome)
#'methylpy.data <- importMethylpy(file, chrom.lengths=arabidopsis_chromosomes)
NULL
#' @describeIn import Import a BSMAP methylation extractor file.
#' @importFrom Biostrings readDNAStringSet vcountPattern reverseComplement
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importBSMAP <- function(file, chrom.lengths=NULL, skip=1, contexts=c(CG='NNCGN', CHG='NNCHG', CHH='NNCHH')) {
## Import data
classes <- c('character', 'numeric', 'character', 'character', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=data.raw$V3, context=factor(NA, levels=names(contexts)), context.full=data.raw$V4)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V7
counts[,"total"] <- data.raw$V8
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
## Rework contexts
ptm <- startTimedMessage("Reworking contexts ...")
if (is.null(contexts)) {
data$context <- factor(data$context.full)
} else {
context.full <- Biostrings::DNAStringSet(data$context.full)
# Reverse complement minus
mask <- strand(data) == '-'
context.full[mask] <- Biostrings::reverseComplement(x = context.full[mask])
for (i1 in 1:length(contexts)) {
ind <- which(Biostrings::vcountPattern(contexts[i1], subject = context.full, fixed = FALSE) == 1)
data$context[ind] <- names(contexts)[i1]
}
}
data$context.full <- NULL
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
return(data)
}
#' @describeIn import Import a Methylpy methylation extractor file.
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importMethylpy <- function(file, chrom.lengths=NULL, skip=1, contexts=c(CG='CGN', CHG='CHG', CHH='CHH')) {
## Import data
classes <- c('character', 'numeric', 'character', 'character', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=data.raw$V3, context=factor(NA, levels=names(contexts)), context.full=data.raw$V4)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V5
counts[,"total"] <- data.raw$V6
data$counts <- counts
data$binomial.test <- data.raw$V7
rm(data.raw)
stopTimedMessage(ptm)
## Rework contexts
ptm <- startTimedMessage("Reworking contexts ...")
if (is.null(contexts)) {
data$context <- factor(data$context.full)
} else {
context.full <- Biostrings::DNAStringSet(data$context.full)
for (i1 in 1:length(contexts)) {
ind <- which(Biostrings::vcountPattern(contexts[i1], subject = context.full, fixed = FALSE) == 1)
data$context[ind] <- names(contexts)[i1]
}
}
data$context.full <- NULL
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
return(data)
}
#' @describeIn import Import a BSSeeker methylation extractor file.
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importBSSeeker <- function(file, chrom.lengths=NULL, skip=0) {
# classes <- c(seqnames='character', nucleotide='character', position='numeric', context='character', context.dinucleotide='character', methylation.level='numeric', counts.methylated='numeric', counts.total='numeric')
classes <- c('character', 'character', 'numeric', 'character', 'character', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V3, end=data.raw$V3), strand=c('C'='+', 'G'='-')[data.raw$V2], context=data.raw$V4)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V7
counts[,"total"] <- data.raw$V8
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
## Make factors
data$context <- factor(data$context)
return(data)
}
#' @describeIn import Import a Bismark methylation extractor file.
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importBismark <- function(file, chrom.lengths=NULL, skip=0) {
# classes <- c(seqnames='character', position='numeric', strand='character', counts.methylated='numeric', counts.total='numeric', context='character', context.trinucleotide='character')
classes <- c('character', 'numeric', 'character', 'numeric', 'numeric', 'character', 'character')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=data.raw$V3, context=data.raw$V6)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V4
counts[,"total"] <- data.raw$V4 + data.raw$V5
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
## Make factors
data$context <- factor(data$context)
return(data)
}
#' Import a Rene methylation extractor file
#'
#' Import a Rene methylation extractor file into a \code{\link[GenomicRanges]{GRanges-class}} object.
#'
#' @param file The file to import.
#' @param chrom.lengths A data.frame with chromosome names in the first, and chromosome lengths in the second column. Only chromosomes named in here will be returned. Alternatively a tab-separated file with such a data.frame (with headers).
#' @param skip The number of lines to skip. Usually 1 if the file contains a header and 0 otherwise.
#' @return A \code{\link{methimputeData}} object.
#'
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @examples
#'## Get an example file in Rene format
#'file <- system.file("extdata","arabidopsis_rene.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'rene.data <- methimpute:::importRene(file, chrom.lengths=arabidopsis_chromosomes)
#'
importRene <- function(file, chrom.lengths=NULL, skip=1) {
classes <- c('character', 'numeric', 'character', 'NULL', 'character', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=c('F'='+', 'R'='-')[data.raw$V3], methylated=data.raw$V10, context=data.raw$V5)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V6
counts[,"total"] <- data.raw$V7
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
seqlevels(data) <- paste0('chr', seqlevels(data))
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
## Make factors
data$context <- factor(data$context)
return(data)
}
#' #' Import a Bismark methylation extractor file
#' #'
#' #' Import a Bismark methylation extractor file into a \code{\link[GenomicRanges]{GRanges-class}} object.
#' #'
#' #' @param files The files to import.
#' #' @param chrom.lengths A data.frame with chromosome names in the first, and chromosome lengths in the second column. Only chromosomes named in here will be returned. Alternatively a tab-separated file with such a data.frame (with headers).
#' #' @param temp.store A folder where to save temporary files on disk. Set to \code{NULL} if no temporary files should be created.
#' #' @return A \code{\link{methimputeData}} object.
#' #'
#' #' @importFrom data.table fread
#' #' @importFrom utils read.table
#' #' @export
#' importBismark <- function(files, chrom.lengths=NULL, temp.store=tempfile("importBismark")) {
#'
#' classes <- c('NULL', 'character', 'character', 'numeric', 'character')
#' datas <- GRangesList()
#'
#' if (!is.null(temp.store)) {
#' ## Create folder for temporary files
#' if (!file.exists(temp.store)) {
#' dir.create(temp.store)
#' }
#' }
#'
#' ### Do each file at a time to avoid memory overflow ###
#' data.dedups <- GRangesList()
#' for (i1 in 1:length(files)) {
#' file <- files[i1]
#' savename <- file.path(temp.store, paste0(basename(file), '.RData'))
#' if (!file.exists(savename)) {
#' ptm <- startTimedMessage("Reading file ", file, " ...")
#' data.raw <- tryCatch({
#' data.raw <- data.table::fread(file, skip=1, sep='\t', colClasses=classes)
#' }, error = function(err) {
#' data.raw <- utils::read.table(file, skip=1, sep='\t', colClasses=classes)
#' })
#' data <- GRanges(seqnames=data.raw$V3, ranges=IRanges(start=data.raw$V4, end=data.raw$V4), strand="*", methylated=factor(data.raw$V2, levels=c("-", "+")), context=data.raw$V5)
#' rm(data.raw)
#' stopTimedMessage(ptm)
#' if (!is.null(temp.store)) {
#' ptm <- startTimedMessage("Saving to file ", savename, " ...")
#' save(data, file=savename)
#' stopTimedMessage(ptm)
#' }
#' } else {
#' ptm <- startTimedMessage("Loading file ", savename, " ...")
#' temp.env <- new.env()
#' data <- get(load(savename, envir=temp.env), envir=temp.env)
#' stopTimedMessage(ptm)
#' }
#'
#' ## Sorting
#' ptm <- startTimedMessage("Sorting ...")
#' data <- sort(data)
#' stopTimedMessage(ptm)
#'
#' ## Get counts at each position
#' ptm <- startTimedMessage("Getting counts ...")
#' data$methylated <- as.logical(as.integer(data$methylated)-1)
#' data$mask.dup <- c(FALSE, as.logical(data@seqnames[-1] == data@seqnames[-length(data)])) & c(FALSE, as.logical(data@ranges@start[-1] == data@ranges@start[-length(data)]))
#' data$unmeth.sum <- rev(cumsum(rev(!data$methylated)))
#' data$meth.sum <- rev(cumsum(rev(data$methylated)))
#' data.dedup <- data[!data$mask.dup]
#' rm(data)
#' data.dedup$counts.unmethylated <- -diff(c(data.dedup$unmeth.sum, 0))
#' data.dedup$counts.methylated <- -diff(c(data.dedup$meth.sum, 0))
#' data.dedup$methylated <- NULL
#' data.dedup$mask.dup <- NULL
#' data.dedup$unmeth.sum <- NULL
#' data.dedup$meth.sum <- NULL
#' data.dedups[[length(data.dedups)+1]] <- data.dedup
#' stopTimedMessage(ptm)
#'
#' ## Merging counts from data objects
#' if (length(data.dedups) == 2) {
#' ptm <- startTimedMessage("Merging counts ...")
#' data <- unlist(data.dedups, use.names=FALSE)
#' data <- sort(data)
#' rm(data.dedups)
#' data$mask.dup <- c(FALSE, as.logical(data@seqnames[-1] == data@seqnames[-length(data)])) & c(FALSE, as.logical(data@ranges@start[-1] == data@ranges@start[-length(data)]))
#' data$unmeth.sum <- rev(cumsum(rev(data$counts.unmethylated)))
#' data$meth.sum <- rev(cumsum(rev(data$counts.methylated)))
#' data.dedup <- data[!data$mask.dup]
#' rm(data)
#' data.dedup$counts.unmethylated <- -diff(c(data.dedup$unmeth.sum, 0))
#' data.dedup$counts.methylated <- -diff(c(data.dedup$meth.sum, 0))
#' data.dedup$mask.dup <- NULL
#' data.dedup$unmeth.sum <- NULL
#' data.dedup$meth.sum <- NULL
#' data.dedups <- GRangesList(data.dedup)
#' stopTimedMessage(ptm)
#' }
#'
#' }
#' data <- data.dedups[[1]]
#' rm(data.dedups, data.dedup)
#'
#' ## Assign seqlengths
#' if (!is.null(chrom.lengths)) {
#' if (is.character(chrom.lengths)) {
#' df <- tryCatch({
#' df <- data.table::fread(chrom.lengths, header=TRUE)
#' }, error = function(err) {
#' df <- utils::read.table(chrom.lengths, header=TRUE)
#' })
#' } else if (is.data.frame(chrom.lengths)) {
#' df <- chrom.lengths
#' }
#' chrom.lengths <- df[,2]
#' names(chrom.lengths) <- df[,1]
#' # Filter by chromosomes supplied in chrom.lengths
#' data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
#' seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
#' }
#'
#' ## Make factors
#' data$context <- factor(data$context)
#'
#' return(data)
#' }
#'
Try the methimpute package in your browser
Any scripts or data that you put into this service are public.
methimpute documentation built on Nov. 8, 2020, 5:47 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 importRene importBismark importBSSeeker importMethylpy importBSMAP
Documented in importBismark importBSMAP importBSSeeker importMethylpy importRene
#' Methimpute data import
#'
#' This page provides an overview of all \pkg{\link{methimpute}} data import functions.
#' @param file The file to import.
#' @param chrom.lengths A data.frame with chromosome names in the first, and chromosome lengths in the second column. Only chromosomes named in here will be returned. Alternatively a tab-separated file with such a data.frame (with headers).
#' @param skip The number of lines to skip. Usually 1 if the file contains a header and 0 otherwise.
#' @param contexts A character vector of the contexts that are to be assigned. Since some programs report 5-letter contexts, this parameter can be used to obtain a reduced number of contexts. Will yield contexts CG, CHG, CHH by default. Set \code{contexts=NULL} to obtain all available contexts.
#' @return A \code{\link{methimputeData}} object.
#' @name import
#' @examples
#'## Get an example file in BSSeeker format
#'file <- system.file("extdata","arabidopsis_bsseeker.txt.gz", package="methimpute")
#'data(arabidopsis_chromosomes)
#'bsseeker.data <- importBSSeeker(file, chrom.lengths=arabidopsis_chromosomes)
#'
#'## Get an example file in Bismark format
#'file <- system.file("extdata","arabidopsis_bismark.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'arabidopsis_chromosomes$chromosome <- sub('chr', '', arabidopsis_chromosomes$chromosome)
#'bismark.data <- importBismark(file, chrom.lengths=arabidopsis_chromosomes)
#'
#'## Get an example file in BSMAP format
#'file <- system.file("extdata","arabidopsis_BSMAP.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'bsmap.data <- importBSMAP(file, chrom.lengths=arabidopsis_chromosomes)
#'
#'## Get an example file in Methylpy format
#'file <- system.file("extdata","arabidopsis_methylpy.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'arabidopsis_chromosomes$chromosome <- sub('chr', '', arabidopsis_chromosomes$chromosome)
#'methylpy.data <- importMethylpy(file, chrom.lengths=arabidopsis_chromosomes)
NULL
#' @describeIn import Import a BSMAP methylation extractor file.
#' @importFrom Biostrings readDNAStringSet vcountPattern reverseComplement
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importBSMAP <- function(file, chrom.lengths=NULL, skip=1, contexts=c(CG='NNCGN', CHG='NNCHG', CHH='NNCHH')) {
## Import data
classes <- c('character', 'numeric', 'character', 'character', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=data.raw$V3, context=factor(NA, levels=names(contexts)), context.full=data.raw$V4)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V7
counts[,"total"] <- data.raw$V8
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
## Rework contexts
ptm <- startTimedMessage("Reworking contexts ...")
if (is.null(contexts)) {
data$context <- factor(data$context.full)
} else {
context.full <- Biostrings::DNAStringSet(data$context.full)
# Reverse complement minus
mask <- strand(data) == '-'
context.full[mask] <- Biostrings::reverseComplement(x = context.full[mask])
for (i1 in 1:length(contexts)) {
ind <- which(Biostrings::vcountPattern(contexts[i1], subject = context.full, fixed = FALSE) == 1)
data$context[ind] <- names(contexts)[i1]
}
}
data$context.full <- NULL
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
return(data)
}
#' @describeIn import Import a Methylpy methylation extractor file.
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importMethylpy <- function(file, chrom.lengths=NULL, skip=1, contexts=c(CG='CGN', CHG='CHG', CHH='CHH')) {
## Import data
classes <- c('character', 'numeric', 'character', 'character', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=data.raw$V3, context=factor(NA, levels=names(contexts)), context.full=data.raw$V4)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V5
counts[,"total"] <- data.raw$V6
data$counts <- counts
data$binomial.test <- data.raw$V7
rm(data.raw)
stopTimedMessage(ptm)
## Rework contexts
ptm <- startTimedMessage("Reworking contexts ...")
if (is.null(contexts)) {
data$context <- factor(data$context.full)
} else {
context.full <- Biostrings::DNAStringSet(data$context.full)
for (i1 in 1:length(contexts)) {
ind <- which(Biostrings::vcountPattern(contexts[i1], subject = context.full, fixed = FALSE) == 1)
data$context[ind] <- names(contexts)[i1]
}
}
data$context.full <- NULL
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
return(data)
}
#' @describeIn import Import a BSSeeker methylation extractor file.
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importBSSeeker <- function(file, chrom.lengths=NULL, skip=0) {
# classes <- c(seqnames='character', nucleotide='character', position='numeric', context='character', context.dinucleotide='character', methylation.level='numeric', counts.methylated='numeric', counts.total='numeric')
classes <- c('character', 'character', 'numeric', 'character', 'character', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V3, end=data.raw$V3), strand=c('C'='+', 'G'='-')[data.raw$V2], context=data.raw$V4)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V7
counts[,"total"] <- data.raw$V8
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
## Make factors
data$context <- factor(data$context)
return(data)
}
#' @describeIn import Import a Bismark methylation extractor file.
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @export
#'
importBismark <- function(file, chrom.lengths=NULL, skip=0) {
# classes <- c(seqnames='character', position='numeric', strand='character', counts.methylated='numeric', counts.total='numeric', context='character', context.trinucleotide='character')
classes <- c('character', 'numeric', 'character', 'numeric', 'numeric', 'character', 'character')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=data.raw$V3, context=data.raw$V6)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V4
counts[,"total"] <- data.raw$V4 + data.raw$V5
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
## Make factors
data$context <- factor(data$context)
return(data)
}
#' Import a Rene methylation extractor file
#'
#' Import a Rene methylation extractor file into a \code{\link[GenomicRanges]{GRanges-class}} object.
#'
#' @param file The file to import.
#' @param chrom.lengths A data.frame with chromosome names in the first, and chromosome lengths in the second column. Only chromosomes named in here will be returned. Alternatively a tab-separated file with such a data.frame (with headers).
#' @param skip The number of lines to skip. Usually 1 if the file contains a header and 0 otherwise.
#' @return A \code{\link{methimputeData}} object.
#'
#' @importFrom data.table fread
#' @importFrom utils read.table
#' @examples
#'## Get an example file in Rene format
#'file <- system.file("extdata","arabidopsis_rene.txt", package="methimpute")
#'data(arabidopsis_chromosomes)
#'rene.data <- methimpute:::importRene(file, chrom.lengths=arabidopsis_chromosomes)
#'
importRene <- function(file, chrom.lengths=NULL, skip=1) {
classes <- c('character', 'numeric', 'character', 'NULL', 'character', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric')
ptm <- startTimedMessage("Reading file ", file, " ...")
data.raw <- tryCatch({
data.raw <- data.table::fread(file, skip=skip, sep='\t', colClasses=classes)
}, error = function(err) {
data.raw <- utils::read.table(file, skip=skip, sep='\t', colClasses=classes)
})
data <- GRanges(seqnames=data.raw$V1, ranges=IRanges(start=data.raw$V2, end=data.raw$V2), strand=c('F'='+', 'R'='-')[data.raw$V3], methylated=data.raw$V10, context=data.raw$V5)
counts <- array(NA, dim=c(length(data), 2), dimnames=list(NULL, c("methylated", "total")))
counts[,"methylated"] <- data.raw$V6
counts[,"total"] <- data.raw$V7
data$counts <- counts
rm(data.raw)
stopTimedMessage(ptm)
seqlevels(data) <- paste0('chr', seqlevels(data))
## Assign seqlengths
if (!is.null(chrom.lengths)) {
if (is.character(chrom.lengths)) {
df <- tryCatch({
df <- data.table::fread(chrom.lengths, header=TRUE)
}, error = function(err) {
df <- utils::read.table(chrom.lengths, header=TRUE)
})
} else if (is.data.frame(chrom.lengths)) {
df <- chrom.lengths
}
chrom.lengths <- df[,2]
names(chrom.lengths) <- df[,1]
# Filter by chromosomes supplied in chrom.lengths
data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
}
## Make factors
data$context <- factor(data$context)
return(data)
}
#' #' Import a Bismark methylation extractor file
#' #'
#' #' Import a Bismark methylation extractor file into a \code{\link[GenomicRanges]{GRanges-class}} object.
#' #'
#' #' @param files The files to import.
#' #' @param chrom.lengths A data.frame with chromosome names in the first, and chromosome lengths in the second column. Only chromosomes named in here will be returned. Alternatively a tab-separated file with such a data.frame (with headers).
#' #' @param temp.store A folder where to save temporary files on disk. Set to \code{NULL} if no temporary files should be created.
#' #' @return A \code{\link{methimputeData}} object.
#' #'
#' #' @importFrom data.table fread
#' #' @importFrom utils read.table
#' #' @export
#' importBismark <- function(files, chrom.lengths=NULL, temp.store=tempfile("importBismark")) {
#'
#' classes <- c('NULL', 'character', 'character', 'numeric', 'character')
#' datas <- GRangesList()
#'
#' if (!is.null(temp.store)) {
#' ## Create folder for temporary files
#' if (!file.exists(temp.store)) {
#' dir.create(temp.store)
#' }
#' }
#'
#' ### Do each file at a time to avoid memory overflow ###
#' data.dedups <- GRangesList()
#' for (i1 in 1:length(files)) {
#' file <- files[i1]
#' savename <- file.path(temp.store, paste0(basename(file), '.RData'))
#' if (!file.exists(savename)) {
#' ptm <- startTimedMessage("Reading file ", file, " ...")
#' data.raw <- tryCatch({
#' data.raw <- data.table::fread(file, skip=1, sep='\t', colClasses=classes)
#' }, error = function(err) {
#' data.raw <- utils::read.table(file, skip=1, sep='\t', colClasses=classes)
#' })
#' data <- GRanges(seqnames=data.raw$V3, ranges=IRanges(start=data.raw$V4, end=data.raw$V4), strand="*", methylated=factor(data.raw$V2, levels=c("-", "+")), context=data.raw$V5)
#' rm(data.raw)
#' stopTimedMessage(ptm)
#' if (!is.null(temp.store)) {
#' ptm <- startTimedMessage("Saving to file ", savename, " ...")
#' save(data, file=savename)
#' stopTimedMessage(ptm)
#' }
#' } else {
#' ptm <- startTimedMessage("Loading file ", savename, " ...")
#' temp.env <- new.env()
#' data <- get(load(savename, envir=temp.env), envir=temp.env)
#' stopTimedMessage(ptm)
#' }
#'
#' ## Sorting
#' ptm <- startTimedMessage("Sorting ...")
#' data <- sort(data)
#' stopTimedMessage(ptm)
#'
#' ## Get counts at each position
#' ptm <- startTimedMessage("Getting counts ...")
#' data$methylated <- as.logical(as.integer(data$methylated)-1)
#' data$mask.dup <- c(FALSE, as.logical(data@seqnames[-1] == data@seqnames[-length(data)])) & c(FALSE, as.logical(data@ranges@start[-1] == data@ranges@start[-length(data)]))
#' data$unmeth.sum <- rev(cumsum(rev(!data$methylated)))
#' data$meth.sum <- rev(cumsum(rev(data$methylated)))
#' data.dedup <- data[!data$mask.dup]
#' rm(data)
#' data.dedup$counts.unmethylated <- -diff(c(data.dedup$unmeth.sum, 0))
#' data.dedup$counts.methylated <- -diff(c(data.dedup$meth.sum, 0))
#' data.dedup$methylated <- NULL
#' data.dedup$mask.dup <- NULL
#' data.dedup$unmeth.sum <- NULL
#' data.dedup$meth.sum <- NULL
#' data.dedups[[length(data.dedups)+1]] <- data.dedup
#' stopTimedMessage(ptm)
#'
#' ## Merging counts from data objects
#' if (length(data.dedups) == 2) {
#' ptm <- startTimedMessage("Merging counts ...")
#' data <- unlist(data.dedups, use.names=FALSE)
#' data <- sort(data)
#' rm(data.dedups)
#' data$mask.dup <- c(FALSE, as.logical(data@seqnames[-1] == data@seqnames[-length(data)])) & c(FALSE, as.logical(data@ranges@start[-1] == data@ranges@start[-length(data)]))
#' data$unmeth.sum <- rev(cumsum(rev(data$counts.unmethylated)))
#' data$meth.sum <- rev(cumsum(rev(data$counts.methylated)))
#' data.dedup <- data[!data$mask.dup]
#' rm(data)
#' data.dedup$counts.unmethylated <- -diff(c(data.dedup$unmeth.sum, 0))
#' data.dedup$counts.methylated <- -diff(c(data.dedup$meth.sum, 0))
#' data.dedup$mask.dup <- NULL
#' data.dedup$unmeth.sum <- NULL
#' data.dedup$meth.sum <- NULL
#' data.dedups <- GRangesList(data.dedup)
#' stopTimedMessage(ptm)
#' }
#'
#' }
#' data <- data.dedups[[1]]
#' rm(data.dedups, data.dedup)
#'
#' ## Assign seqlengths
#' if (!is.null(chrom.lengths)) {
#' if (is.character(chrom.lengths)) {
#' df <- tryCatch({
#' df <- data.table::fread(chrom.lengths, header=TRUE)
#' }, error = function(err) {
#' df <- utils::read.table(chrom.lengths, header=TRUE)
#' })
#' } else if (is.data.frame(chrom.lengths)) {
#' df <- chrom.lengths
#' }
#' chrom.lengths <- df[,2]
#' names(chrom.lengths) <- df[,1]
#' # Filter by chromosomes supplied in chrom.lengths
#' data <- keepSeqlevels(data, seqlevels(data)[seqlevels(data) %in% names(chrom.lengths)])
#' seqlengths(data) <- chrom.lengths[names(seqlengths(data))]
#' }
#'
#' ## Make factors
#' data$context <- factor(data$context)
#'
#' return(data)
#' }
#'
Try the methimpute 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.