#!/usr/bin/env python3
__author__ = 'LieuV'
__version__ = "2.1.1"

import argparse
import re
import getpass
import warnings
import os
import json
import sys
from itertools import groupby
from collections import OrderedDict
import subprocess
from pandas import read_csv
import snakemake
import datetime
import yaml
import input_utils
from sys import stderr
from Bio import SeqIO
# Set paths

ADAPTER_PATH = ""   # Adapter path from trimmomatic.
                    # Should be set during installation

KRAKEN_DB = subprocess.check_output("echo $KRAKEN_DB_PATH",
                                    shell=True).decode("utf-8").strip()
if len(KRAKEN_DB) == 0:
    KRAKEN_DB = " "  # insert space such that snakemake can handle empty value

# Pattern for Illumina readnames
base_pattern = (
    r"(?P<setname>.*)_(?P<lane>L\d{3})_(?P<read>(R1|R2))_(?P<nfiles>\d{3}).*")

# --------------------------------------< Functions >----------------------------------------------------------


def main():
    parser = argparse.ArgumentParser(
        formatter_class=argparse.RawDescriptionHelpFormatter, description="""\
    --------------------------------------------------------------------

                             < QCumber >
              < Quality control and read trimming of NGS data >
           https://gitlab.com/RKIBioinformaticsPipelines/QCumber/

    --------------------------------------------------------------------""",
        epilog=("Example usage: QCumber-2 --input fastq_folder"
                " --reference reference.fasta"))

    # ------------------------------------------------< INPUT >---------------------------------------------------------#

    group_input = parser.add_argument_group("Input")
    group_input.add_argument(
        '--input', '-i', dest='input',
        help=("input sample folder. Illumina filenames should end with"
              "_<lane>_<R1|R2>_number, e.g. Sample_12_345_R1_001.fastq,"
              " to find the right paired set."),
        required=False, nargs="+")
    group_input.add_argument(
        '--read1', '-1', dest='r1', help="Read 1 file", required=False)
    group_input.add_argument(
        '--read2', '-2', dest='r2', help="Read 2 file", required=False)
    group_input.add_argument(
        '--technology', '-T', dest='technology', choices=["Illumina",
                                                          "IonTorrent",
                                                          "PacBio"],
        required=False,
        help=("If not set, automatically determine technology and "
              "search for fastq and bam files. "
              "Set technology to IonTorrent if all files are bam-files,"
              " else set technology to Illumina."))
    group_input.add_argument(
        '--adapter', '-a', dest='adapter',
        choices=['TruSeq2-PE', 'TruSeq2-SE', 'TruSeq3-PE',
                 'TruSeq3-SE', 'TruSeq3-PE-2', 'NexteraPE-PE'],
        help="Adapter name for trimming. Default: all")

    mapping_exclusion = group_input.add_mutually_exclusive_group()
    mapping_exclusion.add_argument(
        '--reference', '-r', dest='reference', required=False,
        help=("Map reads against reference."
              + " Reference needs to be in fasta-format."))
    mapping_exclusion.add_argument(
        '--index', '-I', dest='index', required=False,
        help="Bowtie2 index if available.")

    group_input.add_argument(
        '--kraken_db', '-d', dest='kraken_db',
        help=("Custom Kraken database. Default value is taken from"
              " environment variable KRAKEN_DB_PATH. "
              "Default: %(default)s."),
        required=False, default=KRAKEN_DB)

    group_input.add_argument(
        '--kraken_classified_out', dest='kraken_classified_out',
        help=("Kraken (un)classified-out option."
              " If set, both the --classified-out"
              " and --unclassified-out option are set. "
              "Default: %(default)s."),
        required=False, default=False, action='store_true')

    group_optional = parser.add_argument_group("Optional steps")
    group_optional.add_argument(
        '--sav', '-w', dest='sav', required=False,
        help=("Illumina folder for SAV. Requires RunInfo.xml, RunParameter.xml"
              "and Interop folder."))
    group_optional.add_argument(
        '--trimBetter', choices=["assembly", "mapping", "default"],
        help=("Optimize trimming parameter using 'Per sequence base content'"
              + " from fastqc. Not recommended for amplicons."))
    group_optional.add_argument('--nokraken', '-K', action="store_true")
    group_optional.add_argument('--notrimming', '-Q', action="store_true")

    group_params = parser.add_argument_group("Parameter settings")
    group_params.add_argument(
        '--illuminaclip', '-L', dest='illuminaclip', default="2:30:10",
        help=('Illuminaclip option: '
              '<leading quality>:<trailing quality>:<sliding window>.'
              'Default: %(default)s'))
    group_params.add_argument(
        '--only_trim_adapters', '-A', action='store_true',
        help='If this option is selected, only adapters will be clipped')
    group_params.add_argument(
        '--minlen', '-m', default=50, dest='minlen',
        help=('Minlen parameter for Trimmomatic. Drops read short than minlen.'
              ' Default: %(default)s'),
        type=int)
    group_params.add_argument(
        '--trimOption', '-O', dest="trimOption",
        help=('Additional Trimmomatic input.'
              ' Default (if trimBetter is not set): SLIDINGWINDOW:4:20'),
        type=str)
    group_params.add_argument(
        '--trimBetter_threshold', '-b', dest='trimBetter_threshold',
        help=("Set -trimBetter to use this option.Default setting"
              " for Illumina: 0.15 and for IonTorrent: 0.25."),
        required=False, type=float)
    group_output = parser.add_argument_group("Output")
    group_output.add_argument('--output', '-o', dest='output', default="")
    group_output.add_argument(
        '--rename', '-R', dest="rename", required=False,
        help="TSV File with two columns: <old sample name> <new sample name>")
    group_output.add_argument(
        '--save_mapping', '-S', action="store_true", default=False)

    parser.add_argument('--threads', '-t', dest='threads', default=4, type=int,
                        help="Number of threads. Default:  %(default)s")

    parser.add_argument(
        '--config', '-c', dest='config',
        help=("Configfile to run pipeline. "
              "Additional parameters in the commandline "
              "will override arguments in configfile."
              "If not given and config/config.txt exists in"
              "the directory of the QCumber-2 executable,"
              "that file will be loaded by default."))
    parser.add_argument('--version', '-v',
                        action='version', version='%(prog)s v' + __version__)

    arguments, unknown_args = parser.parse_known_args()
    arguments = vars(arguments)
    if len(sys.argv) == 1:
        parser.print_help()
        sys.exit(1)

    configfile = arguments["config"]
    default_configfile = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            "config",
            "config.txt")
    if not arguments["config"] and os.path.isfile(default_configfile):
        configfile = default_configfile

    if configfile:
        config_args = yaml.load(open(configfile, "r"))
        keep_args = dict()
        for arg in config_args.keys():
            print(arg)
            if (arguments[arg] is None) or arguments[arg] == " ":
                arguments[arg] = config_args[arg]

    arguments["output"] = os.path.abspath(arguments["output"])
    if arguments["only_trim_adapters"]:
        arguments["trimBetter"] = None
    if not os.path.isdir(arguments["output"]):
        os.mkdir(arguments["output"])
    if arguments["reference"]:
        arguments["reference"] = os.path.abspath(arguments["reference"])
    if arguments["sav"]:
        arguments["sav"] = os.path.abspath(arguments["sav"])
    if arguments["rename"]:
        arguments["rename"] = os.path.abspath(arguments["rename"])

    parameter = yaml.load(
        open(os.path.join(
                os.path.dirname(os.path.realpath(__file__)),
                "config",
                "parameter.txt"),
             "r"))
    check_input_validity(arguments)
    # Load adaptive filetypes
    qcumber_path = os.path.dirname(os.path.realpath(__file__))
    sample_file_name = os.path.join(arguments["output"], "samples.yaml")
    with open(os.path.join(qcumber_path,
                           'filenames.yaml'),
              'r') as filetype_h:
        filename_types = yaml.load(filetype_h)
    all_files = []
    for file_or_dir in arguments["input"]:
        if os.path.isdir(file_or_dir):
            for root, dirs, files in os.walk(file_or_dir):
                for file in files:
                    if os.path.getsize(os.path.join(root, file)) != 0:
                        all_files.append(os.path.join(root, file))
        elif os.path.isfile(file_or_dir):
            if os.path.getsize(file_or_dir):
                all_files.append(file_or_dir)
    # Get Parsed Samples fomr input utils module
    formats_found,  discarded = input_utils.parse_sample_info(
            all_files, filename_types, ['pacbio', 'illumina_fastq'])
    try:
        illumina_data = formats_found['illumina_fastq']
        # print(repr(format_known.mfrs).replace('>,', '>,\n'))
    except KeyError:
            exit('No samples found or none met criteria!!\n'
                 'These files were discarded:\n'
                 '%s' % '\n'.join(discarded))
    try:
        pacbio_data = formats_found['pacbio']
        print('looking for pacbio data...')
        pac_samples = pacbio_data.get_samples()
        with open(sample_file_name.replace('.yaml',
                                           '_pacbio.yaml'),
                  'w') as sample_file:
            yaml.dump(pac_samples, sample_file, default_flow_style=False)
        print('looking for illumina data...')
    except KeyError:
        pass
    sample_dict = illumina_data.flatten_naive()
    # flatten naive just drops read info from sampel name
    # if it is paired end data.
    len_known = len(sample_dict)
    if False:
        try:
            salvaged_dict = illumina_data.leftovers.process_leftovers(
                rename=True,
                rename_start_index=len_known+1)
            sample_dict.update(salvaged_dict)
        except input_utils.AmbigiousPairedReadsError as err:
            eprint('Failed parsing files with unrecognized'
                   ' naming convention\n',
                   'Reason:\n', err)
    # Write samples to working directory
    with open(sample_file_name, 'w') as sample_file:
        yaml.dump(sample_dict, sample_file, default_flow_style=False)

    type, samples, joined_samples, name_dict, join_reads = (
        get_input(arguments, parameter))
    get_defaults(arguments, parameter)

    force_run_list = []
    os.makedirs(arguments["output"], exist_ok=True)
    config_file_path = os.path.join(arguments["output"], "config.yaml")
    #if os.path.isfile(config_file_path):
    #    pass
    #else:
    with open(config_file_path, 'w') as config_fh:
        yaml.dump(arguments, config_fh,
                  default_flow_style=False)

    # additional infos
    general_information = OrderedDict()
    general_information["User"] = getpass.getuser()
    general_information["QCumber"] = __version__
    general_information["QCumber_path"] = os.path.dirname(
                                            os.path.realpath(__file__))
    general_information["Execution time"] = datetime.datetime.now().ctime()

    system_info = os.uname()
    general_information["Operating system"] = OrderedDict()
    general_information["Operating system"]["System"] = system_info.sysname
    general_information["Operating system"]["Server"] = system_info.nodename
    general_information["Operating system"]["Operating version"] = (
        system_info.version)
    general_information["Operating system"]["Release"] = system_info.release
    general_information["Operating system"]["Machine"] = system_info.machine
    general_information["Tool versions"] = OrderedDict()
    general_information["Tool versions"]["Python"] = re.sub("\n",
                                                            "", sys.version)
    general_information["Tool versions"]["Snakemake"] = snakemake.__version__
    general_information["Tool versions"]["FastQC"] = (
        get_version("fastqc --version"))

    if not arguments["notrimming"]:
        general_information["Tool versions"]["Trimmomatic"] = (
            get_version("trimmomatic -version", "trimmomatic"))
        if arguments["technology"] == "Illumina" and not arguments["adapter"]:
            general_information["adapter"] = (
                os.path.join(os.path.dirname(os.path.realpath(__file__)),
                             "config", "adapters.fa"))
        elif arguments["technology"] == "Illumina" and arguments["adapter"]:
            general_information["adapter"] = os.path.join(ADAPTER_PATH,
                                                          (arguments["adapter"]
                                                           + ".fa"))
    if arguments["reference"] or arguments["index"]:
        general_information["Tool versions"]["Bowtie2"] = (
            get_version("bowtie2 --version"))
    if not arguments["nokraken"]:
        general_information["Tool versions"]["Kraken"] = (
            get_version("kraken --version"))
    general_information["Sample information"] = OrderedDict()
    general_information["Sample information"]["type"] = type
    general_information["Sample information"]["samples"] = samples
    general_information["Sample information"]["join_reads"] = join_reads
    general_information["Sample information"]["join_lanes"] = joined_samples

    general_information["Sample information"]["rename"] = name_dict
    os.makedirs(os.path.join(arguments["output"],
                             "QCResults", "_data"), exist_ok=True)

    general_information_file = os.path.join(arguments["output"],
                                            "QCResults", "_data",
                                            "general_information.json")
    json.dump(general_information, open(general_information_file, "w"))

    # Fixed:
    # if QCumber is run repeatedly, the rule bowtie_mapping
    # will not be executed.
    # hence it is necessary to force to run the rule everytime
    #  if (arguments["save_mapping"]):
    #      force_run = "--forcerun bowtie_mapping"
    #  else:
    force_run = "--forcerun"
    cmd_string = (
        "snakemake "
        "--configfile {workdir}/config.yaml "
        "--snakefile {snakefile} {additional_commands} "
        "--directory {workdir} "
        "--cores {cores} {targets} {force_run} "
        ).format(
            additional_commands=" ".join(unknown_args),
            snakefile=os.path.join(os.path.dirname(os.path.realpath(__file__)),
                                   "Snakefile"),
            workdir=arguments["output"],
            configfile=general_information_file,
            cores=arguments["threads"],
            targets='',  # "QCResults/batch_report.html",
            force_run=force_run)
    print(cmd_string)
    process = subprocess.Popen(
        ("snakemake "
         "--configfile {workdir}/config.yaml "
         "--snakefile {snakefile} {additional_commands} "
         "--directory {workdir} "
         "--cores {cores} {targets} {force_run} "
         ' -R $(snakemake --list-params-changes '
         #  '--list-input-changes --list-code-changes '
            '--configfile {workdir}/config.yaml '
            '--snakefile {snakefile} {additional_commands} '
            '--directory {workdir} '
            '--cores {cores} {targets} {force_run})'
         ).format(
            additional_commands=" ".join(unknown_args),
            snakefile=os.path.join(os.path.dirname(os.path.realpath(__file__)),
                                   "Snakefile"),
            workdir=arguments["output"],
            configfile=general_information_file,
            cores=arguments["threads"],
            targets='',  # "QCResults/batch_report.html",
            force_run=force_run),
        shell=True)
    process.wait()
    exit(process.returncode)


