Commit 2fdf1691 authored by AITE Meziane's avatar AITE Meziane

removing old scripts

parent 01e2c015
# -*- coding: utf-8 -*-
"""
This file is part of padmet-utils.
padmet-utils is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
padmet-utils is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with padmet-utils. If not, see <http://www.gnu.org/licenses/>.
@author: Meziane AITE, meziane.aite@inria.fr
Description:
1./ extract all data from sgs output
2./ If needed, convert reaction id to metacyc id
usage:
enhanced_sgs_output.py --sgs_output=FILE --padmetRef=FILE --output=FILE [--db=ID] [--mnx=FILE] [-v]
options:
-h --help Show help.
--sgs_output=FILE pathname of a sgs run' result
--db=ID database origin of the reactions. will be converted to metacyc
--padmetRef=FILE pathanme to padmet file corresponding to the database of reference (the repair network)
--mnx=FILE pathanme to metanetx file for reactions (reac_xref.tsv)
--output=FILE pathname to tsv output file
"""
import csv
import re
import docopt
from itertools import izip_longest
def main():
args = docopt.docopt(__doc__)
sgs_output = args["--sgs_output"]
output = args["--output"]
verbose = args["-v"]
output = args["--output"]
verbose = True
sgs_output = "/home/maite/Documents/data/E.faecalis_data/data_askomics/efaecalis_sgs_final.txt"
output = "/home/maite/Documents/data/E.faecalis_data/data_askomics/sgs.txt"
sgs_dict = {}
#1: extract data from shogen output.
count = 1
if verbose: print("Reading sgs output")
with open(sgs_output, 'r') as csvfile:
reader = csv.DictReader(csvfile, delimiter="\t")
for row in reader:
genes, reactions = row["gene_set"], row["reaction_set"]
genes = re.sub("ac:|nc:|na:|\"","",genes).split(" ")
reactions = re.sub("nc:|\"|R_","",reactions).split(" ")
sgs_id = "sgs-"+str(count)+"_"+str(len(genes))
count += 1
sgs_dict[sgs_id] = {"genes": genes, "reactions": reactions,\
"start_position": row["start_position"], "end_position": row["end_position"],\
"length": row["length"], "density": str(round(float(row["density"]),2)),"chr_id": row["chromosome_id"]}
if verbose: print("creating output %s" %output)
with open(output, 'w') as f:
header = ["SGS","contains@gene","concerns@reaction","start_position","end_position","length","density","chromosome_id"]
f.write("\t".join(header)+"\n")
for sgs_id, data in sgs_dict.items():
all_data_ziped = izip_longest(data["genes"],data["reactions"],[data["start_position"]],[data["end_position"]],[data["length"]],[data["density"]],[data["chr_id"]])
for line_data in all_data_ziped:
line = [sgs_id]
line += [x if x else "" for x in line_data]
f.write("\t".join(line)+"\n")
if __name__ == "__main__":
main()
This diff is collapsed.
This diff is collapsed.
# -*- coding: utf-8 -*-
"""
This file is part of padmet-utils.
padmet-utils is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
padmet-utils is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with padmet-utils. If not, see <http://www.gnu.org/licenses/>.
@author: Meziane AITE, meziane.aite@inria.fr
Description:
Run flux balance analyse with cobra package. If the flux is >0. Run also FVA
and return result in standard output
usage:
fba_on_reaction.py --sbml=FILE
option:
-h --help Show help.
--sbml=FILE pathname to the sbml file to test for fba and fva.
"""
from padmet.sbmlPlugin import convert_from_coded_id
from cobra import *
from cobra.io.sbml import create_cobra_model_from_sbml_file
import docopt
def main():
args = docopt.docopt(__doc__)
sbml_file = args["--sbml"]
model=create_cobra_model_from_sbml_file(sbml_file)
solution = model.optimize()
try:
biomassrxn = [rxn for rxn in model.reactions if rxn.objective_coefficient == 1.0][0]
biomassname = biomassrxn.id
except IndexError:
print("Need to set OBJECTIVE COEFFICIENT to '1.0' for the reaction to test")
exit()
print("Testing reaction %s" %biomassname)
print( 'growth rate: %s\nStatus: %s.' %(solution.objective_value, solution.status))
if (solution.objective_value > 1e-5):
#FVA_result = flux_analysis.variability.flux_variability_analysis(model, fraction_of_optimum=1.0)
FVA_result = flux_analysis.variability.flux_variability_analysis(model, fraction_of_optimum=1.0, tolerance_feasibility=1e-6, tolerance_integer=1e-6)
#print(str(FVA_result)+"\n")
essential, alternative, blocked = [[] for i in range(3)]
for k,v in FVA_result.iterrows():
_min = float(v['minimum'])
_max = float(v['maximum'])
if (_min > 0.0001 and _max > 0.0001) or (_min < -0.0001 and _max < -0.0001):
essential.append(k)
elif _min > -0.0001 and _min < 0.0001 and _max > -0.0001 and _max < 0.0001:
blocked.append(k)
else:
alternative.append(k)
print('FVA analysis:')
print('\t'+'Essential reactions: '+str(len(essential)))
print('\t'+'Alternative reactions: '+str(len(alternative)))
print('\t'+'Blocked reactions: '+str(len(blocked)))
print('Essential reactions: '+str(len(essential)))
for r_id in essential:
r_id_decoded = convert_from_coded_id(r_id)[0]
print('\t'+r_id_decoded+' ('+r_id+')')
print('Alternative reactions: '+str(len(alternative)))
for r_id in alternative:
r_id_decoded = convert_from_coded_id(r_id)[0]
print('\t'+r_id_decoded+' ('+r_id+')')
print('Blocked reactions: '+str(len(blocked)))
for r_id in blocked:
r_id_decoded = convert_from_coded_id(r_id)[0]
print('\t'+r_id_decoded+' ('+r_id+')')
else:
# print(metabolites)
# only append to the list compounds that are reactants
bms_reactants = biomassrxn.reactants
biomassrxn.objective_coefficient = 0
# for metabo in metabolites:
# metabolites_list.append(metabo.id)
#print(metabolites_list)
dict_output = {"positive":{},"negative":{}}
for metabolite in bms_reactants:
#lets create a copy of the initial model
model2 = model.copy()
#get the biomass reaction
biomassrxn2 = model2.reactions.get_by_id(biomassname) #R_Ec_biomass_iAF1260_core_59p81M
biomassrxn2.objective_coefficient = 1
# empty the biomass metabolites
for i in range(len(bms_reactants)):
biomassrxn2.reactants.pop(i)
#create a clean dict with only compound of interest and stoichio to -1 (reactant)
metabolitedict = {}
metabolitedict[metabolite]=-1.0
# replace the biomass metabolites with only metabolite1 inside
biomassrxn2.add_metabolites(metabolitedict)
# print(biomassrxn2._metabolites)
#test flux balance analysis
#print(biomassrxn2._metabolites)
solution2 = model2.optimize()
if (solution2.objective_value > 1e-5):
dict_output["positive"][metabolite] = solution2.objective_value
else:
dict_output["negative"][metabolite] = solution2.objective_value
for k,v in dict_output["positive"].items():
print("%s %s positive" %(k,v))
for k,v in dict_output["negative"].items():
print("%s %s negative" %(k,v))
print("%s/%s compounds with positive flux" %(len(dict_output["positive"].keys()), len(bms_reactants)))
print("%s/%s compounds with negative flux" %(len(dict_output["negative"].keys()), len(bms_reactants)))
if __name__ == "__main__":
main()
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
This file is part of padmet-utils.
padmet-utils is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
padmet-utils is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with padmet-utils. If not, see <http://www.gnu.org/licenses/>.
@author: Meziane AITE, meziane.aite@inria.fr
Description:
For a given sbml file, return the DBs of each reaction, based on metanetx file.
usage:
check_db.py --sbml=FILE --mnx=FILE
options:
-h --help Show help.
--sbml=FILE Path to the sbml file to be analysed.
--mnx=FILE path to the metatnetx file for reaction, reac_xref.tsv.
"""
import libsbml
from padmet.utils.sbmlPlugin import convert_from_coded_id
import docopt
def main():
args = docopt.docopt(__doc__)
rxn_mnx = args["--mnx"]
sbml_file = args["--sbml"]
with open(rxn_mnx,'r') as f:
db = dict([tuple(line.split("\t")[0].split(":")[::-1]) for line in f.read().splitlines() if not line.startswith("#")])
reader = libsbml.SBMLReader()
document = reader.readSBML(sbml_file)
for i in range(document.getNumErrors()):
print (document.getError(i).getMessage())
model = document.getModel()
listOfSpecies = model.getListOfSpecies()
listOfReactions = model.getListOfReactions()
nbReactions = len(listOfReactions)
nbSpecies = len(listOfSpecies)
print("nb species: %s" %nbSpecies)
print("nb reactions: %s" %nbReactions)
data = [convert_from_coded_id(r.id)[0] for r in listOfReactions]
result_dict = {"None":[]}
for rxn in data:
if rxn in db.keys():
try:
result_dict[db[rxn]].append(rxn)
except KeyError:
result_dict[db[rxn]] = [rxn]
else:
result_dict["None"].append(rxn)
for k,v in result_dict.iteritems():
print("%s: %s reactions" %(k,len(v)))
for r in v:
print("\t%s" %r)
if __name__ == "__main__":
main()
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 30 16:34:38 2018
@author: maite
toy for ptg tests
Benchmark:
T_lutea
Athalia, Chlamy, e_sill, synecho
Pour chaque model:
1./ Lire output pantograph, récupérer N reactions et faire un sous sbml
2./ Récuéprer la lsit des genes assoc de ces N reactions, lire le fasta Study et faire un sous fasta avec X% de bruit
3./ Lire le sbml du Model, recup la liste des genes assoc et faire
avec du bruit X: % de genes totaux or cela à gader dans le fasta
Tiso: Tiso_gene_9894 / R_ARYLDIALKYL__45__PHOSPHATASE__45__RXN
Ecto: Ec-27_006880
EC-27_006
"""
import libsbml
from Bio import SeqIO
def main():
output_pantograph = "/home/maite/Forge/docker/aureme_workspace/Sjap/networks/output_orthology_based_reconstruction/pantograph/output_pantograph_Nannochloropsis_salina.sbml"
study_fasta = "/home/maite/Forge/docker/aureme_workspace/Sjap/orthology_based_reconstruction/Nannochloropsis_salina/FAA_model.faa"
model_fasta = "/home/maite/Forge/docker/aureme_workspace/Sjap/orthology_based_reconstruction/Nannochloropsis_salina/FAA_model.faa"
reader = libsbml.SBMLReader()
document = reader.readSBML(output_pantograph)
for i in range(document.getNumErrors()):
print (document.getError(i).getMessage())
model = document.getModel()
listOfReactions = model.getListOfReactions()
def sub_sbml(model, list_rxn):
return new_model
def sub_fasta(fasta, list_genes, noise):
return new_fasta
\ No newline at end of file
......@@ -14,3 +14,11 @@
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1540313517 gem_a1562_fibrogem_UserID 1
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1540313517 gem_a1562_fibrogem_UserName Dyliss
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1540313517 gem_a1562_fibrogem_Token 2a50c61b6f36144dfd937110c48df7f5
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 0 gem_a1562_sjapgem__session ku9s1ctm43u42pcs00j5vbe9i2
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1542197507 gem_a1562_sjapgem_UserID 1
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1542197507 gem_a1562_sjapgem_UserName Dyliss
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1542197507 gem_a1562_sjapgem_Token bd9bb61a3f57562a6df9fae873a9e06a
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 0 gem_a1562_ccrgem__session 3o8onclovla4dqtgs69io43nj2
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1542644157 gem_a1562_ccrgem_UserID 1
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1542644157 gem_a1562_ccrgem_UserName Dyliss
#HttpOnly_gem-aureme.irisa.fr FALSE / FALSE 1542644157 gem_a1562_ccrgem_Token 8afc6728636c171388119420d24f5b74
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment