megahit

#!/usr/bin/env python

#-------------------------------------------------------------------------
# MEGAHIT
# Copyright (C) 2014 - 2015 The University of Hong Kong & L3 Bioinformatics Limited
#
# 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/>.
#-------------------------------------------------------------------------

# contact: Dinghua Li <dhli@cs.hku.hk>

from __future__ import print_function

import sys
import getopt
import subprocess
import errno
import os, glob
import shutil
import locale
import signal
import multiprocessing
import logging
import time
import tempfile
from datetime import datetime

megahit_version_str = ""
max_k_allowed = 127
usage_message = '''
Copyright (c) The University of Hong Kong & L3 Bioinformatics Limited
contact: Dinghua Li <dhli@cs.hku.hk>

Usage:
  megahit [options] {-1 <pe1> -2 <pe2> | --12 <pe12> | -r <se>} [-o <out_dir>]

  Input options that can be specified for multiple times (supporting plain text and gz/bz2 extensions)
    -1                       <pe1>          comma-separated list of fasta/q paired-end #1 files, paired with files in <pe2>
    -2                       <pe2>          comma-separated list of fasta/q paired-end #2 files, paired with files in <pe1>
    --12                     <pe12>         comma-separated list of interleaved fasta/q paired-end files
    -r/--read                <se>           comma-separated list of fasta/q single-end files

  Input options that can be specified for at most ONE time (not recommended):
    --input-cmd              <cmd>          command that outputs fasta/q reads to stdout; taken by MEGAHIT as SE reads 

Optional Arguments:
  Basic assembly options:
    --min-count              <int>          minimum multiplicity for filtering (k_min+1)-mers [2]
    --k-list                 <int,int,..>   comma-separated list of kmer size
                                            all must be odd, in the range 15-%d, increment <= 28);
                                            [21,29,39,59,79,99,119,141]

  Another way to set --k-list (overrides --k-list if one of them set):
    --k-min                  <int>          minimum kmer size (<= %d), must be odd number [21]
    --k-max                  <int>          maximum kmer size (<= %d), must be odd number [141]
    --k-step                 <int>          increment of kmer size of each iteration (<= 28), must be even number [12]

  Advanced assembly options:
    --no-mercy                              do not add mercy kmers
    --bubble-level           <int>          intensity of bubble merging (0-2), 0 to disable [2]
    --merge-level            <l,s>          merge complex bubbles of length <= l*kmer_size and similarity >= s [20,0.95]
    --prune-level            <int>          strength of low depth pruning (0-3) [2]
    --prune-depth            <int>          remove unitigs with avg kmer depth less than this value [2]
    --low-local-ratio        <float>        ratio threshold to define low local coverage contigs [0.2]
    --max-tip-len            <int>          remove tips less than this value [2*k]
    --no-local                              disable local assembly
    --kmin-1pass                            use 1pass mode to build SdBG of k_min

  Presets parameters:
    --presets                <str>          override a group of parameters; possible values:
                                            meta-sensitive: '--min-count 1 --k-list 21,29,39,49,...,129,141'
                                            meta-large: '--k-min 27 --k-max 127 --k-step 10'
                                            (large & complex metagenomes, like soil)

  Hardware options:
    -m/--memory              <float>        max memory in byte to be used in SdBG construction
                                            (if set between 0-1, fraction of the machine's total memory) [0.9]
    --mem-flag               <int>          SdBG builder memory mode
                                            0: minimum; 1: moderate; others: use all memory specified by '-m/--memory' [1]
    --use-gpu                               use GPU
    --gpu-mem                <float>        GPU memory in byte to be used. Default: auto detect to use up all free GPU memory. 
    -t/--num-cpu-threads     <int>          number of CPU threads, at least 2 if GPU enabled. [# of logical processors]

  Output options:
    -o/--out-dir             <string>       output directory [./megahit_out]
    --out-prefix             <string>       output prefix (the contig file will be OUT_DIR/OUT_PREFIX.contigs.fa)
    --min-contig-len         <int>          minimum length of contigs to output [200]
    --keep-tmp-files                        keep all temporary files
    --tmp-dir                <string>       set temp directory

Other Arguments:
    --continue                              continue a MEGAHIT run from its last available check point.
                                            please set the output directory correctly when using this option.
    -h/--help                               print the usage message
    -v/--version                            print version
    --verbose                               verbose mode
'''

if sys.version > '3':
    long = int

class Usage(Exception):
    def __init__(self, msg):
        self.msg = msg

class Options():
    def __init__(self):
        self.host_mem = 0.9
        self.gpu_mem = 0
        self.out_dir = "./megahit_out/"
        self.min_contig_len = 200
        self.k_min = 21
        self.k_max = 141
        self.k_step = 10
        self.k_list = [21,29,39,59,79,99,119,141]
        self.auto_k = True
        self.set_list_by_min_max_step = False
        self.min_count = 2
        self.bin_dir = sys.path[0] + "/"
        self.max_tip_len = -1
        self.no_mercy = False
        self.no_local = False
        self.bubble_level = 2
        self.merge_len = 20
        self.merge_similar = 0.95
        self.prune_level = 2
        self.prune_depth = 2
        self.num_cpu_threads = 0
        self.low_local_ratio = 0.2
        self.temp_dir = ""
        self.contig_dir = self.out_dir + "intermediate_contigs/"
        self.keep_tmp_files = False
        self.builder = "megahit_sdbg_build"
        self.use_gpu = False
        self.mem_flag = 1
        self.continue_mode = False;
        self.last_cp = -1;
        self.out_prefix = ""
        self.kmin_1pass = False
        self.verbose = False
        self.pe1 = []
        self.pe2 = []
        self.pe12 = []
        self.se = []
        self.input_cmd = ""
        self.inpipe = dict()
        self.presets = ""
        self.graph_only = False
        self.force_overwrite = False

