cstrittmatter@0: #!/usr/bin/env python3
cstrittmatter@0:
cstrittmatter@0: # -*- coding: utf-8 -*-
cstrittmatter@0:
cstrittmatter@0: """
cstrittmatter@0: rematch.py - Reads mapping against target sequences, checking mapping
cstrittmatter@0: and consensus sequences production
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: Copyright (C) 2019 Miguel Machado
cstrittmatter@0:
cstrittmatter@0: Last modified: August 08, 2019
cstrittmatter@0:
cstrittmatter@0: This program is free software: you can redistribute it and/or modify
cstrittmatter@0: it under the terms of the GNU General Public License as published by
cstrittmatter@0: the Free Software Foundation, either version 3 of the License, or
cstrittmatter@0: (at your option) any later version.
cstrittmatter@0:
cstrittmatter@0: This program is distributed in the hope that it will be useful,
cstrittmatter@0: but WITHOUT ANY WARRANTY; without even the implied warranty of
cstrittmatter@0: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
cstrittmatter@0: GNU General Public License for more details.
cstrittmatter@0:
cstrittmatter@0: You should have received a copy of the GNU General Public License
cstrittmatter@0: along with this program. If not, see .
cstrittmatter@0: """
cstrittmatter@0:
cstrittmatter@0: import os
cstrittmatter@0: import sys
cstrittmatter@0: import time
cstrittmatter@0: import argparse
cstrittmatter@0:
cstrittmatter@0: try:
cstrittmatter@0: from __init__ import __version__
cstrittmatter@0:
cstrittmatter@0: import modules.utils as utils
cstrittmatter@0: import modules.seqFromWebTaxon as seq_from_web_taxon
cstrittmatter@0: import modules.download as download
cstrittmatter@0: import modules.rematch_module as rematch_module
cstrittmatter@0: import modules.checkMLST as check_mlst
cstrittmatter@0: except ImportError:
cstrittmatter@0: from ReMatCh.__init__ import __version__
cstrittmatter@0:
cstrittmatter@0: from ReMatCh.modules import utils as utils
cstrittmatter@0: from ReMatCh.modules import seqFromWebTaxon as seq_from_web_taxon
cstrittmatter@0: from ReMatCh.modules import download as download
cstrittmatter@0: from ReMatCh.modules import rematch_module as rematch_module
cstrittmatter@0: from ReMatCh.modules import checkMLST as check_mlst
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def search_fastq_files(directory):
cstrittmatter@0: files_extensions = ['.fastq.gz', '.fq.gz']
cstrittmatter@0: pair_end_files_separation = [['_R1_001.f', '_R2_001.f'], ['_1.f', '_2.f']]
cstrittmatter@0:
cstrittmatter@0: list_ids = {}
cstrittmatter@0: directories = [d for d in os.listdir(directory) if
cstrittmatter@0: not d.startswith('.') and os.path.isdir(os.path.join(directory, d, ''))]
cstrittmatter@0: for directory_found in directories:
cstrittmatter@0: if directory_found != 'pubmlst':
cstrittmatter@0: directory_path = os.path.join(directory, directory_found, '')
cstrittmatter@0:
cstrittmatter@0: fastq_found = []
cstrittmatter@0: files = [f for f in os.listdir(directory_path) if
cstrittmatter@0: not f.startswith('.') and os.path.isfile(os.path.join(directory_path, f))]
cstrittmatter@0: for file_found in files:
cstrittmatter@0: if file_found.endswith(tuple(files_extensions)):
cstrittmatter@0: fastq_found.append(file_found)
cstrittmatter@0:
cstrittmatter@0: if len(fastq_found) == 1:
cstrittmatter@0: list_ids[directory_found] = [os.path.join(directory_path, f) for f in fastq_found]
cstrittmatter@0: elif len(fastq_found) >= 2:
cstrittmatter@0: file_pair = []
cstrittmatter@0:
cstrittmatter@0: # Search pairs
cstrittmatter@0: for pe_separation in pair_end_files_separation:
cstrittmatter@0: for fastq in fastq_found:
cstrittmatter@0: if pe_separation[0] in fastq or pe_separation[1] in fastq:
cstrittmatter@0: file_pair.append(fastq)
cstrittmatter@0:
cstrittmatter@0: if len(file_pair) == 2:
cstrittmatter@0: break
cstrittmatter@0: else:
cstrittmatter@0: file_pair = []
cstrittmatter@0:
cstrittmatter@0: # Search single
cstrittmatter@0: if len(file_pair) == 0:
cstrittmatter@0: for pe_separation in pair_end_files_separation:
cstrittmatter@0: for fastq in fastq_found:
cstrittmatter@0: if pe_separation[0] not in fastq or pe_separation[1] not in fastq:
cstrittmatter@0: file_pair.append(fastq)
cstrittmatter@0:
cstrittmatter@0: if len(file_pair) >= 1:
cstrittmatter@0: file_pair = file_pair[0]
cstrittmatter@0:
cstrittmatter@0: if len(file_pair) >= 1:
cstrittmatter@0: list_ids[directory_found] = [os.path.join(directory_path, f) for f in file_pair]
cstrittmatter@0:
cstrittmatter@0: return list_ids
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def get_list_ids_from_file(file_list_ids):
cstrittmatter@0: list_ids = []
cstrittmatter@0:
cstrittmatter@0: with open(file_list_ids, 'rtU') as lines:
cstrittmatter@0: for line in lines:
cstrittmatter@0: line = line.rstrip('\r\n')
cstrittmatter@0: if len(line) > 0:
cstrittmatter@0: list_ids.append(line)
cstrittmatter@0:
cstrittmatter@0: if len(list_ids) == 0:
cstrittmatter@0: sys.exit('No runIDs were found in ' + file_list_ids)
cstrittmatter@0:
cstrittmatter@0: return list_ids
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def get_taxon_run_ids(taxon_name, outputfile):
cstrittmatter@0: seq_from_web_taxon.run_seq_from_web_taxon(taxon_name, outputfile, True, True, True, False)
cstrittmatter@0:
cstrittmatter@0: run_ids = []
cstrittmatter@0: with open(outputfile, 'rtU') as reader:
cstrittmatter@0: for line in reader:
cstrittmatter@0: line = line.rstrip('\r\n')
cstrittmatter@0: if len(line) > 0:
cstrittmatter@0: if not line.startswith('#'):
cstrittmatter@0: line = line.split('\t')
cstrittmatter@0: run_ids.append(line[0])
cstrittmatter@0:
cstrittmatter@0: return run_ids
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def get_list_ids(workdir, file_list_ids, taxon_name):
cstrittmatter@0: searched_fastq_files = False
cstrittmatter@0: list_ids = []
cstrittmatter@0: if file_list_ids is None and taxon_name is None:
cstrittmatter@0: list_ids = search_fastq_files(workdir)
cstrittmatter@0: searched_fastq_files = True
cstrittmatter@0: elif file_list_ids is not None:
cstrittmatter@0: list_ids = get_list_ids_from_file(os.path.abspath(file_list_ids))
cstrittmatter@0: elif taxon_name is not None and file_list_ids is None:
cstrittmatter@0: list_ids = get_taxon_run_ids(taxon_name, os.path.join(workdir, 'IDs_list.seqFromWebTaxon.tab'))
cstrittmatter@0:
cstrittmatter@0: if len(list_ids) == 0:
cstrittmatter@0: sys.exit('No IDs were found')
cstrittmatter@0: return list_ids, searched_fastq_files
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def format_gene_info(gene_specific_info, minimum_gene_coverage, minimum_gene_identity, reported_data_type, summary,
cstrittmatter@0: sample, genes_present):
cstrittmatter@0: info = None
cstrittmatter@0: if gene_specific_info['gene_coverage'] >= minimum_gene_coverage and \
cstrittmatter@0: gene_specific_info['gene_identity'] >= minimum_gene_identity:
cstrittmatter@0: if summary and sample not in genes_present:
cstrittmatter@0: genes_present[sample] = {}
cstrittmatter@0:
cstrittmatter@0: if gene_specific_info['gene_number_positions_multiple_alleles'] == 0:
cstrittmatter@0: s = str(gene_specific_info[reported_data_type])
cstrittmatter@0: info = str(s)
cstrittmatter@0: if summary:
cstrittmatter@0: genes_present[sample][gene_specific_info['header']] = str(s)
cstrittmatter@0: else:
cstrittmatter@0: s = 'multiAlleles_' + str(gene_specific_info[reported_data_type])
cstrittmatter@0: info = str(s)
cstrittmatter@0: if summary:
cstrittmatter@0: genes_present[sample][gene_specific_info['header']] = str(s)
cstrittmatter@0: else:
cstrittmatter@0: info = 'absent_' + str(gene_specific_info[reported_data_type])
cstrittmatter@0:
cstrittmatter@0: return info, genes_present
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def write_data_by_gene(gene_list_reference, minimum_gene_coverage, sample, data_by_gene, outdir, time_str, run_times,
cstrittmatter@0: minimum_gene_identity, reported_data_type, summary, genes_present):
cstrittmatter@0: combined_report = \
cstrittmatter@0: os.path.join(outdir,
cstrittmatter@0: 'combined_report.data_by_gene.' + run_times + '.' + reported_data_type + '.' + time_str + '.tab')
cstrittmatter@0:
cstrittmatter@0: if reported_data_type == 'coverage_depth':
cstrittmatter@0: reported_data_type = 'gene_mean_read_coverage'
cstrittmatter@0: elif reported_data_type == 'sequence_coverage':
cstrittmatter@0: reported_data_type = 'gene_coverage'
cstrittmatter@0:
cstrittmatter@0: combined_report_exist = os.path.isfile(combined_report)
cstrittmatter@0: with open(combined_report, 'at') as writer:
cstrittmatter@0: seq_list = sorted(gene_list_reference.keys())
cstrittmatter@0: if not combined_report_exist:
cstrittmatter@0: writer.write('#sample' + '\t' + '\t'.join([gene_list_reference[seq] for seq in seq_list]) + '\n')
cstrittmatter@0:
cstrittmatter@0: results = {}
cstrittmatter@0: headers = []
cstrittmatter@0: for i in data_by_gene:
cstrittmatter@0: results[data_by_gene[i]['header']], genes_present = format_gene_info(data_by_gene[i], minimum_gene_coverage,
cstrittmatter@0: minimum_gene_identity,
cstrittmatter@0: reported_data_type, summary, sample,
cstrittmatter@0: genes_present)
cstrittmatter@0: headers.append(data_by_gene[i]['header'])
cstrittmatter@0:
cstrittmatter@0: if len(headers) != gene_list_reference:
cstrittmatter@0: for gene in gene_list_reference:
cstrittmatter@0: if gene not in headers:
cstrittmatter@0: results[gene] = 'NA'
cstrittmatter@0:
cstrittmatter@0: writer.write(sample + '\t' + '\t'.join([results[seq] for seq in seq_list]) + '\n')
cstrittmatter@0:
cstrittmatter@0: return genes_present
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def write_sample_report(sample, outdir, time_str, file_size, run_successfully_fastq, run_successfully_rematch_first,
cstrittmatter@0: run_successfully_rematch_second, time_taken_fastq, time_taken_rematch_first,
cstrittmatter@0: time_taken_rematch_second, time_taken_sample, sequencing_information, sample_data_general_first,
cstrittmatter@0: sample_data_general_second, fastq_used):
cstrittmatter@0: sample_report = os.path.join(outdir, 'sample_report.' + time_str + '.tab')
cstrittmatter@0: report_exist = os.path.isfile(sample_report)
cstrittmatter@0:
cstrittmatter@0: header_general = ['sample', 'sample_run_successfully', 'sample_run_time', 'files_size', 'download_run_successfully',
cstrittmatter@0: 'download_run_time', 'rematch_run_successfully_first', 'rematch_run_time_first',
cstrittmatter@0: 'rematch_run_successfully_second', 'rematch_run_time_second']
cstrittmatter@0: header_data_general = ['number_absent_genes', 'number_genes_multiple_alleles', 'mean_sample_coverage']
cstrittmatter@0: header_sequencing = ['run_accession', 'instrument_platform', 'instrument_model', 'library_layout', 'library_source',
cstrittmatter@0: 'extra_run_accession', 'nominal_length', 'read_count', 'base_count', 'date_download']
cstrittmatter@0:
cstrittmatter@0: with open(sample_report, 'at') as writer:
cstrittmatter@0: if not report_exist:
cstrittmatter@0: writer.write('#' + '\t'.join(header_general) + '\t' + '_first\t'.join(header_data_general) + '_first\t' +
cstrittmatter@0: '_second\t'.join(header_data_general) + '_second\t' + '\t'.join(header_sequencing) + '\t' +
cstrittmatter@0: 'fastq_used' + '\n')
cstrittmatter@0:
cstrittmatter@0: writer.write('\t'.join([sample,
cstrittmatter@0: str(all([run_successfully_fastq is not False,
cstrittmatter@0: run_successfully_rematch_first is not False,
cstrittmatter@0: run_successfully_rematch_second is not False])),
cstrittmatter@0: str(time_taken_sample),
cstrittmatter@0: str(file_size),
cstrittmatter@0: str(run_successfully_fastq),
cstrittmatter@0: str(time_taken_fastq),
cstrittmatter@0: str(run_successfully_rematch_first),
cstrittmatter@0: str(time_taken_rematch_first),
cstrittmatter@0: str(run_successfully_rematch_second),
cstrittmatter@0: str(time_taken_rematch_second)]) +
cstrittmatter@0: '\t' + '\t'.join([str(sample_data_general_first[i]) for i in header_data_general]) +
cstrittmatter@0: '\t' + '\t'.join([str(sample_data_general_second[i]) for i in header_data_general]) +
cstrittmatter@0: '\t' + '\t'.join([str(sequencing_information[i]) for i in header_sequencing]) +
cstrittmatter@0: '\t' + ','.join(fastq_used) + '\n')
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def concatenate_extra_seq_2_consensus(consensus_sequence, reference_sequence, extra_seq_length, outdir):
cstrittmatter@0: reference_dict, ignore, ignore = rematch_module.get_sequence_information(reference_sequence, extra_seq_length)
cstrittmatter@0: consensus_dict, genes, ignore = rematch_module.get_sequence_information(consensus_sequence, 0)
cstrittmatter@0: number_consensus_with_sequences = 0
cstrittmatter@0: for k, values_consensus in list(consensus_dict.items()):
cstrittmatter@0: for values_reference in list(reference_dict.values()):
cstrittmatter@0: if values_reference['header'] == values_consensus['header']:
cstrittmatter@0: if len(set(consensus_dict[k]['sequence'])) > 1:
cstrittmatter@0: number_consensus_with_sequences += 1
cstrittmatter@0: if extra_seq_length <= len(values_reference['sequence']):
cstrittmatter@0: right_extra_seq = \
cstrittmatter@0: '' if extra_seq_length == 0 else values_reference['sequence'][-extra_seq_length:]
cstrittmatter@0: consensus_dict[k]['sequence'] = \
cstrittmatter@0: values_reference['sequence'][:extra_seq_length] + \
cstrittmatter@0: consensus_dict[k]['sequence'] + \
cstrittmatter@0: right_extra_seq
cstrittmatter@0: consensus_dict[k]['length'] += extra_seq_length + len(right_extra_seq)
cstrittmatter@0:
cstrittmatter@0: consensus_concatenated = os.path.join(outdir, 'consensus_concatenated_extraSeq.fasta')
cstrittmatter@0: with open(consensus_concatenated, 'wt') as writer:
cstrittmatter@0: for i in consensus_dict:
cstrittmatter@0: writer.write('>' + consensus_dict[i]['header'] + '\n')
cstrittmatter@0: fasta_sequence_lines = rematch_module.chunkstring(consensus_dict[i]['sequence'], 80)
cstrittmatter@0: for line in fasta_sequence_lines:
cstrittmatter@0: writer.write(line + '\n')
cstrittmatter@0:
cstrittmatter@0: return consensus_concatenated, genes, consensus_dict, number_consensus_with_sequences
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def clean_headers_reference_file(reference_file, outdir, extra_seq):
cstrittmatter@0: problematic_characters = ["|", " ", ",", ".", "(", ")", "'", "/", ":"]
cstrittmatter@0: print('Checking if reference sequences contain ' + str(problematic_characters) + '\n')
cstrittmatter@0: # headers_changed = False
cstrittmatter@0: new_reference_file = str(reference_file)
cstrittmatter@0: sequences, genes, headers_changed = rematch_module.get_sequence_information(reference_file, extra_seq)
cstrittmatter@0: if headers_changed:
cstrittmatter@0: print('At least one of the those characters was found. Replacing those with _' + '\n')
cstrittmatter@0: new_reference_file = \
cstrittmatter@0: os.path.join(outdir, os.path.splitext(os.path.basename(reference_file))[0] + '.headers_renamed.fasta')
cstrittmatter@0: with open(new_reference_file, 'wt') as writer:
cstrittmatter@0: for i in sequences:
cstrittmatter@0: writer.write('>' + sequences[i]['header'] + '\n')
cstrittmatter@0: fasta_sequence_lines = rematch_module.chunkstring(sequences[i]['sequence'], 80)
cstrittmatter@0: for line in fasta_sequence_lines:
cstrittmatter@0: writer.write(line + '\n')
cstrittmatter@0: return new_reference_file, genes, sequences
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def write_mlst_report(sample, run_times, consensus_type, st, alleles_profile, loci_order, outdir, time_str):
cstrittmatter@0: mlst_report = os.path.join(outdir, 'mlst_report.' + time_str + '.tab')
cstrittmatter@0: mlst_report_exist = os.path.isfile(mlst_report)
cstrittmatter@0: with open(mlst_report, 'at') as writer:
cstrittmatter@0: if not mlst_report_exist:
cstrittmatter@0: writer.write('\t'.join(['#sample', 'ReMatCh_run', 'consensus_type', 'ST'] + loci_order) + '\n')
cstrittmatter@0: writer.write('\t'.join([sample, run_times, consensus_type, str(st)] + alleles_profile.split(',')) + '\n')
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def run_get_st(sample, mlst_dicts, consensus_sequences, mlst_consensus, run_times, outdir, time_str):
cstrittmatter@0: if mlst_consensus == 'all':
cstrittmatter@0: for consensus_type in consensus_sequences:
cstrittmatter@0: print('Searching MLST for ' + consensus_type + ' consensus')
cstrittmatter@0: st, alleles_profile = check_mlst.get_st(mlst_dicts, consensus_sequences[consensus_type])
cstrittmatter@0: write_mlst_report(sample, run_times, consensus_type, st, alleles_profile, mlst_dicts[2], outdir, time_str)
cstrittmatter@0: print('ST found: ' + str(st) + ' (' + alleles_profile + ')')
cstrittmatter@0: else:
cstrittmatter@0: st, alleles_profile = check_mlst.get_st(mlst_dicts, consensus_sequences[mlst_consensus])
cstrittmatter@0: write_mlst_report(sample, run_times, mlst_consensus, st, alleles_profile, mlst_dicts[2], outdir, time_str)
cstrittmatter@0: print('ST found for ' + mlst_consensus + ' consensus: ' + str(st) + ' (' + alleles_profile + ')')
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def write_summary_report(outdir, reported_data_type, time_str, gene_list_reference, genes_present):
cstrittmatter@0: with open(os.path.join(outdir,
cstrittmatter@0: 'summary.{reported_data_type}.{time_str}.tab'.format(reported_data_type=reported_data_type,
cstrittmatter@0: time_str=time_str)), 'wt') as writer:
cstrittmatter@0: seq_list = []
cstrittmatter@0: for info in list(genes_present.values()):
cstrittmatter@0: seq_list.extend(list(info.keys()))
cstrittmatter@0: seq_list = list(set(seq_list))
cstrittmatter@0: writer.write('#sample' + '\t' + '\t'.join([gene_list_reference[seq] for seq in sorted(seq_list)]) + '\n')
cstrittmatter@0: for sample, info in list(genes_present.items()):
cstrittmatter@0: data = []
cstrittmatter@0: for seq in sorted(seq_list):
cstrittmatter@0: if seq in info:
cstrittmatter@0: data.append(info[seq])
cstrittmatter@0: else:
cstrittmatter@0: data.append('NF')
cstrittmatter@0: writer.write(sample + '\t' + '\t'.join(data) + '\n')
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def run_rematch(args):
cstrittmatter@0: workdir = os.path.abspath(args.workdir)
cstrittmatter@0: if not os.path.isdir(workdir):
cstrittmatter@0: os.makedirs(workdir)
cstrittmatter@0:
cstrittmatter@0: aspera_key = os.path.abspath(args.asperaKey.name) if args.asperaKey is not None else None
cstrittmatter@0:
cstrittmatter@0: # Start logger
cstrittmatter@0: logfile, time_str = utils.start_logger(workdir)
cstrittmatter@0:
cstrittmatter@0: # Get general information
cstrittmatter@0: script_path = utils.general_information(logfile, __version__, workdir, time_str, args.doNotUseProvidedSoftware,
cstrittmatter@0: aspera_key, args.downloadCramBam, args.SRA, args.SRAopt)
cstrittmatter@0:
cstrittmatter@0: # Set list_ids
cstrittmatter@0: list_ids, searched_fastq_files = get_list_ids(workdir, args.listIDs.name if args.listIDs is not None else None,
cstrittmatter@0: args.taxon)
cstrittmatter@0:
cstrittmatter@0: mlst_sequences = None
cstrittmatter@0: mlst_dicts = None
cstrittmatter@0: if args.mlst is not None:
cstrittmatter@0: time_taken_pub_mlst, mlst_dicts, mlst_sequences = check_mlst.download_pub_mlst_xml(args.mlst,
cstrittmatter@0: args.mlstSchemaNumber,
cstrittmatter@0: workdir)
cstrittmatter@0: args.softClip_recodeRun = 'first'
cstrittmatter@0:
cstrittmatter@0: if args.reference is None:
cstrittmatter@0: if args.mlst is not None:
cstrittmatter@0: reference_file = check_mlst.check_existing_schema(args.mlst, args.mlstSchemaNumber, script_path)
cstrittmatter@0: args.extraSeq = 200
cstrittmatter@0: if reference_file is None:
cstrittmatter@0: print('It was not found provided MLST scheme sequences for ' + args.mlst)
cstrittmatter@0: print('Trying to obtain reference MLST sequences from PubMLST')
cstrittmatter@0: if len(mlst_sequences) > 0:
cstrittmatter@0: reference_file = check_mlst.write_mlst_reference(args.mlst, mlst_sequences, workdir, time_str)
cstrittmatter@0: args.extraSeq = 0
cstrittmatter@0: else:
cstrittmatter@0: sys.exit('It was not possible to download MLST sequences from PubMLST!')
cstrittmatter@0: else:
cstrittmatter@0: print('Using provided scheme as referece: ' + reference_file)
cstrittmatter@0: else:
cstrittmatter@0: sys.exit('Need to provide at least one of the following options: "--reference" and "--mlst"')
cstrittmatter@0: else:
cstrittmatter@0: reference_file = os.path.abspath(args.reference.name)
cstrittmatter@0:
cstrittmatter@0: # Run ReMatCh for each sample
cstrittmatter@0: print('\n' + 'STARTING ReMatCh' + '\n')
cstrittmatter@0:
cstrittmatter@0: # Clean sequences headers
cstrittmatter@0: reference_file, gene_list_reference, reference_dict = clean_headers_reference_file(reference_file, workdir,
cstrittmatter@0: args.extraSeq)
cstrittmatter@0:
cstrittmatter@0: if args.mlst is not None:
cstrittmatter@0: problem_genes = False
cstrittmatter@0: for header in mlst_sequences:
cstrittmatter@0: if header not in gene_list_reference:
cstrittmatter@0: print('MLST gene {header} not found between reference sequences'.format(header=header))
cstrittmatter@0: problem_genes = True
cstrittmatter@0: if problem_genes:
cstrittmatter@0: sys.exit('Missing MLST genes from reference sequences (at least sequences names do not match)!')
cstrittmatter@0:
cstrittmatter@0: if len(gene_list_reference) == 0:
cstrittmatter@0: sys.exit('No sequences left')
cstrittmatter@0:
cstrittmatter@0: # To use in combined report
cstrittmatter@0:
cstrittmatter@0: number_samples_successfully = 0
cstrittmatter@0: genes_present_coverage_depth = {}
cstrittmatter@0: genes_present_sequence_coverage = {}
cstrittmatter@0: for sample in list_ids:
cstrittmatter@0: sample_start_time = time.time()
cstrittmatter@0: print('\n\n' + 'Sample ID: ' + sample)
cstrittmatter@0:
cstrittmatter@0: # Create sample outdir
cstrittmatter@0: sample_outdir = os.path.join(workdir, sample, '')
cstrittmatter@0: if not os.path.isdir(sample_outdir):
cstrittmatter@0: os.mkdir(sample_outdir)
cstrittmatter@0:
cstrittmatter@0: run_successfully_fastq = None
cstrittmatter@0: time_taken_fastq = 0
cstrittmatter@0: sequencing_information = {'run_accession': None, 'instrument_platform': None, 'instrument_model': None,
cstrittmatter@0: 'library_layout': None, 'library_source': None, 'extra_run_accession': None,
cstrittmatter@0: 'nominal_length': None, 'read_count': None, 'base_count': None, 'date_download': None}
cstrittmatter@0: if not searched_fastq_files:
cstrittmatter@0: # Download Files
cstrittmatter@0: time_taken_fastq, run_successfully_fastq, fastq_files, sequencing_information = \
cstrittmatter@0: download.run_download(sample, args.downloadLibrariesType, aspera_key, sample_outdir,
cstrittmatter@0: args.downloadCramBam, args.threads, args.downloadInstrumentPlatform, args.SRA,
cstrittmatter@0: args.SRAopt)
cstrittmatter@0: else:
cstrittmatter@0: fastq_files = list_ids[sample]
cstrittmatter@0:
cstrittmatter@0: file_size = None
cstrittmatter@0:
cstrittmatter@0: run_successfully_rematch_first = None
cstrittmatter@0: run_successfully_rematch_second = None
cstrittmatter@0: time_taken_rematch_first = 0
cstrittmatter@0: time_taken_rematch_second = 0
cstrittmatter@0: sample_data_general_first = None
cstrittmatter@0: sample_data_general_second = None
cstrittmatter@0: if run_successfully_fastq is not False:
cstrittmatter@0: file_size = sum(os.path.getsize(fastq) for fastq in fastq_files)
cstrittmatter@0: # Run ReMatCh
cstrittmatter@0: time_taken_rematch_first, run_successfully_rematch_first, data_by_gene, sample_data_general_first, \
cstrittmatter@0: consensus_files, consensus_sequences = \
cstrittmatter@0: rematch_module.run_rematch_module(sample, fastq_files, reference_file, args.threads, sample_outdir,
cstrittmatter@0: args.extraSeq, args.minCovPresence, args.minCovCall,
cstrittmatter@0: args.minFrequencyDominantAllele, args.minGeneCoverage,
cstrittmatter@0: args.debug, args.numMapLoc, args.minGeneIdentity,
cstrittmatter@0: 'first', args.softClip_baseQuality, args.softClip_recodeRun,
cstrittmatter@0: reference_dict, args.softClip_cigarFlagRecode,
cstrittmatter@0: args.bowtieAlgo, args.bowtieOPT,
cstrittmatter@0: gene_list_reference, args.notWriteConsensus, clean_run=True)
cstrittmatter@0: if run_successfully_rematch_first:
cstrittmatter@0: if args.mlst is not None and (args.mlstRun == 'first' or args.mlstRun == 'all'):
cstrittmatter@0: run_get_st(sample, mlst_dicts, consensus_sequences, args.mlstConsensus, 'first', workdir, time_str)
cstrittmatter@0: genes_present_coverage_depth = write_data_by_gene(gene_list_reference, args.minGeneCoverage, sample,
cstrittmatter@0: data_by_gene, workdir, time_str, 'first_run',
cstrittmatter@0: args.minGeneIdentity, 'coverage_depth', args.summary,
cstrittmatter@0: genes_present_coverage_depth)
cstrittmatter@0: if args.reportSequenceCoverage:
cstrittmatter@0: genes_present_sequence_coverage = write_data_by_gene(gene_list_reference, args.minGeneCoverage,
cstrittmatter@0: sample, data_by_gene, workdir, time_str,
cstrittmatter@0: 'first_run', args.minGeneIdentity,
cstrittmatter@0: 'sequence_coverage', args.summary,
cstrittmatter@0: genes_present_sequence_coverage)
cstrittmatter@0: if args.doubleRun:
cstrittmatter@0: rematch_second_outdir = os.path.join(sample_outdir, 'rematch_second_run', '')
cstrittmatter@0: if not os.path.isdir(rematch_second_outdir):
cstrittmatter@0: os.mkdir(rematch_second_outdir)
cstrittmatter@0: consensus_concatenated_fasta, consensus_concatenated_gene_list, consensus_concatenated_dict, \
cstrittmatter@0: number_consensus_with_sequences = \
cstrittmatter@0: concatenate_extra_seq_2_consensus(consensus_files['noMatter'], reference_file, args.extraSeq,
cstrittmatter@0: rematch_second_outdir)
cstrittmatter@0: if len(consensus_concatenated_gene_list) > 0:
cstrittmatter@0: if args.mlst is None or \
cstrittmatter@0: (args.mlst is not None and number_consensus_with_sequences == len(gene_list_reference)):
cstrittmatter@0: time_taken_rematch_second, run_successfully_rematch_second, data_by_gene, \
cstrittmatter@0: sample_data_general_second, consensus_files, consensus_sequences = \
cstrittmatter@0: rematch_module.run_rematch_module(sample, fastq_files, consensus_concatenated_fasta,
cstrittmatter@0: args.threads, rematch_second_outdir, args.extraSeq,
cstrittmatter@0: args.minCovPresence, args.minCovCall,
cstrittmatter@0: args.minFrequencyDominantAllele, args.minGeneCoverage,
cstrittmatter@0: args.debug, args.numMapLoc,
cstrittmatter@0: args.minGeneIdentity, 'second',
cstrittmatter@0: args.softClip_baseQuality, args.softClip_recodeRun,
cstrittmatter@0: consensus_concatenated_dict,
cstrittmatter@0: args.softClip_cigarFlagRecode,
cstrittmatter@0: args.bowtieAlgo, args.bowtieOPT,
cstrittmatter@0: gene_list_reference, args.notWriteConsensus,
cstrittmatter@0: clean_run=True)
cstrittmatter@0: if not args.debug:
cstrittmatter@0: os.remove(consensus_concatenated_fasta)
cstrittmatter@0: if run_successfully_rematch_second:
cstrittmatter@0: if args.mlst is not None and (args.mlstRun == 'second' or args.mlstRun == 'all'):
cstrittmatter@0: run_get_st(sample, mlst_dicts, consensus_sequences, args.mlstConsensus, 'second',
cstrittmatter@0: workdir, time_str)
cstrittmatter@0: _ = write_data_by_gene(gene_list_reference, args.minGeneCoverage, sample, data_by_gene,
cstrittmatter@0: workdir, time_str, 'second_run', args.minGeneIdentity,
cstrittmatter@0: 'coverage_depth', False, {})
cstrittmatter@0: if args.reportSequenceCoverage:
cstrittmatter@0: _ = write_data_by_gene(gene_list_reference, args.minGeneCoverage, sample,
cstrittmatter@0: data_by_gene, workdir, time_str, 'second_run',
cstrittmatter@0: args.minGeneIdentity, 'sequence_coverage', False, {})
cstrittmatter@0: else:
cstrittmatter@0: print('Some sequences missing after ReMatCh module first run. Second run will not be'
cstrittmatter@0: ' performed')
cstrittmatter@0: if os.path.isfile(consensus_concatenated_fasta):
cstrittmatter@0: os.remove(consensus_concatenated_fasta)
cstrittmatter@0: if os.path.isdir(rematch_second_outdir):
cstrittmatter@0: utils.remove_directory(rematch_second_outdir)
cstrittmatter@0: else:
cstrittmatter@0: print('No sequences left after ReMatCh module first run. Second run will not be performed')
cstrittmatter@0: if os.path.isfile(consensus_concatenated_fasta):
cstrittmatter@0: os.remove(consensus_concatenated_fasta)
cstrittmatter@0: if os.path.isdir(rematch_second_outdir):
cstrittmatter@0: utils.remove_directory(rematch_second_outdir)
cstrittmatter@0:
cstrittmatter@0: if not searched_fastq_files and not args.keepDownloadedFastq and fastq_files is not None:
cstrittmatter@0: for fastq in fastq_files:
cstrittmatter@0: if os.path.isfile(fastq):
cstrittmatter@0: os.remove(fastq)
cstrittmatter@0:
cstrittmatter@0: time_taken = utils.run_time(sample_start_time)
cstrittmatter@0:
cstrittmatter@0: write_sample_report(sample, workdir, time_str, file_size, run_successfully_fastq,
cstrittmatter@0: run_successfully_rematch_first, run_successfully_rematch_second, time_taken_fastq,
cstrittmatter@0: time_taken_rematch_first, time_taken_rematch_second, time_taken, sequencing_information,
cstrittmatter@0: sample_data_general_first if run_successfully_rematch_first else
cstrittmatter@0: {'number_absent_genes': None, 'number_genes_multiple_alleles': None,
cstrittmatter@0: 'mean_sample_coverage': None},
cstrittmatter@0: sample_data_general_second if run_successfully_rematch_second else
cstrittmatter@0: {'number_absent_genes': None, 'number_genes_multiple_alleles': None,
cstrittmatter@0: 'mean_sample_coverage': None},
cstrittmatter@0: fastq_files if fastq_files is not None else '')
cstrittmatter@0:
cstrittmatter@0: if all([run_successfully_fastq is not False,
cstrittmatter@0: run_successfully_rematch_first is not False,
cstrittmatter@0: run_successfully_rematch_second is not False]):
cstrittmatter@0: number_samples_successfully += 1
cstrittmatter@0:
cstrittmatter@0: if args.summary:
cstrittmatter@0: write_summary_report(workdir, 'coverage_depth', time_str, gene_list_reference, genes_present_coverage_depth)
cstrittmatter@0: if args.reportSequenceCoverage:
cstrittmatter@0: write_summary_report(workdir, 'sequence_coverage', time_str, gene_list_reference,
cstrittmatter@0: genes_present_sequence_coverage)
cstrittmatter@0:
cstrittmatter@0: return number_samples_successfully, len(list_ids)
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: def main():
cstrittmatter@0: if sys.version_info[0] < 3:
cstrittmatter@0: sys.exit('Must be using Python 3. Try calling "python3 rematch.py"')
cstrittmatter@0:
cstrittmatter@0: parser = argparse.ArgumentParser(prog='rematch.py',
cstrittmatter@0: description='Reads mapping against target sequences, checking mapping and'
cstrittmatter@0: ' consensus sequences production',
cstrittmatter@0: formatter_class=argparse.ArgumentDefaultsHelpFormatter)
cstrittmatter@0: parser.add_argument('--version', help='Version information', action='version',
cstrittmatter@0: version='{prog} v{version}'.format(prog=parser.prog, version=__version__))
cstrittmatter@0:
cstrittmatter@0: parser_optional_general = parser.add_argument_group('General facultative options')
cstrittmatter@0: parser_optional_general.add_argument('-r', '--reference', type=argparse.FileType('r'),
cstrittmatter@0: metavar='/path/to/reference_sequence.fasta',
cstrittmatter@0: help='Fasta file containing reference sequences', required=False)
cstrittmatter@0: parser_optional_general.add_argument('-w', '--workdir', type=str, metavar='/path/to/workdir/directory/',
cstrittmatter@0: help='Path to the directory where ReMatCh will run and produce the outputs'
cstrittmatter@0: ' with reads (ended with fastq.gz/fq.gz and, in case of PE data, pair-end'
cstrittmatter@0: ' direction coded as _R1_001 / _R2_001 or _1 / _2) already'
cstrittmatter@0: ' present (organized in sample folders) or to be downloaded',
cstrittmatter@0: required=False, default='.')
cstrittmatter@0: parser_optional_general.add_argument('-j', '--threads', type=int, metavar='N', help='Number of threads to use',
cstrittmatter@0: required=False, default=1)
cstrittmatter@0: parser_optional_general.add_argument('--mlst', type=str, metavar='"Streptococcus agalactiae"',
cstrittmatter@0: help='Species name (same as in PubMLST) to be used in MLST'
cstrittmatter@0: ' determination. ReMatCh will use Bowtie2 very-sensitive-local mapping'
cstrittmatter@0: ' parameters and will recode the soft clip CIGAR flags of the first run',
cstrittmatter@0: required=False)
cstrittmatter@0: parser_optional_general.add_argument('--doNotUseProvidedSoftware', action='store_true',
cstrittmatter@0: help='Tells ReMatCh to not use Bowtie2, Samtools and Bcftools that are'
cstrittmatter@0: ' provided with it')
cstrittmatter@0:
cstrittmatter@0: parser_optional_download_exclusive = parser.add_mutually_exclusive_group()
cstrittmatter@0: parser_optional_download_exclusive.add_argument('-l', '--listIDs', type=argparse.FileType('r'),
cstrittmatter@0: metavar='/path/to/list_IDs.txt',
cstrittmatter@0: help='Path to list containing the IDs to be'
cstrittmatter@0: ' downloaded (one per line)', required=False)
cstrittmatter@0: parser_optional_download_exclusive.add_argument('-t', '--taxon', type=str, metavar='"Streptococcus agalactiae"',
cstrittmatter@0: help='Taxon name for which ReMatCh will download fastq files',
cstrittmatter@0: required=False)
cstrittmatter@0:
cstrittmatter@0: parser_optional_rematch = parser.add_argument_group('ReMatCh module facultative options')
cstrittmatter@0: parser_optional_rematch.add_argument('--extraSeq', type=int, metavar='N',
cstrittmatter@0: help='Sequence length added to both ends of target sequences (usefull to'
cstrittmatter@0: ' improve reads mapping to the target one) that will be trimmed in'
cstrittmatter@0: ' ReMatCh outputs', required=False, default=0)
cstrittmatter@0: parser_optional_rematch.add_argument('--minCovPresence', type=int, metavar='N',
cstrittmatter@0: help='Reference position minimum coverage depth to consider the position to be'
cstrittmatter@0: ' present in the sample', required=False, default=5)
cstrittmatter@0: parser_optional_rematch.add_argument('--minCovCall', type=int, metavar='N',
cstrittmatter@0: help='Reference position minimum coverage depth to perform a base call. Lower'
cstrittmatter@0: ' coverage will be coded as N', required=False, default=10)
cstrittmatter@0: parser_optional_rematch.add_argument('--minFrequencyDominantAllele', type=float, metavar='0.6',
cstrittmatter@0: help='Minimum relative frequency of the dominant allele coverage depth (value'
cstrittmatter@0: ' between [0, 1]). Positions with lower values will be considered as'
cstrittmatter@0: ' having multiple alleles (and will be coded as N)', required=False,
cstrittmatter@0: default=0.6)
cstrittmatter@0: parser_optional_rematch.add_argument('--minGeneCoverage', type=int, metavar='N',
cstrittmatter@0: help='Minimum percentage of target reference gene sequence covered'
cstrittmatter@0: ' by --minCovPresence to consider a gene to be present (value'
cstrittmatter@0: ' between [0, 100])', required=False, default=70)
cstrittmatter@0: parser_optional_rematch.add_argument('--minGeneIdentity', type=int, metavar='N',
cstrittmatter@0: help='Minimum percentage of identity of reference gene sequence covered'
cstrittmatter@0: ' by --minCovCall to consider a gene to be present (value'
cstrittmatter@0: ' between [0, 100]). One INDEL will be considered as one difference',
cstrittmatter@0: required=False, default=80)
cstrittmatter@0: parser_optional_rematch.add_argument('--numMapLoc', type=int, metavar='N', help=argparse.SUPPRESS, required=False,
cstrittmatter@0: default=1)
cstrittmatter@0: # parser_optional_rematch.add_argument('--numMapLoc', type=int, metavar='N', help='Maximum number of locations to which a read can map (sometimes useful when mapping against similar sequences)', required=False, default=1)
cstrittmatter@0: parser_optional_rematch.add_argument('--doubleRun', action='store_true',
cstrittmatter@0: help='Tells ReMatCh to run a second time using as reference the noMatter'
cstrittmatter@0: ' consensus sequence produced in the first run. This will improve'
cstrittmatter@0: ' consensus sequence determination for sequences with high percentage of'
cstrittmatter@0: ' target reference gene sequence covered')
cstrittmatter@0: parser_optional_rematch.add_argument('--reportSequenceCoverage', action='store_true',
cstrittmatter@0: help='Produce an extra combined_report.data_by_gene with the sequence coverage'
cstrittmatter@0: ' instead of coverage depth')
cstrittmatter@0: parser_optional_rematch.add_argument('--summary', action='store_true',
cstrittmatter@0: help='Produce extra report files containing only sequences present in at least'
cstrittmatter@0: ' one sample (usefull when using a large number of reference'
cstrittmatter@0: ' sequences, and only for first run)')
cstrittmatter@0: parser_optional_rematch.add_argument('--notWriteConsensus', action='store_true',
cstrittmatter@0: help='Do not write consensus sequences')
cstrittmatter@0: parser_optional_rematch.add_argument('--bowtieAlgo', type=str, metavar='"--very-sensitive-local"',
cstrittmatter@0: help='Bowtie2 alignment mode. It can be an end-to-end alignment (unclipped'
cstrittmatter@0: ' alignment) or local alignment (soft clipped alignment). Also, can'
cstrittmatter@0: ' choose between fast or sensitive alignments. Please check Bowtie2'
cstrittmatter@0: ' manual for extra'
cstrittmatter@0: ' information: http://bowtie-bio.sourceforge.net/bowtie2/index.shtml .'
cstrittmatter@0: ' This option should be provided between quotes and starting with'
cstrittmatter@0: ' an empty space (like --bowtieAlgo " --very-fast") or using equal'
cstrittmatter@0: ' sign (like --bowtieAlgo="--very-fast")',
cstrittmatter@0: required=False, default='--very-sensitive-local')
cstrittmatter@0: parser_optional_rematch.add_argument('--bowtieOPT', type=str, metavar='"--no-mixed"',
cstrittmatter@0: help='Extra Bowtie2 options. This option should be provided between quotes and'
cstrittmatter@0: ' starting with an empty space (like --bowtieOPT " --no-mixed") or using'
cstrittmatter@0: ' equal sign (like --bowtieOPT="--no-mixed")',
cstrittmatter@0: required=False)
cstrittmatter@0: parser_optional_rematch.add_argument('--debug', action='store_true',
cstrittmatter@0: help='DeBug Mode: do not remove temporary files')
cstrittmatter@0:
cstrittmatter@0: parser_optional_mlst = parser.add_argument_group('MLST facultative options')
cstrittmatter@0: parser_optional_rematch.add_argument('--mlstReference', action='store_true',
cstrittmatter@0: help='If the curated scheme for MLST alleles is available, tells ReMatCh to'
cstrittmatter@0: ' use these as reference (force Bowtie2 to run with very-sensitive-local'
cstrittmatter@0: ' parameters, and sets --extraSeq to 200), otherwise ReMatCh uses the'
cstrittmatter@0: ' first alleles of each MLST gene fragment in PubMLST as reference'
cstrittmatter@0: ' sequences (force Bowtie2 to run with very-sensitive-local parameters,'
cstrittmatter@0: ' and sets --extraSeq to 0)')
cstrittmatter@0: parser_optional_mlst.add_argument('--mlstSchemaNumber', type=int, metavar='N',
cstrittmatter@0: help='Number of the species PubMLST schema to be used in case of multiple schemes'
cstrittmatter@0: ' available (by default will use the first schema)', required=False)
cstrittmatter@0: parser_optional_mlst.add_argument('--mlstConsensus', choices=['noMatter', 'correct', 'alignment', 'all'], type=str,
cstrittmatter@0: metavar='noMatter',
cstrittmatter@0: help='Consensus sequence to be used in MLST'
cstrittmatter@0: ' determination (available options: %(choices)s)', required=False,
cstrittmatter@0: default='noMatter')
cstrittmatter@0: parser_optional_mlst.add_argument('--mlstRun', choices=['first', 'second', 'all'], type=str, metavar='first',
cstrittmatter@0: help='ReMatCh run outputs to be used in MLST determination (available'
cstrittmatter@0: ' options: %(choices)s)', required=False, default='all')
cstrittmatter@0:
cstrittmatter@0: parser_optional_download = parser.add_argument_group('Download facultative options')
cstrittmatter@0: parser_optional_download.add_argument('-a', '--asperaKey', type=argparse.FileType('r'),
cstrittmatter@0: metavar='/path/to/asperaweb_id_dsa.openssh',
cstrittmatter@0: help='Tells ReMatCh to download fastq files from ENA using Aspera'
cstrittmatter@0: ' Connect. With this option, the path to Private-key file'
cstrittmatter@0: ' asperaweb_id_dsa.openssh must be provided (normaly found in'
cstrittmatter@0: ' ~/.aspera/connect/etc/asperaweb_id_dsa.openssh).', required=False)
cstrittmatter@0: parser_optional_download.add_argument('-k', '--keepDownloadedFastq', action='store_true',
cstrittmatter@0: help='Tells ReMatCh to keep the fastq files downloaded')
cstrittmatter@0: parser_optional_download.add_argument('--downloadLibrariesType', type=str, metavar='PAIRED',
cstrittmatter@0: help='Tells ReMatCh to download files with specific library'
cstrittmatter@0: ' layout (available options: %(choices)s)',
cstrittmatter@0: choices=['PAIRED', 'SINGLE', 'BOTH'], required=False, default='BOTH')
cstrittmatter@0: parser_optional_download.add_argument('--downloadInstrumentPlatform', type=str, metavar='ILLUMINA',
cstrittmatter@0: help='Tells ReMatCh to download files with specific library layout (available'
cstrittmatter@0: ' options: %(choices)s)', choices=['ILLUMINA', 'ALL'], required=False,
cstrittmatter@0: default='ILLUMINA')
cstrittmatter@0: parser_optional_download.add_argument('--downloadCramBam', action='store_true',
cstrittmatter@0: help='Tells ReMatCh to also download cram/bam files and convert them to fastq'
cstrittmatter@0: ' files')
cstrittmatter@0:
cstrittmatter@0: parser_optional_sra = parser.add_mutually_exclusive_group()
cstrittmatter@0: parser_optional_sra.add_argument('--SRA', action='store_true',
cstrittmatter@0: help='Tells getSeqENA.py to download reads in fastq format only from NCBI SRA'
cstrittmatter@0: ' database (not recommended)')
cstrittmatter@0: parser_optional_sra.add_argument('--SRAopt', action='store_true',
cstrittmatter@0: help='Tells getSeqENA.py to download reads from NCBI SRA if the download from ENA'
cstrittmatter@0: ' fails')
cstrittmatter@0:
cstrittmatter@0: parser_optional_soft_clip = parser.add_argument_group('Soft clip facultative options')
cstrittmatter@0: parser_optional_soft_clip.add_argument('--softClip_baseQuality', type=int, metavar='N',
cstrittmatter@0: help='Base quality phred score in reads soft clipped regions',
cstrittmatter@0: required=False,
cstrittmatter@0: default=7)
cstrittmatter@0: parser_optional_soft_clip.add_argument('--softClip_recodeRun', type=str, metavar='first',
cstrittmatter@0: help='ReMatCh run to recode soft clipped regions (available'
cstrittmatter@0: ' options: %(choices)s)', choices=['first', 'second', 'both', 'none'],
cstrittmatter@0: required=False, default='none')
cstrittmatter@0: parser_optional_soft_clip.add_argument('--softClip_cigarFlagRecode', type=str, metavar='M',
cstrittmatter@0: help='CIGAR flag to recode CIGAR soft clip (available options: %(choices)s)',
cstrittmatter@0: choices=['M', 'I', 'X'], required=False, default='X')
cstrittmatter@0:
cstrittmatter@0: args = parser.parse_args()
cstrittmatter@0:
cstrittmatter@0: msg = []
cstrittmatter@0: if args.reference is None and not args.mlstReference:
cstrittmatter@0: msg.append('At least --reference or --mlstReference should be provided')
cstrittmatter@0: elif args.reference is not None and args.mlstReference:
cstrittmatter@0: msg.append('Only --reference or --mlstReference should be provided')
cstrittmatter@0: else:
cstrittmatter@0: if args.mlstReference:
cstrittmatter@0: if args.mlst is None:
cstrittmatter@0: msg.append('Please provide species name using --mlst')
cstrittmatter@0: if args.minFrequencyDominantAllele < 0 or args.minFrequencyDominantAllele > 1:
cstrittmatter@0: msg.append('--minFrequencyDominantAllele should be a value between [0, 1]')
cstrittmatter@0: if args.minGeneCoverage < 0 or args.minGeneCoverage > 100:
cstrittmatter@0: msg.append('--minGeneCoverage should be a value between [0, 100]')
cstrittmatter@0: if args.minGeneIdentity < 0 or args.minGeneIdentity > 100:
cstrittmatter@0: msg.append('--minGeneIdentity should be a value between [0, 100]')
cstrittmatter@0: if args.notWriteConsensus and args.doubleRun:
cstrittmatter@0: msg.append('--notWriteConsensus and --doubleRun cannot be used together.'
cstrittmatter@0: ' Maybe you only want to use --doubleRun')
cstrittmatter@0:
cstrittmatter@0: if len(msg) > 0:
cstrittmatter@0: argparse.ArgumentParser.error('\n'.join(msg))
cstrittmatter@0:
cstrittmatter@0: start_time = time.time()
cstrittmatter@0:
cstrittmatter@0: number_samples_successfully, samples_total_number = run_rematch(args)
cstrittmatter@0:
cstrittmatter@0: print('\n' + 'END ReMatCh')
cstrittmatter@0: print('\n' +
cstrittmatter@0: str(number_samples_successfully) + ' samples out of ' + str(samples_total_number) + ' run successfully')
cstrittmatter@0: time_taken = utils.run_time(start_time)
cstrittmatter@0: del time_taken
cstrittmatter@0:
cstrittmatter@0: if number_samples_successfully == 0:
cstrittmatter@0: sys.exit('No samples run successfully!')
cstrittmatter@0:
cstrittmatter@0:
cstrittmatter@0: if __name__ == "__main__":
cstrittmatter@0: main()