Main convenience wrapper for using usearch to filter/cluster seqs

    The complete 'usearch_qf' process is a multistep process with many calls
    to usearch with various parameters.  It is likely to change from the
    original implementation.  A lot.

    fasta_filepath = fasta filepath to filtering/clustering (e.g., output
     seqs.fna file from split_libraries.py)
    refseqs_fp = fasta filepath for ref-based otu picking.
    output_dir = directory to store the otu mapping file, as well logs and
     the intermediate files created if save_intermediate_files is True.
    percent_ID = percent ID for clustering sequences.
    percent_ID_err = percent ID for filtering out chimeras
    minsize = Minimum size of cluster for retention after chimera removal.
    abundance_skew = threshold setting for chimera removal with de novo
     chimera detection.
    db_filepath = filepath of reference fasta sequence set for ref based
     chimera detection.
    rev = search plus and minus strands of sequences, used in ref based chimera
     detection.
    label_prefix = optional prefix added to filtered fasta file.
    label_suffix = optional suffix added to filtered fasta file.
    retain_label_as_comment = option to add usearch generated label to
     enumerated fasta labels.
    count_start = integer to begin counting at for sequence enumeration.
    perc_id_blast = percent identity setting for using blast algorithm to
     assign original sequence labels to filtered fasta.
    global_alignment = Setting for assignment of original seq labels to filtered
     seqs.
    sizein = not defined in usearch helpstring
    sizeout = not defined in usearch helpstring
    w = Word length for U-sorting
    slots = Size of compressed index table. Should be prime, e.g. 40000003.
     Should also specify --w, typical is --w 16 or --w 32.
    maxrejects = Max rejected targets, 0=ignore, default 32.
    save_intermediate_files = retain all the intermediate files created during
     this process.
    minlen = (not specified in usearch helpstring), but seems like a good bet
     that this refers to the minimum length of the sequences for dereplication.
    HALT_EXEC = used to debug app controller problems.
    de_novo_chimera_detection = If True, will detect chimeras de novo
    reference_chimera_detection = If True, will detect chimeras ref based
    cluster_size_filtering = If True, will filter OTUs according to seq counts.
    remove_usearch_logs = If True, will not call the --log function for each
     usearch call.
    usersort = Used for specifying custom sorting (i.e., non-length based
     sorting) with usearch/uclust.
    suppress_new_clusters = with reference based OTU picking, if enabled,
     will prevent new clusters that do not match the reference from being
     clustered.
    chimeras_retention = accepts either 'intersection' or 'union'.  Will test
     for chimeras against the full input error clustered sequence set, and
     retain sequences flagged as non-chimeras by either (union) or
     only those flagged as non-chimeras by both (intersection).

