INNUca.py

#!/usr/bin/env python2

# -*- coding: utf-8 -*-


"""
INNUca - Reads Control and Assembly
INNUca.py - INNUENDO quality control of reads, de novo assembly and contigs
quality assessment, and possible contamination detection
<https://github.com/B-UMMI/INNUca>

Copyright (C) 2017 Miguel Machado <mpmachado@medicina.ulisboa.pt>

Last modified: September 24, 2018

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

import modules.utils as utils
from modules.kraken import run_for_innuca as kraken
from modules.kraken import get_kraken_version as kraken_version
import modules.fastQintegrity as fastQintegrity
import modules.estimated_coverage as coverage
import modules.fastqc as fastqc
import modules.trimmomatic as trimmomatic
import modules.pear as pear
import modules.spades as spades
import modules.assembly_mapping as assembly_mapping
import modules.pilon as pilon
import modules.mlst as mlst
import modules.trueCoverage_rematch as trueCoverage
import modules.combine_reports as combine_reports
import time
import os
import sys


def get_trueCoverage_config(skipTrueCoverage, trueConfigFile, speciesExpected, script_path):
    trueCoverage_config = None
    if not skipTrueCoverage:
        trueCoverage_reference = None
        trueCoverage_config = None

        if trueConfigFile is None:
            print 'No trueCoverage_ReMatCh config file was provided. Search for default files'
            trueCoverage_config_file, trueCoverage_reference = trueCoverage.check_existing_default_config(speciesExpected, script_path)
        else:
            trueCoverage_config_file = trueConfigFile

        if trueCoverage_config_file is not None:
            trueCoverage_config = trueCoverage.parse_config(trueCoverage_config_file)
        if trueConfigFile is None and trueCoverage_config is not None:
            trueCoverage_config['reference_file'] = trueCoverage_reference

        if trueCoverage_config is not None:
            print 'The following trueCoverage_ReMatCh config file will be used: ' + trueCoverage_config_file
            print 'The following trueCoverage_ReMatCh reference file will be used: ' + trueCoverage_config['reference_file'] + '\n'
        else:
            print 'No trueCoverage_ReMatCh config file was found'
    return trueCoverage_config


def include_rematch_dependencies_path(doNotUseProvidedSoftware):
    command = ['which', 'rematch.py']
    run_successfully, stdout, stderr = utils.runCommandPopenCommunicate(command, False, None, False)
    if run_successfully:
        rematch = stdout.splitlines()[0]
        path_variable = os.environ['PATH']
        script_folder = os.path.dirname(rematch)
        if not doNotUseProvidedSoftware:
            bcftools = os.path.join(script_folder, 'src', 'bcftools-1.3.1', 'bin')
            os.environ['PATH'] = str(':'.join([bcftools, path_variable]))


version_kraken_global = None


def main():
    version = '4.0'
    args = utils.parseArguments(version)

    general_start_time = time.time()
    time_str = time.strftime("%Y%m%d-%H%M%S")

    # Check if output directory exists
    outdir = os.path.abspath(os.path.join(args.outdir, ''))
    if not os.path.isdir(outdir):
        os.makedirs(outdir)

    # Start logger
    if not args.noLog:
        sys.stdout = utils.Logger(outdir, time_str)

    print '\n' + '==========> INNUca.py <=========='
    print '\n' + 'Program start: ' + time.ctime()

    # Tells where the logfile will be stored
    if not args.noLog:
        print '\n' + 'LOGFILE:'
        print sys.stdout.getLogFile()

    # Print command
    print '\n' + 'COMMAND:'
    script_path = os.path.abspath(sys.argv[0])
    print sys.executable + ' ' + script_path + ' ' + ' '.join(sys.argv[1:])

    # Print directory where programme was lunch
    print '\n' + 'PRESENT DIRECTORY:'
    print os.getcwd()

    # Print program version
    print '\n' + 'VERSION INNUca.py:'
    utils.script_version_git(version=version, current_directory=os.getcwd(), script_path=script_path,
                             no_git_info=args.noGitInfo)

    # Get CPU information
    utils.get_cpu_information(outdir, time_str)

    # Get trueCoverage_ReMatCh settings
    trueCoverage_config = get_trueCoverage_config(args.skipTrueCoverage, args.trueConfigFile.name if args.trueConfigFile is not None else None, args.speciesExpected, script_path)

    # Check programms
    programs_version_dictionary = {}
    programs_version_dictionary['gunzip'] = ['--version', '>=', '1.6']

    # Java check first for java dependents check next
    if not (args.skipFastQC and args.skipTrimmomatic and (args.skipPilon or args.skipSPAdes)):
        # programs_version_dictionary['java'] = ['-version', '>=', '1.8']
        programs_version_dictionary['java'] = [None, '>=', '1.8']  # For OpenJDK compatibility
    missingPrograms, programs_version_dictionary = utils.checkPrograms(programs_version_dictionary)
    if len(missingPrograms) > 0:
        sys.exit('\n' + 'Errors:' + '\n' + '\n'.join(missingPrograms))

    if args.runKraken:
        global version_kraken_global
        version_kraken_global = kraken_version()
        if version_kraken_global == 2:
            programs_version_dictionary['kraken2'] = ['--version', '>=', '2.0.6']
        else:
            programs_version_dictionary['kraken'] = ['--version', '>=', '0.10.6']
            programs_version_dictionary['kraken-repor'] = ['--version', '>=', '0.10.6']
    if not args.skipTrueCoverage or trueCoverage_config is not None:
        include_rematch_dependencies_path(args.doNotUseProvidedSoftware)
        programs_version_dictionary['rematch.py'] = ['--version', '>=', '3.2']
        programs_version_dictionary['bcftools'] = ['--version', '==', '1.3.1']
    if not (args.skipTrueCoverage and ((args.skipAssemblyMapping and args.skipPilon) or args.skipSPAdes)):
        programs_version_dictionary['bowtie2'] = ['--version', '>=', '2.2.9']
        programs_version_dictionary['samtools'] = ['--version', '==', '1.3.1']
    if not args.skipFastQC:
        programs_version_dictionary['fastqc'] = ['--version', '==', '0.11.5']
    if not args.skipTrimmomatic:
        programs_version_dictionary['trimmomatic-0.36.jar'] = ['-version', '==', '0.36']
    if args.runPear:
        programs_version_dictionary['pear'] = ['--version', '>=', '0.9.10']
    if not args.skipSPAdes:
        programs_version_dictionary['spades.py'] = ['--version', '>=', '3.9.0']
    if not (args.skipPilon or args.skipSPAdes):
        programs_version_dictionary['pilon-1.18.jar'] = ['--version', '==', '1.18']
    if not (args.skipMLST or args.skipSPAdes):
        programs_version_dictionary['mlst'] = ['--version', '>=', '2.4']

    # Set and print PATH variable
    utils.setPATHvariable(args, script_path)

    missingPrograms, programs_version_dictionary = utils.checkPrograms(programs_version_dictionary)
    if len(missingPrograms) > 0:
        sys.exit('\n' + 'Errors:' + '\n' + '\n'.join(missingPrograms))

    # .jar paths
    jar_path_trimmomatic = None
    if not args.skipTrimmomatic:
        jar_path_trimmomatic = programs_version_dictionary['trimmomatic-0.36.jar'][3]

    jar_path_pilon = None
    if not args.skipPilon and not args.skipSPAdes:
        jar_path_pilon = programs_version_dictionary['pilon-1.18.jar'][3]

    rematch_script = None
    # ReMatCh path
    if not args.skipTrueCoverage:
        rematch_script = programs_version_dictionary['rematch.py'][3]

    # pairEnd_filesSeparation_list = args.pairEnd_filesSeparation
    pairEnd_filesSeparation_list = None
    samples, inputDirectory, removeCreatedSamplesDirectories, indir_same_outdir = get_samples(args.inputDirectory, args.fastq, outdir, pairEnd_filesSeparation_list)

    # Start running the analysis
    print '\n' + 'RUNNING INNUca.py'

    # Prepare run report file
    samples_report_path = os.path.join(outdir, 'samples_report.' + time_str + '.tab')
    utils.start_sample_report_file(samples_report_path)

    number_samples_successfully = 0
    number_samples_pass = 0
    number_samples_warning = 0

    # Get MLST scheme to use
    scheme = 'unknown'
    species_genus, mlst_scheme_genus = None, None
    if not args.skipMLST and not args.skipSPAdes:
        scheme, species_genus, mlst_scheme_genus = mlst.getScheme(args.speciesExpected)
        # Print path to blastn
        mlst.getBlastPath()

    # Memory
    available_memory_GB = utils.get_free_memory() / (1024.0 ** 2)
    # Determine SPAdes maximum memory
    spadesMaxMemory = None
    if not args.skipSPAdes:
        print ''
        spadesMaxMemory = spades.define_memory(args.spadesMaxMemory, args.threads, available_memory_GB)
    # Determine .jar maximum memory
    jarMaxMemory = 'off'
    if not (args.skipTrimmomatic and (args.skipSPAdes or args.skipPilon)):
        print ''
        jarMaxMemory = utils.define_jar_max_memory(args.jarMaxMemory, args.threads, available_memory_GB)

    # Run INNUca for each sample
    sample_report_json = {}
    for sample in samples:
        sample_start_time = time.time()

        print '\n' + 'Sample: ' + sample + '\n'

        # Create sample outdir
        sample_outdir = os.path.abspath(os.path.join(outdir, sample, ''))
        if not os.path.isdir(sample_outdir):
            os.makedirs(sample_outdir)

        # Get fastq files
        fastq_files = utils.searchFastqFiles(os.path.join(inputDirectory, sample, ''), pairEnd_filesSeparation_list, False)
        if len(fastq_files) == 1:
            print 'Only one fastq file was found: ' + str(fastq_files)
            print 'Pair-End sequencing is required. Moving to the next sample'
            continue
        elif len(fastq_files) == 0:
            print 'No compressed fastq files were found. Continue to the next sample'
            continue

        print 'The following files will be used:'
        print str(fastq_files) + '\n'

        # Run INNUca.py analysis
        run_successfully, pass_qc, run_report = run_innuca(sample, sample_outdir, fastq_files, args, script_path, scheme, spadesMaxMemory, jar_path_trimmomatic, jar_path_pilon, jarMaxMemory, trueCoverage_config, rematch_script, species_genus, mlst_scheme_genus)

        # Save sample fail report
        utils.write_fail_report(os.path.join(sample_outdir, 'fail_report.txt'), run_report)

        # Save warning report
        write_warning_report(os.path.join(sample_outdir, 'warning_report.txt'), run_report)

        # Get raw reads files size
        fileSize = sum(os.path.getsize(fastq) for fastq in fastq_files)

        # Remove sample directory if it was created during the process
        if removeCreatedSamplesDirectories and not indir_same_outdir:
            utils.removeDirectory(os.path.join(inputDirectory, sample, ''))

        # Remove reads folder if --fastq was provided
        if args.fastq is not None:
            utils.removeDirectory(os.path.join(outdir, 'reads', ''))

        print 'END ' + sample + ' analysis'
        time_taken = utils.runTime(sample_start_time)

        # Save run report
        warning, json_pass_qc = utils.write_sample_report(samples_report_path, sample, run_successfully, pass_qc, time_taken, fileSize, run_report)

        # Save runs statistics
        if run_successfully:
            number_samples_successfully += 1
        if pass_qc:
            if warning:
                number_samples_warning += 1
            else:
                number_samples_pass += 1

        sample_report_json[sample] = {'run_successfully': run_successfully, 'pass_qc': json_pass_qc, 'modules_run_report': run_report}

    # Save combine_samples_reports
    combine_reports.combine_reports(outdir, outdir, args.json, time_str, len(samples))

    # Save sample_report in json
    if args.json:
        import json
        with open(os.path.join(outdir, 'samples_report.' + time_str + '.json'), 'wt') as writer:
            json.dump(sample_report_json, writer)

    # Remove temporary folder with symlink to fastq files in case of --fastq use
    if args.inputDirectory is None and args.fastq is not None:
        utils.removeDirectory(os.path.join(inputDirectory, ''))

    # Run report
    print '\n' + 'END INNUca.py'
    print '\n' + 'Pipeline problems: {not_run_successfully} samples'.format(not_run_successfully=(len(samples) - number_samples_successfully))
    print '\n' + 'FAIL: {number_samples_fail} samples'.format(number_samples_fail=(len(samples) - number_samples_pass - number_samples_warning))
    print '\n' + 'WARNING: {number_samples_warning} samples'.format(number_samples_warning=number_samples_warning)
    print '\n' + 'PASS: {number_samples_pass} samples'.format(number_samples_pass=number_samples_pass)
    time_taken = utils.runTime(general_start_time)
    del time_taken

    # Check whether INNUca.py run at least one sample successfully
    if number_samples_successfully == 0:
        sys.exit('No samples run successfully!')


def write_warning_report(warning_report_path, run_report):
    with open(warning_report_path, 'wt') as writer_warningReport:
        warnings = []
        for step in ('reads_Kraken', 'first_FastQC', 'Trimmomatic', 'second_FastQC', 'Pear', 'SPAdes',
                     'Assembly_Mapping', 'MLST', 'reads_Kraken'):
            if len(run_report[step][4]) > 0:
                if step == 'first_FastQC' and \
                        run_report['second_FastQC'][1] is not False and \
                        len(run_report['second_FastQC'][4]) == 0:
                    continue
                else:
                    if len(run_report[step][4]) == 0:
                        continue
                    else:
                        warnings.append('#' + step)
                        for key, warning_reasons in run_report[step][4].items():
                            warnings.append('>' + str(key))
                            if isinstance(warning_reasons, (list, tuple)):
                                for reasons in warning_reasons:
                                    warnings.append(str(reasons))
                            else:
                                warnings.append(str(warning_reasons))
        writer_warningReport.write('\n'.join(warnings))


def get_sample_args_fastq(fastq_files_list, outdir, pairEnd_filesSeparation_list):
    new_indir = os.path.join(outdir, 'reads', '')
    utils.removeDirectory(new_indir)
    os.mkdir(new_indir)
    samples = []
    for fastq in fastq_files_list:
        fastq_link = os.path.join(new_indir, os.path.basename(fastq))
        os.symlink(fastq, fastq_link)
    samples, removeCreatedSamplesDirectories, indir_same_outdir = utils.checkSetInputDirectory(new_indir, outdir, pairEnd_filesSeparation_list)
    return new_indir, samples, removeCreatedSamplesDirectories, indir_same_outdir


def get_samples(args_inputDirectory, args_fastq, outdir, pairEnd_filesSeparation_list):
    if args_fastq is None:
        # Check if input directory exists with fastq files and store samples name that have fastq files
        inputDirectory = os.path.abspath(os.path.join(args_inputDirectory, ''))
        print ''
        samples, removeCreatedSamplesDirectories, indir_same_outdir = utils.checkSetInputDirectory(inputDirectory, outdir, pairEnd_filesSeparation_list)
    elif args_inputDirectory is None:
        fastq_files = [os.path.abspath(fastq.name) for fastq in args_fastq]
        if fastq_files[0] == fastq_files[1]:
            sys.exit('Same fastq file provided twice')
        inputDirectory, samples, removeCreatedSamplesDirectories, indir_same_outdir = get_sample_args_fastq(fastq_files, outdir, pairEnd_filesSeparation_list)

    return samples, inputDirectory, removeCreatedSamplesDirectories, indir_same_outdir


def run_innuca(sample_name, outdir, fastq_files, args, script_path, scheme, spades_max_memory, jar_path_trimmomatic,
               jar_path_pilon, jar_max_memory, true_coverage_config, rematch_script, species_genus, mlst_scheme_genus):
    threads = args.threads
    adapters_fasta = args.adapters
    if adapters_fasta is not None:
        adapters_fasta = os.path.abspath(adapters_fasta.name)
    genome_size = args.genomeSizeExpectedMb
    # run_successfully, pass_qc, time_taken, failing, warning, file_size
    skipped = [None, None, 0, {'sample': 'Skipped'}, {}, 'NA']
    not_run = [None, None, 0, {'sample': 'Not run'}, {}, 'NA']

    runs = {}

    # Run FastQ integrity check
    not_corruption_found, pass_qc, time_taken, failing, fastq_encoding, min_reads_length, max_reads_length = \
        fastQintegrity.runFastQintegrity(fastq_files, threads, outdir)
    runs['FastQ_Integrity'] = [not_corruption_found, pass_qc, time_taken, failing, {}, 'NA']

    pear_folder = None
    trimmomatic_folder = None
    if not_corruption_found:
        # Run Kraken
        # most_abundant_taxon_percent = None
        if args.runKraken:
            print('\n'
                  '--runKraken set. Running Kraken for reads')
            run_successfully, pass_qc, time_taken, failing, warning, _ = \
                kraken(species=args.speciesExpected, files_to_classify=fastq_files, kraken_db=args.krakenDB,
                       files_type='fastq', outdir=outdir, version_kraken=version_kraken_global,
                       db_mem=args.krakenMemory, quick=args.krakenQuick, min_percent_covered=args.krakenMinCov,
                       max_unclassified_frag=args.krakenMaxUnclass, min_base_quality=args.krakenMinQual,
                       threads=threads)
            runs['reads_Kraken'] = [run_successfully, pass_qc, time_taken, failing, warning, 'NA']
        else:
            runs['reads_Kraken'] = skipped

        if (run_successfully and not pass_qc) and not args.krakenProceed and not args.krakenIgnoreQC:
            print('\n'
                  'This sample does not pass Kraken module QA/QC. It will not proceed with INNUca pipeline')
        else:
            # Run first Estimated Coverage
            run_successfully_estimated_coverage = False
            estimated_coverage = None
            run_successfully_true_coverage = False
            pass_qc_true_coverage = False
            if not args.skipEstimatedCoverage:
                # Check whether the Estimated Coverage output is already present
                report_file = os.path.join(outdir, 'coverage_report.txt')
                if os.path.isfile(report_file):
                    os.remove(report_file)
                # Run getEstimatedCoverage
                run_successfully_estimated_coverage, pass_qc, time_taken, failing, estimated_coverage = \
                    coverage.getEstimatedCoverage(fastq_files, genome_size, outdir, threads,
                                                  args.estimatedMinimumCoverage)
                runs['first_Coverage'] = [run_successfully_estimated_coverage, pass_qc, time_taken, failing, {}, 'NA']
            else:
                print '--skipEstimatedCoverage set. Skipping First Estimated Coverage analysis'
                runs['first_Coverage'] = skipped

            trimmomatic_run_successfully = False

            # # Correct first estimation coverage with Kraken percentage
            # # Does not seem to be a good idea (at least for Streptococcus agalactiae)
            # if args.runKraken and \
            #         (runs['Kraken'][0] and runs['Kraken'][1]) and \
            #         most_abundant_taxon_percent is not None and \
            #         estimated_coverage is not None:
            #     new_estimation = estimated_coverage * (most_abundant_taxon_percent / 100)
            #     print('\n'
            #           'Correct estimated coverage ({estimated}x) with Kraken taxon percentage'
            #           ' coverage ({percent}%): {new_estimation}x'.format(estimated=estimated_coverage,
            #                                                              percent=most_abundant_taxon_percent,
            #                                                              new_estimation=new_estimation))
            #     estimated_coverage = new_estimation

            if args.skipEstimatedCoverage or (run_successfully_estimated_coverage and
                                              not estimated_coverage < args.estimatedMinimumCoverage):
                if not args.skipTrueCoverage and true_coverage_config is not None:
                    # Run True Coverage
                    run_successfully_true_coverage, pass_qc_true_coverage, time_taken, failing = \
                        trueCoverage.runTrueCoverage(sample_name, fastq_files, true_coverage_config['reference_file'],
                                                     threads, outdir,
                                                     true_coverage_config['length_extra_seq'],
                                                     true_coverage_config['minimum_depth_presence'],
                                                     true_coverage_config['minimum_depth_call'],
                                                     true_coverage_config['minimum_depth_frequency_dominant_allele'],
                                                     true_coverage_config['minimum_gene_coverage'], False,
                                                     true_coverage_config['minimum_gene_identity'],
                                                     true_coverage_config, rematch_script)
                    runs['trueCoverage_ReMatCh'] = [run_successfully_true_coverage, pass_qc_true_coverage, time_taken,
                                                    failing, {}, 'NA']
                else:
                    print '\n' + '--skipTrueCoverage set. Skipping True coverage analysis'
                    runs['trueCoverage_ReMatCh'] = skipped

                if args.skipTrueCoverage or true_coverage_config is None or (run_successfully_true_coverage and
                                                                             pass_qc_true_coverage):
                    # Run first FastQC
                    nts2clip_based_nts_content = None
                    if not args.skipFastQC:
                        run_successfully, pass_qc, time_taken, failing, warning, maximum_reads_length, \
                            nts2clip_based_nts_content = fastqc.runFastQCanalysis(outdir, threads, adapters_fasta,
                                                                                  fastq_files, args.fastQCkeepFiles,
                                                                                  'first_run')
                        runs['first_FastQC'] = [run_successfully, pass_qc, time_taken, failing, warning, 'NA']
                    else:
                        print '--skipFastQC set. Skipping First FastQC analysis'
                        runs['first_FastQC'] = skipped

                    # Run Trimmomatic
                    not_empty_fastq = True
                    if not args.skipTrimmomatic:
                        run_successfully, not_empty_fastq, time_taken, failing, paired_reads, trimmomatic_folder, \
                            file_size, warning = trimmomatic.runTrimmomatic(jar_path_trimmomatic, sample_name, outdir,
                                                                            threads, adapters_fasta, script_path,
                                                                            args.doNotSearchAdapters, fastq_files,
                                                                            max_reads_length, args.doNotTrimCrops,
                                                                            args.trimCrop, args.trimHeadCrop,
                                                                            args.trimLeading, args.trimTrailing,
                                                                            args.trimSlidingWindow, args.trimMinLength,
                                                                            nts2clip_based_nts_content, jar_max_memory,
                                                                            fastq_encoding)
                        runs['Trimmomatic'] = [run_successfully, None, time_taken, failing, warning, file_size]
                        trimmomatic_run_successfully = run_successfully

                        if run_successfully and not_empty_fastq:
                            fastq_files = paired_reads
                            min_reads_length = args.trimMinLength

                            # Run second Estimated Coverage
                            if not args.skipEstimatedCoverage:
                                run_successfully_estimated_coverage, pass_qc, time_run, failing, estimated_coverage = \
                                    coverage.getEstimatedCoverage(fastq_files, genome_size, outdir, threads,
                                                                  args.estimatedMinimumCoverage)
                                runs['second_Coverage'] = [run_successfully_estimated_coverage, pass_qc, time_run,
                                                           failing, {}, 'NA']
                            else:
                                print '--skipEstimatedCoverage set. Skipping Second Estimated Coverage analysis'
                                runs['second_Coverage'] = skipped

                            if args.skipEstimatedCoverage or (run_successfully_estimated_coverage and
                                                              not estimated_coverage < args.estimatedMinimumCoverage):
                                # Run second FastQC
                                if not args.skipFastQC:
                                    run_successfully, pass_qc, time_taken, failing, warning, maximum_reads_length, \
                                        nts2clip_based_nts_content = fastqc.runFastQCanalysis(outdir, threads,
                                                                                              adapters_fasta,
                                                                                              fastq_files,
                                                                                              args.fastQCkeepFiles,
                                                                                              'second_run')
                                    runs['second_FastQC'] = [run_successfully, pass_qc, time_taken, failing, warning,
                                                             'NA']
                                    if run_successfully:
                                        max_reads_length = maximum_reads_length
                                else:
                                    print '--skipFastQC set. Skipping Second FastQC analysis'
                                    runs['second_FastQC'] = skipped
                            else:
                                print('\n'
                                      'Estimated coverage is too lower (< {estimatedMinimumCoverage}x). This sample'
                                      ' will not proceed with INNUca'
                                      ' pipeline'.format(estimatedMinimumCoverage=args.estimatedMinimumCoverage))
                        else:
                            print('Trimmomatic did not run successfully or return zero reads! Skipping Second Estimated'
                                  ' Coverage analysis and FastQC analysis')
                            runs['second_Coverage'] = skipped
                            runs['second_FastQC'] = skipped
                    else:
                        print('--skipTrimmomatic set. Skipping Trimmomatic, but also Second FastQC analysis and Second'
                              ' Estimated Coverage analysis')
                        runs['Trimmomatic'] = skipped
                        runs['second_Coverage'] = skipped
                        runs['second_FastQC'] = skipped

                    if not args.skipFastQC and \
                            (runs['second_FastQC'][1] or (runs['second_FastQC'][1] is None and
                                                          runs['first_FastQC'][1])) is False and \
                            not not_empty_fastq and not args.fastQCproceed:
                        print('\n'
                              'This sample does not pass FastQC module QA/QC. It will not proceed with INNUca pipeline')
                else:
                    print('\n'
                          'This sample does not pass True Coverage module QA/QC. This sample will not proceed with'
                          ' INNUca pipeline')
            else:
                print('\n'
                      'Estimated coverage is too lower (< {estimatedMinimumCoverage}x). This sample will not proceed'
                      ' with INNUca pipeline'.format(estimatedMinimumCoverage=args.estimatedMinimumCoverage))

        continue_second_part = False
        if not args.runKraken or \
                (runs['reads_Kraken'][0] is True and runs['reads_Kraken'][1] is True) or \
                args.krakenProceed or \
                args.krakenIgnoreQC:
            if args.skipEstimatedCoverage or (run_successfully_estimated_coverage and
                                              not estimated_coverage < args.estimatedMinimumCoverage):
                if args.skipTrueCoverage or true_coverage_config is None or (run_successfully_true_coverage and
                                                                             pass_qc_true_coverage):
                    if args.skipFastQC or (runs['second_FastQC'][1] or
                                           (runs['second_FastQC'][1] is None and
                                            runs['first_FastQC'][1])) is not False or \
                            args.fastQCproceed:
                        continue_second_part = True

        if continue_second_part:
            unassembled_pe_reads = None
            assembled_se_reads = None
            # Run Pear
            if args.runPear:
                print '--runPear set. Running Pear'
                pear_min_overlap = pear.determine_minimum_overlap(args.pearMinOverlap, min_reads_length,
                                                                  max_reads_length)
                run_successfully, pass_qc, time_taken, failing, unassembled_pe_reads, assembled_se_reads, \
                    pear_folder, warning = pear.runPear(fastq_files, threads, outdir, sample_name,
                                                        fastq_encoding, trimmomatic_run_successfully,
                                                        pear_min_overlap)
                runs['Pear'] = [run_successfully, pass_qc, time_taken, failing, warning, 'NA']
            else:
                runs['Pear'] = skipped

            # Run SPAdes
            if not args.skipSPAdes:
                run_successfully, pass_qc, time_taken, failing, contigs_spades, warning = \
                    spades.runSpades(sample_name, outdir, threads,
                                     unassembled_pe_reads if unassembled_pe_reads is not None else fastq_files,
                                     args.spadesNotUseCareful, spades_max_memory,
                                     args.spadesMinCoverageAssembly, args.spadesMinContigsLength, genome_size,
                                     args.spadesKmers, max_reads_length, args.spadesDefaultKmers,
                                     args.spadesMinKmerCovContigs, assembled_se_reads, args.saveExcludedContigs,
                                     args.maxNumberContigs)
                runs['SPAdes'] = [run_successfully, pass_qc, time_taken, failing, warning, 'NA']

                if run_successfully:
                    contigs = contigs_spades

                    # Run Assembly Mapping check
                    bam_file = None
                    original_bam = None
                    assembly_mapping_folder = None
                    possible_assemblies_bam_remove = {}
                    if not args.skipAssemblyMapping:
                        run_successfully, pass_qc, time_taken, failing, assembly_filtered, bam_file, \
                            assembly_mapping_folder, warning, original_bam = \
                            assembly_mapping.run_assembly_mapping(fastq_files=fastq_files, reference_file=contigs,
                                                                  outdir=outdir, estimated_genome_size_mb=genome_size,
                                                                  max_number_contigs=args.maxNumberContigs,
                                                                  save_excluded_contigs=args.saveExcludedContigs,
                                                                  min_coverage_assembly=args.assemblyMinCoverageContigs,
                                                                  keep_bam=args.keepBAM, threads=threads)
                        runs['Assembly_Mapping'] = [run_successfully, pass_qc, time_taken, failing, warning,
                                                    'NA']

                        if run_successfully:
                            # Assembly to remove
                            if not args.keepIntermediateAssemblies:
                                if os.path.isfile(contigs_spades) and \
                                        assembly_filtered is not None and \
                                        assembly_filtered != contigs_spades:
                                    if not args.keepBAM:
                                        os.remove(contigs_spades)
                                    else:
                                        possible_assemblies_bam_remove['assembly_mapping'] = contigs_spades

                            if assembly_filtered is not None and \
                                    assembly_filtered != contigs_spades and \
                                    os.path.isfile(assembly_filtered):
                                contigs = assembly_filtered
                    else:
                        print '--skipAssemblyMapping set. Skipping Assembly Mapping check'
                        runs['Assembly_Mapping'] = skipped

                    # Run Pilon
                    pilon_new_bam = False
                    pilon_bam = None
                    if not args.skipPilon:
                        run_successfully, _, time_taken, failing, assembly_polished, pilon_folder, pilon_new_bam, \
                            pilon_bam = pilon.run_pilon(jar_path_pilon=jar_path_pilon, assembly=contigs,
                                                        fastq_files=fastq_files, outdir=outdir,
                                                        jar_max_memory=jar_max_memory, alignment_file=bam_file,
                                                        keep_bam=args.keepBAM, threads=threads)
                        runs['Pilon'] = [run_successfully, None, time_taken, failing, {}, 'NA']

                        if run_successfully:
                            if not args.keepIntermediateAssemblies:
                                if os.path.isfile(contigs) and \
                                        assembly_polished is not None and \
                                        os.path.isfile(assembly_polished):
                                    if not args.keepBAM:
                                        os.remove(contigs)
                                    else:
                                        if not pilon_new_bam:
                                            possible_assemblies_bam_remove['pilon'] = contigs

                            if assembly_polished is not None and \
                                    os.path.isfile(assembly_polished):
                                contigs = assembly_polished

                        if not args.pilonKeepFiles and os.path.isdir(pilon_folder):
                            utils.removeDirectory(pilon_folder)

                    else:
                        print '--skipPilon set. Skipping Pilon correction'
                        runs['Pilon'] = skipped

                    if not args.keepBAM:
                        if bam_file is not None:
                            if os.path.isfile(bam_file):
                                os.remove(bam_file)
                            if os.path.isfile(bam_file + '.bai'):
                                os.remove(bam_file + '.bai')

                        if original_bam is not None and os.path.isfile(original_bam):
                            os.remove(original_bam)

                        if pilon_bam is not None and os.path.isfile(pilon_bam):
                            os.remove(pilon_bam)

                        if 'assembly_mapping' in possible_assemblies_bam_remove and \
                                os.path.isfile(possible_assemblies_bam_remove['assembly_mapping']):
                            os.remove(possible_assemblies_bam_remove['assembly_mapping'])
                        if 'pilon' in possible_assemblies_bam_remove and \
                                os.path.isfile(possible_assemblies_bam_remove['pilon']):
                            os.remove(possible_assemblies_bam_remove['pilon'])
                    else:
                        if pilon_new_bam:
                            if bam_file is not None:
                                if os.path.isfile(bam_file):
                                    os.remove(bam_file)
                                if os.path.isfile(bam_file + '.bai'):
                                    os.remove(bam_file + '.bai')

                            if original_bam is not None and os.path.isfile(original_bam):
                                os.remove(original_bam)

                            if 'assembly_mapping' in possible_assemblies_bam_remove and \
                                    os.path.isfile(possible_assemblies_bam_remove['assembly_mapping']):
                                os.remove(possible_assemblies_bam_remove['assembly_mapping'])
                        else:
                            if original_bam is not None and os.path.isfile(original_bam) and \
                                    bam_file is not None and os.path.isfile(bam_file):
                                os.remove(bam_file)
                            if 'pilon' in possible_assemblies_bam_remove and \
                                    os.path.isfile(possible_assemblies_bam_remove['pilon']):
                                os.remove(possible_assemblies_bam_remove['pilon'])

                    if not args.skipAssemblyMapping:
                        utils.removeDirectory(assembly_mapping_folder)

                    print '\n' + 'Final assembly: ' + contigs
                    with open(os.path.join(outdir, 'final_assembly.txt'), 'wt') as writer:
                        writer.write(contigs + '\n')

                    # Run MLST
                    if not args.skipMLST:
                        run_successfully, pass_qc, time_taken, failing, warning = \
                            mlst.runMlst(contigs, scheme, outdir, species_genus, mlst_scheme_genus)
                        runs['MLST'] = [run_successfully, pass_qc, time_taken, failing, warning, 'NA']
                    else:
                        print '--skipMLST set. Skipping MLST analysis'
                        runs['MLST'] = skipped

                    # Run Kraken
                    if args.runKraken:
                        print('\n'
                              '--runKraken set. Running Kraken for assembly')
                        run_successfully, pass_qc, time_taken, failing, warning, _ = \
                            kraken(species=args.speciesExpected, files_to_classify=[contigs], kraken_db=args.krakenDB,
                                   files_type='fasta', outdir=outdir, version_kraken=version_kraken_global,
                                   db_mem=args.krakenMemory, quick=args.krakenQuick,
                                   min_percent_covered=args.krakenMinCov,
                                   max_unclassified_frag=args.krakenMaxUnclass, min_base_quality=args.krakenMinQual,
                                   threads=threads)
                        runs['assembly_Kraken'] = [run_successfully, pass_qc, time_taken, failing, warning, 'NA']
                    else:
                        runs['assembly_Kraken'] = skipped

                else:
                    print('SPAdes did not run successfully! Skipping Pilon correction, Assembly Mapping check'
                          ' and MLST analysis')
            else:
                print('--skipSPAdes set. Skipping SPAdes, Pilon correction, Assembly Mapping check and MLST'
                      ' analysis')
                runs['SPAdes'] = skipped
                runs['Assembly_Mapping'] = skipped
                runs['Pilon'] = skipped
                runs['MLST'] = skipped
    else:
        print 'Moving to the next sample'

    for step in ('reads_Kraken', 'first_Coverage', 'trueCoverage_ReMatCh', 'first_FastQC', 'Trimmomatic',
                 'second_Coverage', 'second_FastQC', 'Pear', 'SPAdes', 'Assembly_Mapping', 'Pilon', 'MLST',
                 'assembly_Kraken'):
        if step not in runs:
            runs[step] = not_run

    # Remove Pear directory
    if not args.pearKeepFiles and pear_folder is not None:
        utils.removeDirectory(pear_folder)
    # Remove Trimmomatic directory with cleaned reads
    if not args.trimKeepFiles and trimmomatic_folder is not None:
        utils.removeDirectory(trimmomatic_folder)

    # Check run
    run_successfully = all(runs[step][0] or runs[step][0] is None for step in runs)

    pass_fastq_integrity = runs['FastQ_Integrity'][0]
    pass_reads_kraken = runs['reads_Kraken'][1] is not False or args.krakenIgnoreQC
    pass_cov = (runs['second_Coverage'][1] or (runs['second_Coverage'][1] is None and
                                               runs['first_Coverage'][1])) is not False
    pass_true_cov = runs['trueCoverage_ReMatCh'][1] is not False
    pass_fastqc = (runs['second_FastQC'][1] or (runs['second_FastQC'][1] is None and
                                                runs['first_FastQC'][1])) is not False
    # pass_trimmomatic = runs['Trimmomatic'][1] is not False
    # pass_pear = runs['Pear'][1] is not False
    # pass_spades = runs['SPAdes'][1] is not False or runs['Assembly_Mapping'][1] is True
    pass_spades = runs['SPAdes'][1] is not False
    pass_assembly_mapping = runs['Assembly_Mapping'][1] is not False
    pass_pilon = runs['Pilon'][0] is not False
    pass_mlst = runs['MLST'][1] is not False or args.mlstIgnoreQC
    pass_assembly_kraken = runs['assembly_Kraken'][1] is not False or args.krakenIgnoreQC
    pass_qc = all([pass_fastq_integrity, pass_reads_kraken, pass_cov, pass_true_cov, pass_fastqc, pass_spades,
                   pass_assembly_mapping, pass_pilon, pass_mlst, pass_assembly_kraken])

    return run_successfully, pass_qc, runs


if __name__ == "__main__":
    main()