adesomics.cli

  1import argparse
  2from pathlib import Path
  3from xml.parsers.expat import model
  4import numpy as np
  5import pandas as pd
  6from adesomics.io.count_matrix import CountMatrix
  7from adesomics.mapping.connection import Connector
  8from adesomics.dge.statistics import DifferentialExpression
  9from adesomics.population.population_model import PopulationSet, _save
 10from adesomics.io.metadata import Metadata
 11from adesomics.io.gsmm import Model
 12
 13
 14def main():
 15    
 16    parser = argparse.ArgumentParser(
 17        prog="adesomics",
 18        description="adesomics: RNA-seq to FBA pipeline toolkit",
 19    )
 20    parser.add_argument(
 21        "--countmatrix",
 22        metavar="PATH",
 23        help="Path to a count matrix file from featureCounts ().txt/.tsv)",
 24    )
 25    parser.add_argument(
 26        "--model",
 27        metavar="PATH",
 28        help="Path to the metabolic model file (.xml or .SBML)",
 29    )
 30    parser.add_argument(
 31        "--contrast",
 32        nargs=2,
 33        metavar=("TREATED", "REFERENCE")
 34    )
 35    parser.add_argument(
 36        "--stats"
 37    )
 38    parser.add_argument(
 39        "--pop",
 40        action="store_true",
 41        help="Computes the growth rate into a population model and saves respective plots."
 42    )
 43    parser.add_argument(
 44        "--label",
 45        nargs="+",
 46        metavar="LABEL",
 47        help="One label per sample (in count-matrix column order), used in the prints instead of the sample IDs.",
 48    )
 49
 50    src = parser.add_mutually_exclusive_group(required=False)
 51    src.add_argument("--sra", metavar="CSV",
 52                    help="SraRunTable.csv from NCBI Run Selector")
 53    args = parser.parse_args()
 54
 55    if args.countmatrix is None:
 56        path_str = input("Path to count matrix file: ").strip()
 57    else:
 58        path_str = args.countmatrix
 59
 60    path = Path(path_str)
 61
 62    if args.model is None:
 63        model_path_str = input("Path to metabolic model file: ").strip()
 64    else:
 65        model_path_str = args.model
 66
 67    model_path = Path(model_path_str)
 68
 69    print(f"Loading count matrix from '{path}' ...")
 70    try:
 71        cm = CountMatrix.from_file(path)
 72        samples = list(cm.sample_names)
 73        if args.label is not None and len(args.label) != len(samples):
 74            print(f"Error: --label expects {len(samples)} label(s) "
 75                f"(one per sample), got {len(args.label)}.")
 76            raise SystemExit(1)
 77    except FileNotFoundError:
 78        print(f"Error: File not found: '{path}'")
 79        raise SystemExit(1)
 80    except ValueError as e:
 81        print(f"Error: {e}")
 82        raise SystemExit(1)
 83
 84    if args.sra:
 85        md = Metadata.from_sra_run_table(args.sra)
 86        md_samples = list(md.align(cm.counts).groups)
 87    else: 
 88        md = None
 89        md_samples = cm.sample_names
 90
 91    if args.stats:
 92        de = DifferentialExpression(cm, md).fit()
 93        
 94    if args.contrast:
 95        a, b = args.contrast
 96        de = DifferentialExpression(cm, md).fit()
 97        
 98        de.results(a, b).to_csv(f"results/de_{a}_vs_{b}.csv")
 99
100        ax = de.volcano(a, b)
101        _save(ax.figure, f"volcano_{a}_vs_{b}.png", "results")
102
103    print(f"Loading metabolic model from '{model_path}' ...")
104
105    try:
106        model = Model.from_file(model_path)
107    except FileNotFoundError:
108        print(f"Error: File not found: '{model_path}'")
109        raise SystemExit(1)
110    except ValueError as e:
111        print(f"Error: {e}")
112        raise SystemExit(1)
113
114    print(cm.summary())
115    print(model.summary())
116    print("---- Mapping genes between count matrix and metabolic model ----")
117
118    conn  = Connector(model, cm)
119    sol = conn.compare(cm.sample_names)
120
121    print("---- Successfull optimization ----")        
122
123    if args.label is not None:
124        md_samples = args.label
125    for sample, md_sample in zip(cm.sample_names, md_samples):
126        print(f"growth rate mue for {md_sample}: {sol[sample].objective_value}")
127
128    if args.pop:
129        pops = PopulationSet.from_connector(conn, md)
130        ax = pops.plot_growth_curves(t_max=6)  
131        pops.plot_growth_rates() 
132
133    print("-----------------------------------------")
134    print("|   A   D   E   S   O   M   I   C   S   |")
135    print("-----------------------------------------")