def _bowtie(workflow, conf):
    for target in conf.sample_targets:
        bowtie = attach_back(workflow,
            ShellCommand(
                "{tool} -p {param[threads]} -S -m {param[max_align]} \
                {param[genome_index]} {input[fastq]} {output[sam]} 2> {output[bowtie_summary]}",
                input={"genome_dir": os.path.dirname(conf.get_path("lib", "genome_index")),
                       "fastq": target + ".fastq"},
                output={"sam": target + ".sam",
                        "bowtie_summary": target + "_bowtie_summary.txt", },
                tool="bowtie",
                param={"threads": 4,
                       "max_align": 1,
                       "genome_index": conf.get_path("lib", "genome_index")}))

        bowtie.update(param=conf.items("bowtie"))
    __sam2bam(workflow, conf)

    ## using bowtie standard error output
    attach_back(workflow,
        PythonCommand(stat_bowtie,
            input={"bowtie_summaries": [t + "_bowtie_summary.txt" for t in conf.sample_targets],
                   "db": ChiLinQC_db,
                   "template": rlang_template},
            output={"json": conf.json_prefix + "_bowtie.json",
                    "R": conf.prefix + "_bowtie.R",
                    "pdf": conf.prefix + "_bowtie.pdf"},
            param={"sams": [t + ".sam" for t in conf.sample_targets], }))

def sample_bam_stat(workflow, conf, tex):
    """ sample non chrm bam to 15M for NSC and PBC
    sample non chrm bam to 5M for spot
    """
    for i, target in enumerate(conf.treatment_targets):
        ## for PE, use name sorted in order to calculate PBC
        input_bam = target + "_name_sorted.bam" if conf.pe else target + "_final_nochrm.bam"
        attach_back(workflow, ShellCommand(
            "{tool} {input[namesorted]} {param[run_spp]} {output[bamstat]} {output[sppstat]}  {param[pe]} {output[pbc]}",
            tool = "eap_dnase_stats",
            input = {"namesorted": input_bam},
            output = {"bamstat": target + "_bam_stat.qc",  ## 15M
                      "sppstat": target + "_spp.qc",
                      "pbc": target + "_final_nochrm_15M_pbc.qc"},
            param = {"pe": "pe" if conf.pe else "se",
                     "run_spp": conf.get("tool", "spp")}))

        if not "macs" in conf.get("tool", "peak_calling"):

            attach_back(workflow, ShellCommand(
                "{tool} {input[bamwithoutchrm]} {param[genome]} {param[readsize]} {output[spot]} {param[hotspot_dir]} {param[hotspot_output]} {param[hotspot_tmp]} {param[spot_tmp]}",
                tool = "dac_spot", ## 5M
                input = {"bamwithoutchrm": target + "_final_nochrm.bam"},
                output = {"spot": target + "_spot_nochrm_5M.qc"},

                param = {"genome": conf.species,
                         "spot_tmp": conf.hotspot_reps_tmp_prefix[i] + "_final_nochrm.bam.5000000.spot.out",
                         "readsize": conf.readsize,
                         "hotspot_dir": conf.get("tool", "peak_calling"),
                         "hotspot_output": target + "_hotspot",
                         "hotspot_tmp": target + "_hotspot_tmp"}))

def fragment(workflow, conf):
    ## this is done after FRiP
    if conf.get("tool", "macs2"):
       macs2_bin = conf.get("tool", "macs2")
    else:
        macs2_bin = "macs2"
    for target in conf.treatment_targets:
        fragment_size = attach_back(workflow, ShellCommand(
            "{tool} predictd -i {input[bam]} --rfile {param[prefix]} -g {param[species]}",
            tool = macs2_bin,
            input = {"bam": target + ".bam"},
            output = {"R": target + "_model.R"},
            param = {"prefix": target + "_model.R",
                     "species": 'hs'}))
        fragment_size.update(param = conf.items("macs2"))
        ## except two few peaks for modeling
        fragment_size.allow_fail = True
        fragment_size.allow_dangling = True
        
    ## extract standard deviation from MACS2 model.R,
    ## use m, p, and pileup value for standard deviation; mean fragment size is provided (choose the one with highest correlation)
    frag_qc = attach_back(workflow, PythonCommand(
        stat_frag_std,
        input = {"r": [target + "_model.R" for target in conf.treatment_targets]},
        output = {"json": conf.json_prefix + "_frag.json", "r": [ target + "_frag_sd.R" for target in conf.treatment_targets ]},
        param = {"samples": conf.treatment_bases,
                 "frag_tool": "BAMSE"},
        name = "macs2 model R script parser"))
    frag_qc.allow_fail = True
    frag_qc.allow_dangling = True