opt = Options()
cp = 0

def log_file_name():
    if opt.out_prefix == "":
        return opt.out_dir + "log"
    else:
        return opt.out_dir + opt.out_prefix + ".log"

def opt_file_name():
    return opt.out_dir + "opts.txt"

def make_out_dir():
    if not os.path.exists(opt.out_dir):
        os.mkdir(opt.out_dir)

    if not opt.continue_mode:
        if opt.temp_dir == opt.out_dir + "tmp/":
            if not os.path.exists(opt.temp_dir):
                os.mkdir(opt.temp_dir)
        else:
            opt.temp_dir = tempfile.mkdtemp(dir = opt.temp_dir, prefix = 'megahit_tmp_') + "/"

    if not os.path.exists(opt.contig_dir):
        os.mkdir(opt.contig_dir)

def parse_opt(argv):
    try:
        opts, args = getopt.getopt(argv, "hm:o:r:t:v1:2:l:f", 
                                    ["help",
                                     "read=",
                                     "12=",
                                     "input-cmd=",
                                     "memory=",
                                     "out-dir=",
                                     "min-contig-len=",
                                     "use-gpu",
                                     "num-cpu-threads=",
                                     "gpu-mem=",
                                     "kmin-1pass",
                                     "k-min=",
                                     "k-max=",
                                     "k-step=",
                                     "k-list=",
                                     "num-cpu-threads=",
                                     "min-count=",
                                     "no-mercy",
                                     "no-local",
                                     "max-tip-len=",
                                     "bubble-level=",
                                     "prune-level=",
                                     "prune-depth=",
                                     "merge-level=",
                                     "low-local-ratio=",
                                     "keep-tmp-files",
                                     "tmp-dir=",
                                     "mem-flag=",
                                     "continue",
                                     "version",
                                     "out-prefix=",
                                     "verbose",
                                     "presets=",
                                     "force"
                                     "graph-only",
                                     "max-read-len=",
                                     "no-low-local",
                                     "cpu-only"])
    except getopt.error as msg:
        raise Usage(megahit_version_str + '\n' + str(msg))
    if len(opts) == 0:
        raise Usage(megahit_version_str + '\n' + usage_message)

    global opt
    need_continue = False

    for option, value in opts:
        if option in ("-h", "--help"):
            print(megahit_version_str + '\n' + usage_message)
            exit(0)
        elif option in ("-o", "--out-dir"):
            if opt.continue_mode == 0:
                opt.out_dir = value + "/"
        elif option in ("-m", "--memory"):
            opt.host_mem = float(value)
        elif option == "--gpu-mem":
            opt.gpu_mem = long(float(value))
        elif option == "--min-contig-len":
            opt.min_contig_len = int(value)
        elif option in ("-t", "--num-cpu-threads"):
            opt.num_cpu_threads = int(value)
        elif option == "--kmin-1pass":
            opt.kmin_1pass = True
        elif option == "--k-min":
            opt.k_min = int(value)
            opt.set_list_by_min_max_step = True
            opt.auto_k = False
        elif option == "--k-max":
            opt.k_max = int(value)
            opt.set_list_by_min_max_step = True
            opt.auto_k = False
        elif option == "--k-step":
            opt.k_step = int(value)
            opt.set_list_by_min_max_step = True
            opt.auto_k = False
        elif option == "--k-list":
            opt.k_list = list(map(int, value.split(",")))
            opt.k_list.sort()
            opt.auto_k = False
            opt.set_list_by_min_max_step = False
        elif option == "--min-count":
            opt.min_count = int(value)
        elif option == "--max-tip-len":
            opt.max_tip_len = int(value)
        elif option == "--merge-level":
            (opt.merge_len, opt.merge_similar) = map(float, value.split(","))
            opt.merge_len = int(opt.merge_len)
        elif option == "--prune-level":
            opt.prune_level = int(value)
        elif option == "--prune-depth":
            opt.prune_depth = float(value)
        elif option == "--bubble-level":
            opt.bubble_level = int(value)
        elif option == "--no-mercy":
            opt.no_mercy = True
        elif option == "--no-local":
            opt.no_local = True
        elif option == "--low-local-ratio":
            opt.low_local_ratio = float(value)
        elif option == "--keep-tmp-files":
            opt.keep_tmp_files = True
        elif option == "--use-gpu":
            opt.use_gpu = True
            opt.builder = "megahit_sdbg_build_gpu"
        elif option == "--mem-flag":
            opt.mem_flag = int(value)
        elif option in ("-v", "--version"):
            print(megahit_version_str)
            exit(0)
        elif option == "--verbose":
            opt.verbose = True
        elif option == "--continue":
            if opt.continue_mode == 0: # avoid check again again again...
                need_continue = True
        elif option == "--out-prefix":
            opt.out_prefix = value
        elif option == "--tmp-dir":
            if opt.continue_mode == 0:
                opt.temp_dir = value + "/"
        elif option in ("--cpu-only", "-l", "--max-read-len", "--no-low-local"):
            continue # historical options, just ignore
        elif option in ("-r", "--read"):
            opt.se += value.split(",")
        elif option == "-1":
            opt.pe1 += value.split(",")
        elif option == "-2":
            opt.pe2 += value.split(",")
        elif option == "--12":
            opt.pe12 += value.split(",")
        elif option == "--input-cmd":
            opt.input_cmd = value
        elif option == "--presets":
            opt.presets = value
        elif option == "--graph-only":
            opt.graph_only = True;
        elif option in ("-f", "--force"):
            opt.force_overwrite = True

        else:
            raise Usage("Invalid option %s", option)

    if opt.temp_dir == "":
        opt.temp_dir = opt.out_dir + "tmp/"
    opt.contig_dir = opt.out_dir + "intermediate_contigs/"

    if need_continue:
        prepare_continue()
    elif opt.continue_mode == 0 and not opt.force_overwrite and os.path.exists(opt.out_dir):
        raise Usage("Output directory " + opt.out_dir + " already exists, please change the parameter -o to another value to avoid overwriting.")

