<commit_before># ----------------------------------------------------------------------------
# Copyright (c) 2015, The Deblur Development Team.
#
# Distributed under the terms of the BSD 3-clause License.
#
# The full license is in the file LICENSE, distributed with this software.
# ----------------------------------------------------------------------------

from os.path import splitext, join, basename, isfile, split
from datetime import datetime
from os import stat
from glob import glob
import logging
import re
import scipy
import numpy as np
import subprocess
import time
import warnings

from skbio.parse.sequences import parse_fasta
from biom.table import Table
from biom.util import biom_open, HAVE_H5PY

from deblur.deblurring import deblur


def trim_seqs(input_seqs, trim_len):
    """Trim FASTA sequences to specified length.

    Parameters
    ----------
    input_seqs : iterable of (str, str)
        The list of input sequences in (label, sequence) format
    trim_len : int
        Sequence trimming length

    Returns
    -------
    Generator of (str, str)
        The trimmed sequences in (label, sequence) format
    """

    for label, seq in input_seqs:
        if len(seq) >= trim_len:
            yield label, seq[:trim_len]


def dereplicate_seqs(seqs_fp,
                     output_fp,
                     min_size=2,
                     use_log=False,
                     threads=1):
    """Dereplicate FASTA sequences and remove singletons using VSEARCH.

    Parameters
    ----------
    seqs_fp : string
        filepath to FASTA sequence file
    output_fp : string
        file path to dereplicated sequences (FASTA format)
    min_size : integer, optional
        discard sequences with an abundance value smaller
        than integer
    use_log: boolean, optional
        save the vsearch logfile as well (to output_fp.log)
        default=False
    threads : int, optional
        number of threads to use (0 for all available)
    """
    logger = logging.getLogger(__name__)
    logger.info('dereplicate seqs file %s' % seqs_fp)

    log_name = "%s.log" % output_fp

    params = ['vsearch', '--derep_fulllength', seqs_fp,
              '--output', output_fp, '--sizeout',
              '--fasta_width', '0', '--minuniquesize', str(min_size),
              '--quiet', '--threads', str(threads)]
    if use_log:
        params.extend(['--log', log_name])
    sout, serr, res = _system_call(params)
    if not res == 0:
        logger.error('Problem running vsearch dereplication on file %s' %
                     seqs_fp)
        logger.debug('parameters used:\n%s' % params)
        logger.debug('stdout: %s' % sout)
        logger.debug('stderr: %s' % serr)
        return


def build_index_sortmerna(ref_fp, working_dir):
    """Build a SortMeRNA index for all reference databases.

    Parameters
    ----------
    ref_fp: tuple
        filepaths to FASTA reference databases
    working_dir: string
        working directory path where to store the indexed database

    Returns
    -------
    all_db: tuple
        filepaths to SortMeRNA indexed reference databases
    """
    logger = logging.getLogger(__name__)
    logger.info('build_index_sortmerna files %s to'
                ' dir %s' % (ref_fp, working_dir))
    all_db = []
    for db in ref_fp:
        fasta_dir, fasta_filename = split(db)
        index_basename = splitext(fasta_filename)[0]
        db_output = join(working_dir, index_basename)
        logger.debug('processing file %s into location %s' % (db, db_output))
        params = ['indexdb_rna', '--ref', '%s,%s' %
                  (db, db_output), '--tmpdir', working_dir]
        sout, serr, res = _system_call(params)
        if not res == 0:
            logger.error('Problem running indexdb_rna on file %s to dir %s. '
                         'database not indexed' % (db, db_output))
            logger.debug('stdout: %s' % sout)
            logger.debug('stderr: %s' % serr)
            logger.critical('execution halted')
            raise RuntimeError('Cannot index database file %s' % db)
        logger.debug('file %s indexed' % db)
        all_db.append(db_output)
    return all_db


