<commit_before>import slice_miseq
import sam_handler
import os, sys
import subprocess
import Utility
import logging
import argparse
import pool_traceback

LOGGER = logging.getLogger(__name__)
LOGGER.setLevel(logging.DEBUG)
console_handler = logging.StreamHandler(sys.stdout)
formatter = logging.Formatter('%(asctime)s - [%(levelname)s] [%(name)s] [%(process)d] %(message)s')
console_handler.setFormatter(formatter)
LOGGER.addHandler(console_handler)

SELECTION_BF = "QuickSelectionDetection.bf"
HYPHY_EXE = "HYPHYMP"
HYPHY_BASEDIR = "/usr/local/lib/hyphy/TemplateBatchFiles/"

FASTTREE_EXE = "FastTreeMP"
ENV_OMP_NUM_THREADS = 'OMP_NUM_THREADS'

NUC_PER_CODON = 3

CMDLINE_OPTIONS = [
    "sam=",
    "ref=",
    "ref_len=",
    "out_dir=",
    "map_q_thresh=",
    "read_q_thresh=",
    "max_N_thresh=",
    "window_size=",
    "window_breadth_thresh=",
    "window_depth_thresh=",
    "start_nucpos=",
    "end_nucpos=",
    "pvalue=",
    "threads_per_window=",
    "concurrent_windows=",
    "dnds_tsv=",
    "hyphy_exe",
    "hyphy_basedir",
    "fastree_exe"]



def create_full_msa_fasta(sam_filename, out_dir, ref, ref_len, mapping_cutoff, read_qual_cutoff, max_prop_N):
    # Create a pseudo multiple-sequence aligned fasta file
    # TODO:  handle indels in multiple sequence align file
    sam_filename_nopath = os.path.split(sam_filename)[1]
    sam_filename_prefix = os.path.splitext(sam_filename_nopath)[0]
    msa_filename_prefix = out_dir + os.sep + sam_filename_prefix + "." + ref + ".msa"
    msa_fasta_filename = msa_filename_prefix + ".fasta"

    LOGGER.debug("Start Full MSA-Fasta from SAM for ref " + ref)
    if not os.path.exists(msa_fasta_filename) or os.path.getsize(msa_fasta_filename) <= 0:
        sam_handler.create_msa_fasta_from_sam(sam_filename=sam_filename, ref=ref, ref_len=ref_len,
                                              out_fasta_filename=msa_fasta_filename, mapping_cutoff=mapping_cutoff,
                                              read_qual_cutoff=read_qual_cutoff, max_prop_N=max_prop_N)
        LOGGER.debug("Done Full MSA-Fasta from SAM for ref " + ref)
    else:
        LOGGER.warn("Found existing Full MSA-Fasta from SAM for ref " + ref + ".  Not regenerating")

    return msa_fasta_filename


def tabulate_dnds(dnds_tsv_dir, ref, ref_nuc_len, pvalue_thresh, output_dnds_tsv_filename, comments):
    LOGGER.debug("Start Ave Dn/DS for all windows for ref " + ref + " " + output_dnds_tsv_filename)
    seq_dnds_info = slice_miseq.get_seq_dnds_info(dnds_tsv_dir=dnds_tsv_dir, pvalue_thresh=pvalue_thresh, ref=ref,
                                                    ref_codon_len=ref_nuc_len/NUC_PER_CODON)
    with open(output_dnds_tsv_filename, 'w') as dnds_fh:
        dnds_fh.write("# " + comments + "\n")
        dnds_fh.write("Ref\tSite\tdNdS\n")
        for site in range(seq_dnds_info.get_seq_len()):
            site_dnds = seq_dnds_info.get_site_ave_dnds(site)
            dnds_fh.write(ref + "\t" + str(site) + "\t" + str(site_dnds) + "\n")

    LOGGER.debug("Done Ave Dn/DS for all windows  for ref " + ref + ".  Wrote to " + output_dnds_tsv_filename)
    return seq_dnds_info


def eval_windows(ref, ref_len, sam_filename, out_dir, output_dnds_tsv_filename,
                 mapping_cutoff, read_qual_cutoff, max_prop_N,
                 start_nucpos, end_nucpos,
                 window_depth_thresh, window_breadth_thresh, windowsize=300,
                 pvalue=0.05, threads=1, ):
    """
    Slides window along genome.
    Creates the multiple sequence aligned fasta files for the window.
    Feeds the multiple-sequence aligned fasta files to fasttree2 to create a tree.
    Feeds the tree into HyPhy to obtain dn/ds values.
    """
    # TODO:  how to get ORFs?
    # TODO:  do not use kalign2 from PATH
    # TODO:  do not hyphy from PATH

    # Create a pseudo multiple-sequence aligned fasta file
    # TODO:  handle indels in multiple sequence align file
    msa_fasta_filename = create_full_msa_fasta(sam_filename=sam_filename, out_dir=out_dir, ref=ref, ref_len=ref_len,
                          mapping_cutoff=mapping_cutoff, read_qual_cutoff=read_qual_cutoff, max_prop_N=max_prop_N)

    # All nucleotide positions are 1-based<commit_msg>updated sliding_window_tree to output fuller output about substitutions.  Should put this into one function.<commit_after>import slice_miseq
import sam_handler
import os, sys
import subprocess
import Utility
import logging
import argparse
import pool_traceback

LOGGER = logging.getLogger(__name__)
LOGGER.setLevel(logging.DEBUG)
console_handler = logging.StreamHandler(sys.stdout)
formatter = logging.Formatter('%(asctime)s - [%(levelname)s] [%(name)s] [%(process)d] %(message)s')
console_handler.setFormatter(formatter)
LOGGER.addHandler(console_handler)