def check_opt():
    global opt
    if opt.host_mem <= 0:
        raise Usage("Please specify a positive number for -m flag.")
    elif opt.host_mem < 1:
        total_mem = detect_available_mem()
        opt.host_mem = long(total_mem * opt.host_mem)
        if total_mem <= 0:
            raise Usage("Failed to detect available memory. Please specify the value in bytes using -m flag.")
        else:
            print(str(round(total_mem/(1024**3),3)) + "Gb memory in total.", file=sys.stderr)
            print("Using: " + str(round(float(opt.host_mem)/(1024**3),3)) + "Gb.", file=sys.stderr)
    else:
        opt.host_mem = long(opt.host_mem)

    # set mode
    if opt.presets != "":
        opt.auto_k = True;

        if opt.presets == "meta-sensitive":
            opt.min_count = 1
            opt.k_list = [21,29,39,49,59,69,79,89,99,109,119,129,141]
            opt.set_list_by_min_max_step = False
        elif opt.presets == "meta-large":
            opt.min_count = 1
            opt.k_min = 27
            opt.k_max = 127
            opt.k_step = 10
            opt.set_list_by_min_max_step = True
        else:
            raise Usage("Invalid preset: " + opt.presets)

    if opt.set_list_by_min_max_step:
        if opt.k_step % 2 == 1:
            raise Usage("k-step must be even number!")
        if opt.k_min > opt.k_max:
            raise Usage("Error: k_min > k_max!")

        opt.k_list = list()
        k = opt.k_min
        while k < opt.k_max:
            opt.k_list.append(k)
            k = k + opt.k_step
        opt.k_list.append(opt.k_max)

    if len(opt.k_list) == 0:
        raise Usage("k list should not be empty!")

    if opt.k_list[0] < 15 or opt.k_list[-1] > max_k_allowed:
        raise Usage("All k's should be in range [15, %d]" % max_k_allowed)

    for k in opt.k_list:
        if k % 2 == 0:
            raise Usage("All k must be odd number!")

    for i in range(1, len(opt.k_list)):
        if opt.k_list[i] - opt.k_list[i-1] > 28:
            raise Usage("k-step/adjacent k difference must be <= 28")

    opt.k_min, opt.k_max = opt.k_list[0], opt.k_list[-1]

    if opt.use_gpu == 0:
        opt.gpu_mem = 0
    if opt.k_max < opt.k_min:
        raise Usage("k_min should be no larger than k_max.")
    if opt.min_count <= 0:
        raise Usage("min_count must be greater than 0.")
    elif opt.min_count == 1:
        opt.kmin_1pass = True
        opt.no_mercy = True
    if opt.prune_level < 0 or opt.prune_level > 3:
        raise Usage("prune level must be in 0-3.")
    if opt.merge_len < 0:
        raise Usage("merge_level: length must be >= 0")
    if opt.merge_similar < 0 or opt.merge_similar > 1:
        raise Usage("merge_level: similarity must be in [0, 1]")
    if opt.low_local_ratio <= 0 or opt.low_local_ratio > 0.5:
        raise Usage("low_local_ratio should be in (0, 0.5].")
    if opt.num_cpu_threads > multiprocessing.cpu_count():
        print("Maximum number of available CPU thread is %d." % multiprocessing.cpu_count(), file=sys.stderr);
        print("Number of thread is reset to the %d." % multiprocessing.cpu_count(), file=sys.stderr);
        opt.num_cpu_threads = multiprocessing.cpu_count()
    if opt.num_cpu_threads == 1 and opt.use_gpu:
        opt.num_cpu_threads = 2
        print("GPU enabled. Number of thread is reset to the 2.", file=sys.stderr);
    if opt.num_cpu_threads == 0:
        opt.num_cpu_threads = multiprocessing.cpu_count()
    if opt.prune_depth < 0 and opt.prune_level < 3:
        opt.prune_depth = opt.min_count
    if opt.bubble_level < 0:
        print("Reset bubble level to 0.", file=sys.stderr)
        opt.bubble_level = 0
    if opt.bubble_level > 2:
        print("Reset bubble level to 2.", file=sys.stderr)
        opt.bubble_level = 2

def check_reads():
    # reads
    global opt
    global cp
    if (not opt.continue_mode) or (cp > opt.last_cp):
        if len(opt.pe1) != len(opt.pe2):
            raise Usage("Number of paired-end files not match!")
        for r in opt.pe1 + opt.pe2 + opt.se + opt.pe12:
            if not os.path.exists(r):
                raise Usage("Cannot find file " + r)

        if opt.input_cmd == "" and len(opt.pe1 + opt.pe2 + opt.se + opt.pe12) == 0:
            raise Usage("No input files or input command!")

    write_cp()