def remove_artifacts_seqs(seqs_fp,
                          ref_fp,
                          working_dir,
                          ref_db_fp,
                          negate=False,
                          threads=1,
                          verbose=False,
                          sim_thresh=None,
                          coverage_thresh=None):
    """Remove artifacts from FASTA file using SortMeRNA.

    Parameters
    ----------
    seqs_fp: string
        file path to FASTA input sequence file
    ref_fp: tuple
        file path(s) to FASTA database file
    working_dir: string
        working directory path
    ref_db_fp: tuple
        file path(s) to indexed FASTA database
    negate: boolean, optional
        if True, discard all input sequences aligning
        to reference database
    threads: integer, optional
        number of threads to use for SortMeRNA
    verbose: boolean, optional
        If true, output SortMeRNA errors
    sim_thresh: float, optional
        The minimal similarity threshold (between 0 and 1)
        for keeping the sequence
        if None, the default values used are 0.65 for negate=False,
        0.95 for negate=True
    coverage_thresh: float, optional
        The minimal coverage threshold (between 0 and 1)
        for alignments for keeping the sequence
        if Nonr, the default values used are 0.3 for negate=False,
        0.95 for negate=True
    """
    logger = logging.getLogger(__name__)
    logger.info('remove_artifacts_seqs file %s' % seqs_fp)

    if stat(seqs_fp).st_size == 0:
        logger.warn('file %s has size 0, continuing' % seqs_fp, UserWarning)
        return

    if coverage_thresh is None:
        if negate:
            coverage_thresh = 0.95
        else:
            coverage_thresh = 0.3

    if sim_thresh is None:
        if negate:
            sim_thresh = 0.95
        else:
            sim_thresh = 0.65

    output_fp = join(working_dir,
                     "%s.no_artifacts" % basename(seqs_fp))
    blast_output = join(working_dir,
                        '%s.sortmerna' % basename(seqs_fp))
    aligned_seq_ids = set()
    for i, db in enumerate(ref_fp):
        logger.debug('running on ref_fp %s working dir %s refdb_fp %s seqs %s'
                     % (db, working_dir, ref_db_fp[i], seqs_fp))
        # run SortMeRNA
        params = ['sortmerna', '--reads', seqs_fp, '--ref', '%s,%s' %
                  (db, ref_db_fp[i]),
                  '--aligned', blast_output, '--blast', '3', '--best', '1',
                  '--print_all_reads', '-v', '-e', '10']

        sout, serr, res = _system_call(params)
        if not res == 0:
            logger.error('sortmerna error on file %s' % seqs_fp)
            logger.error('stdout : %s' % sout)
            logger.error('stderr : %s' % serr)
            return output_fp

        with open('%s.blast' % blast_output, 'r') as bfl:
            for line in bfl:
                line = line.strip().split('\t')
                # if * means no match
                if line[1] == '*':
                    continue
                # check if % identity[2] and coverage[13] are large enough
                # note e-value is [10]
                if negate:
                    if (float(line[2]) >= sim_thresh*100) and \
                       (float(line[13]) >= coverage_thresh*100):
                        aligned_seq_ids.add(line[0])
                else:
                    if float(line[10]) <= 10:
                        aligned_seq_ids.add(line[0])

    if negate:
        def op(x): return x not in aligned_seq_ids
    else:
        def op(x): return x in aligned_seq_ids

    # if negate = False, only output sequences
    # matching to at least one of the databases
    totalseqs = 0
    okseqs = 0
    badseqs = 0
    with open(seqs_fp, 'U') as seqs_f, open(output_fp, 'w') as out_f:
        for label, seq in parse_fasta(seqs_f):
            totalseqs += 1
            label = label.split()[0]
            if op(label):
                out_f.write(">%s\n%s\n" % (label, seq))
                okseqs += 1
            else:
                badseqs += 1
    logger.info('total sequences %d, passing sequences %d, '
                'failing sequences %d' % (totalseqs, okseqs, badseqs))
    return output_fp


