gen2Rule: get reaction expression state from gene expression

Description Usage Arguments Value Author(s) See Also Examples

Description

function to get state of expression of GPR rules from gene expression data

# rules in brief: #1-Complexes: average, 2-isoenzymes: sum #3-multifunctioning: divide by count

N.B.: GPR rules should be in form Sum-of-products (AND to OR)

Usage

1
gene2Rule(model, geneExpr, selected_rxns = NULL)

Arguments

model

An object of class modelorg.

geneExpr

a data frame: GeneID, expr_val for each gene.

selected_rxns

optional parameter to select only a set of reactions.

Value

return list with main slot:

ruleExpr: rxn_id,expr_val

Author(s)

Abdelmoneim Amer Desouki

See Also

modelorg, FECorr

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
## Not run: 
##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (model, geneExpr, selected_rxns = NULL) 
{
    if (length(selected_rxns) == 0) {
        ugpr = as.vector(unique(gpr(model)))
    }
    else {
        ugpr = as.vector(unique(gpr(model)[selected_rxns]))
    }
    ugpr = ugpr[ugpr != ""]
    gprExpr = NULL
    for (v_rule in ugpr) {
        rl = gsub("\)", " ) ", v_rule)
        rl = gsub("\(", " ( ", rl)
        pr = lapply(strsplit(unlist(strsplit(rl, " or ")), " and "), 
            function(x) gsub("[() ]", "", x))
        expr_val = 0
        for (p in 1:length(pr)) {
            gene_ind = match(pr[[p]], geneExpr$geneID)
            if (length(gene_ind) < length(pr[[p]])) {
                warning(sprintf("Rule %s containing gene names not in geneID list, 
				term no: %d term: %s ", 
                  v_rule, p, pr[[p]][1]))
            }
            else {
                expr_val = expr_val + mean(geneExpr[gene_ind, 
                  "expr_val"])
            }
        }
        cnt = sum(gpr(model) == v_rule)
        gprExpr = rbind(gprExpr, cbind(rxn_id = react_id(model)[gpr(model) == 
            v_rule], expr_val = expr_val/cnt, gpr = v_rule, cnt = rep(cnt, 
            cnt)))
    }
    return(gprExpr)
  }

## End(Not run)

sybilEFBA documentation built on May 29, 2017, 9:35 a.m.