-
Notifications
You must be signed in to change notification settings - Fork 0
/
ess_genes_script.r
25 lines (20 loc) · 860 Bytes
/
ess_genes_script.r
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
##jmd
##1.23.12
##ko_all.r
source('/msc/neurospora/FBA/farm/farm_header.r')
essential.genes.thresh <- 10**-6
#each gene's ko'd rxns, which are in model
ko.rxns <- NCUvector2rxns(ncu.v=genes, gpr=gpr, gpr.name='model GPR')
#intitialize matrix of genes that individually KO >= 1 rxn
g.mat <- data.frame(genes=names(ko.rxns), growth=NA)
rownames(g.mat) <- names(genes)[match(g.mat$genes, genes)]
#loop
g.mat$growth <- sapply(ko.rxns, FUN=function(x){
fba.tmp <- fba.na(a=s.al, ko=x, fba.ub=fva.ub, quiet=TRUE, control=list(trace=0))
ret <- signif(fba.tmp$obj, 3)
})
ess.ind <- which(g.mat$growth<=essential.genes.thresh)
nrow(essential.genes <- g.mat[ess.ind,])
(essential.genes <- essential.genes[order(rownames(essential.genes)),])
write.csv(essential.genes, 'pred_essentials.csv', quote=FALSE)
write.csv(g.mat, 'singleKO_growth.csv', quote=FALSE)