def usearch_qf(
    fasta_filepath,
    refseqs_fp=None,
    output_dir=None,
    percent_id=0.97,
    percent_id_err=0.97,
    minsize=4,
    abundance_skew=2.0,
    db_filepath=None,
    rev=False,
    label_prefix="",
    label_suffix="",
    retain_label_as_comment=False,
    count_start=0,
    perc_id_blast=0.97,
    save_intermediate_files=False,
    HALT_EXEC=False,
    global_alignment=True,
    sizein=True,
    sizeout=True,
    w=64,
    slots=16769023,
    maxrejects=64,
    minlen=64,
    de_novo_chimera_detection=True,
    derep_fullseq=False,
    reference_chimera_detection=True,
    cluster_size_filtering=True,
    remove_usearch_logs=False,
    usersort=True,
    suppress_new_clusters=False,
    chimeras_retention="union",
    verbose=False
):
    """ Main convenience wrapper for using usearch to filter/cluster seqs

    The complete 'usearch_qf' process is a multistep process with many calls
    to usearch with various parameters.  It is likely to change from the
    original implementation.  A lot.

    fasta_filepath = fasta filepath to filtering/clustering (e.g., output
     seqs.fna file from split_libraries.py)
    refseqs_fp = fasta filepath for ref-based otu picking.
    output_dir = directory to store the otu mapping file, as well logs and
     the intermediate files created if save_intermediate_files is True.
    percent_ID = percent ID for clustering sequences.
    percent_ID_err = percent ID for filtering out chimeras
    minsize = Minimum size of cluster for retention after chimera removal.
    abundance_skew = threshold setting for chimera removal with de novo
     chimera detection.
    db_filepath = filepath of reference fasta sequence set for ref based
     chimera detection.
    rev = search plus and minus strands of sequences, used in ref based chimera
     detection.
    label_prefix = optional prefix added to filtered fasta file.
    label_suffix = optional suffix added to filtered fasta file.
    retain_label_as_comment = option to add usearch generated label to
     enumerated fasta labels.
    count_start = integer to begin counting at for sequence enumeration.
    perc_id_blast = percent identity setting for using blast algorithm to
     assign original sequence labels to filtered fasta.
    global_alignment = Setting for assignment of original seq labels to filtered
     seqs.
    sizein = not defined in usearch helpstring
    sizeout = not defined in usearch helpstring
    w = Word length for U-sorting
    slots = Size of compressed index table. Should be prime, e.g. 40000003.
     Should also specify --w, typical is --w 16 or --w 32.
    maxrejects = Max rejected targets, 0=ignore, default 32.
    save_intermediate_files = retain all the intermediate files created during
     this process.
    minlen = (not specified in usearch helpstring), but seems like a good bet
     that this refers to the minimum length of the sequences for dereplication.
    HALT_EXEC = used to debug app controller problems.
    de_novo_chimera_detection = If True, will detect chimeras de novo
    reference_chimera_detection = If True, will detect chimeras ref based
    cluster_size_filtering = If True, will filter OTUs according to seq counts.
    remove_usearch_logs = If True, will not call the --log function for each
     usearch call.
    usersort = Used for specifying custom sorting (i.e., non-length based
     sorting) with usearch/uclust.
    suppress_new_clusters = with reference based OTU picking, if enabled,
     will prevent new clusters that do not match the reference from being
     clustered.
    chimeras_retention = accepts either 'intersection' or 'union'.  Will test
     for chimeras against the full input error clustered sequence set, and
     retain sequences flagged as non-chimeras by either (union) or
     only those flagged as non-chimeras by both (intersection).
    """

    # Save a list of intermediate filepaths in case they are to be removed.
    intermediate_files = []

    # Need absolute paths to avoid problems with app controller
    if output_dir:
        output_dir = abspath(output_dir) + '/'

    fasta_filepath = abspath(fasta_filepath)

    try:

        if verbose:
            print "Sorting sequences by length..."
        # Sort seqs by length
        app_result, output_filepath_len_sorted =\
            usearch_fasta_sort_from_filepath(fasta_filepath, output_filepath=
                                             join(
                                                 output_dir,
                                                 'len_sorted.fasta'),
                                             save_intermediate_files=save_intermediate_files,
                                             remove_usearch_logs=remove_usearch_logs,
                                             working_dir=output_dir, HALT_EXEC=HALT_EXEC)

        intermediate_files.append(output_filepath_len_sorted)

        if verbose:
            print "Dereplicating sequences..."
        # Dereplicate sequences
        app_result, output_filepath_dereplicated =\
            usearch_dereplicate_exact_subseqs(output_filepath_len_sorted,
                                              output_filepath=join(
                                                  output_dir,
                                                  'dereplicated_seqs.fasta'),
                                              minlen=minlen, w=w, slots=slots, sizeout=sizeout,
                                              maxrejects=maxrejects, save_intermediate_files=save_intermediate_files,
                                              remove_usearch_logs=remove_usearch_logs,
                                              working_dir=output_dir, HALT_EXEC=HALT_EXEC)

        intermediate_files.append(output_filepath_dereplicated)

        if verbose:
            print "Sorting by abundance..."
        # Sort by abundance, initially no filter based on seqs/otu
        app_result, output_fp =\
            usearch_sort_by_abundance(output_filepath_dereplicated,
                                      output_filepath=join(
                                          output_dir,
                                          'abundance_sorted.fasta'),
                                      usersort=True, sizein=sizein, sizeout=sizeout, minsize=0,
                                      remove_usearch_logs=remove_usearch_logs, working_dir=output_dir,
                                      HALT_EXEC=HALT_EXEC)

        intermediate_files.append(output_fp)

        if verbose:
            print "Clustering sequences for error correction..."

        # Create .uc file of clusters file, to identify original sequences
        # later
        output_uc_filepath = output_dir + 'err_corrected_clusters.uc'

        app_result, error_clustered_output_fp =\
            usearch_cluster_error_correction(output_fp,
                                             output_filepath=join(output_dir,
                                                                  'clustered_error_corrected.fasta'),
                                             output_uc_filepath=output_uc_filepath,
                                             usersort=True, percent_id_err=percent_id_err, sizein=sizein,
                                             sizeout=sizeout, w=w, slots=slots, maxrejects=maxrejects,
                                             remove_usearch_logs=remove_usearch_logs,
                                             save_intermediate_files=save_intermediate_files,
                                             working_dir=output_dir, HALT_EXEC=HALT_EXEC)

        intermediate_files.append(error_clustered_output_fp)
        intermediate_files.append(output_uc_filepath)

        # Series of conditional tests, using generic 'output_fp' name so the
        # conditional filtering, if any/all are selected, do not matter.
        if de_novo_chimera_detection:

            if verbose:
                print "Performing de novo chimera detection..."
            app_result, output_fp_de_novo_nonchimeras =\
                usearch_chimera_filter_de_novo(error_clustered_output_fp,
                                               abundance_skew=abundance_skew, output_chimera_filepath=
                                               join(
                                                   output_dir,
                                                   'de_novo_chimeras.fasta'),
                                               output_non_chimera_filepath=join(
                                                   output_dir,
                                                   'de_novo_non_chimeras.fasta'), usersort=True,
                                               save_intermediate_files=save_intermediate_files,
                                               remove_usearch_logs=remove_usearch_logs, working_dir=output_dir,
                                               HALT_EXEC=HALT_EXEC)

            intermediate_files.append(output_fp_de_novo_nonchimeras)

            output_fp = output_fp_de_novo_nonchimeras

        if reference_chimera_detection:
            if verbose:
                print "Performing reference based chimera detection..."

            app_result, output_fp_ref_nonchimeras =\
                usearch_chimera_filter_ref_based(error_clustered_output_fp,
                                                 db_filepath=db_filepath, output_chimera_filepath=
                                                 join(
                                                     output_dir,
                                                     'reference_chimeras.fasta'),
                                                 output_non_chimera_filepath=
                                                 join(output_dir, 'reference_non_chimeras.fasta'), usersort=True,
                                                 save_intermediate_files=save_intermediate_files, rev=rev,
                                                 remove_usearch_logs=remove_usearch_logs, working_dir=output_dir,
                                                 HALT_EXEC=HALT_EXEC)

            intermediate_files.append(output_fp_ref_nonchimeras)

            output_fp = output_fp_ref_nonchimeras

        # get intersection or union if both ref and de novo chimera detection
        if de_novo_chimera_detection and reference_chimera_detection:
            if verbose:
                print "Finding %s of non-chimeras..." % chimeras_retention
            output_fp = get_retained_chimeras(
                output_fp_de_novo_nonchimeras, output_fp_ref_nonchimeras,
                output_combined_fp=
                join(output_dir, 'combined_non_chimeras.fasta'),
                chimeras_retention=chimeras_retention)

            intermediate_files.append(output_fp)

        if cluster_size_filtering:
            # Test for empty filepath following filters, raise error if all seqs
            # have been removed
            if verbose:
                print "Filtering by cluster size..."
            # chimera detection was not performed, use output file of step 4 as input
            # to filtering by cluster size
            if not (reference_chimera_detection and de_novo_chimera_detection):
                output_fp = error_clustered_output_fp
            app_result, output_fp =\
                usearch_sort_by_abundance(output_fp, output_filepath=
                                          join(output_dir, 'abundance_sorted_minsize_' + str(minsize) +
                                               '.fasta'),
                                          minsize=minsize, sizein=sizein, sizeout=sizeout,
                                          remove_usearch_logs=remove_usearch_logs, working_dir=output_dir,
                                          HALT_EXEC=HALT_EXEC)

            intermediate_files.append(output_fp)

        # cluster seqs
        # Should we add in option to use alternative OTU picking here?
        # Seems like it will be a bit of a mess...maybe after we determine
        # if usearch_qf should become standard.
        if refseqs_fp:
            if verbose:
                print "Clustering against reference sequences..."
            app_result, output_filepath =\
                usearch_cluster_seqs_ref(output_fp, output_filepath=
                                         join(
                                             output_dir,
                                             'ref_clustered_seqs.uc'),
                                         percent_id=percent_id, sizein=sizein,
                                         sizeout=sizeout, w=w, slots=slots, maxrejects=maxrejects,
                                         save_intermediate_files=save_intermediate_files,
                                         remove_usearch_logs=remove_usearch_logs,
                                         suppress_new_clusters=suppress_new_clusters, refseqs_fp=refseqs_fp,
                                         output_dir=output_dir, working_dir=output_dir, rev=rev,
                                         HALT_EXEC=HALT_EXEC
                                         )

        else:
            if verbose:
                print "De novo clustering sequences..."
            app_result, output_filepath =\
                usearch_cluster_seqs(output_fp, output_filepath=
                                     join(output_dir, 'clustered_seqs.fasta'),
                                     percent_id=percent_id, sizein=sizein,
                                     sizeout=sizeout, w=w, slots=slots, maxrejects=maxrejects,
                                     save_intermediate_files=save_intermediate_files,
                                     remove_usearch_logs=remove_usearch_logs, working_dir=output_dir,
                                     HALT_EXEC=HALT_EXEC)

        intermediate_files.append(output_filepath)

        # Enumerate the OTUs in the clusters
        if not suppress_new_clusters:
            if verbose:
                print "Enumerating OTUs..."
            output_filepath =\
                enumerate_otus(output_filepath, output_filepath=
                               join(output_dir, 'enumerated_otus.fasta'),
                               label_prefix=label_prefix,
                               label_suffix=label_suffix, count_start=count_start,
                               retain_label_as_comment=retain_label_as_comment)

            intermediate_files.append(output_filepath)

        # Get original sequence label identities
        if verbose:
            print "Assigning sequences to clusters..."
        app_result, clusters_file = assign_reads_to_otus(fasta_filepath,
                                                         filtered_fasta=output_filepath, output_filepath=join(
                                                             output_dir,
                                                             'assign_reads_to_otus.uc'), perc_id_blast=percent_id,
                                                         global_alignment=global_alignment,
                                                         remove_usearch_logs=remove_usearch_logs, working_dir=output_dir,
                                                         HALT_EXEC=HALT_EXEC)

        intermediate_files.append(clusters_file)

    except ApplicationError:
        raise ApplicationError('Error running usearch. Possible causes are '
                               'unsupported version (current supported version is usearch ' +
                               'v5.2.236) is installed or improperly formatted input file was ' +
                               'provided')
    except ApplicationNotFoundError:
        remove_files(files_to_remove)
        raise ApplicationNotFoundError('usearch not found, is it properly ' +
                                       'installed?')

    # Get dict of clusters, list of failures
    # Set OTU ID field to 9 for the case of closed reference OTU picking
    if suppress_new_clusters:
        otu_id_field = 9
    else:
        otu_id_field = 1
    clusters, failures = clusters_from_blast_uc_file(open(clusters_file, "U"),
                                                     otu_id_field)

    # Remove temp files unless user specifies output filepath
    if not save_intermediate_files:
        remove_files(intermediate_files)

    return clusters, failures