def multiple_sequence_alignment(seqs_fp, threads=1):
    """Perform multiple sequence alignment on FASTA file using MAFFT.

    Parameters
    ----------
    seqs_fp: string
        filepath to FASTA file for multiple sequence alignment
    threads: integer, optional
        number of threads to use. 0 to use all threads

    Returns
    -------
    msa_fp : str
        name of output alignment file or None if error encountered
    """
    logger = logging.getLogger(__name__)
    logger.info('multiple_sequence_alignment seqs file %s' % seqs_fp)

    # for mafft we use -1 to denote all threads and not 0
    if threads == 0:
        threads = -1

    if stat(seqs_fp).st_size == 0:
        logger.warning('msa failed. file %s has no reads' % seqs_fp)
        return None
    msa_fp = seqs_fp + '.msa'
    params = ['mafft', '--quiet', '--parttree', '--auto',
              '--thread', str(threads), seqs_fp]
    sout, serr, res = _system_call(params, stdoutfilename=msa_fp)
    if not res == 0:
        logger.info('msa failed for file %s (maybe only 1 read?)' % seqs_fp)
        logger.debug('stderr : %s' % serr)
        return None
    return msa_fp


def remove_chimeras_denovo_from_seqs(seqs_fp, working_dir, threads=1):
    """Remove chimeras de novo using UCHIME (VSEARCH implementation).

    Parameters
    ----------
    seqs_fp: string
        file path to FASTA input sequence file
    output_fp: string
        file path to store chimera-free results
    threads : int
        number of threads (0 for all cores)

    Returns
    -------
    output_fp
        the chimera removed fasta file name
    """
    logger = logging.getLogger(__name__)
    logger.info('remove_chimeras_denovo_from_seqs seqs file %s'
                'to working dir %s' % (seqs_fp, working_dir))

    output_fp = join(
        working_dir, "%s.no_chimeras" % basename(seqs_fp))

    # we use the parameters dn=0.000001, xn=1000, minh=10000000
    # so 1 mismatch in the A/B region will cancel it being labeled as chimera
    # and ~3 unique reads in each region will make it a chimera if
    # no mismatches
    params = ['vsearch', '--uchime_denovo', seqs_fp,
              '--nonchimeras', output_fp,
              '-dn', '0.000001', '-xn', '1000',
              '-minh', '10000000', '--mindiffs', '5',
              '--fasta_width', '0', '--threads', str(threads)]
    sout, serr, res = _system_call(params)
    if not res == 0:
        logger.error('problem with chimera removal for file %s' % seqs_fp)
        logger.debug('stdout : %s' % sout)
        logger.debug('stderr : %s' % serr)
    return output_fp


def split_sequence_file_on_sample_ids_to_files(seqs,
                                               outdir):
    """Split FASTA file on sample IDs.

    Parameters
    ----------
    seqs: file handler
        file handler to demultiplexed FASTA file
    outdir: string
        dirpath to output split FASTA files
    """
    logger = logging.getLogger(__name__)
    logger.info('split_sequence_file_on_sample_ids_to_files'
                ' for file %s into dir %s' % (seqs, outdir))

    outputs = {}
    for bits in parse_fasta(seqs):
        sample = bits[0].split('_', 1)[0]
        if sample not in outputs:
            outputs[sample] = open(join(outdir, sample + '.fasta'), 'w')
        outputs[sample].write(">%s\n%s\n" % (bits[0], bits[1]))
    for sample in outputs:
        outputs[sample].close()
    logger.info('split to %d files' % len(outputs))


def write_biom_table(table, biom_fp):
    """Write BIOM table to file.

    Parameters
    ----------
    table: biom.table
        an instance of a BIOM table
    biom_fp: string
        filepath to output BIOM table
    """
    logger = logging.getLogger(__name__)<commit_msg>fix split problem for samples with _ in name<commit_after># ----------------------------------------------------------------------------