def detect_available_mem():
    mem = long()
    if sys.platform.find("linux") != -1:
        try:
            mem = long(float(os.popen("free").readlines()[1].split()[1]) * 1024)
        except IndexError:
            mem = 0
    elif sys.platform.find("darwin") != -1:
        try:
            mem = long(float(os.popen("sysctl hw.memsize").readlines()[0].split()[1]))
        except IndexError:
            mem = 0
    else:
        mem = 0
    return mem

def write_opt(argv):
    with open(opt_file_name(), "w") as f:
        print("\n".join(argv), file=f)
        print("MEGAHIT_TEMP_DIR: " + opt.temp_dir, file=f)
    f.close()

def prepare_continue():
    global opt # out_dir is already set
    if not os.path.exists(opt_file_name()):
        print("Cannot find " + opt.out_dir + "opts.txt", file=sys.stderr)
        print("Please check whether the output directory is correctly set by \"-o\"", file=sys.stderr)
        print("Now switching to normal mode.", file=sys.stderr)
        return

    print("Continue mode activated. Ignore all options other than -o/--out-dir.", file=sys.stderr)

    with open(opt_file_name(), "r") as f:
        argv = []
        for line in f:
            if line.startswith("MEGAHIT_TEMP_DIR: "):
                temp_dir = line.strip()[18:]
            else:
                argv.append(line.strip())
        print("Continue with options: " + " ".join(argv), file=sys.stderr)
        t_dir = opt.out_dir
        opt = Options()
        opt.out_dir = t_dir
        opt.temp_dir = temp_dir;
        opt.continue_mode = True # avoid dead loop
        parse_opt(argv)
    f.close()

    opt.last_cp = -1
    if os.path.exists(opt.temp_dir + "cp.txt"):
        with open(opt.temp_dir + "cp.txt", "r") as cpf:
            for line in cpf:
                a = line.strip().split()
                if len(a) == 2 and a[1] == "done":
                    opt.last_cp = int(a[0])
        cpf.close()
    print("Continue from check point " + str(opt.last_cp), file=sys.stderr)

def check_bin():
    for subprogram in ["megahit_asm_core", "megahit_toolkit"]:
        if not os.path.exists(opt.bin_dir + subprogram):
            raise Usage("Cannot find sub-program \"" + subprogram + "\", please recompile.")

def get_version():
    global megahit_version_str
    global usage_message
    global max_k_allowed
    megahit_version_str = "MEGAHIT " + \
                          subprocess.Popen([opt.bin_dir + "megahit_asm_core", "dumpversion"],
                                           stdout=subprocess.PIPE).communicate()[0].rstrip().decode('utf-8')

    max_k_allowed = int(subprocess.Popen([opt.bin_dir + "megahit_asm_core", "kmax"],
                                           stdout=subprocess.PIPE).communicate()[0].rstrip().decode('utf-8'))

    usage_message = usage_message % (max_k_allowed, max_k_allowed, max_k_allowed)


def check_builder():
    if not os.path.exists(opt.bin_dir + opt.builder):
        usg = megahit_version_str + '\n' + "Cannot find sub-program \"%s\", please recompile." % opt.builder
        if opt.use_gpu == 0:
            usg += "\nOr if you want to use the GPU version, please run MEGAHIT with \"--use-gpu\""
        raise Usage(usg)

def graph_prefix(kmer_k):
    if not os.path.exists(opt.temp_dir + "k" + str(kmer_k)):
        os.mkdir(opt.temp_dir + "k" + str(kmer_k))
    return opt.temp_dir + "k" + str(kmer_k) + "/" + str(kmer_k)

def contig_prefix(kmer_k):
    return opt.contig_dir + "k" + str(kmer_k)

def delect_file_if_exist(file_name):
    if os.path.exists(file_name):
        os.remove(file_name)

def delete_tmp_after_build(kmer_k):
    for i in range(0, opt.num_cpu_threads):
        delect_file_if_exist(graph_prefix(kmer_k) + ".edges." + str(i))
    for i in range(0, 64):
        delect_file_if_exist(graph_prefix(kmer_k) + ".mercy_cand." + str(i))
    for i in range(0, opt.num_cpu_threads - 1):
        delect_file_if_exist(graph_prefix(kmer_k) + ".mercy." + str(i))
    delect_file_if_exist(graph_prefix(kmer_k) + ".cand")

def delete_tmp_after_assemble(kmer_k):
    for extension in ["w", "last", "isd", "dn", "f", "mul", "mul2"]:
        delect_file_if_exist(graph_prefix(kmer_k) + "." + extension)
    for i in range(0, opt.num_cpu_threads):
        delect_file_if_exist(graph_prefix(kmer_k) + ".sdbg." + str(i))

def delete_tmp_after_iter(kmer_k):
    delect_file_if_exist(graph_prefix(kmer_k) + ".rr.pb")

def write_cp():
    global cp
    cpf = open(opt.temp_dir + "cp.txt", "a")
    print(str(cp) + "\t" + "done", file=cpf);
    cp = cp + 1
    cpf.close()

def inpipe_cmd(file_name):
    if file_name.endswith('.gz'):
        return "gzip -cd '" + file_name + "'"
    elif file_name.endswith('.bz2'):
        return "bzip2 -cd '" + file_name + "'"
    else:
        return "cat '" +  file_name + "'"