def get_basename(abs_name):
    return os.path.basename(os.path.splitext(
        os.path.splitext(os.path.basename(abs_name))[0])[0])


def getFilenameWithoutExtension(string, getBase=False):
    if getBase:
        string = os.path.basename(string)
    string = os.path.splitext(string)[0]
    i = 0
    while os.path.splitext(string)[-1] in [
            ".gz", ".gzip", ".zip", ".bz", ".fasta", ".fastq", ".bam"]:
        string = os.path.splitext(string)[0]
    return string


def get_setname(filename, base=True, grouping=True):
    """ Get Setname

    Args:
        filename (obj::`str`): filename;

    Kwargs:
        base (bool): transform filename to basename (default = True)
        grouping (bool): if not grouping, only return basename of file
            (default = True)
    Returns:
    """
    if not grouping:
        return get_basename(filename)
    try:
        if base:
            filename = get_basename(filename)
        sep = iter([";", "\t"])
        if arguments["rename"]:
            rename_dict = read_csv(arguments["rename"],
                                   index_col=0, header=None)

            while (rename_dict.columns.__len__() == 0
                   or sep.__length_hint__() != 0):
                rename_dict = read_csv(arguments["rename"], index_col=0,
                                       header=None, sep=next(sep))
            if rename_dict.columns.__len__() == 0:
                warnings.warn(
                    ("Problems reading rename-file %s. "
                     "Valid delimiters are"
                     " ';', ',' or '\\t'.") % arguments["rename"],
                    UserWarning, stacklevel=2)

            new_name = [(x, rename_dict.ix[x][1]) for x
                        in rename_dict.index
                        if filename.startswith(x)]
            if len(new_name) != 1:
                warnings.warn(
                    ("Could not find unique renames."
                     " Found %s for %s") % (new_name, filename),
                    UserWarning, stacklevel=2)
                filename = getFilenameWithoutExtension(filename, True)
            else:
                filename = filename.replace(new_name[0][0], new_name[0][1])
    except:
        # print("Couldnt rename sample files.")
        pass
    try:
        paired_reads_pattern = base_pattern
        setname_pattern = re.search(paired_reads_pattern,
                                    os.path.basename(filename))
        if setname_pattern:
            return setname_pattern.group("setname")
        else:
            return filename
    except:
        print("Problems getting samplenames: %s" % filename)
        return filename