# Copyright (c) 2015, The Deblur Development Team.
#
# Distributed under the terms of the BSD 3-clause License.
#
# The full license is in the file LICENSE, distributed with this software.
# ----------------------------------------------------------------------------

from os.path import splitext, join, basename, isfile, split
from datetime import datetime
from os import stat
from glob import glob
import logging
import re
import scipy
import numpy as np
import subprocess
import time
import warnings

from skbio.parse.sequences import parse_fasta
from biom.table import Table
from biom.util import biom_open, HAVE_H5PY

from deblur.deblurring import deblur


def trim_seqs(input_seqs, trim_len):
    """Trim FASTA sequences to specified length.

    Parameters
    ----------
    input_seqs : iterable of (str, str)
        The list of input sequences in (label, sequence) format
    trim_len : int
        Sequence trimming length

    Returns
    -------
    Generator of (str, str)
        The trimmed sequences in (label, sequence) format
    """

    for label, seq in input_seqs:
        if len(seq) >= trim_len:
            yield label, seq[:trim_len]


def dereplicate_seqs(seqs_fp,
                     output_fp,
                     min_size=2,
                     use_log=False,
                     threads=1):
    """Dereplicate FASTA sequences and remove singletons using VSEARCH.

    Parameters
    ----------
    seqs_fp : string
        filepath to FASTA sequence file
    output_fp : string
        file path to dereplicated sequences (FASTA format)
    min_size : integer, optional
        discard sequences with an abundance value smaller
        than integer
    use_log: boolean, optional
        save the vsearch logfile as well (to output_fp.log)
        default=False
    threads : int, optional
        number of threads to use (0 for all available)
    """
    logger = logging.getLogger(__name__)
    logger.info('dereplicate seqs file %s' % seqs_fp)

    log_name = "%s.log" % output_fp

    params = ['vsearch', '--derep_fulllength', seqs_fp,
              '--output', output_fp, '--sizeout',
              '--fasta_width', '0', '--minuniquesize', str(min_size),
              '--quiet', '--threads', str(threads)]
    if use_log:
        params.extend(['--log', log_name])
    sout, serr, res = _system_call(params)
    if not res == 0:
        logger.error('Problem running vsearch dereplication on file %s' %
                     seqs_fp)
        logger.debug('parameters used:\n%s' % params)
        logger.debug('stdout: %s' % sout)
        logger.debug('stderr: %s' % serr)
        return


def build_index_sortmerna(ref_fp, working_dir):
    """Build a SortMeRNA index for all reference databases.

    Parameters
    ----------
    ref_fp: tuple
        filepaths to FASTA reference databases
    working_dir: string
        working directory path where to store the indexed database

    Returns
    -------
    all_db: tuple
        filepaths to SortMeRNA indexed reference databases
    """
    logger = logging.getLogger(__name__)
    logger.info('build_index_sortmerna files %s to'
                ' dir %s' % (ref_fp, working_dir))
    all_db = []
    for db in ref_fp:
        fasta_dir, fasta_filename = split(db)
        index_basename = splitext(fasta_filename)[0]
        db_output = join(working_dir, index_basename)
        logger.debug('processing file %s into location %s' % (db, db_output))
        params = ['indexdb_rna', '--ref', '%s,%s' %
                  (db, db_output), '--tmpdir', working_dir]
        sout, serr, res = _system_call(params)
        if not res == 0:
            logger.error('Problem running indexdb_rna on file %s to dir %s. '
                         'database not indexed' % (db, db_output))
            logger.debug('stdout: %s' % sout)
            logger.debug('stderr: %s' % serr)
            logger.critical('execution halted')
            raise RuntimeError('Cannot index database file %s' % db)
        logger.debug('file %s indexed' % db)
        all_db.append(db_output)
    return all_db