def write_lib():
    global opt
    opt.lib = opt.temp_dir + "reads.lib"
    lib = open(opt.lib, "w")
    for i in range(0, len(opt.pe12)):
        print(opt.pe12[i], file=lib)

        if inpipe_cmd(opt.pe12[i]) != "":
            print("interleaved " + opt.temp_dir + "inpipe.pe12." + str(i), file=lib)
        else:
            print("interleaved " + opt.pe12[i], file=lib)

    for i in range(0, len(opt.pe1)):
        if inpipe_cmd(opt.pe1[i]) != "":
            f1 = opt.temp_dir + "inpipe.pe1." + str(i)
        else:
            f1 = opt.pe1[i]

        if inpipe_cmd(opt.pe2[i]) != "":
            f2 = opt.temp_dir + "inpipe.pe2." + str(i)
        else:
            f2 = opt.pe2[i]

        print(','.join([opt.pe1[i], opt.pe2[i]]), file=lib)
        print("pe " + f1 + " " + f2, file=lib)

    for i in range(0, len(opt.se)):
        print(opt.se[i], file=lib)

        if inpipe_cmd(opt.se[i]) != "":
            print("se " + opt.temp_dir + "inpipe.se." + str(i), file=lib)
        else:
            print("se " + opt.se[i], file=lib)

    if opt.input_cmd != "":
        print('\"' + opt.input_cmd + '\"', file=lib)
        print("se " + "-", file=lib)

    lib.close()

def build_lib():
    global cp
    if (not opt.continue_mode) or (cp > opt.last_cp):
        build_lib_cmd = [opt.bin_dir + "megahit_asm_core", "buildlib",
                         opt.lib,
                         opt.lib]

        fifos = list()
        pipes = list()
        try:
            # create inpipe

            for i in range(0, len(opt.pe12)):
                if inpipe_cmd(opt.pe12[i]) != "":
                    delect_file_if_exist(opt.temp_dir + "inpipe.pe12." + str(i))
                    os.mkfifo(opt.temp_dir + "inpipe.pe12." + str(i))
                    fifos.append(opt.temp_dir + "inpipe.pe12." + str(i))

            for i in range(0, len(opt.pe1)):
                if inpipe_cmd(opt.pe1[i]) != "":
                    delect_file_if_exist(opt.temp_dir + "inpipe.pe1." + str(i))
                    os.mkfifo(opt.temp_dir + "inpipe.pe1." + str(i))
                    fifos.append(opt.temp_dir + "inpipe.pe1." + str(i))
                
                if inpipe_cmd(opt.pe2[i]) != "":
                    delect_file_if_exist(opt.temp_dir + "inpipe.pe2." + str(i))
                    os.mkfifo(opt.temp_dir + "inpipe.pe2." + str(i))
                    fifos.append(opt.temp_dir + "inpipe.pe2." + str(i))

            for i in range(0, len(opt.se)):
                if inpipe_cmd(opt.se[i]) != "":
                    delect_file_if_exist(opt.temp_dir + "inpipe.se." + str(i))
                    os.mkfifo(opt.temp_dir + "inpipe.se." + str(i))
                    fifos.append(opt.temp_dir + "inpipe.se." + str(i))

            logging.info("--- [%s] Converting reads to binaries ---" % datetime.now().strftime("%c"))
            logging.debug("%s" % (" ").join(build_lib_cmd))

            if opt.input_cmd != "":
                logging.debug("input cmd: " + opt.input_cmd)
                input_thread = subprocess.Popen(opt.input_cmd, shell = True, stdout = subprocess.PIPE)
                p = subprocess.Popen(build_lib_cmd, stdin = input_thread.stdout, stdout = subprocess.PIPE, stderr=subprocess.PIPE)
            else:
                p = subprocess.Popen(build_lib_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)

            # output to inpipe

            for i in range(0, len(opt.pe12)):
                if inpipe_cmd(opt.pe12[i]) != "":
                    ip_thread12 = subprocess.Popen(inpipe_cmd(opt.pe12[i]) + " > " + opt.temp_dir + "inpipe.pe12." + str(i), shell = True, preexec_fn = os.setsid)
                    pipes.append(ip_thread12)

            for i in range(0, len(opt.pe1)):
                if inpipe_cmd(opt.pe1[i]) != "":
                    ip_thread1 = subprocess.Popen(inpipe_cmd(opt.pe1[i]) + " > " + opt.temp_dir + "inpipe.pe1." + str(i), shell = True, preexec_fn = os.setsid)
                    pipes.append(ip_thread1)
                
                if inpipe_cmd(opt.pe2[i]) != "":
                    ip_thread2 = subprocess.Popen(inpipe_cmd(opt.pe2[i]) + " > " + opt.temp_dir + "inpipe.pe2." + str(i), shell = True, preexec_fn = os.setsid)
                    pipes.append(ip_thread2)

            for i in range(0, len(opt.se)):
                if inpipe_cmd(opt.se[i]) != "":
                    ip_thread_se = subprocess.Popen(inpipe_cmd(opt.se[i]) + " > " + opt.temp_dir + "inpipe.se." + str(i), shell = True, preexec_fn = os.setsid)
                    pipes.append(ip_thread_se)
            
            while True:
                line = p.stderr.readline().rstrip()
                if not line:
                    break;
                logging.info(line)

            ret_code = p.wait()

            if ret_code != 0:
                logging.error("Error occurs when running \"megahit_asm_core buildlib\"; please refer to %s for detail" % log_file_name())
                logging.error("[Exit code %d]" % ret_code)
                exit(ret_code)

            while len(pipes) > 0:
                pp = pipes.pop()
                pp_ret = pp.wait()
                if pp_ret != 0:
                    logging.error("Error occurs when reading inputs")
                    exit(pp_ret)

        except OSError as o:
            if o.errno == errno.ENOTDIR or o.errno == errno.ENOENT:
                logging.error("Error: sub-program megahit_asm_core not found, please recompile MEGAHIT")
            exit(1)
        except KeyboardInterrupt:
            p.terminate()
            exit(1)

        finally:
            for p in pipes:
                os.killpg(p.pid, signal.SIGTERM)
            for f in fifos:
                delect_file_if_exist(f)

    write_cp()

