Description Details Author(s) See Also Examples
The package sybilccFBA
implements some methods to get cost constrained fluxes.
It is required to supply the molecular weights. It can be calculated from genome data using function calc_MW
.
Also requires kinetic data along with the model.
Package: | sybilccFBA |
Type: | Package |
Version: | 0.0.1 |
Date: | 2013-06-03 |
License: | GPL Version 3 |
LazyLoad: | yes |
Depends: | sybil , methods |
Abdelmoneim Amer Desouki
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 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 | ## Not run:
data(Ec_core)
model=Ec_core
genedef=read.csv(paste0(path.package("sybilccFBA"), '/extdata/Ec_core_genedef.csv'),
stringsAsFactors=FALSE)
mw=data.frame(gene=genedef[,'gene'],mw=genedef[,'mw'],stringsAsFactors=FALSE)
mw[mw[,1]=='s0001','mw']=0.001#spontenious
kl=read.csv(stringsAsFactors=FALSE,paste0(path.package("sybilccFBA"),
'/extdata/','allKcats_upd34_dd_h.csv'))
kl=kl[!is.na(kl[,'ijo_id']),]
kcat=data.frame(rxn_id=kl[,'ijo_id'],val=kl[,'kcat_max'],dirxn=kl[,'dirxn'],src=kl[,'src'],
stringsAsFactors=FALSE)
kcat=kcat[kcat[,'rxn_id']%in% react_id(model),]
kcat[(is.na(kcat[,'src'])),'src']='Max'
mod2=mod2irrev(model)
uppbnd(mod2)[react_id(mod2)=="EX_o2(e)_b"]=1000
lowbnd(mod2)[react_id(mod2)=="ATPM"]=0
uppbnd(mod2)[react_id(mod2)=="ATPM"]=0
nr=react_num(mod2)
medianVal=median(kcat[,'val'])
# sum(is.na(genedef[,'mw']))
C=0.13#as not all genes exist
sim_name=paste0('Ec_core_med',round(medianVal,2),'_C',100*C)
cpx_stoich =read.csv(paste0(path.package("sybilccFBA"),
'/extdata/','cpx_stoich_me.csv'),stringsAsFactors=FALSE)
CSList=c("R_EX_glc_e__b","R_EX_glyc_e__b","R_EX_ac_e__b","R_EX_fru_e__b",
"R_EX_pyr_e__b","R_EX_gal_e__b",
"R_EX_lac_L_e__b","R_EX_malt_e__b","R_EX_mal_L_e__b","R_EX_fum_e__b",
"R_EX_xyl_D_e__b","R_EX_man_e__b","R_EX_tre_e__b",
"R_EX_mnl_e__b","R_EX_g6p_e__b","R_EX_succ_e__b","R_EX_gam_e__b",
"R_EX_sbt_D_e__b","R_EX_glcn_e__b",
"R_EX_rib_D_e__b","R_EX_gsn_e__b","R_EX_ala_L_e__b",
"R_EX_akg_e__b","R_EX_acgam_e__b")
msrd=c(0.66,0.47,0.29,0.54,0.41,0.24,0.41,0.52,0.55,0.47,0.51,0.35,0.48,0.61,
0.78,0.50,0.40,0.48,0.68,0.41,0.37,0.24,0.24,0.61)
CA=c(6,3,2,6,3,6,3,12,4,4,5,6,12,6,6,4,6,6,6,5,10,3,5,8)
CSList=substring(gsub('_e_','(e)',CSList),3)
react_name(mod2)[react_id(mod2) %in% CSList]
msrd=msrd[CSList %in% react_id(mod2) ]
CA=CA[CSList %in% react_id(mod2)]
CSList=CSList[CSList %in% react_id(mod2)]
uppbnd(mod2)[react_id(mod2) %in% CSList]=0
mod2R=mod2
##---------------------
# xt=FALSE
st=TRUE
selected_rxns=react_id(model)[gpr(model)!=""]
if(st){
mrMat = getccFBA_mat(model,mod2,kcat,MW=mw,verbose=2,RHS=C,cpx_stoich=cpx_stoich,
medval=3600*medianVal,selected_rxns=selected_rxns)
}else{#nost
mrMat = getccFBA_mat(model,mod2,kcat,MW=mw,verbose=2,RHS=C,cpx_stoich=NULL,
medval=3600*medianVal,selected_rxns=selected_rxns)
}
mrMat_st = mrMat
r1i=getGpr1iso(model,MW=mw,cpx_stoich=cpx_stoich)
gpr1iso=r1i$gpr1iso;
mod_cpx_mw=r1i$mod_cpx_mw; rgst=r1i$rgst;
sum(is.na(gpr1iso[,'stoich']))
nOR = sapply(gprRules(model),function(x) nchar(x)-nchar(gsub('|','',x,fixed=TRUE)))
nAND = sapply(gprRules(model),function(x) nchar(x)-nchar(gsub('&','',x,fixed=TRUE)))
rgm=mrMat$rxnGeneCol
write.csv(file='rgm_aa.csv',cbind(rgm,nOR[rgm[,'revRxn']],gpr(model)[rgm[,'revRxn']]))
rgm=rgm[!is.na(rgm[,'gCol']) & rgm[,'stoich']!=0,]# only used genes
cnames=mrMat$cnames
gpr1iso=mrMat$gpr1iso
bm_rxn=which(obj_coef(mod2)!=0)
all_flx=NULL
all_flx_mr=NULL
all_flx_org=NULL
Kcatorg=kcat
solver='glpkAPI'
solverParm=NA
for(cs in 1:length(CSList)){
print(CSList[cs])
mod2=mod2R
uppbnd(mod2)[react_id(mod2) %in% CSList]=0
uppbnd(mod2)[react_id(mod2)==CSList[cs]]=1000
prob=sysBiolAlg(mod2,LHS=mrMat$LHS,rlb=mrMat$rlb,rub=mrMat$rub,rtype=mrMat$rtype,
lb=mrMat$clb,ub=mrMat$cub,obj=mrMat$obj_cf,lpdir='max',cnames=mrMat$cnames,solver=solver,
algorithm="ccFBA",solverParm=solverParm,
pname=sprintf("ccFBA mfe %s: %s,cpxst,kcatupd25",solver,mod_name(model)),
writeProbToFileName=sprintf('EC_ccFBA_Mat_%s.lp',solver));
sol=optimizeProb(prob);
str(sol)
print(getMeanStatus(sol$stat,solver))
sol_mr=cfba_moment_mr(model,mod2,kcat,MW=mw,verbose=1,RHS=C,solver="glpkAPI",
medval=3600*medianVal)
sol_org=cfba_moment(model,mod2,kcat,MW=mw,verbose=1,RHS=C,solver="glpkAPI",
medval=3600*medianVal)
### preparing output -------------------
all_flx=rbind(all_flx,data.frame(stringsAsFactors=FALSE,cs,csname=CSList[cs],
rxn_id=react_id(mod2),flx=sol$fluxes[1:nr]))
all_flx_mr=rbind(all_flx_mr,data.frame(stringsAsFactors=FALSE,cs,csname=CSList[cs],
rxn_id=react_id(mod2),flx=sol_mr$sol$fluxes[1:nr]))
all_flx_org=rbind(all_flx_org,data.frame(stringsAsFactors=FALSE,cs,
csname=CSList[cs],rxn_id=react_id(mod2),flx=sol_org$sol$fluxes[1:nr]))
}
upt=all_flx[all_flx[,'csname']==all_flx[,'rxn_id'],]
#bm=all_flx[all_flx[,'obj']!=0,]
bm=all_flx[react_id(mod2)[obj_coef(mod2)!=0]==all_flx[,'rxn_id'],]
bm_mr=all_flx_mr[react_id(mod2)[obj_coef(mod2)!=0]==all_flx_mr[,'rxn_id'],]
bm_org=all_flx_org[react_id(mod2)[obj_coef(mod2)!=0]==all_flx_org[,'rxn_id'],]
cor.test(bm[,'flx'],msrd,method='spearman')
cor.test(bm_mr[,'flx'],msrd,method='spearman')
cor.test(bm_org[,'flx'],msrd,method='spearman')
bm=cbind(bm,msrd)
#Wong
cor(as.numeric(bm[,"flx"]),(as.numeric(bm[,"flx"])/as.numeric(upt[,"flx"]))^0.5 )
\dontrun{
plot(msrd,bm[,'flx'],ylab="pred bm",xlim=c(0,0.8),ylim=c(0,max(bm_org[,'flx'])),
pch=19,col='red',
main="Predicting growth rate using Ec_core")
abline(a=0,b=1,col="grey",lwd=2,lty=2)
points(msrd,bm_org[,'flx'],ylab="pred bm",xlim=c(0,0.8),ylim=c(0,max(bm_org[,'flx'])),
pch=19,col='green')
points(msrd,bm_mr[,'flx'],ylab="pred bm",xlim=c(0,0.8),ylim=c(0,max(bm_org[,'flx'])),
pch=19,col='darkblue')
legend(legend=c('MOMENT','MOMENT_mr','ccFBA st'),col=c('green','darkblue','red'),
x=0,y=1,cex=0.75,pch=19)
################################top 5 genes####################
# selected_genes =c('b2323','b2472','b3774','b2779','b0431') # excl 'b0241', too many rxns 248
selected_genes =c('b1773','b0118','b1779') # excl 'b0241', too many rxns 248
sel_rxns =react_id(mod2)[(rowSums(rxnGeneMat(mod2)[,allGenes(mod2)%in% selected_genes])>0)]
sel_rxns=c(sel_rxns,CSList[cs],react_id(mod2)[obj_coef(mod2)!=0])
write.csv(file=sprintf('iAF_top4g_flx_le%s.csv',sim_name),
all_flx[all_flx[,'rxn_id']%in% sel_rxns,])
}
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.