def get_defaults(arguments, parameter):
    if arguments["only_trim_adapters"]:
        return
    if arguments["trimBetter"] == "assembly":
        # if arguments["forAssembly"]:
        if not arguments["trimBetter_threshold"]:
            arguments["trimBetter_threshold"] = (
                parameter["forAssembly." + arguments['technology']]
                         ['trimBetter_threshold'])
        if not arguments["trimOption"]:
            arguments["trimOption"] = (
                    parameter["forAssembly." + arguments['technology']]
                             ["trimOption"])

    elif arguments["trimBetter"] == "mapping":
        # elif arguments["forMapping"]:
        if not arguments["trimBetter_threshold"]:
            arguments["trimBetter_threshold"] = (
                parameter["forMapping." + arguments['technology']]
                         ["trimBetter_threshold"])
        if not arguments["trimOption"]:
            arguments["trimOption"] = (
                parameter["forMapping." + arguments['technology']]
                         ["trimOption"])
    elif arguments["trimBetter"] == "default":
        arguments["trimBetter_threshold"] = (parameter["Trimmomatic"]
                                                      ["trimBetter_threshold"])

    if not arguments["trimOption"]:
        arguments["trimOption"] = parameter["Trimmomatic"]["trimOption"]

    if arguments["trimBetter_threshold"]:
        arguments["trimBetter_threshold"] = (parameter["Trimmomatic"]
                                                      ["trimBetter_threshold"])