SELECTION_BF = "QuickSelectionDetection.bf"
HYPHY_EXE = "HYPHYMP"
HYPHY_BASEDIR = "/usr/local/lib/hyphy/TemplateBatchFiles/"

FASTTREE_EXE = "FastTreeMP"
ENV_OMP_NUM_THREADS = 'OMP_NUM_THREADS'

NUC_PER_CODON = 3

CMDLINE_OPTIONS = [
    "sam=",
    "ref=",
    "ref_len=",
    "out_dir=",
    "map_q_thresh=",
    "read_q_thresh=",
    "max_N_thresh=",
    "window_size=",
    "window_breadth_thresh=",
    "window_depth_thresh=",
    "start_nucpos=",
    "end_nucpos=",
    "pvalue=",
    "threads_per_window=",
    "concurrent_windows=",
    "dnds_tsv=",
    "hyphy_exe",
    "hyphy_basedir",
    "fastree_exe"]



def create_full_msa_fasta(sam_filename, out_dir, ref, ref_len, mapping_cutoff, read_qual_cutoff, max_prop_N):
    # Create a pseudo multiple-sequence aligned fasta file
    # TODO:  handle indels in multiple sequence align file
    sam_filename_nopath = os.path.split(sam_filename)[1]
    sam_filename_prefix = os.path.splitext(sam_filename_nopath)[0]
    msa_filename_prefix = out_dir + os.sep + sam_filename_prefix + "." + ref + ".msa"
    msa_fasta_filename = msa_filename_prefix + ".fasta"

    LOGGER.debug("Start Full MSA-Fasta from SAM for ref " + ref)
    if not os.path.exists(msa_fasta_filename) or os.path.getsize(msa_fasta_filename) <= 0:
        sam_handler.create_msa_fasta_from_sam(sam_filename=sam_filename, ref=ref, ref_len=ref_len,
                                              out_fasta_filename=msa_fasta_filename, mapping_cutoff=mapping_cutoff,
                                              read_qual_cutoff=read_qual_cutoff, max_prop_N=max_prop_N)
        LOGGER.debug("Done Full MSA-Fasta from SAM for ref " + ref)
    else:
        LOGGER.warn("Found existing Full MSA-Fasta from SAM for ref " + ref + ".  Not regenerating")

    return msa_fasta_filename


def tabulate_dnds(dnds_tsv_dir, ref, ref_nuc_len, pvalue_thresh, output_dnds_tsv_filename, comments):
    LOGGER.debug("Start Ave Dn/DS for all windows for ref " + ref + " " + output_dnds_tsv_filename)
    seq_dnds_info = slice_miseq.get_seq_dnds_info(dnds_tsv_dir=dnds_tsv_dir, pvalue_thresh=pvalue_thresh, ref=ref,
                                                    ref_codon_len=ref_nuc_len/NUC_PER_CODON)
    with open(output_dnds_tsv_filename, 'w') as dnds_fh:
        dnds_fh.write("# " + comments + "\n")
        dnds_fh.write("Ref\tSite\tdNdS\tWindows\tCodons\tNonSyn\tSyn\tSubst\n")
        for site in range(1, seq_dnds_info.get_seq_len() + 1):
            site_dnds = seq_dnds_info.get_site_ave_dnds(site_1based=site)
            window = seq_dnds_info.get_site_window_cov(site_1based=site)
            reads = seq_dnds_info.get_site_ave_read_cov(site_1based=site)
            nonsyn = seq_dnds_info.get_site_ave_nonsyn_subs(site_1based=site)
            syn = seq_dnds_info.get_site_ave_syn_subs(site_1based=site)
            subs = seq_dnds_info.get_site_ave_subs(site_1based=site)
            line = "\t".join((ref, str(site), str(site_dnds), str(window), str(reads),  str(nonsyn), str(syn), str(subs)))
            dnds_fh.write(line + "\n")

    LOGGER.debug("Done Ave Dn/DS for all windows  for ref " + ref + ".  Wrote to " + output_dnds_tsv_filename)
    return seq_dnds_info


def eval_windows(ref, ref_len, sam_filename, out_dir, output_dnds_tsv_filename,
                 mapping_cutoff, read_qual_cutoff, max_prop_N,
                 start_nucpos, end_nucpos,
                 window_depth_thresh, window_breadth_thresh, windowsize=300,
                 pvalue=0.05, threads=1, ):
    """
    Slides window along genome.
    Creates the multiple sequence aligned fasta files for the window.
    Feeds the multiple-sequence aligned fasta files to fasttree2 to create a tree.
    Feeds the tree into HyPhy to obtain dn/ds values.
    """
    # TODO:  how to get ORFs?
    # TODO:  do not use kalign2 from PATH
    # TODO:  do not hyphy from PATH

    # Create a pseudo multiple-sequence aligned fasta file
    # TODO:  handle indels in multiple sequence align file
    msa_fasta_filename = create_full_msa_fasta(sam_filename=sam_filename, out_dir=out_dir, ref=ref, ref_len=ref_len,
                          mapping_cutoff=mapping_cutoff, read_qual_cutoff=read_qual_cutoff, max_prop_N=max_prop_N)

    # All nucleotide positions are 1-based
{'fraction_non_alphanumeric': 0.0961283185840708, 'fraction_numerical': 0.003982300884955752, 'mean_word_length': 3.484623015873016, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2449526', 'n_tokens_mistral': 3309, 'n_tokens_neox': 3205, 'n_words': 619}