|
| 1 | +# Extraction of features from Recon data frames with reaction annotation |
| 2 | + |
| 3 | +# Mapping of reactions to gene identifiers ------- |
| 4 | + |
| 5 | +#' Extract gene identifiers from reaction annotations. |
| 6 | +#' |
| 7 | +#' @description |
| 8 | +#' The functions extract Entrez ID gene identifiers mapped to all or |
| 9 | +#' user-specified reactions present in a data frame with reaction annotation |
| 10 | +#' and gene - reaction association rules. |
| 11 | +#' Entrez IDs of genes associated with reactions are listed, and character |
| 12 | +#' strings with the gene - reaction association rules are translated to R |
| 13 | +#' expressions. |
| 14 | +#' `as_reactDB()` offers a simplified tools for generation of |
| 15 | +#' \code{\link{reactDB}} data frames, |
| 16 | +#' while `extract_genes()` allows for selection of reactions of interest and |
| 17 | +#' parsing error diagnostics. |
| 18 | +#' |
| 19 | +#' @details |
| 20 | +#' Association rules with unrecognized Entrez IDs are removed from the output |
| 21 | +#' and a parsing warning is raised. |
| 22 | +#' The gene association rules have to operate with "bare" Entrez IDs without |
| 23 | +#' the version information (i.e. numbers after a dot). |
| 24 | +#' |
| 25 | +#' @return a data frame of class \code{\link{reactDB}} containing reaction |
| 26 | +#' IDs (`id`), reaction names (`name`), subsystem information (`subsystem`), |
| 27 | +#' list of gene identifiers (`entrez_id`) and a list of expressions to be |
| 28 | +#' evaluated by at calculating reaction activity estimates. |
| 29 | +#' |
| 30 | +#' @param x a data frame with the following obligatory columns: `id` with reaction |
| 31 | +#' identifiers, `name` with reaction names, `subsystem` with subsystem assignment, |
| 32 | +#' and `gene_association` with character strings with gene association rules. |
| 33 | +#' @param react_id BiGG reaction ID, with or without the leading 'R_' string. |
| 34 | +#' Defaults to `NULL`, which means that all reactions are mapped to genes. |
| 35 | +#' @param inspect_errors logical. If `TRUE`, the function returns parsing errors. |
| 36 | +#' @param ... additional arguments passed to `extract_genes()`. |
| 37 | +#' |
| 38 | +#' @export |
| 39 | + |
| 40 | + extract_genes <- function(x, |
| 41 | + react_id = NULL, |
| 42 | + inspect_errors = FALSE) { |
| 43 | + |
| 44 | + ## entry control ---------- |
| 45 | + |
| 46 | + if(!is.data.frame(x)) stop("`x` has to be a data frame.", call. = FALSE) |
| 47 | + |
| 48 | + fix_cols <- c("id", "name", "subsystem", "gene_association") |
| 49 | + |
| 50 | + missing_cols <- setdiff(fix_cols, names(x)) |
| 51 | + |
| 52 | + if(length(missing_cols) > 0) { |
| 53 | + |
| 54 | + stop(paste("The following obligatory columns are missing from `x`:", |
| 55 | + paste(missing_cols, collapse = ", ")), |
| 56 | + call. = FALSE) |
| 57 | + |
| 58 | + } |
| 59 | + |
| 60 | + x <- x[, fix_cols] |
| 61 | + |
| 62 | + class_check <- map_lgl(x, is.character) |
| 63 | + |
| 64 | + if(any(!class_check)) { |
| 65 | + |
| 66 | + stop(paste("The following obligatory column in `x` must be of character type:", |
| 67 | + paste(fix_cols[!class_check]), collapse = ", "), |
| 68 | + call. = FALSE) |
| 69 | + |
| 70 | + } |
| 71 | + |
| 72 | + if(any(stri_detect(na.omit(x[["gene_association"]]), |
| 73 | + regex = "\\d+\\.\\d+"))) { |
| 74 | + |
| 75 | + stop(paste("Gene association rules contain identifiers with", |
| 76 | + "version information, i.e. number after a dot.", |
| 77 | + "Please remove it to proceed."), |
| 78 | + call. = FALSE) |
| 79 | + |
| 80 | + } |
| 81 | + |
| 82 | + if(!is.null(react_id)) { |
| 83 | + |
| 84 | + stopifnot(is.character(react_id)) |
| 85 | + |
| 86 | + react_id <- |
| 87 | + ifelse(!stri_detect(react_id, regex = '^R_'), |
| 88 | + paste0('R_', react_id), react_id) |
| 89 | + |
| 90 | + x <- filter(x, .data[["id"]] %in% react_id) |
| 91 | + |
| 92 | + if(nrow(x) == 0) { |
| 93 | + |
| 94 | + stop("No reactions to precess after filtering with `react_id`.", |
| 95 | + call. = FALSE) |
| 96 | + |
| 97 | + } |
| 98 | + |
| 99 | + } |
| 100 | + |
| 101 | + stopifnot(is.logical(inspect_errors)) |
| 102 | + inspect_errors <- inspect_errors[1] |
| 103 | + |
| 104 | + ## filtering the reaction annotation df --------- |
| 105 | + |
| 106 | + ### getting rid of empty rules and space-only rules in the annotation |
| 107 | + ### data frame |
| 108 | + |
| 109 | + proc_df <- |
| 110 | + filter(x[, c("id", "gene_association")], |
| 111 | + !is.na(.data[["id"]]), |
| 112 | + .data[["id"]] != "", |
| 113 | + !is.na(.data[["gene_association"]]), |
| 114 | + .data[["gene_association"]] != "", |
| 115 | + !stri_detect(.data[["gene_association"]], |
| 116 | + regex = "^\\s+$")) |
| 117 | + |
| 118 | + ### removal of leading and trailing spaces in the gene association rules |
| 119 | + |
| 120 | + proc_df[["gene_association"]] <- |
| 121 | + stri_replace_all(proc_df[["gene_association"]], |
| 122 | + regex = "^\\s+", |
| 123 | + replacement = "") |
| 124 | + |
| 125 | + proc_df[["gene_association"]] <- |
| 126 | + stri_replace_all(proc_df[["gene_association"]], |
| 127 | + regex = "\\s+$", |
| 128 | + replacement = "") |
| 129 | + |
| 130 | + ## extraction of the Entrez IDs -------- |
| 131 | + |
| 132 | + proc_df[["entrez_id"]] <- |
| 133 | + map(proc_df[["gene_association"]], |
| 134 | + stri_extract_all, regex = "\\d+") |
| 135 | + |
| 136 | + proc_df[["entrez_id"]] <- map(proc_df[["entrez_id"]], |
| 137 | + unlist) |
| 138 | + |
| 139 | + proc_df[["entrez_id"]] <- map(proc_df[["entrez_id"]], |
| 140 | + unique) |
| 141 | + |
| 142 | + ## translation of gene assignment rules to R expressions ------ |
| 143 | + |
| 144 | + exp_lst <- |
| 145 | + map(set_names(proc_df[["gene_association"]], |
| 146 | + proc_df[["id"]]), |
| 147 | + escape_numbers) |
| 148 | + |
| 149 | + exp_lst <- |
| 150 | + map(exp_lst, |
| 151 | + stri_replace_all, |
| 152 | + regex = "and|AND", |
| 153 | + replacement = "%AND%") |
| 154 | + |
| 155 | + exp_lst <- |
| 156 | + map(exp_lst, |
| 157 | + stri_replace_all, |
| 158 | + regex = "or|OR", |
| 159 | + replacement = "%OR%") |
| 160 | + |
| 161 | + exp_lst <- map(exp_lst, safely(str2lang)) |
| 162 | + |
| 163 | + parse_errors <- compact(map(exp_lst, ~.x$error)) |
| 164 | + |
| 165 | + if(length(parse_errors) > 0) { |
| 166 | + |
| 167 | + warning(paste('There were', length(parse_errors), 'parsing errors.'), |
| 168 | + call. = FALSE) |
| 169 | + |
| 170 | + } |
| 171 | + |
| 172 | + if(inspect_errors) return(parse_errors) |
| 173 | + |
| 174 | + exp_lst <- compact(map(exp_lst, ~.x$result)) |
| 175 | + |
| 176 | + ## output ------- |
| 177 | + |
| 178 | + id <- NULL |
| 179 | + exprs <- NULL |
| 180 | + |
| 181 | + map_df <- tibble(id = names(exp_lst), |
| 182 | + exprs = exp_lst) |
| 183 | + |
| 184 | + out_df <- left_join(proc_df, map_df, by = "id") |
| 185 | + |
| 186 | + out_tbl <- left_join(x[, c("id", "name", "subsystem")], |
| 187 | + out_df, |
| 188 | + by = "id") |
| 189 | + |
| 190 | + return(reactDB(out_tbl)) |
| 191 | + |
| 192 | + } |
| 193 | + |
| 194 | +#' @rdname extract_genes |
| 195 | +#' @export |
| 196 | + |
| 197 | + as_reactDB <- function(x, ...) extract_genes(x, ...) |
| 198 | + |
| 199 | +# END ------ |
0 commit comments