def stat_bedAnnotate(workflow, conf, has_dhs, has_velcro):
    """ Describe peaks' distribution
    # collect meta gene distribution info
    """
    collect_meta2 = attach_back(workflow, PythonCommand(
        json_meta2,
        input={"meta": conf.prefix + ".meta"},
        output={"json": conf.json_prefix + "_meta.json"},
        param={"id": conf.id},
        name="bedAnnotate summary"))
    collect_meta2.allow_fail = True
    collect_meta2.allow_dangling = True

    if has_dhs:
        collect_dhs = attach_back(workflow, PythonCommand(
            json_dhs,
            input={"dhs": conf.prefix + ".dhs",
                   "top_peaks": 5000},
            output={"json": conf.json_prefix + "_dhs.json"},
            name="DHS summary"))
        collect_dhs.allow_dangling = True
        collect_dhs.allow_fail = True

    if has_velcro:
        collect_velcro = attach_back(workflow, PythonCommand(
            json_velcro,
            input={"velcro": conf.prefix + ".velcro",
                   "top_peaks": 5000},
            output={"json": conf.json_prefix + "_velcro.json"},
            name="Velcro summary"))
        collect_velcro.allow_fail = True
        collect_velcro.allow_dangling = True

def tex_fastqc(workflow, conf):
    quality = attach_back(workflow,
        PythonCommand(
            load_latex,
            input={"json": conf.json_prefix + "_fastqc.json",
                   "template": resource_filename("chilin2.modules.fastqc", "fastqc.tex"),
                   "pdf": conf.prefix + "_raw_sequence_qc.pdf"},
            output={"latex": conf.latex_prefix + "_fastqc.tex"}))

    quality.allow_fail = True
    quality.allow_dangling = True

    #these are name, png pairings
    if not conf.pe:
        gccontent_graphs = [(nm.replace("_"," "),
                             os.path.join(conf.target_dir, "%s_100k_fastqc" % nm,
                                          "Images","per_sequence_gc_content.png"))\
                                for nm in conf.sample_bases]
    else:
        gccontent_graphs = [(nm.replace("_"," "),
                             os.path.join(conf.target_dir, "%spair1_100k_fastqc" % nm,
                                          "Images","per_sequence_gc_content.png"))\
                                for nm in conf.sample_bases]
    gc = attach_back(workflow,
        PythonCommand(
            load_gc_latex,
            input={"template": resource_filename("chilin2.modules.fastqc", "fastqc_gc.tex"),
                   "gccontent_graphs":gccontent_graphs },
            output={"latex": conf.latex_prefix + "_fastqc_gc.tex"}))

    gc.allow_fail = True
    gc.allow_dangling = True

def tex_conserv(workflow, conf):
    attach_back(workflow,
        PythonCommand(
            latex_conservation,
            input={"template": resource_filename("chilin2.modules.conservation", "conservation.tex")},
            output={"latex": conf.latex_prefix + "_conserv.tex"},
            param = {"prefix": conf.prefix}))

def _bowtie_latex(workflow, conf):
    attach_back(workflow,
        PythonCommand(
            latex_bowtie,
            input={"json": conf.json_prefix + "_bowtie.json",
                   "template": latex_template},
            output={"latex": conf.latex_prefix + "_bowtie.latex"}))