def remove_artifacts_seqs(seqs_fp,
                          ref_fp,
                          working_dir,
                          ref_db_fp,
                          negate=False,
                          threads=1,
                          verbose=False,
                          sim_thresh=None,
                          coverage_thresh=None):
    """Remove artifacts from FASTA file using SortMeRNA.

    Parameters
    ----------
    seqs_fp: string
        file path to FASTA input sequence file
    ref_fp: tuple
        file path(s) to FASTA database file
    working_dir: string
        working directory path
    ref_db_fp: tuple
        file path(s) to indexed FASTA database
    negate: boolean, optional
        if True, discard all input sequences aligning
        to reference database
    threads: integer, optional
        number of threads to use for SortMeRNA
    verbose: boolean, optional
        If true, output SortMeRNA errors
    sim_thresh: float, optional
        The minimal similarity threshold (between 0 and 1)
        for keeping the sequence
        if None, the default values used are 0.65 for negate=False,
        0.95 for negate=True
    coverage_thresh: float, optional
        The minimal coverage threshold (between 0 and 1)
        for alignments for keeping the sequence
        if Nonr, the default values used are 0.3 for negate=False,
        0.95 for negate=True
    """
    logger = logging.getLogger(__name__)
    logger.info('remove_artifacts_seqs file %s' % seqs_fp)

    if stat(seqs_fp).st_size == 0:
        logger.warn('file %s has size 0, continuing' % seqs_fp, UserWarning)
        return

    if coverage_thresh is None:
        if negate:
            coverage_thresh = 0.95
        else:
            coverage_thresh = 0.3

    if sim_thresh is None:
        if negate:
            sim_thresh = 0.95
        else:
            sim_thresh = 0.65

    output_fp = join(working_dir,
                     "%s.no_artifacts" % basename(seqs_fp))
    blast_output = join(working_dir,
                        '%s.sortmerna' % basename(seqs_fp))
    aligned_seq_ids = set()
    for i, db in enumerate(ref_fp):
        logger.debug('running on ref_fp %s working dir %s refdb_fp %s seqs %s'
                     % (db, working_dir, ref_db_fp[i], seqs_fp))
        # run SortMeRNA
        params = ['sortmerna', '--reads', seqs_fp, '--ref', '%s,%s' %
                  (db, ref_db_fp[i]),
                  '--aligned', blast_output, '--blast', '3', '--best', '1',
                  '--print_all_reads', '-v', '-e', '10']

        sout, serr, res = _system_call(params)
        if not res == 0:
            logger.error('sortmerna error on file %s' % seqs_fp)
            logger.error('stdout : %s' % sout)
            logger.error('stderr : %s' % serr)
            return output_fp

        with open('%s.blast' % blast_output, 'r') as bfl:
            for line in bfl:
                line = line.strip().split('\t')
                # if * means no match
                if line[1] == '*':
                    continue
                # check if % identity[2] and coverage[13] are large enough
                # note e-value is [10]
                if negate:
                    if (float(line[2]) >= sim_thresh*100) and \
                       (float(line[13]) >= coverage_thresh*100):
                        aligned_seq_ids.add(line[0])
                else:
                    if float(line[10]) <= 10:
                        aligned_seq_ids.add(line[0])

    if negate:
        def op(x): return x not in aligned_seq_ids
    else:
        def op(x): return x in aligned_seq_ids

    # if negate = False, only output sequences
    # matching to at least one of the databases
    totalseqs = 0
    okseqs = 0
    badseqs = 0
    with open(seqs_fp, 'U') as seqs_f, open(output_fp, 'w') as out_f:
        for label, seq in parse_fasta(seqs_f):
            totalseqs += 1
            label = label.split()[0]
            if op(label):
                out_f.write(">%s\n%s\n" % (label, seq))
                okseqs += 1
            else:
                badseqs += 1
    logger.info('total sequences %d, passing sequences %d, '
                'failing sequences %d' % (totalseqs, okseqs, badseqs))
    return output_fp