def get_version(cmd, jar=None):
    try:
        return re.match(r"(?P<version>.*)\n",
                        subprocess.check_output(
                            cmd, shell=True, stderr=subprocess.PIPE
                        ).decode("utf-8")).group("version")
    except:
        try:
            try:
                return subprocess.check_output(
                        cmd, shell=True,
                        stderr=subprocess.PIPE
                    ).decode("utf-8")
            except:
                from debpackageinfo import get_upstream_version
                return get_upstream_version(jar)
        except:
            return "NaN"


def get_input(arguments, parameter):
    """
    Get Read Input

    Checks wether input reads are from Ion Torrent or Illumina, by checking
    file extensions of bam and/or fastq. Detects presence of read pair in input
    and returns type of reads used (single end / paired end reads)

    Args:
        None

    Returns:
        (obj::`str`): type - PE or SE  (Paired Ends / Single Ends)
        (obj::`collections.OrderedDict`): sample_dict
        (obj::`collections.OrderedDict`): join_lanes
        (obj::`dict`): name_dict
        (obj::`dict`): join_reads
    """
    bam_ext = [x.strip(" ") for x in parameter["Fileextension"]["bam"]]
    fastq_ext = [x.strip(" ") for x in parameter["Fileextension"]["fastq"]]

    sample_dict = OrderedDict()
    all_files = []
    name_dict = {}
    join_reads = {}
    type = "PE"
    join_lanes = OrderedDict()

    if arguments["r1"]:
        assert os.path.getsize(arguments["r1"]) != 0, (
               "File %s is empty." % arguments["r1"])

        if any([arguments["r1"].endswith(ext) for ext in bam_ext]):
            arguments["technology"] = "IonTorrent"
        else:
            arguments["technology"] = "Illumina"
        if arguments["r2"]:
            assert os.path.getsize(arguments["r2"]) != 0, (
                   "File %s is empty." % arguments["r2"])

            sample_dict[get_setname(arguments["r1"])] = (
                [os.path.abspath(arguments["r1"]),
                 os.path.abspath(arguments["r2"])])
            name_dict[get_basename(arguments["r1"])] = (
                get_setname(arguments["r1"]) + "_R1")
            name_dict[get_basename(arguments["r2"])] = (
                get_setname(arguments["r2"]) + "_R2")
        else:
            type = "SE"
            sample_dict[get_setname(arguments["r1"])] = (
                [os.path.abspath(arguments["r1"])])
            name_dict[get_basename(arguments["r1"])] = (
                get_setname(arguments["r1"]))
    else:
        if os.path.isdir(arguments["input"][0]):
            for root, dirs, files in os.walk(arguments["input"][0]):
                for file in files:
                    if any([file.endswith(ext)
                            for ext in fastq_ext + bam_ext]):
                        if os.path.getsize(os.path.join(root, file)) != 0:
                            all_files.append(os.path.join(root, file))
                        else:
                            warnings.warn("Skip empty file %s" % file,
                                          stacklevel=2)
        else:
            all_files = arguments["input"]
            assert all_files, (
                   ("Check input again. "
                    "No files found for pattern %s ") % arguments["input"])

        if len([x for x
                in all_files
                if any([ext in x for ext in bam_ext])
                ]
               ) == len(all_files):
            arguments["technology"] = "IonTorrent"
        else:
            arguments["technology"] = "Illumina"

        if (len(all_files) == 0):
            sys.exit(str(arguments["input"])
                     + " does not contain fastq or bam files.")

        # find read pairs
        all_files = sorted(list(all_files))

        if all([re.search(base_pattern, x) for x in all_files]):
            for setname, files in groupby(all_files,
                                          key=lambda x: re.search(
                                              base_pattern,
                                              x
                                          ).group("setname")):
                read_pairs = dict()
                setname = get_setname(setname)
                for lane, lane_file in groupby(list(files),
                                               key=lambda x: re.search(
                                                                 base_pattern,
                                                                 x
                                                             ).group("lane")):
                    read_pairs[lane] = []
                    for readgroup, readfiles in groupby(
                            list(lane_file),
                            key=lambda x: re.search(base_pattern, x
                                                    ).group("read")):
                        readfiles = list(readfiles)
                        if len(readfiles) != 0:
                            if len(readfiles) > 1:
                                concat_reads = (
                                        "QCResults/tmp/join_reads/"
                                        + "_".join(
                                            re.search(base_pattern,
                                                      os.path.basename(
                                                          readfiles[0])
                                                      ).groups()[:-2])
                                        + "_000.fastq.gz")
                                join_reads[concat_reads] = [os.path.abspath(x)
                                                            for x in readfiles]
                                readfiles = concat_reads

                            else:
                                readfiles = os.path.abspath(readfiles[0])
                            read_pairs[lane].append(readfiles)
                # Multiple lanes
                if len(read_pairs) > 1:
                    join_lanes[setname] = []
                for key in sorted(read_pairs.keys()):
                    samplename = setname + "_" + key
                    sample_dict[samplename] = read_pairs[key]
                    if len(read_pairs[key]) == 2:  # type == "PE":
                        name_dict[get_basename(read_pairs[key][0])] = (
                            get_setname(samplename) + "_R1")
                        name_dict[get_basename(read_pairs[key][1])] = (
                            get_setname(samplename) + "_R2")
                    else:
                        type = "SE"
                        name_dict[get_basename(read_pairs[key][0])] = (
                            get_setname(samplename, grouping=False))
                    if len(read_pairs) > 1:
                        join_lanes[setname].append(samplename)
        else:  # treat each file as sample
            print("Treat files as single end")
            sample_dict = OrderedDict(
                    [[get_setname(x, grouping=False), [os.path.abspath(x)]]
                     for x in all_files])
            name_dict = dict(
                    [get_basename(x), get_setname(x, grouping=False)]
                    for x in all_files)
            type = "SE"

    return type, sample_dict, join_lanes, name_dict, join_reads