def _macs2_cor_latex(workflow, conf):
    attach_back(workflow,
        PythonCommand(
            latex_cor,
            input={"json": conf.json_prefix + "_cor.json",
                   "template": latex_template},
            output={"latex": conf.latex_prefix + "_cor.latex"}))

def tex_bwa(workflow, conf):
    attach_back(workflow,
                PythonCommand(
                    long_tex,
                    input = {"template": resource_filename("chilin2.modules.bwa", "bwa.tex"),
                             "figure": conf.prefix + "_bwa_compare.pdf"},
                    output = {"latex": conf.latex_prefix + "_map.tex"}))

def _begin_latex(workflow, conf):
    attach_back(workflow,
        PythonCommand(
            latex_start,
            input={"template": latex_template},
            output={"latex": conf.latex_prefix + "_start.latex"},
            param={"id": conf.id}))

def prepare_clean_up(workflow, conf):
    """
    package all the necessary results and delete temporary files
    """
    p_list = ['*.bam', '*.xls', '*_summits.bed', '*_peaks.bed', '*.bw',
              '*.png', '*.pdf', '*.R', '*.zip', '*cor*', 'json', "*summary*",
              "*seqpos","*fastqc", '*latex', "*.conf"]

    p_pattern = [os.path.join(conf.target_dir, p) for p in p_list]

    final_dir = conf.target_dir + '/dataset_' + conf.id
    attach_back(workflow,
        ShellCommand("if [ ! -d '{output}' ]; then mkdir -p {output}; fi",
            output=final_dir))

    for pf in p_pattern:
        if not glob(pf):
            print(pf)
            continue
        move = attach_back(workflow,
            ShellCommand('mv {param[preserve_files]} {output[dir]} \n# Pattern: {param[p_pattern]}',
                output={"dir": final_dir},
                param={"preserve_files": " ".join(glob(pf)),
                       "p_pattern": pf}, ))
        move.allow_fail = True

def _seqpos_latex(workflow, conf):
    attach_back(workflow,
        PythonCommand(
            latex_seqpos,
            input={"json": conf.json_prefix + "_seqpos.json",
                   "template": latex_template},
            output={"latex": conf.latex_prefix + "_seqpos.latex"}))

def _summary_table_latex(workflow, conf):
    attach_back(workflow,
        PythonCommand(
            latex_summary_table,
            input={"template": latex_template},
            output={"latex": conf.latex_prefix + "_summary_table.latex"},
            param={"conf": conf}))

def contamination_check(workflow, conf):
    """
    bowtie mapping back to different species
    """
    if conf.items("contamination"):
        for target in conf.sample_targets:
            for species in dict(conf.items("contamination")):
                index = conf.get("contamination", species)
                if conf.mapper == "bwa":
                    output = target + species + ".sam"
                    if conf.pe:
                        outsai = [target + species + "pair1.sai", target + species + "pair2.sai"]
                        targets = [ target + "pair1", target + "pair2" ]
                    else:
                        outsai = target + species + ".sai"
                        targets = target
                    bwa(workflow, conf, targets, output, outsai, index)
                elif conf.mapper == "bowtie":
                    output = target + species + ".sam"
                    bowtie(workflow, conf, target, output, index)
                elif conf.mapper == "star":
                    output = target + species + "Aligned.out.sam"
                    star(workflow, conf, target, output, index)

                sam2bam = attach_back(workflow,  ## use mapping quality 1 defined by samtools official FAQ
                            ShellCommand(
                                """
                                {tool} view -bS -t {param[genome]} -q {param[mapq]} {input[sam]} > {param[tmp_bam]} && {tool} sort -m {param[max_mem]} {param[tmp_bam]} {param[output_prefix]}
                                """,
                                tool="samtools",
                                input={"sam": output},
                                output={"bam":target + species + ".bam"},
                                param={"tmp_bam": target + species + ".tmp.bam", "output_prefix": target + species,
                                       "mapq": 1,
                                       "genome": conf.get(conf.get("basics", "species"), "chrom_len"),
                                       "max_mem": 4000000000},
                                name = "filtering mapping and convert")) # Use 5G memory as default

                sam2bam.update(param=conf.items("sam2bam"))
                sam2bam.allow_dangling = True
                sam2bam.allow_fail = True

                rem = attach_back(workflow, ShellCommand(
                """
                {tool} view -Sc {input[sam]} > {output[total]}
                {tool} flagstat {input[bam]} > {output[stat]}
                """,
                tool = "samtools",
                input = {"bam": target + species + ".bam",
                         "sam": output},
                output = {"stat": target + species + "_mapped." + conf.mapper,
                          "total": target + species + "_total." + conf.mapper},
                name = "contamination calculation"))
                rem.allow_fail = True
                rem.allow_dangling = True

        ## QC part
        stat_contamination(workflow, conf)
        if conf.long:
            tex_contamination(workflow, conf)

