Bioinformatics-package-installation-and-sample-processing/Gapseq_testing at main · suetli19/Bioinformatics-package-installation-and-sample-processing · GitHub

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
#this bash script intends to run gapseq of an example genome from genbank to validate its functionality. Assume that gapseq has been installed in ubuntu.

#change directory to gapseq
cd ./gapseq

#Set path
export PATH=/home/suetli19/gapseq:$PATH

#check path
which gapseq

#download f.prausnitzii from genbank to local folder - DS483503.1 Faecalibacterium prausnitzii M21/2 Scfld_02_24 genomic scaffold, whole genome shotgun sequence (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000154385.1/)

#zip the fasta file if needed (optional)
#tar -czvf bifido.fna.gz ./bifido.fna

#run gapseq #the fasta file is renamed to ASM15438_genomic.fna
# 1. gapseq find -p all is to predict functional pathways
gapseq find -p all -b 200 -m Bacteria ASM15438_genomic.fna

# 2. gapseq find-transport is to predict transporters
gapseq find-transport -b 200 ASM15438_genomic.fna

#3. Build draft model based on pathway, reactions and transporter prediction
gapseq draft -r ASM15438_genomic-all-Reactions.tbl -t ASM15438_genomic-Transporter.tbl -p ASM15438_genomic-all-Pathways.tbl -c ASM15438_genomic.fna -u 200 -l 100

#4. Gapfilling - if there is no available 'medium.csv' from gapseq dat/media folder, then we can predict the growth medium from the results generated by gapseq find and gapseq draft.
#a. The predicted growth medium can directly be used also as gapfilling medium in gapseq fill
gapseq medium -m ASM15438_genomic-draft.RDS -p ASM15438_genomic-all-Pathways.tbl

#b. Now, we can run gapfilling using the predicted medium.csv
gapseq fill -m ASM15438_genomic-draft.RDS -n ASM15438_genomic-medium.csv -c ASM15438_genomic-rxnWeights.RDS -g ASM15438_genomic-rxnXgenes.RDS -b 100


#5. Flux-Balance-Analysis (FBA) and Flux-Variability-Analysis (FVA) prediction of metabolic by-products #Start R in terminal by just typing "R"
#define a function, that automatically performs FBA with the minimalization of total flux (MTF) as secondary objective and FVA for all exchange reactions. The function also summarizes the results in a sorted data.table.

getMetaboliteProduction <- function(mod) {
  require(sybil)
  require(data.table)

  # MTF
  sol.mtf <- optimizeProb(mod, algorithm = "mtf")
  dt.mtf  <- data.table(ex = mod@react_id,
                        mtf.flux = sol.mtf@fluxdist@fluxes[1:mod@react_num])
  dt.mtf.tmp <- copy(dt.mtf[grepl("^EX_cpd[0-9]+_e0", ex)])

  # FVA
  sol.fv <- fluxVar(mod, react = mod@react_id[grep("^EX_cpd[0-9]+_e0", mod@react_id)])

  dt <- data.table(ex       = rep(mod@react_id[grep("^EX_cpd[0-9]+_e0", mod@react_id)],2),
                   rxn.name = rep(mod@react_name[grep("^EX_cpd[0-9]+_e0", mod@react_id)],2),
                   dir      = c(rep("l",length(grep("^EX_cpd[0-9]+_e0", mod@react_id))),
                                rep("u",length(grep("^EX_cpd[0-9]+_e0", mod@react_id)))),
                   fv       = sol.fv@lp_obj)
  dt <- dcast(dt, ex + rxn.name ~ dir, value.var = "fv")[(u>1e-6 & l >= 0)]

  dt <- merge(dt, dt.mtf, by = "ex")

  return(dt[order(-mtf.flux)])
}

#load sybil library
library(sybil)

#load the RDS file of the previous gapfilling output/constructed model for the genome/bacteria
fp <- readRDS("ASM15438_genomic.RDS")

#calculate metabolite production and show top 10 predicted metabolic by-product
getMetaboliteProduction(fp)[1:10]