1+ #!/usr/bin/env python
2+
3+ ###############################################################################
4+ # #
5+ # This program is free software: you can redistribute it and/or modify #
6+ # it under the terms of the GNU General Public License as published by #
7+ # the Free Software Foundation, either version 3 of the License, or #
8+ # (at your option) any later version. #
9+ # #
10+ # This program is distributed in the hope that it will be useful, #
11+ # but WITHOUT ANY WARRANTY; without even the implied warranty of #
12+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
13+ # GNU General Public License for more details. #
14+ # #
15+ # You should have received a copy of the GNU General Public License #
16+ # along with this program. If not, see <http://www.gnu.org/licenses/>. #
17+ # #
18+ ###############################################################################
19+
20+ __prog_name__ = 'generate_training_set.py'
21+ __prog_desc__ = ('Generate the a new training set based on a gtdb taxonomy. Parse the metadata file from gtdb,'
22+ 'select the representative genomes and up to 99 additional genomes per species. For this set of genomes '
23+ 'The script will calculate the feature values - coding density, GC content, Trp and Gly ratios - and save the resulting table for training the classifiers.'
24+ 'The script will aslo create a file use to pick the ground truth Translation table.' )
25+
26+ __author__ = 'Pierre Chaumeil'
27+ __copyright__ = 'Copyright 2026'
28+ __credits__ = ['Pierre Chaumeil' ]
29+ __license__ = 'GPL3'
30+ __version__ = '0.0.3'
31+ __maintainer__ = 'Pierre Chaumeil'
32+ __email__ = 'uqpchaum@uq.edu.au'
33+ __status__ = 'Development'
34+
35+ import argparse
36+ import os
37+ import sys
38+
39+ import joblib
40+ import numpy as np
41+ import pandas as pd
42+
43+ # Sklearn Imports
44+ from sklearn .ensemble import AdaBoostClassifier , RandomForestClassifier
45+ from sklearn .metrics import balanced_accuracy_score
46+ from sklearn .model_selection import StratifiedKFold , train_test_split
47+ from sklearn .preprocessing import MinMaxScaler , LabelEncoder
48+ from sklearn .tree import DecisionTreeClassifier
49+ from sklearn .neighbors import KNeighborsClassifier
50+ from sklearn .pipeline import Pipeline
51+ from sklearn .utils import compute_class_weight
52+ from lightgbm import LGBMClassifier
53+
54+ from mlxtend .feature_selection import SequentialFeatureSelector
55+ from script_logger import CustomLogger
56+
57+ # Import classifier configurations
58+ from classifier_parameters import classifier_configs
59+
60+ # prevents warnings from being displayed.
61+ import warnings
62+
63+ warnings .simplefilter ("ignore" )
64+
65+
66+ class TrainingSetGenerator (object ):
67+ """
68+ A class to merge data, perform feature engineering, and train multiple
69+ classifiers with pipeline-based feature scaling and model persistence.
70+ """
71+
72+ def __init__ (self , feature_file , output_dir , threads_count , seed = None ):
73+
74+ self .feature_file = feature_file
75+ self .output_dir = output_dir
76+ self .thread_count = threads_count
77+ self .seed = seed
78+
79+ # Initialize the logger
80+ logger_instance = CustomLogger (output_dir , __prog_name__ )
81+ self .logger = logger_instance .get_logger ()
82+
83+ # make sure the output directory exists
84+ if not os .path .exists (self .output_dir ):
85+ os .makedirs (self .output_dir )
86+
87+ # Map for instantiating raw models
88+ self .model_mapping = {
89+ "KNeighbors" : KNeighborsClassifier (n_jobs = self .thread_count ),
90+ "RandomForest" : RandomForestClassifier (random_state = self .seed , n_jobs = self .thread_count ),
91+ "AdaBoost" : AdaBoostClassifier (random_state = self .seed ),
92+ "LGBM" : LGBMClassifier (random_state = self .seed , n_jobs = self .thread_count , verbose = - 1 ),
93+ "DecisionTree" : DecisionTreeClassifier (random_state = self .seed )
94+ }
95+
96+ def load_and_preprocess_data (self ):
97+ """Loads the 4 input TSVs, merges them, and engineers features."""
98+ self .logger .info ("Loading and merging data files..." )
99+
100+ try :
101+ df_f = pd .read_csv (self .feature_file , sep = '\t ' )
102+ except FileNotFoundError as e :
103+ self .logger .error (f"Error loading files: { e } " )
104+ raise
105+
106+ df_f ['Genome' ] = df_f ['Genome' ].astype (str ).str .strip ()
107+
108+
109+ # Filter Target
110+ valid_targets = ['11' , '4' , '25' ]
111+ df = df_f [df_f ['Ground_truth' ].astype (str ).isin (valid_targets )].copy ()
112+
113+ self .logger .info (f"Merged dataset shape after target filtering: { df .shape } " )
114+
115+ # Feature Engineering
116+ df ['Coding_density_4' ] = df ['Coding_density_4' ]
117+ df ['Coding_density_11' ] = df ['Coding_density_11' ]
118+ df ['Density_Diff' ] = df ['Coding_density_4' ] - df ['Coding_density_11' ]
119+
120+ df ['Trp_ratio' ] = np .log ((df ['tt4_uga_count' ] + 1 ) / (df ['tt4_ugg_count' ] + 1 )).clip (- 6.0 , 5.0 )
121+ df ['Trp_magnitude' ] = np .log (df ['tt4_uga_count' ] + df ['tt4_ugg_count' ] + 1 )
122+ df ['Gly_ratio' ] = np .log ((df ['tt4_uga_count' ] + 1 ) / (df ['tt4_std_gly_count' ] + 1 )).clip (- 10.0 , 0.0 )
123+
124+ #df['feature_uga_density'] = df['tt4_uga_count'] / (df['Coding_density_4'] + 0.001)
125+ df ['UGG_density' ] = df ['tt4_ugg_count' ] / df ['tt4_std_gly_count' ]
126+ #df['feature_uga_vs_gly_pool'] = df['tt4_uga_count'] / (df['tt4_std_gly_count'] + df['tt4_uga_count'] + 1)
127+
128+ self .feature_cols = [
129+ 'Coding_density_4' , 'Coding_density_11' , 'Density_Diff' , 'GC' ,
130+ 'Trp_ratio' , 'Trp_magnitude' , 'Gly_ratio' ,
131+ 'UGG_density'
132+ ]
133+
134+ # Handle NaNs from mathematical operations
135+ if df [self .feature_cols ].isnull ().any ().any ():
136+ self .logger .warning ("NaN values generated during feature engineering. Filling with 0." )
137+ df [self .feature_cols ] = df [self .feature_cols ].fillna (0 )
138+
139+ return df
140+
141+ def run (self , split_data = False ):
142+ """Execute the full training pipeline across multiple classifiers."""
143+ self .logger .info (f'{ __prog_name__ } { " " .join (sys .argv [1 :])} ' )
144+
145+ seed_to_use = self .seed if self .seed is not None else np .random .randint (0 , 1000 )
146+ self .logger .info (f"We use seed: { seed_to_use } " )
147+
148+ # 1. Load & Preprocess
149+ df = self .load_and_preprocess_data ()
150+
151+
152+
153+
154+
155+ # 2. Prepare Data
156+ X = df [self .feature_cols ]
157+ labels_raw = df ['Ground_truth' ].astype (str )
158+
159+ # lets show the first few rows of the dataframe to check if everything is correct
160+ self .logger .info ("First few rows of the preprocessed dataframe:" )
161+ self .logger .info (X .head ().to_string ())
162+
163+ # Encode multi-class labels safely
164+ le = LabelEncoder ()
165+ y = le .fit_transform (labels_raw )
166+
167+ # Save LabelEncoder to decode predictions later
168+ joblib .dump (le , os .path .join (self .output_dir , 'label_encoder.pkl' ))
169+
170+ # count class distribution
171+ class_counts = np .bincount (y )
172+ self .logger .info (f"Class counts: { class_counts } " )
173+
174+ class_weights = compute_class_weight (class_weight = 'balanced' , classes = np .unique (y ), y = y )
175+ weights_dict = dict (zip (le .classes_ , class_weights ))
176+ self .logger .info (f"Class weights mapping: { weights_dict } " )
177+
178+ if split_data :
179+ self .logger .info ("Splitting the data into training and validation sets" )
180+ X_train , X_val , y_train , y_val = train_test_split (X , y , test_size = 0.2 , stratify = y ,
181+ random_state = seed_to_use )
182+ self .logger .info (f"Training data shape: { X_train .shape } " )
183+ else :
184+ self .logger .info ("Using the full training data" )
185+ X_train , X_val , y_train , y_val = X , X , y , y
186+
187+ # 3. Iterate Over Configurations
188+ for config in classifier_configs :
189+ model_name = config ["name" ]
190+ short_name = config ["short_name" ]
191+ best_params = config ["best_params" ]
192+
193+ self .logger .info (f"--- Training { model_name } ---" )
194+
195+ base_model = self .model_mapping .get (model_name )
196+ if base_model is None :
197+ self .logger .warning (f"Model { model_name } not found in mapping. Skipping." )
198+ continue
199+
200+ # Build Pipeline (MinMaxScaler for distance-based models)
201+ if model_name in ['KNeighbors' ]:
202+ pipe = Pipeline ([
203+ ('scaler' , MinMaxScaler (feature_range = (0 , 1 ))),
204+ ('model' , base_model )
205+ ])
206+ else :
207+ pipe = Pipeline ([
208+ ('model' , base_model )
209+ ])
210+
211+ # Apply best parameters dynamically
212+ pipe .set_params (** best_params )
213+
214+ X_train_current = X_train .copy ()
215+ X_val_current = X_val .copy ()
216+
217+ # Fit the pipeline
218+ pipe .fit (X_train_current , y_train )
219+
220+ # Save the fitted pipeline
221+ model_out_path = os .path .join (self .output_dir , f'{ short_name } _multi_class.pkl' )
222+ try :
223+ joblib .dump (pipe , model_out_path )
224+ self .logger .info (f"Saved { model_name } pipeline to { model_out_path } " )
225+ except Exception as e :
226+ self .logger .error (f"Error saving { model_name } : { e } " )
227+ continue
228+
229+ # Validation
230+ predictions = pipe .predict (X_val_current )
231+ balanced_accuracy = balanced_accuracy_score (y_val , predictions )
232+ self .logger .info (f"Balanced accuracy for { model_name } : { balanced_accuracy :.4f} \n " )
233+
234+
235+ if __name__ == "__main__" :
236+ print (__prog_name__ + ' v' + __version__ + ': ' + __prog_desc__ )
237+ print (' by ' + __author__ + ' (' + __email__ + ')' + '\n ' )
238+
239+ parser = argparse .ArgumentParser (formatter_class = argparse .ArgumentDefaultsHelpFormatter )
240+
241+ # Input files matching your execution loop script
242+ parser .add_argument ('--metadata_file_bac' , required = True , help = 'GTDB metadata file for bacteria (tsv).' )
243+ parser .add_argument ('--metadata_file_arc' , required = True , help = 'GTDB metadata file for archaea (tsv).' )
244+ parser .add_argument ('--genome_path_file' , required = True , help = 'File containing paths to genome assemblies (tsv).' )
245+ parser .add_argument ('--output_dir' , '-od' , required = True , help = 'Directory to save the features and ground_truth tables.' )
246+ args = parser .parse_args ()
247+
248+ try :
249+ classifier = TrainingSetGenerator (
250+ feature_file = args .features ,
251+ output_dir = args .od ,
252+ threads_count = args .threads ,
253+ seed = args .seed
254+ )
255+ classifier .run (args .split_data )
256+ except Exception as e :
257+ print (f"ERROR: { str (e )} " , file = sys .stderr )
258+ sys .exit (1 )
0 commit comments