def reg_potential(workflow, conf):
    """

    """
    get_top_peaks = attach_back(workflow,
                                ShellCommand(
                                    "{tool} -n {param[peaks]} {input} | cut -f 1,2,3,4,9> {output}",
                                    tool="head",
                                    input=conf.prefix + "_sort_peaks.narrowPeak" if conf.get("macs2", "type") in ["both", "narrow"] else conf.prefix + "_b_sort_peaks.broadPeak",
                                    output=conf.prefix + "_peaks_top_reg.bed",
                                    param={"peaks": 10000},
                                    name="top summits for regpotential"))
    get_top_peaks.update(param=conf.items("reg_potential"))

    reg = attach_back(workflow,
                      ShellCommand(
                          "{tool} -t {input[peaks]} -g {param[geneTable]} -n {param[prefix]} -d {param[dist]}",
                          tool = "RegPotential.py",
                          input = {"peaks": conf.prefix + "_peaks_top_reg.bed"},
                          output = {"potential": conf.prefix + "_gene_score.txt"},
                          param = {"geneTable": conf.get_path(conf.get("basics", "species"), "geneTable"),
                                   "tool": resource_filename("chilin2.modules", "regulatory/RegPotential.py"),
                                   "prefix": conf.prefix,
                                   "dist": 100000},
                          name = "Regulatory Potential"))
    reg.update(param=conf.items("reg_potential"))

def read_quality(workflow, conf, tex):
    if conf.pe:
        for raw, target in conf.treatment_pairs_pe:
            attach_back(workflow,
                        ShellCommand(
                            "{tool} {input[fastq][0]} {input[fastq][1]} {output[stat][0]} {output[stat][1]}",
                            tool = "dac_pe_read_quality",
                            input = {"fastq": raw},
                            output = {"stat": [ i + "_read_quality.qc" for i in target ]}))

        # attach_back(workflow, PythonCommand(stat_fastqStat,
        #                                     input = {"seq": [ [ p + "_100k.seq" for p in target ] for target in conf.treatment_pair_data ]},
        #                                     output = {"json": conf.json_prefix + "_seq_quality.json"},
        #                                     param = {"samples": conf.treatment_bases, "seq_type": conf.pe}))
        # attach_back(workflow, PythonCommand(
        #     seq_quality_doc,
        #     input = {"tex": tex, "json": conf.json_prefix + "_seq_quality.json"},
        #     output = {"seq": conf.latex_prefix + "seq_quality.tex", "len": conf.latex_prefix + "len.tex"},
        #     param = {"seq_type": conf.seq_type, "reps": len(conf.treatment_pairs),
        #              "pe_samples": conf.treatment_bases}))
    else:
        for raw, target in conf.treatment_pairs:
            sample_fq = {"stat": target + "_read_quality.qc"}
            attach_back(workflow,
                        ShellCommand(
                            "{tool} {input} {output[stat]}",
                            tool = "dac_se_read_quality",
                            input = raw,
                            output = sample_fq,
                            name = "100k read sequence quality and sequence length"))

