Mercurial > repos > bimib > cobraxy
comparison COBRAxy/metabolic_model_setting.py @ 490:c6ea189ea7e9 draft
Uploaded
| author | francesco_lapi |
|---|---|
| date | Mon, 29 Sep 2025 15:13:21 +0000 |
| parents | 5b625d91bc7f |
| children |
comparison
equal
deleted
inserted
replaced
| 489:97eea560a10f | 490:c6ea189ea7e9 |
|---|---|
| 14 import pandas as pd | 14 import pandas as pd |
| 15 import utils.general_utils as utils | 15 import utils.general_utils as utils |
| 16 from typing import Optional, Tuple, List | 16 from typing import Optional, Tuple, List |
| 17 import utils.model_utils as modelUtils | 17 import utils.model_utils as modelUtils |
| 18 import logging | 18 import logging |
| 19 from pathlib import Path | |
| 20 | |
| 19 | 21 |
| 20 ARGS : argparse.Namespace | 22 ARGS : argparse.Namespace |
| 21 def process_args(args: List[str] = None) -> argparse.Namespace: | 23 def process_args(args: List[str] = None) -> argparse.Namespace: |
| 22 """ | 24 """ |
| 23 Parse command-line arguments for metabolic_model_setting. | 25 Parse command-line arguments for metabolic_model_setting. |
| 145 try: | 147 try: |
| 146 os.makedirs(os.path.dirname(path) or ".", exist_ok=True) | 148 os.makedirs(os.path.dirname(path) or ".", exist_ok=True) |
| 147 df.to_csv(path, sep="\t", index=False) | 149 df.to_csv(path, sep="\t", index=False) |
| 148 except Exception as e: | 150 except Exception as e: |
| 149 raise utils.DataErr(path, f"failed writing tabular output: {e}") | 151 raise utils.DataErr(path, f"failed writing tabular output: {e}") |
| 152 | |
| 153 def is_placeholder(gid) -> bool: | |
| 154 """Return True if the gene id looks like a placeholder (e.g., 0/NA/NAN/empty).""" | |
| 155 if gid is None: | |
| 156 return True | |
| 157 s = str(gid).strip().lower() | |
| 158 return s in {"0", "", "na", "nan"} # lowercase for simple matching | |
| 159 | |
| 160 def sample_valid_gene_ids(genes, limit=10): | |
| 161 """Yield up to `limit` valid gene IDs, skipping placeholders (e.g., the first 0 in RECON).""" | |
| 162 out = [] | |
| 163 for g in genes: | |
| 164 gid = getattr(g, "id", getattr(g, "gene_id", g)) | |
| 165 if not is_placeholder(gid): | |
| 166 out.append(str(gid)) | |
| 167 if len(out) >= limit: | |
| 168 break | |
| 169 return out | |
| 150 | 170 |
| 151 | 171 |
| 152 ###############################- ENTRY POINT -################################ | 172 ###############################- ENTRY POINT -################################ |
| 153 def main(args:List[str] = None) -> None: | 173 def main(args:List[str] = None) -> None: |
| 154 """ | 174 """ |
| 198 # Apply selected medium uptake bounds (negative for uptake) | 218 # Apply selected medium uptake bounds (negative for uptake) |
| 199 for reaction, value in medium.items(): | 219 for reaction, value in medium.items(): |
| 200 if value is not None: | 220 if value is not None: |
| 201 model.reactions.get_by_id(reaction).lower_bound = -float(value) | 221 model.reactions.get_by_id(reaction).lower_bound = -float(value) |
| 202 | 222 |
| 223 # Initialize translation_issues dictionary | |
| 224 translation_issues = {} | |
| 225 | |
| 203 if (ARGS.name == "Recon" or ARGS.name == "ENGRO2") and ARGS.gene_format != "Default": | 226 if (ARGS.name == "Recon" or ARGS.name == "ENGRO2") and ARGS.gene_format != "Default": |
| 204 logging.basicConfig(level=logging.INFO) | 227 logging.basicConfig(level=logging.INFO) |
| 205 logger = logging.getLogger(__name__) | 228 logger = logging.getLogger(__name__) |
| 206 | 229 |
| 207 model = modelUtils.translate_model_genes( | 230 model, translation_issues = modelUtils.translate_model_genes( |
| 208 model=model, | 231 model=model, |
| 209 mapping_df= pd.read_csv(ARGS.tool_dir + "/local/mappings/genes_human.csv", dtype={'entrez_id': str}), | 232 mapping_df= pd.read_csv(ARGS.tool_dir + "/local/mappings/genes_human.csv", dtype={'entrez_id': str}), |
| 210 target_nomenclature=ARGS.gene_format, | 233 target_nomenclature=ARGS.gene_format, |
| 211 source_nomenclature='HGNC_symbol', | 234 source_nomenclature='HGNC_symbol', |
| 212 logger=logger | 235 logger=logger |
| 213 ) | 236 ) |
| 214 | 237 |
| 238 if ARGS.name == "Custom_model" and ARGS.gene_format != "Default": | |
| 239 logging.basicConfig(level=logging.INFO) | |
| 240 logger = logging.getLogger(__name__) | |
| 241 | |
| 242 tmp_check = [] | |
| 243 for g in model.genes[1:5]: # check first 3 genes only | |
| 244 tmp_check.append(modelUtils.gene_type(g.id, "Custom_model")) | |
| 245 | |
| 246 if len(set(tmp_check)) > 1: | |
| 247 raise utils.DataErr("Custom_model", "The custom model contains genes with mixed or unrecognized nomenclature. Please ensure all genes use the same recognized nomenclature before applying gene_format conversion.") | |
| 248 else: | |
| 249 source_nomenclature = tmp_check[0] | |
| 250 | |
| 251 if source_nomenclature != ARGS.gene_format: | |
| 252 model, translation_issues = modelUtils.translate_model_genes( | |
| 253 model=model, | |
| 254 mapping_df= pd.read_csv(ARGS.tool_dir + "/local/mappings/genes_human.csv", dtype={'entrez_id': str}), | |
| 255 target_nomenclature=ARGS.gene_format, | |
| 256 source_nomenclature=source_nomenclature, | |
| 257 logger=logger | |
| 258 ) | |
| 259 | |
| 260 | |
| 261 | |
| 262 | |
| 263 if ARGS.name == "Custom_model" and ARGS.gene_format != "Default": | |
| 264 logger = logging.getLogger(__name__) | |
| 265 | |
| 266 # Take a small, clean sample of gene IDs (skipping placeholders like 0) | |
| 267 ids_sample = sample_valid_gene_ids(model.genes, limit=10) | |
| 268 if not ids_sample: | |
| 269 raise utils.DataErr( | |
| 270 "Custom_model", | |
| 271 "No valid gene IDs found (many may be placeholders like 0)." | |
| 272 ) | |
| 273 | |
| 274 # Detect source nomenclature on the sample | |
| 275 types = [] | |
| 276 for gid in ids_sample: | |
| 277 try: | |
| 278 t = modelUtils.gene_type(gid, "Custom_model") | |
| 279 except Exception as e: | |
| 280 # Keep it simple: skip problematic IDs | |
| 281 logger.debug(f"gene_type failed for {gid}: {e}") | |
| 282 t = None | |
| 283 if t: | |
| 284 types.append(t) | |
| 285 | |
| 286 if not types: | |
| 287 raise utils.DataErr( | |
| 288 "Custom_model", | |
| 289 "Could not detect a known gene nomenclature from the sample." | |
| 290 ) | |
| 291 | |
| 292 unique_types = set(types) | |
| 293 if len(unique_types) > 1: | |
| 294 raise utils.DataErr( | |
| 295 "Custom_model", | |
| 296 "Mixed or inconsistent gene nomenclatures detected. " | |
| 297 "Please unify them before converting." | |
| 298 ) | |
| 299 | |
| 300 source_nomenclature = types[0] | |
| 301 | |
| 302 # Convert only if needed | |
| 303 if source_nomenclature != ARGS.gene_format: | |
| 304 model, translation_issues = modelUtils.translate_model_genes( | |
| 305 model=model, | |
| 306 mapping_df= pd.read_csv(ARGS.tool_dir + "/local/mappings/genes_human.csv", dtype={'entrez_id': str}), | |
| 307 target_nomenclature=ARGS.gene_format, | |
| 308 source_nomenclature=source_nomenclature, | |
| 309 logger=logger | |
| 310 ) | |
| 311 | |
| 215 # generate data | 312 # generate data |
| 216 rules = modelUtils.generate_rules(model, asParsed = False) | 313 rules = modelUtils.generate_rules(model, asParsed = False) |
| 217 reactions = modelUtils.generate_reactions(model, asParsed = False) | 314 reactions = modelUtils.generate_reactions(model, asParsed = False) |
| 218 bounds = modelUtils.generate_bounds(model) | 315 bounds = modelUtils.generate_bounds(model) |
| 219 medium = modelUtils.get_medium(model) | 316 medium = modelUtils.get_medium(model) |
| 223 compartments = modelUtils.generate_compartments(model) | 320 compartments = modelUtils.generate_compartments(model) |
| 224 | 321 |
| 225 df_rules = pd.DataFrame(list(rules.items()), columns = ["ReactionID", "GPR"]) | 322 df_rules = pd.DataFrame(list(rules.items()), columns = ["ReactionID", "GPR"]) |
| 226 df_reactions = pd.DataFrame(list(reactions.items()), columns = ["ReactionID", "Formula"]) | 323 df_reactions = pd.DataFrame(list(reactions.items()), columns = ["ReactionID", "Formula"]) |
| 227 | 324 |
| 325 # Create DataFrame for translation issues | |
| 326 df_translation_issues = pd.DataFrame([ | |
| 327 {"ReactionID": rxn_id, "TranslationIssues": issues} | |
| 328 for rxn_id, issues in translation_issues.items() | |
| 329 ]) | |
| 330 | |
| 228 df_bounds = bounds.reset_index().rename(columns = {"index": "ReactionID"}) | 331 df_bounds = bounds.reset_index().rename(columns = {"index": "ReactionID"}) |
| 229 df_medium = medium.rename(columns = {"reaction": "ReactionID"}) | 332 df_medium = medium.rename(columns = {"reaction": "ReactionID"}) |
| 230 df_medium["InMedium"] = True | 333 df_medium["InMedium"] = True |
| 231 | 334 |
| 232 merged = df_reactions.merge(df_rules, on = "ReactionID", how = "outer") | 335 merged = df_reactions.merge(df_rules, on = "ReactionID", how = "outer") |
| 233 merged = merged.merge(df_bounds, on = "ReactionID", how = "outer") | 336 merged = merged.merge(df_bounds, on = "ReactionID", how = "outer") |
| 234 merged = merged.merge(objective_function, on = "ReactionID", how = "outer") | 337 merged = merged.merge(objective_function, on = "ReactionID", how = "outer") |
| 235 if ARGS.name == "ENGRO2": | 338 if ARGS.name == "ENGRO2": |
| 236 merged = merged.merge(compartments, on = "ReactionID", how = "outer") | 339 merged = merged.merge(compartments, on = "ReactionID", how = "outer") |
| 237 merged = merged.merge(df_medium, on = "ReactionID", how = "left") | 340 merged = merged.merge(df_medium, on = "ReactionID", how = "left") |
| 341 | |
| 342 # Add translation issues column | |
| 343 if not df_translation_issues.empty: | |
| 344 merged = merged.merge(df_translation_issues, on = "ReactionID", how = "left") | |
| 345 merged["TranslationIssues"] = merged["TranslationIssues"].fillna("") | |
| 346 else: | |
| 347 # Add empty TranslationIssues column if no issues found | |
| 348 #merged["TranslationIssues"] = "" | |
| 349 pass | |
| 238 | 350 |
| 239 merged["InMedium"] = merged["InMedium"].fillna(False) | 351 merged["InMedium"] = merged["InMedium"].fillna(False) |
| 240 | 352 |
| 241 merged = merged.sort_values(by = "InMedium", ascending = False) | 353 merged = merged.sort_values(by = "InMedium", ascending = False) |
| 242 | 354 |