def check_fasta(fastafile, filetype='File'):
    seq_record = SeqIO.parse(fastafile, "fasta")
    try:
        for seq in seq_record:
            break # break after first sequence
    except FileNotFoundError as err:
        print("%s %s does not exist" % (filetype, fastafile), file=stderr)
        return False
    try:
        s=str(seq) # we somehow need to touch the sequence element to trigger an error or not
        return True
    except:
        print('%s %s does not contain valid fasta data' % (filetype, fastafile), file=stderr)
        return False


def check_input_validity(arguments):
    if arguments["reference"]:
        ref_file = arguments["reference"]
        seq_record = ""
        if not check_fasta(arguments["reference"],"Reference"):
            sys.exit(1)
        try:
            if ref_file[-2:] == "gz":
                with subprocess.Popen(["gzip", "-cd", ref_file],
                                      stdout=subprocess.PIPE) as gz_proc:
                    seq_record = gz_proc.stdout.readline().decode()
                    gz_proc.terminate()
            else:
                with open(arguments["reference"], "r") as ref_fh:
                    seq_record = ref_fh.readline()
            assert seq_record.startswith(">"), ("Error: Reference"
                                                "file is not valid.")
        except AssertionError as e:
            sys.exit(e)
        except FileNotFoundError:
            sys.exit("Error: Reference file not found")

    #
    # Check validity of Kraken DB
    if not arguments["nokraken"]:
        if not os.path.exists(arguments["kraken_db"]):
            sys.exit("ERROR: %s does not exist.i"
                     " Enter a valid database"
                     " for kraken" % arguments["kraken_db"])
        else:
            if "database.kdb" not in os.listdir(arguments["kraken_db"]):
                sys.exit("ERROR: database "
                         + arguments["kraken_db"]
                         + " does not contain necessary file database.kdb")
    #
    # Check input
    if arguments["input"]:
        if not all([os.path.exists(x[0]) for x in arguments["input"]]):
            sys.exit(str(arguments["input"]) + " does not exist.")
    else:
        if not arguments["r1"]:
            sys.exit("Pleaser enter an input file (--input or -1/-2)")
        if not os.path.isfile(arguments["r1"]):
            sys.exit(arguments["r1"]
                     + " does not exist. Input file required."
                     + " Use option -input or -1 / -2.")

        if arguments["r2"]:
            if not os.path.isfile(arguments["r2"]):
                sys.exit(arguments["r2"] + " does not exist.")

    if arguments["trimBetter_threshold"] and not arguments["trimBetter"]:
        sys.exit("--trimBetter must be set to use --trimbetter_threshold."
                 " Add --trimBetter to your command "
                 "or remove --trimbetter_threshold.")

# --------------------------------------< main >-------------------------------


if __name__ == "__main__":
    main()