def _macs2_cor(workflow, conf):
    cor_on_bw = attach_back(workflow,
        ShellCommand(
            template=
            """{tool} \
            -s {param[wig_correlation_step]}  \
            --min-score {param[wig_correlation_min]} --max-score {param[wig_correlation_max]} \
            -r {output[R]} {param[bw]}  {param[rep]} && \
            mv {output[R]}.pdf {output[pdf]}""",
            tool="bigwig_correlation.py",
            input=[target + "_treat.bw" for target in conf.treatment_targets],
            output={"R": conf.prefix + "_cor.R", "pdf": conf.prefix + "_cor.pdf"},
            param={"wig_correlation_method": "mean",
                   "wig_correlation_min": 2,
                   "wig_correlation_max": 50,
                   "wig_correlation_step": 10,},
            name="cor_on_bw"))
    cor_on_bw.param["bw"] = " ".join(cor_on_bw.input)
    cor_on_bw.param["rep"] = " ".join([" -l replicate_%s" % (x + 1) for x in range(len(conf.treatment_pairs))])
    cor_on_bw.update(param=conf.items("correlation"))
    cor_on_bw.allow_fail = True

    attach_back(workflow,
        PythonCommand(
            stat_cor,
            input={"correlation_R": conf.prefix + "_cor.R",
                   "cor_pdf": conf.prefix + "_cor.pdf"},
            output={"json": conf.json_prefix + "_cor.json"}))

def stat_fastqc(workflow, conf):  # collect raw reads quality and GC contents
    """
    long: generate long pages or not
    """
    sums = []
    for raw, target in conf.sample_pairs:
        if conf.pe:
            sums.append(target[0] + "_100k_fastqc/fastqc_data.txt")
        else:
            sums.append(target + "_100k_fastqc/fastqc_data.txt")

    attach_back(workflow,
                PythonCommand(
                    json_fastqc,
                    input={"fastqc_summaries": sums},
                    output={"json": conf.json_prefix + "_fastqc.json"},
                    param={"ids": conf.sample_bases,
                           "id": conf.id},
                    name = "collect fastqc results"))
    
    if conf.long:  ## prepare long document images and tex
        attach_back(workflow,
        PythonCommand(fastqc_detailed_figure,
                      input = {"dbaccessor": resource_filename("chilin2.modules.dbaccessor", "ChiLinQC.db"),
                               "template": resource_filename("chilin2.modules.summary", "R_culmulative_plot.R"), 
                               "json": conf.json_prefix + "_fastqc.json"},
                      output = {"R": conf.prefix + "_raw_sequence_qc.R",
                                "pdf": conf.prefix + "_raw_sequence_qc.pdf"},
                      param={"ids": conf.sample_bases}))

def _raw_QC_latex(workflow, conf):
    attach_back(workflow,
        PythonCommand(
            latex_fastqc,
            input={"json": conf.json_prefix + "_fastqc.json",
                   "template": latex_template},
            output={"latex": conf.latex_prefix + "_fastqc.latex"}))

def _bwa(workflow, conf):
    """
    incorpate ENCODE ChIP-seq alignment parameters
    """
    for raw, target in conf.treatment_pairs:
        param = {"threads": conf.threads,
                 "index":conf.get(conf.species, "genome_index"),
                 "prefix": target + "_raw_sorted",
                 "qc2": target + "_rawbam_stats.qc"}

        if conf.pe:
            bwa = attach_back(workflow, ShellCommand(
                "{tool} {param[threads]} {param[index]} {input[fastq][0]} {input[fastq][1]} {output[bam]} {output[qc]} {param[prefix]} {param[qc2]}",
                tool = "eap_run_bwa_pe",
                input = {"fastq": raw},
                output = {"bam": target + "_raw_sorted.bam", "qc": target + "_rawbam.qc"},
                param = param,
                name = "pair end mapping"))
        else:
            bwa = attach_back(workflow, ShellCommand(
                "{tool} {param[threads]} {param[index]} {input[fastq]} {output[bam]} {output[qc]} {param[prefix]} {param[qc2]}",
                tool = "eap_run_bwa_se",
                input = {"fastq": raw},
                output = {"bam": target + "_raw_sorted.bam", "qc": target + "_rawbam.qc"},
                param = param,
                name = "single end mapping"))
        bwa.update(param = conf.items("bwa"))