def get_max_read_len():
    ret = 0
    with open(opt.lib + ".lib_info") as read_info:
        for info in read_info.readlines()[2::2]:
            ret = max(ret, int(info.split()[2]))
    read_info.close()
    return ret


def set_max_k_by_lib():
    if not opt.auto_k or len(opt.k_list) == 1:
        return False

    max_read_len = get_max_read_len()
    for i in range(0, len(opt.k_list)):
        if opt.k_list[i] >= max_read_len + 20:
            opt.k_max = max_read_len + 19 - max_read_len % 2
            if i > 0 and opt.k_list[i-1] >= opt.k_max - 2:
                opt.k_list = opt.k_list[:i]
                opt.k_max = opt.k_list[-1]
            else:
                opt.k_list = opt.k_list[:(i+1)]
                opt.k_list[-1] = opt.k_max
            return True


def build_first_graph():
    global cp
    phase1_out_threads = max(1, int(opt.num_cpu_threads / 3))
    if (not opt.continue_mode) or (cp > opt.last_cp):
        count_opt = ["-k", str(opt.k_min),
                     "-m", str(opt.min_count),
                     "--host_mem", str(opt.host_mem),
                     "--mem_flag", str(opt.mem_flag),
                     "--gpu_mem", str(opt.gpu_mem),
                     "--output_prefix", graph_prefix(opt.k_min),
                     "--num_cpu_threads", str(opt.num_cpu_threads),
                     "--num_output_threads", str(phase1_out_threads),
                     "--read_lib_file", opt.lib]

        cmd = []
        if opt.kmin_1pass:
            cmd = [opt.bin_dir + opt.builder, "read2sdbg"] + count_opt
            if not opt.no_mercy:
                cmd.append("--need_mercy")
        else:
            cmd = [opt.bin_dir + opt.builder, "count"] + count_opt

        try:
            if opt.kmin_1pass:
                logging.info("--- [%s] Extracting solid (k+1)-mers and building sdbg for k = %d ---" % (datetime.now().strftime("%c"), opt.k_min))
            else:
                logging.info("--- [%s] Extracting solid (k+1)-mers for k = %d ---" % (datetime.now().strftime("%c"), opt.k_min))

            logging.debug("cmd: %s" % (" ").join(cmd))
            p = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)


            while True:
                line = p.stderr.readline().rstrip()
                if not line:
                    break;
                logging.debug(line)

            ret_code = p.wait()

            if ret_code != 0:
                logging.error("Error occurs when running \"sdbg_builder count/read2sdbg\", please refer to %s for detail" % log_file_name())
                logging.error("[Exit code %d] " % ret_code)
                exit(ret_code)

        except OSError as o:
            if o.errno == errno.ENOTDIR or o.errno == errno.ENOENT:
                logging.error("Error: sub-program sdbg_builder not found, please recompile MEGAHIT")
            exit(1)
        except KeyboardInterrupt:
            p.terminate()
            exit(1)

    write_cp()
    if not opt.kmin_1pass:
        build_graph(opt.k_min, 0, phase1_out_threads)
    elif not opt.keep_tmp_files:
        delete_tmp_after_build(opt.k_min)

