text: """Machine Learning Quality Assessment of NGS Data
  
  Main script for seqQscorer. "python seqQscorer.py --help" will display
  a formatted help text on the console. A comprehensive description is provided
  in the GitHub README that includes examples as well. In short, seqQscorer
  uses dedicated machine learning algorithms to train classification models
  to perform automatic NGS quality control.
  
  date:	2020-11-26
  author:	Steffen Albrecht
  
  """
  
  from sys import *
  import os
  import pickle
  import pandas as pd
  import numpy as np
  import json
  import random
  import argparse
  
  import warnings
  warnings.filterwarnings("ignore")
  
  # import project utils
  import utils.Exceptions as myExceptions
  import utils.utils as utils
  import utils.parser as parser
  
  # parse command line arguments
  script_dir = './'
  if argv[0].find('/') >= 0:
  	script_dir = argv[0][: - argv[0][::-1].find('/')]
  utils_dir = '%sutils/'%(script_dir)
  
  argsParser = argparse.ArgumentParser(description='seqQscorer - A machine learning application for quality assessment of NGS data')
  argsParser.add_argument('--indir', '-i', type=str, required=True, help='Input directory containing the feature set files. The feature set files are perfectly fomated by the script "deriveFeatures.py": the file names (until the ".") define the sample ID while the file endings define the corresponding feature set RAW, MAP, LOC, and TSS. By default seqQscorer applies the machine learning model to all samples from the given directory within milliseconds. However, it can be restricted to one sample using --sampleID.')
  argsParser.add_argument('--species', '-s', type=str, default='generic', 
  						choices=['generic','human', 'mouse'],  help='Species specifying the model used.')
  argsParser.add_argument('--assay', '-a', type=str, default='generic', 
  						choices=['generic','ChIP-seq','DNase-seq','RNA-seq'], help='Assay specifying the model used.')
  argsParser.add_argument('--runtype', '-r', type=str, default='generic',
  						choices=['generic','single-end','paired-end'], help='Run-Type specifying the model used.')
  argsParser.add_argument('--model', '-m', type=str, default=None, help='Path to a serialized model, trained on own data. If used, the parameters --species, --assay, and --runtype have no impact on the classification model.')
  argsParser.add_argument('--noRAW', action='store_true', help='Ignore all RAW features.')
  argsParser.add_argument('--noMAP', action='store_true', help='Ignore all MAP features.')
  argsParser.add_argument('--noLOC', action='store_true', help='Ignore all LOC features.')
  argsParser.add_argument('--noTSS', action='store_true', help='Ignore all TSS features.')
  argsParser.add_argument('--noFS', action='store_true', help='Switch off feature selection. (has only an impact if the best performance was achieved with chi2 or RFE)')
  argsParser.add_argument('--bestCalib', action='store_true', help='Classifier setting is used that achieved the lowest brier score, hence the best calibration of the probabilities.')
  argsParser.add_argument('--peaktype', '-pt', type=str, default=None, choices=['narrow','broad'], help='Optionally specify the peak-type for ChIP-seq data.')
  argsParser.add_argument('--probOut', '-po', type=str, default=None,
  						help='To specify an output file for the probabilities. Output will be tab-separated.')
  argsParser.add_argument('--compOut', '-co', type=str, default=None,
  						help='To specify an out file for the comprehensive output. Output will be kind of tab-separated.')
  argsParser.add_argument('--inputOut', '-io', type=str, default=None,
  						help='To specify an out file that will contain the parsed input. Output will be tab-separated.')
  argsParser.add_argument('--noVerbose', '-nv', action='store_true', help='Turn off verboseness, without being quiet.')
  argsParser.add_argument('--seed', '-rs', type=int, default=1, help='Some classifiers apply randomization. Use --seed to make results reproducible. By default the seed 1 is used, set it to -1 if using a seed is not desired. For K-nearest neighbor and Naive Bayes the seed has no impact.')
  argsParser.add_argument('--sampleID', '-id', type=str, default=None,
  						help='Restrict application of seqQscorer to only one sample defined by the ID.')
  
  args = argsParser.parse_args()
  
  if not os.path.isdir(args.indir):
  	raise myExceptions.WrongFeatureInputException(
  						'"%s" is not a directory'%(args.indir))
  
  feature_sets = ['RAW','MAP','LOC','TSS']
  
  # restrict feature sets used according to given optional parameters
  if args.noRAW:
  	feature_sets.remove('RAW')
  if args.noMAP:
  	feature_sets.remove('MAP')
  if args.noLOC:
  	feature_sets.remove('LOC')
  if args.noTSS:
  	feature_sets.remove('TSS')
  
  model_sel_metric = 'brier' if args.bestCalib else 'auROC'
  fs_suffix = '_noFS' if args.noFS else ''
  
  # initiate the classification model and other data needed
  species, assay, run_type = args.species, args.assay, args.runtype
  
  if args.peaktype != None:
  	if args.assay != 'ChIP-seq' or args.runtype != 'single-end':
  		raise myExceptions.WrongSettingException(
  			'Peak-type specification can only be used for single-ended ChIP-seq.')
  	assay = args.peaktype + assay
  
  best_clf, feature_selection, selection, parameters, auROC, brier = utils.get_best_classifier(utils_dir, species, assay, 
  															run_type, feature_sets, fs_suffix, model_sel_metric)
  
  if best_clf == None and args.model == None:
  	message = '''\nPlease check the given setting:
  	assay:\t\t%s\n\tspecies:\t%s\n\trun-type:\t%s\n'''%(assay, species, run_type)
  	message += '\tfeature sets:\t%s\n'%('-'.join(feature_sets))
  	message += 'A specialized model for this setting is not available,\n'
  	message += 'the generic model is used to proceed.\n'
  	print(message)
  	species, assay, run_type = 'generic', 'generic', 'generic'
  
  application_case = '%s_%s_%s_%s'%(species, assay, run_type, '-'.join(feature_sets))
  application_case += '_%s%s'%(model_sel_metric, fs_suffix)
  model_file_path = '%smodels/%s_%d.model'%(script_dir, application_case, args.seed)
  
  if args.model != None:
  	print('An external model is provided.')
  	model_file_path = args.model
  	try: 
  		pickle.load(open(model_file_path, 'rb'))
  	except:
  		raise myExceptions.IncorrectModelException(
  			'The provided model from file "%s" could not be loaded.'%(model_file_path))
  
  # load median values organized by subset, needed to impute missing values
  medians = pickle.load(open('%sutils/medians.dict'%(script_dir), 'rb'))
  medians = medians[species][assay][run_type]
  
  # parse given input files
  input_data, feature_columns = parser.generate_input_data(args.indir, feature_sets, run_type, medians, args.noVerbose, args.sampleID)
  
  # if the particula model is used for the very first time, it is trained and serialized
  best_clf, feature_selection, selection, parameters, auROC, brier = utils.get_best_classifier(utils_dir, 
  																	 species, assay, run_type, feature_sets, 
  																	 fs_suffix, model_sel_metric)
  
  if not os.path.exists(model_file_path) and args.model == None:
  	print('\nThe required model was not used so far.')
  	print('It needs to be trained and serialized...')
  	
  	clf = utils.get_clf_algos()[best_clf]
  	if not best_clf in ['GNB','KNN']:
  		if args.seed != -1:
  			parameters['random_state'] = args.seed
  	
  	clf_setup = clf.set_params(**parameters)
  	
  	data_file_path = '%sutils/datasets/%s_%s_%s.tsv'%(script_dir, assay, species, run_type)
  	train_data = pd.read_csv(data_file_path, sep='\t')
  	
  	y = np.array(train_data['status'])
  	train_data = train_data[feature_columns]
  	if selection != None:
  		train_data = train_data.loc[:,selection]
  	X = np.array(train_data)
  	
  	model = clf_setup.fit(X,y)
  	if args.seed != -1:
  		pickle.dump(model, open(model_file_path, 'wb'))
  		print('... training and serialization is done!')
  		print('The model is instantly available from now!')
  	else:
  		print('... training is done, but only reproducible models are serialized.')
  		print('Because no seed was used the model was not serialized.')
  
  # load the model
  model = pickle.load(open(model_file_path, 'rb'))
  
  # prepare input data format
  input_values = input_data[feature_columns]
  if selection != None:
  	input_values = input_values.loc[:,selection]
  
  # apply model on given samples to get the probabilities
  probabilities = model.predict_proba(np.array(input_values))
  fileIDs = list(input_data['sampleID'])
  fileID_score = list(zip(fileIDs, [prob[1] for prob in probabilities]))
  
  if args.model == None:
  	print('\nThe best predictive performance was achived by %s'%(utils.clf_full_names(best_clf)))
  	print('%s feature selection is applied (using %s%s of the features)'%(
  		tuple(feature_selection.split('-') + ['%'])))
  	print('Within the cross-validated gird search this model achived:')
  	print('\tauROC: %s'%(auROC))
  	print('\tBrier: %s'%(brier))
  	print('\nThe model used, achieved these measures for different decision thresholds applied\non the probabilities within the grid-search (ten-fold cross-validation):')
  	table = utils.read_in_measure_table(utils_dir, species, assay, run_type, 
  										feature_sets, fs_suffix, model_sel_metric)
  	utils.print_nice_table(table)
  	print('')
  
  # print the scores to the console
  probas_str = ''
  for fileID, score in sorted(fileID_score, key=lambda x: x[1]):
  	probas_str += '%s\t%f\n'%(fileID, score)
  	print(fileID, '%.3f '%(score), '(probability for being of low quality)', sep='\t')
  
  # write probabilities to file if a file-path is given
  if args.probOut != None:
  	try:
  		open(args.probOut, 'w').write(probas_str)
  	except:
  		raise myExceptions.WrongOutputFileException(
  			'Unable to write the probabilities to file!')
  
  # write comprehensive output to file if a file-path is given
  if args.compOut != None:
  	comp_out  = 'Model trained by: %s\n'%(utils.clf_full_names(best_clf))
  	comp_out += '%s feature selection applied\n'%(feature_selection.split('-')[0])
  	comp_out += '%s %s of the features are used\n'%(feature_selection.split('-')[1], '%')
  	comp_out += 'auROC: %s\n'%(auROC)
  	comp_out += 'Brier: %s\n\n'%(brier)
  	comp_out += 'Metric table:\n'
  	for row in table:
  		comp_out += '\t'.join(row) + '\n'
  	comp_out += '\n' + probas_str
  	
  	try:
  		open(args.compOut, 'w').write(comp_out)
  	except:
  		raise myExceptions.WrongOutputFileException(
  			'Unable to write the comprehensive output to file!')
  
  # write the parsed input into a file, if a pth is given
  if args.inputOut != None:
  	try:
  		input_data.to_csv(args.inputOut, sep='\t', index=False)
  	except:
  		raise myExceptions.WrongOutputFileException(
  			'Unable to write the parsed input to file')
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  