def multiple_sequence_alignment(seqs_fp, threads=1):
    """Perform multiple sequence alignment on FASTA file using MAFFT.

    Parameters
    ----------
    seqs_fp: string
        filepath to FASTA file for multiple sequence alignment
    threads: integer, optional
        number of threads to use. 0 to use all threads

    Returns
    -------
    msa_fp : str
        name of output alignment file or None if error encountered
    """
    logger = logging.getLogger(__name__)
    logger.info('multiple_sequence_alignment seqs file %s' % seqs_fp)

    # for mafft we use -1 to denote all threads and not 0
    if threads == 0:
        threads = -1

    if stat(seqs_fp).st_size == 0:
        logger.warning('msa failed. file %s has no reads' % seqs_fp)
        return None
    msa_fp = seqs_fp + '.msa'
    params = ['mafft', '--quiet', '--parttree', '--auto',
              '--thread', str(threads), seqs_fp]
    sout, serr, res = _system_call(params, stdoutfilename=msa_fp)
    if not res == 0:
        logger.info('msa failed for file %s (maybe only 1 read?)' % seqs_fp)
        logger.debug('stderr : %s' % serr)
        return None
    return msa_fp


def remove_chimeras_denovo_from_seqs(seqs_fp, working_dir, threads=1):
    """Remove chimeras de novo using UCHIME (VSEARCH implementation).

    Parameters
    ----------
    seqs_fp: string
        file path to FASTA input sequence file
    output_fp: string
        file path to store chimera-free results
    threads : int
        number of threads (0 for all cores)

    Returns
    -------
    output_fp
        the chimera removed fasta file name
    """
    logger = logging.getLogger(__name__)
    logger.info('remove_chimeras_denovo_from_seqs seqs file %s'
                'to working dir %s' % (seqs_fp, working_dir))

    output_fp = join(
        working_dir, "%s.no_chimeras" % basename(seqs_fp))

    # we use the parameters dn=0.000001, xn=1000, minh=10000000
    # so 1 mismatch in the A/B region will cancel it being labeled as chimera
    # and ~3 unique reads in each region will make it a chimera if
    # no mismatches
    params = ['vsearch', '--uchime_denovo', seqs_fp,
              '--nonchimeras', output_fp,
              '-dn', '0.000001', '-xn', '1000',
              '-minh', '10000000', '--mindiffs', '5',
              '--fasta_width', '0', '--threads', str(threads)]
    sout, serr, res = _system_call(params)
    if not res == 0:
        logger.error('problem with chimera removal for file %s' % seqs_fp)
        logger.debug('stdout : %s' % sout)
        logger.debug('stderr : %s' % serr)
    return output_fp


def split_sequence_file_on_sample_ids_to_files(seqs,
                                               outdir):
    """Split FASTA file on sample IDs.

    Parameters
    ----------
    seqs: file handler
        file handler to demultiplexed FASTA file
    outdir: string
        dirpath to output split FASTA files
    """
    logger = logging.getLogger(__name__)
    logger.info('split_sequence_file_on_sample_ids_to_files'
                ' for file %s into dir %s' % (seqs, outdir))

    outputs = {}
    for bits in parse_fasta(seqs):
        sample = bits[0].rsplit('_', 1)[0]
        if sample not in outputs:
            outputs[sample] = open(join(outdir, sample + '.fasta'), 'w')
        outputs[sample].write(">%s\n%s\n" % (bits[0], bits[1]))
    for sample in outputs:
        outputs[sample].close()
    logger.info('split to %d files' % len(outputs))


def write_biom_table(table, biom_fp):
    """Write BIOM table to file.

    Parameters
    ----------
    table: biom.table
        an instance of a BIOM table
    biom_fp: string
        filepath to output BIOM table
    """
    logger = logging.getLogger(__name__)
{'fraction_non_alphanumeric': 0.11488143908421913, 'fraction_numerical': 0.010956663941128373, 'mean_word_length': 1.9340290272280196, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5066535', 'n_tokens_mistral': 7556, 'n_tokens_neox': 7119, 'n_words': 2224}