Commit f5e1e18b authored by Eric Coissac's avatar Eric Coissac

debug the following

parent 6e719b68
This diff is collapsed.
......@@ -16,7 +16,8 @@
#' @export
#'
#' @examples
plot_read_x_motus <- function(metabar,mapping = NULL,
plot_read_x_hill <- function(metabar, q = 0,
mapping = NULL,
data = NULL, stat = "identity",
position = "identity", ...,
na.rm = FALSE,
......@@ -25,6 +26,7 @@ plot_read_x_motus <- function(metabar,mapping = NULL,
fcall <- match.call()
fcall$metabar=NULL
fcall$q = NULL
fcall[[1]] <- quote(geom_point)
gp <- eval(fcall)
......@@ -33,15 +35,150 @@ plot_read_x_motus <- function(metabar,mapping = NULL,
cats <- categories_colname(metabar)
left_join(metabar %>% samples_read_count(),
metabar %>% samples_motu_count(),
metabar %>% samples_hill(q = q),
by = sids) %>%
left_join(samples(metabar), by = sids) %>%
filter(read_count > 0 & hill > 0) %>%
ggplot(aes(x=read_count,
y=motu_count,
y=hill,
col=get(cats))) +
gp +
scale_x_log10() +
scale_y_log10() +
xlab("Read counts") +
ylab("MOTU counts")
}
\ No newline at end of file
ylab(glue::glue("Diversity (q={q})"))
}
#' @export
plot_read_x_pcr_tags <- function(metabar,
only,
library_colname,
only_forward_tags,
only_reverse_tags,
ncol_library,
mapping = NULL,
data = NULL, stat = "identity",
position = "identity", ...,
na.rm = FALSE,
show.legend = NA,
inherit.aes = TRUE) {
fcall <- match.call()
fcall$metabar=NULL
fcall$only = NULL
fcall$library_colname = NULL
fcall$only_forward_tags = NULL
fcall$only_reverse_tags = NULL
fcall$ncol_library = NULL
fcall[[1]] <- quote(geom_point)
forward_factor <- factor
reverse_factor <- factor
gp <- eval(fcall)
data <- if (missing(only))
metabar$data
else {
if (is_character(only))
only = tibble(motu = only)
else if (is_robimotu(only)) {
only = ids(only)
colnames(only) <- "motu"
}
else if (is.data.frame(only) && is_character(only[[1]]))
only = tibble(motu = only[[1]])
else
abort("I cannot guess the only format. Must be a character or a robimutu")
suppressWarnings(only %>% left_join(metabar$data, by = "motu"))
}
bys = c(id = samples_ids_colname(metabar))
selector = c(id = "id",
forward_tag = forward_tags_colname(metabar),
reverse_tag = reverse_tags_colname(metabar),
read_count = "read_count",
category = categories_colname(metabar))
if (!missing(library_colname))
selector = c(selector,library = library_colname)
data %>%
group_by(sample) %>% rename(id=sample) %>%
summarise(read_count=sum(count)) %>%
left_join(samples(metabar), by = bys) %>%
select(selector) -> data
if (!missing(only_forward_tags)) {
data <- data %>% filter(forward_tag %in% only_forward_tags)
forward_factor <- function(x) factor(x,levels=only_forward_tags)
}
if (!missing(only_reverse_tags)) {
data <- data %>% filter(reverse_tag %in% only_reverse_tags)
reverse_factor <- function(x) factor(x,levels=only_reverse_tags)
}
ggplot(data,aes(x = forward_factor(forward_tag),
y = reverse_factor(reverse_tag),
col = category,
cex = sqrt(read_count))) + gp -> gg
if (!missing(library_colname)) {
if (missing(ncol_library))
ncol_library <- floor(sqrt(length(unique(data$library))))
gg <- gg +
facet_wrap(~ library, ncol = ncol_library)
}
gg +
xlab("Forward tags") +
ylab("Reverse tags") +
theme(axis.text.x = element_text(size=rel(0.5),angle = 90,family="mono"),
axis.text.y = element_text(size=rel(0.5),family="mono"))
}
#' @export
plot_tag_jumps <- function(metabar,
maximum = FALSE,
mapping = NULL,
data = NULL, stat = "identity",
position = "identity", ...,
na.rm = FALSE,
show.legend = NA,
inherit.aes = TRUE) {
fcall <- match.call()
fcall$metabar = NULL
fcall$maximum = NULL
fcall[[1]] <- quote(geom_histogram)
gh <- eval(fcall)
metabar$data %>%
left_join(samples(metabar),
by = c(sample="id")) %>%
filter(category == "sequencing_blk") -> data
if (maximum)
data %>%
group_by(motu) %>%
summarise(count = max(count)) -> data
data %>% ggplot(aes(x=count)) + gh -> gg
if(maximum)
gg +
xlab("Maximum read count") +
ylab("MOTUs")
else
gg +
xlab("Read count") +
ylab("MOTU occurrences")
}
#' Title
#'
#' @param data
#'
#' @return
#' @export
#'
#' @examples
seek_contaminents <- function(data) {
UseMethod("seek_contaminents",data)
}
#' @rdname seek_contaminents
#' @export
seek_contaminents.robimetabar <- function(data) {
sid_name <- samples_ids_colname(data)
mid_name <- motus_ids_colname(data)
cat_name <- categories_colname(data)
suppressWarnings(
data %>%
motus_sample_max() %>%
left_join(motus(data) %>% select(mid_name,taxid),
by = mid_name) %>%
left_join(samples(data) %>% select(sid_name,cat_name),
by = c(sample_max=sid_name))) %>%
filter(category != "sample" & category != "positive_ctrl") %>%
arrange(category,desc(sample_count_max)) %>%
as_robimotu()
}
......@@ -22,7 +22,7 @@ ecotag_best_identity.robimotu <- function(object,taxonomy) {
cbind(ids(object)) %>%
gather(key = "db",value = "identity", -all_of(mids)) %>%
group_by_at(mids) %>%
mutate(rank__ = rank(identity)) %>%
mutate(rank__ = order(identity,decreasing = TRUE)) %>%
filter(rank__ <= 1) %>%
select(-rank__) %>% left_join(
object$best_match %>%
......
......@@ -254,7 +254,7 @@ left_parse_sample_ids.robimetabar <- function(object,parts,n,pattern="_",
replace = TRUE,
prefix="Left") {
call <- match.call()
call[[1]] <- quote(left_parse_sample_ids.robisample)
call[[1]] <- quote(left_parse_sample_ids.robidata)
call[[2]] <- quote(object$samples)
eval(call) -> i
i
......@@ -270,7 +270,7 @@ right_parse_sample_ids.robimetabar <- function(object,parts,n,pattern="_",
replace = TRUE,
prefix="Right") {
call <- match.call()
call[[1]] <- quote(right_parse_sample_ids.robisample)
call[[1]] <- quote(right_parse_sample_ids.robidata)
call[[2]] <- quote(object$samples)
eval(call) -> i
i
......
......@@ -26,6 +26,7 @@ read_metabar <- function(file,
obiclean_prefix = "obiclean",
obiclean_status = paste0(obiclean_prefix,"_status"),
obiclean_cluster = paste0(obiclean_prefix,"_cluster"),
obiclean_count = paste0(obiclean_prefix,"_count"),
sumaclust_prefix = "cluster",
main_taxid = "taxid",
other_taxid = c("family","genus","order","species","taxid_by_db"),
......@@ -38,7 +39,7 @@ read_metabar <- function(file,
merged_read_count <- sprintf("merged_%s",read_count)
features$sequence <- mfasta$sequence
features$id <- robiuniqueid(mfasta$id)
features$id <- as_robiuniqueid(mfasta$id)
data <- feature_as_three_columns(features[[merged_read_count]],
motus_ids = features$id) %>%
......@@ -69,12 +70,23 @@ read_metabar <- function(file,
cd <- colnames(data)
cd[cd == "value"] <- obiclean_cluster
colnames(data) <- cd
data[[obiclean_cluster]] = factor(data[[obiclean_cluster]],
levels = levels(features$id))
features[[obiclean_cluster]] <- NULL
}
if (obiclean_count %in% colnames(features)) {
data <- data %>%
full_join(features[[obiclean_count]] %>%
feature_as_three_columns(motus_ids = features$id),
by = c("motu", "sample"))
cd <- colnames(data)
cd[cd == "value"] <- obiclean_count
colnames(data) <- cd
features[[obiclean_count]] <- NULL
}
features %>% lapply(function(x) {
if (is(x,"list")) feature_as_matrix(x) %>% as_tibble()
else x}) %>% as_tibble() -> features
......@@ -98,7 +110,7 @@ read_metabar <- function(file,
}
if (main_taxid %in% colnames(features)) {
features[[main_taxid]] <- robitaxid_master(features[[main_taxid]])
features[[main_taxid]] <- as_robitaxid_master(features[[main_taxid]])
}
for (c in other_taxid)
......@@ -106,9 +118,9 @@ read_metabar <- function(file,
fs <- features[[c]]
if (is_tibble(fs))
for (j in seq_along(fs))
fs[[j]] <- robitaxid(fs[[j]])
fs[[j]] <- as_robitaxid(fs[[j]])
else
fs <- robitaxid(fs)
fs <- as_robitaxid(fs)
features[[c]] <- fs
}
......@@ -117,7 +129,7 @@ read_metabar <- function(file,
if (!is.null(ngsfilter)) {
m <- match(usampleids,ids(ngsfilter) %>% pull(),NA)
missing <- usampleids[which(is.na(m))]
robiassert_arg(length(missing)==0,
robiassert_arg(length(missing) == 0,
"ngsfilter",
"The provided ngsfilter doesn't define all the samples. ({paste(missing,collapse=', ')}) are missing ")
samples <- robisample(ngsfilter)
......@@ -126,7 +138,7 @@ read_metabar <- function(file,
}
robimetabar(motus = robimotu(features %>% select(id,everything())),
samples = robisample(unique(data$sample)),
samples = samples,
data = data,
verbose = verbose)
}
......@@ -39,13 +39,13 @@ read_ngsfilter <- function(filename) {
separate(tags,
into = c("forward_tag","reverse_tag"),
sep=":", remove = TRUE, fill = "right") %>%
mutate(id = robiuniqueid(id),
mutate(id = as_robiuniqueid(id),
reverse_tag = ifelse(is.na(reverse_tag),forward_tag,reverse_tag),
full_length = parse_logical(full_length)) %>%
mutate(forward_tag = robitag_forward(ifelse(forward_tag=="-",
mutate(forward_tag = as_robitag_forward(ifelse(forward_tag=="-",
NA,
forward_tag)),
reverse_tag = robitag_reverse(ifelse(reverse_tag=="-",
reverse_tag = as_robitag_reverse(ifelse(reverse_tag=="-",
NA,
reverse_tag))) -> left
......
......@@ -39,7 +39,8 @@ validate_object.robiatomic <- function(object, verbose=is_robi_verbose()) {
#' @export
new_robiatomic <- function(object,...,class=NULL) {
vc <- new_vctr(x, class = c(class, "robiatomic","robiobject"))
vc <- new_vctr(object,
class = c(class, "robiatomic","robiobject"))
validate_object(vc)
}
......@@ -87,52 +88,7 @@ as_robiatomic.robiatomic <- function(object) {
#' @export
c.robiatomic <- function(...) {
#fcall = match.call()
d <- list(...)
all_robiatomic <- all(sapply(d,is_robiatomic))
#fcall[[1]] = quote(base::c)
result <- NextMethod()
if (all_robiatomic) {
classes <- lapply(d, class)
classes <- lapply(classes, rev)
mclass <- min(sapply(classes,length))
common_class <- classes[[1]][(which(sapply(1:mclass, function(i)
length(unique(sapply(classes,function(x) x[[i]]))) == 1
)) %>% rev())[1]]
fcall = call(paste0("new_",common_class),quote(result))
result <- tryCatch(eval(fcall),
error = function(e) robiatomic(result))
}
result
vec_c(...)
}
#' @export
rep.robiatomic <- function(x, ...) {
constructor <- paste0("new_",class(x)[1])
result <- NextMethod()
eval(call(constructor,quote(result)))
}
#' @export
`[.robiatomic` <- function(x, ..., drop = FALSE) {
constructor <- paste0("new_",class(x)[1])
result <- NextMethod()
eval(call(constructor,quote(result)))
}
#' @export
`[<-.robiatomic` <- function(x, ..., value) {
constructor <- paste0("new_",class(x)[1])
result <- NextMethod()
eval(call(constructor,quote(result)))
}
#' @include robiobject.R
#' @include robiatomic.R
#' @import tibble
#' @import dplyr
#' @import vctrs
#'
NULL
.robi_all_categories <- c("sample", "positive_ctrl", "extraction_blk", "pcr_blk", "sequencing_blk")
......@@ -28,7 +30,7 @@ robicategory_mapping <- function(sample="s",positive_ctrl="PCRPos",
m <- as.character(c(sample, positive_ctrl, extraction_blk, pcr_blk, sequencing_blk))
mapping <- function(data) (
robicategory(.robi_all_categories[match(as.character(data),m,nomatch = NA)],
new_robicategory(.robi_all_categories[match(as.character(data),m,nomatch = NA)],
verbose = verbose)
)
}
......@@ -61,20 +63,35 @@ is_robicategory <- function(object) {
#' new_robiuniqueid(c("A", "A", "B", "C"))
#' @author Eric Coissac <eric.coissac@metabarcoding.org>
#' @export
new_robicategory <- function(object, ..., class = character()) {
new_robicategory <- function(object, ...,
verbose = is_robi_verbose(),
class = character()) {
object <- vec_cast(object,character())
fac <- match(object, .robi_all_categories)
if (!is.factor(object)) {
stop("robicategory can only be based on q factor",
call. = FALSE
)
if (verbose && any(is.na(fac))) {
warning(ngettext(
sum(is.na(fac)),
sprintf(
"The element %s doesn't belong one of the following categories",
(object[is.na(fac)])[1],
substr(paste(deparse(.robi_all_categories), collapse = ""), 2, 1000)
),
sprintf(
"The elements %s don't belong one of the following categories",
substr(paste(deparse(unique(object[is.na(fac)])), collapse = ""), 2, 1000),
substr(paste(deparse(.robi_all_categories), collapse = ""), 2, 1000)
),
))
}
validate_object(
new_robiatomic(object,
...,
class = c(class, "robicategory")
)
)
fac <- new_factor(fac,levels = .robi_all_categories,
class = c(class, "robicategory",
"robiatomic","robiobject"))
validate_object(fac)
}
validate_object.robicategory <- function(object) {
......@@ -111,32 +128,103 @@ validate_object.robicategory <- function(object) {
#' ))
#' @author Eric Coissac <eric.coissac@metabarcoding.org>
#' @export
robicategory <- function(object = character(0), verbose = is_robi_verbose()) {
robicategory <- function(x = 0) {
new_robicategory(rep("sample",x))
}
object <- as.character(object)
fac <- match(object, .robi_all_categories)
#' @export
vec_ptype2.robicategory <- function(x, y, ...) {
UseMethod("vec_ptype2.robicategory", y)
}
#' @expmethod vec_ptype2.robicategory default
#' @export
vec_ptype2.robicategory.default <- function(x, y, ..., x_arg = "x", y_arg = "y") {
vec_default_ptype2(x, y, x_arg = x_arg, y_arg = y_arg)
}
#' @method vec_ptype2.robicategory robicategory
#' @export
vec_ptype2.robicategory.robicategory <- function(x, y, ...) robicategory()
#' @method vec_ptype2.robicategory integer
#' @export
vec_ptype2.robicategory.integer <- function(x, y, ...) integer()
#' @method vec_ptype2.integer robicategory
#' @export
vec_ptype2.integer.robicategory <- function(x, y, ...) integer()
#' @method vec_ptype2.robicategory double
#' @export
vec_ptype2.robicategory.double <- function(x, y, ...) double()
#' @method vec_ptype2.double robicategory
#' @export
vec_ptype2.double.robicategory <- function(x, y, ...) double()
#' @method vec_ptype2.robicategory character
#' @export
vec_ptype2.robicategory.character <- function(x, y, ...) character()
#' @method vec_ptype2.character robicategory
#' @export
vec_ptype2.character.robicategory <- function(x, y, ...) character()
#' @method vec_ptype2.robicategory factor
#' @export
vec_ptype2.robicategory.factor <- function(x, y, ...) factor()
#' @method vec_ptype2.factor robicategory
#' @export
vec_ptype2.factor.robicategory <- function(x, y, ...) factor()
#' @export
vec_cast.robicategory <- function(x, to, ...)
UseMethod("vec_cast.robicategory")
#' @method vec_cast.robicategory default
#' @export
vec_cast.robicategory.default <- function(x, to, ...)
vec_default_cast(x, to)
#' @method vec_cast.robicategory robicategory
#' @export
vec_cast.robicategory.robicategory <- function(x, to, ...) x
#' @method vec_cast.integer robicategory
#' @export
vec_cast.integer.robicategory <- function(x, to, ...)
vec_data(x)
#' @method vec_cast.double robicategory
#' @export
vec_cast.double.robicategory <- function(x, to, ...)
vec_cast(vec_data(x),double())
#' @method vec_cast.factor robicategory
#' @export
vec_cast.factor.robicategory <- function(x, to, ...)
new_factor(vec_data(x),levels = levels(x))
#' @method vec_cast.robicategory factor
#' @export
vec_cast.robicategory.factor <- function(x, to, ...)
new_robicategory(vec_cast(x,character()))
#' @method vec_cast.robicategory character
#' @export
vec_cast.robicategory.character <- function(x, to, ...)
new_robicategory(x)
#' @method vec_cast.character robicategory
#' @export
vec_cast.character.robicategory <- function(x, to, ...)
levels(x)[vec_data(x)]
if (verbose && any(is.na(fac))) {
warning(ngettext(
sum(is.na(fac)),
sprintf(
"The element %s doesn't belong one of the following categories",
(object[is.na(fac)])[1],
substr(paste(deparse(.robi_all_categories), collapse = ""), 2, 1000)
),
sprintf(
"The elements %s don't belong one of the following categories",
substr(paste(deparse(unique(object[is.na(fac)])), collapse = ""), 2, 1000),
substr(paste(deparse(.robi_all_categories), collapse = ""), 2, 1000)
),
))
}
class(fac) <- "factor"
attr(fac, "levels") <- .robi_all_categories
new_robicategory(fac)
}
......@@ -157,38 +245,11 @@ as_robicategory <- function(object) {
UseMethod("as_robicategory", object)
}
#' @rdname as_robicategory
#' @export
as_robicategory.character <- function(object) {
new_robicategory(factor(object,
levels = .robi_all_categories))
}
#' @rdname as_robicategory
#' @export
as_robicategory.factor <- function(object) {
categories <- c("sample", "positive_ctrl", "extraction_blk", "pcr_blk", "sequencing_blk")
if (identical(levels(object), .robi_all_categories)) {
new_robicategory(object)
} else {
as_robicategory.character(as.character(object))
}
}
#' @rdname as_robicategory
#' @export
as_robicategory.default <- function(object) {
as_robicategory.character(as.character(object))
}
#' @rdname as_robicategory
#' @export
as_robicategory.robicategory <- function(object) {
classes <- class(object)
rclass <- which(classes == 'robicategory')
class(object) <- classes[rclass:length(classes)]
object
vec_cast(object,robicategory())
}
......@@ -218,33 +279,17 @@ as_robicategory.robicategory <- function(object) {
#' @author Eric Coissac <eric.coissac@metabarcoding.org>
#' @export
`levels <- .robicategory` <- function(x, value) {
stop("levels of a robicategory cannot be change")
}
robiassert_arg(identical(value,.robi_all_categories),
"value",
"levels of a robicategory cannot be change")
x
leves}
#' Formats a \code{robicategory} instance for printing
#'
#' Actually call the format method for factor
#'
#' @param x a \code{robicategory} to format
#' @param ... further arguments passed to or from other methods.
#'
#'
#' @examples
#'
#' x <- robicategory(rep(
#' c("sample", "positive_ctrl", "pcr_blk"),
#' c(3, 2, 1)
#' ))
#'
#' format(x)
#' @author Eric Coissac <eric.coissac@metabarcoding.org>
#' @export
format.robicategory <- format.factor
#' @author Eric Coissac <eric.coissac@metabarcoding.org>
#' @export
print.robicategory <- print.factor
vec_ptype_abbr.robicategory <- function(x, ...) {
"robicategory"
}
#' @importFrom pillar type_sum
#' @export
......
#' @include robiobject.R
#' @import tibble
#' @import dplyr
#' @import vctrs
NULL
#' Test if an object belongs \code{robidata} class
......@@ -192,11 +194,10 @@ rbind.robidata <- function(..., deparse.level = 1) {
}
# # @export
# rbind.robidata <- function(..., deparse.level = 1,
# make.row.names = TRUE,
# stringsAsFactors = default.stringsAsFactors(),
# factor.exclude = NA) {
#' @export
rbind.robidata <- function(..., deparse.level = 1) {
vec_rbind(...)
}
#
# fcall <- match.call()
# fcall[[1]] <- quote(base::rbind.data.frame)
......
......@@ -603,5 +603,8 @@ samples_ids.robimetabar <- function(object) {
samples(object) %>% ids()
}
#' @export
clone.robimetabar <- function(object) {
object$clone()
}
......@@ -74,7 +74,7 @@ build_motu_names <- function(n,
snames <- sprintf("%s%0*d",prefixes,ndigits,ns)