def build_graph(kmer_k, kmer_from, num_edge_files):
    global cp
    if (not opt.continue_mode) or (cp > opt.last_cp):
        build_comm_opt = ["--host_mem", str(opt.host_mem),
                             "--mem_flag", str(opt.mem_flag),
                             "--gpu_mem", str(opt.gpu_mem),
                             "--output_prefix", graph_prefix(kmer_k),
                             "--num_cpu_threads", str(opt.num_cpu_threads),
                             "-k", str(kmer_k), 
                             "--kmer_from", str(kmer_from),
                             "--num_edge_files", str(num_edge_files)]

        build_cmd = [opt.bin_dir + opt.builder, "seq2sdbg"] + build_comm_opt

        file_size = 0

        if (os.path.exists(graph_prefix(kmer_k) + ".edges.0")):
            build_cmd += ["--input_prefix", graph_prefix(kmer_k)]
            file_size += os.path.getsize(graph_prefix(kmer_k) + ".edges.0")
            tid = 1
            while (os.path.exists(graph_prefix(kmer_k) + ".edges." + str(tid))):
                file_size += os.path.getsize(graph_prefix(kmer_k) + ".edges." + str(tid))
                tid += 1

        if (os.path.exists(contig_prefix(kmer_from) + ".addi.fa")):
            build_cmd += ["--addi_contig", contig_prefix(kmer_from) + ".addi.fa"]
            file_size += os.path.getsize(contig_prefix(kmer_from) + ".addi.fa")

        if (os.path.exists(contig_prefix(kmer_from) + ".local.fa")):
            build_cmd += ["--local_contig", contig_prefix(kmer_from) + ".local.fa"]
            file_size += os.path.getsize(contig_prefix(kmer_from) + ".local.fa")

        if (os.path.exists(contig_prefix(kmer_from) + ".contigs.fa")):
            build_cmd += ["--contig", contig_prefix(kmer_from) + ".contigs.fa"]
            build_cmd += ["--bubble", contig_prefix(kmer_from) + ".bubble_seq.fa"]

        if file_size == 0:
            return False # not build

        if not opt.no_mercy and kmer_k == opt.k_min:
            build_cmd.append("--need_mercy")

        try:
            logging.info("--- [%s] Building graph for k = %d ---" % (datetime.now().strftime("%c"), kmer_k))
            logging.debug("%s" % (" ").join(build_cmd))

            p = subprocess.Popen(build_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
            
            while True:
                line = p.stderr.readline().rstrip()
                if not line:
                    break;
                logging.debug(line)

            ret_code = p.wait()

            if ret_code != 0:
                logging.error("Error occurs when running \"builder build\" for k = %d; please refer to %s for detail" % (kmer_k, log_file_name()))
                logging.error("[Exit code %d]" % ret_code)
                exit(ret_code)

        except OSError as o:
            if o.errno == errno.ENOTDIR or o.errno == errno.ENOENT:
                logging.error("Error: sub-program builder not found, please recompile MEGAHIT")
            exit(1)
        except KeyboardInterrupt:
            p.terminate()
            exit(1)

    write_cp()
    if not opt.keep_tmp_files:
        delete_tmp_after_build(kmer_k)
    return True

def iterate(cur_k, step):
    global cp
    if (not opt.continue_mode) or (cp > opt.last_cp):
        next_k = cur_k + step

        iterate_cmd = [opt.bin_dir + "megahit_asm_core", "iterate",
                       "-c", contig_prefix(cur_k) + ".contigs.fa",
                       "-b", contig_prefix(cur_k) + ".bubble_seq.fa",
                       "-t", str(opt.num_cpu_threads),
                       "-k", str(cur_k),
                       "-s", str(step),
                       "-o", graph_prefix(next_k),
                       "-r", opt.lib + ".bin",
                       "-f", "binary"]

        try:
            logging.info("--- [%s] Extracting iterative edges from k = %d to %d ---" % (datetime.now().strftime("%c"), cur_k, next_k))
            logging.debug("cmd: %s" % (" ").join(iterate_cmd))

            p = subprocess.Popen(iterate_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
            
            while True:
                line = p.stderr.readline().rstrip()
                if not line:
                    break;
                logging.debug(line)

            ret_code = p.wait()

            if ret_code != 0:
                logging.error("Error occurs when running iterator for k = %d to k = %d, please refer to %s for detail" % (cur_k, next_k, log_file_name()))
                logging.error("[Exit code %d]" % ret_code)
                exit(ret_code)

        except OSError as o:
            if o.errno == errno.ENOTDIR or o.errno == errno.ENOENT:
                logging.error("Error: sub-program megahit_iter not found, please recompile MEGAHIT")
            exit(1)

    write_cp()
    if not opt.keep_tmp_files:
        delete_tmp_after_iter(cur_k)

def assemble(cur_k):
    global cp
    if (not opt.continue_mode) or (cp > opt.last_cp):
        min_standalone = max(min(opt.k_max * 3 - 1, opt.min_contig_len * 1.5), opt.min_contig_len)
        if (opt.max_tip_len >= 0):
            min_standalone = max(opt.max_tip_len + opt.k_max - 1, opt.min_contig_len)

        assembly_cmd = [opt.bin_dir + "megahit_asm_core", "assemble",
                        "-s", graph_prefix(cur_k),
                        "-o", contig_prefix(cur_k),
                        "-t", str(opt.num_cpu_threads),
                        "--min_standalone", str(min_standalone),
                        "--prune_level", str(opt.prune_level),
                        "--merge_len", str(int(opt.merge_len)),
                        "--merge_similar", str(opt.merge_similar),
                        "--low_local_ratio", str(opt.low_local_ratio),
                        "--min_depth", str(opt.prune_depth),
                        "--bubble_level", str(opt.bubble_level)]

        if opt.max_tip_len == -1 and cur_k * 3 - 1 > opt.min_contig_len * 1.5:
            assembly_cmd += ["--max_tip_len", str(max(1, opt.min_contig_len * 1.5 + 1 - cur_k))]
        else:
            assembly_cmd += ["--max_tip_len", str(opt.max_tip_len)]

        # if cur_k < (opt.k_min + opt.k_max) / 2 and cur_k < 60:
        if cur_k < opt.k_max:
            assembly_cmd.append("--careful_bubble")

        if cur_k == opt.k_max:
            assembly_cmd.append("--is_final_round")

        if opt.no_local:
            assembly_cmd.append("--output_standalone")
            
        try:
            logging.info("--- [%s] Assembling contigs from SdBG for k = %d ---" % (datetime.now().strftime("%c"), cur_k))
            logging.debug("cmd: %s" % (" ").join(assembly_cmd))

            p = subprocess.Popen(assembly_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
            
            while True:
                line = p.stderr.readline().rstrip()
                if not line:
                    break;
                logging.debug(line)

            ret_code = p.wait()

            if ret_code != 0:
                logging.error("Error occurs when assembling contigs for k = %d, please refer to %s for detail" % (cur_k, log_file_name()))
                logging.error("[Exit code %d]" % ret_code)
                exit(ret_code)
            
        except OSError as o:
            if o.errno == errno.ENOTDIR or o.errno == errno.ENOENT:
                logging.error("Error: sub-program megahit_assemble not found, please recompile MEGAHIT")
            exit(1)
        except KeyboardInterrupt:
            p.terminate()
            exit(1)

    write_cp()
    if (not opt.keep_tmp_files) and (cur_k != opt.k_max):
        delete_tmp_after_assemble(cur_k)

def local_assemble(cur_k, kmer_to):
    global cp
    if (not opt.continue_mode) or (cp > opt.last_cp):
        la_cmd = [opt.bin_dir + "megahit_asm_core", "local",
                  "-c", contig_prefix(cur_k) + ".contigs.fa",
                  "-l", opt.lib,
                  "-t", str(opt.num_cpu_threads),
                  "-o", contig_prefix(cur_k) + ".local.fa",
                  # "--sparsity", "1",
                  "--kmax", str(kmer_to)]
        try:
            logging.info("--- [%s] Local assembling k = %d ---" % (datetime.now().strftime("%c"), cur_k))
            logging.debug("cmd: %s" % (" ").join(la_cmd))

            p = subprocess.Popen(la_cmd, stderr=subprocess.PIPE)
            
            while True:
                line = p.stderr.readline().rstrip()
                if not line:
                    break;
                logging.debug(line)

            ret_code = p.wait()

            if ret_code != 0:
                logging.error("Error occurs when running local assembly for k = %d to k = %d, please refer to %s for detail" % (cur_k, min(cur_k + opt.k_step, opt.k_max), log_file_name()))
                logging.error("[Exit code %d]" % ret_code)
                exit(ret_code)

        except OSError as o:
            if o.errno == errno.ENOTDIR or o.errno == errno.ENOENT:
                logging.error("Error: sub-program megahit_local_asm not found, please recompile MEGAHIT")
            exit(1)
        except KeyboardInterrupt:
            p.terminate()
            exit(1)

    write_cp()

def merge_final():
    global cp
    if (not opt.continue_mode) or (cp > opt.last_cp):
        logging.info("--- [%s] Merging to output final contigs ---" % (datetime.now().strftime("%c")))
        final_contig_name = opt.out_dir + "final.contigs.fa"
        if opt.out_prefix != "":
            final_contig_name = opt.out_dir + opt.out_prefix + ".contigs.fa"

        with open(final_contig_name, "w") as final_contigs:
            merge_cmd = "cat " + opt.contig_dir + "*.final.contigs.fa " + contig_prefix(opt.k_max) + ".contigs.fa | " + \
                        opt.bin_dir + "megahit_toolkit filterbylen " + str(opt.min_contig_len)
            p = subprocess.Popen(merge_cmd, shell = True, stdout = final_contigs, stderr = subprocess.PIPE)

            while True:
                line = p.stderr.readline().rstrip().decode('utf-8')
                if not line:
                    break;
                logging.info("--- [STAT] " + line)

            ret_code = p.wait()

        if ret_code != 0:
            logging.error("Error occurs when merging final contigs, please refer to %s for detail" % log_file_name())
            logging.error("[Exit code %d]" % ret_code)
            exit(ret_code)

    write_cp()

def main(argv = None):
    if argv is None:
        argv = sys.argv

    try:
        start_time = time.time()

        check_bin()
        get_version()
        parse_opt(argv[1:])
        check_opt()
        check_builder()
        make_out_dir()
        check_reads()

        logging.basicConfig(level = logging.NOTSET,
                            format = '%(message)s',
                            filename = log_file_name(),
                            filemode = 'a')

        console = logging.StreamHandler()
        console.setLevel(logging.INFO)
        if opt.verbose:
            console.setLevel(logging.NOTSET)

        formatter = logging.Formatter('%(message)s')
        console.setFormatter(formatter)
        logging.getLogger('').addHandler(console)

        logging.info(megahit_version_str)
        logging.info("--- [%s] Start assembly. Number of CPU threads %d ---" % (datetime.now().strftime("%c"), opt.num_cpu_threads))
        logging.info("--- [%s] Available memory: %d, used: %d" % (datetime.now().strftime("%c"), detect_available_mem(), opt.host_mem))

        if not opt.continue_mode:
            write_opt(argv[1:]) # for --continue

        write_lib()
        build_lib()

        if set_max_k_by_lib():
            logging.info("--- [%s] k-max reset to: %d ---" % (datetime.now().strftime("%c"), opt.k_max))

        logging.info("--- [%s] k list: %s ---" % (datetime.now().strftime("%c"), ','.join(map(str, opt.k_list))))

        build_first_graph()

        if not opt.graph_only:
            assemble(opt.k_min)

            cur_k = opt.k_min
            next_k_idx = 0

            while cur_k < opt.k_max:
                next_k_idx += 1
                next_k = opt.k_list[next_k_idx]
                k_step = next_k - cur_k

                if not opt.no_local:
                    local_assemble(cur_k, next_k)

                iterate(cur_k, k_step)

                if not build_graph(next_k, cur_k, 1):
                    opt.k_max = cur_k # for merging cur_k.contigs.fa to final
                    break

                assemble(next_k)
                cur_k = next_k
            # end while

            merge_final()

            if not opt.keep_tmp_files and os.path.exists(opt.temp_dir):
                shutil.rmtree(opt.temp_dir)

            open(opt.out_dir + "done", "w").close()

        logging.info("--- [%s] ALL DONE. Time elapsed: %f seconds ---" % (datetime.now().strftime("%c"), time.time() - start_time))

    except Usage as err:
        print(sys.argv[0].split("/")[-1] + ": " + str(err.msg), file=sys.stderr)
        return 2

if __name__ == "__main__":
    sys.exit(main())