count.py

import logging
import os
import re
from urllib.parse import urlparse

from .config import get_bustools_binary_path, get_kallisto_binary_path
from .constants import (
    ADATA_PREFIX,
    BUS_CDNA_PREFIX,
    BUS_FILENAME,
    BUS_FILTERED_FILENAME,
    BUS_FILTERED_SUFFIX,
    BUS_INTRON_PREFIX,
    BUS_S_FILENAME,
    BUS_SC_FILENAME,
    BUS_UNFILTERED_FILENAME,
    BUS_UNFILTERED_SUFFIX,
    COUNTS_PREFIX,
    ECMAP_FILENAME,
    FILTER_WHITELIST_FILENAME,
    FILTERED_COUNTS_DIR,
    INSPECT_FILENAME,
    TCC_PREFIX,
    TXNAMES_FILENAME,
    UNFILTERED_COUNTS_DIR,
    WHITELIST_FILENAME,
)
from .utils import (
    copy_whitelist,
    import_matrix_as_anndata,
    import_tcc_matrix_as_anndata,
    make_directory,
    overlay_anndatas,
    run_executable,
    stream_file,
    sum_anndatas,
    whitelist_provided,
)
from .validate import validate_files

logger = logging.getLogger(__name__)

INSPECT_PARSER = re.compile(r'^.*?(?P<count>[0-9]+)')


@validate_files()
def kallisto_bus(fastqs, index_path, technology, out_dir, threads=8):
    """Runs `kallisto bus`.

    :param fastqs: list of FASTQ file paths
    :type fastqs: list
    :param index_path: path to kallisto index
    :type index_path: str
    :param technology: single-cell technology used
    :type technology: str
    :param out_dir: path to output directory
    :type out_dir: str
    :param threads: number of threads to use, defaults to `8`
    :type threads: int, optional

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    logger.info('Generating BUS file from')
    for fastq in fastqs:
        logger.info((' ' * 8) + fastq)
    command = [get_kallisto_binary_path(), 'bus']
    command += ['-i', index_path]
    command += ['-o', out_dir]
    command += ['-x', technology]
    command += ['-t', threads]
    command += fastqs
    run_executable(command)

    bus_path = os.path.join(out_dir, BUS_FILENAME)
    return {
        'bus': bus_path,
        'ecmap': os.path.join(out_dir, ECMAP_FILENAME),
        'txnames': os.path.join(out_dir, TXNAMES_FILENAME),
    }


@validate_files()
def bustools_sort(bus_path, out_path, temp_dir='tmp', threads=8, memory='4G'):
    """Runs `bustools sort`.

    :param bus_path: path to BUS file to sort
    :type bus_path: str
    :param out_dir: path to output directory
    :type out_dir: str
    :param temp_dir: path to temporary directory, defaults to `tmp`
    :type temp_dir: str, optional
    :param threads: number of threads to use, defaults to `8`
    :type threads: int, optional
    :param memory: amount of memory to use, defaults to `4G`
    :type memory: str, optional

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    logger.info('Sorting BUS file {} to {}'.format(bus_path, out_path))
    command = [get_bustools_binary_path(), 'sort']
    command += ['-o', out_path]
    command += ['-T', temp_dir]
    command += ['-t', threads]
    command += ['-m', memory]
    command += [bus_path]
    run_executable(command)
    return {'bus': out_path}


@validate_files(pre=False)
def bustools_inspect(bus_path, out_path, whitelist_path, ecmap_path):
    """Runs `bustools inspect`.

    :param bus_path: path to BUS file to sort
    :type bus_path: str
    :param out_path: path to output inspect JSON file
    :type out_path: str
    :param whitelist_path: path to whitelist
    :type whitelist_path: str
    :param ecmap_path: path to ecmap file, as generated by `kallisto bus`
    :type ecmap_path: str

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    logger.info('Inspecting BUS file {}'.format(bus_path))
    command = [get_bustools_binary_path(), 'inspect']
    command += ['-o', out_path]
    command += ['-w', whitelist_path]
    command += ['-e', ecmap_path]
    command += [bus_path]
    run_executable(command)
    return {'inspect': out_path}


@validate_files(pre=False)
def bustools_correct(bus_path, out_path, whitelist_path):
    """Runs `bustools correct`.

    :param bus_path: path to BUS file to correct
    :type bus_path: str
    :param out_path: path to output corrected BUS file
    :type out_path: str
    :param whitelist_path: path to whitelist
    :type whitelist_path: str

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    logger.info(
        'Correcting BUS records in {} to {} with whitelist {}'.format(
            bus_path, out_path, whitelist_path
        )
    )
    command = [get_bustools_binary_path(), 'correct']
    command += ['-o', out_path]
    command += ['-w', whitelist_path]
    command += [bus_path]
    run_executable(command)
    return {'bus': out_path}


@validate_files(pre=False)
def bustools_count(
    bus_path, out_prefix, t2g_path, ecmap_path, txnames_path, tcc=False
):
    """Runs `bustools count`.

    :param bus_path: path to BUS file to correct
    :type bus_path: str
    :param out_prefix: prefix of the output files to generate
    :type out_prefix: str
    :param t2g_path: path to output transcript-to-gene mapping
    :type t2g_path: str
    :param ecmap_path: path to ecmap file, as generated by `kallisto bus`
    :type ecmap_path: str
    :param txnames_path: path to transcript names file, as generated by `kallisto bus`
    :type txnames_path: str
    :param tcc: whether to generate a TCC matrix instead of a gene count matrix,
                defaults to `False`
    :type tcc: bool, optional

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    logger.info(
        'Generating count matrix {} from BUS file {}'.format(
            out_prefix, bus_path
        )
    )
    command = [get_bustools_binary_path(), 'count']
    command += ['-o', out_prefix]
    command += ['-g', t2g_path]
    command += ['-e', ecmap_path]
    command += ['-t', txnames_path]
    if not tcc:
        command += ['--genecounts']
    command += [bus_path]
    run_executable(command)
    return {
        'mtx':
            '{}.mtx'.format(out_prefix),
        'ec' if tcc else 'genes':
            '{}.ec.txt'.format(out_prefix)
            if tcc else '{}.genes.txt'.format(out_prefix),
        'barcodes':
            '{}.barcodes.txt'.format(out_prefix),
    }


@validate_files(pre=False)
def bustools_capture(
    bus_path,
    out_path,
    capture_path,
    ecmap_path,
    txnames_path,
    capture_type='transcripts'
):
    """Runs `bustools capture`.

    :param bus_path: path to BUS file to capture
    :type bus_path: str
    :param out_path: path to BUS file to generate
    :type out_path: str
    :param capture_path: path transcripts-to-capture list
    :type capture_path: str
    :param ecmap_path: path to ecmap file, as generated by `kallisto bus`
    :type ecmap_path: str
    :param txnames_path: path to transcript names file, as generated by `kallisto bus`
    :type txnames_path: str
    :param capture_type: the type of information in the capture list.
                      can be one of `transcripts`, `umis`, `barcode`.
    :type capture_type: str

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    logger.info(
        'Capturing records from BUS file {} to {} with capture list {}'.format(
            bus_path, out_path, capture_path
        )
    )
    command = [get_bustools_binary_path(), 'capture']
    command += ['-o', out_path]
    command += ['-c', capture_path]
    command += ['-e', ecmap_path]
    command += ['-t', txnames_path]
    command += ['--complement']
    command += ['--{}'.format(capture_type)]
    command += [bus_path]
    run_executable(command)
    return {'bus': out_path}


@validate_files(pre=False)
def bustools_whitelist(bus_path, out_path):
    """Runs `bustools whitelist`.

    :param bus_path: path to BUS file generate the whitelist from
    :type bus_path: str
    :param out_path: path to output whitelist
    :type out_path: str

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    logger.info(
        'Generating whitelist {} from BUS file {}'.format(out_path, bus_path)
    )
    command = [
        get_bustools_binary_path(), 'whitelist', '-o', out_path, bus_path
    ]
    run_executable(command)
    return {'whitelist': out_path}


def convert_matrix(
    counts_dir,
    matrix_path,
    barcodes_path,
    genes_path=None,
    ec_path=None,
    txnames_path=None,
    loom=False,
    h5ad=False,
    tcc=False,
    threads=8,
):
    """Convert a gene count or TCC matrix to loom or h5ad.

    :param counts_dir: path to counts directory
    :type counts_dir: str
    :param matrix_path: path to matrix
    :type matrix_path: str
    :param barcodes_path: list of paths to barcodes.txt
    :type barcodes_path: str
    :param genes_path: path to genes.txt, defaults to `None`
    :type genes_path: str, optional
    :param ec_path: path to ec.txt, defaults to `None`
    :type ec_path: str, optional
    :param txnames_path: path to transcripts.txt, defaults to `None`
    :type txnames_path: str, optional
    :param loom: whether to generate loom file, defaults to `False`
    :type loom: bool, optional
    :param h5ad: whether to generate h5ad file, defaults to `False`
    :type h5ad: bool, optional
    :param tcc: whether the matrix is a TCC matrix, defaults to `False`
    :type tcc: bool, optional
    :param threads: number of threads to use, defaults to `8`
    :type threads: int, optional

    :return: dictionary of generated files
    :rtype: dict
    """
    results = {}
    logger.info('Reading matrix {}'.format(matrix_path))
    adata = import_tcc_matrix_as_anndata(
        matrix_path, barcodes_path, ec_path, txnames_path, threads=threads
    ) if tcc else import_matrix_as_anndata(
        matrix_path, barcodes_path, genes_path
    )
    if loom:
        loom_path = os.path.join(counts_dir, '{}.loom'.format(ADATA_PREFIX))
        logger.info('Writing matrix to loom {}'.format(loom_path))
        adata.write_loom(loom_path)
        results.update({'loom': loom_path})
    if h5ad:
        h5ad_path = os.path.join(counts_dir, '{}.h5ad'.format(ADATA_PREFIX))
        logger.info('Writing matrix to h5ad {}'.format(h5ad_path))
        adata.write(h5ad_path)
        results.update({'h5ad': h5ad_path})
    return results


def convert_matrices(
    counts_dir,
    matrix_paths,
    barcodes_paths,
    genes_paths=None,
    ec_paths=None,
    txnames_path=None,
    loom=False,
    h5ad=False,
    nucleus=False,
    tcc=False,
    threads=8,
):
    """Convert a gene count or TCC matrix to loom or h5ad.

    :param counts_dir: path to counts directory
    :type counts_dir: str
    :param matrix_paths: list of paths to matrices
    :type matrix_paths: str
    :param barcodes_paths: list of paths to barcodes.txt
    :type barcodes_paths: str
    :param genes_paths: list of paths to genes.txt, defaults to `None`
    :type genes_paths: str, optional
    :param ec_paths: list of path to ec.txt, defaults to `None`
    :type ec_paths: str, optional
    :param txnames_path: list of paths to transcripts.txt, defaults to `None`
    :type txnames_path: str, optional
    :param loom: whether to generate loom file, defaults to `False`
    :type loom: bool, optional
    :param h5ad: whether to generate h5ad file, defaults to `False`
    :type h5ad: bool, optional
    :param nucleus: whether the matrices contain single nucleus counts, defaults to `False`
    :type nucleus: bool, optional
    :param tcc: whether the matrix is a TCC matrix, defaults to `False`
    :type tcc: bool, optional
    :param threads: number of threads to use, defaults to `8`
    :type threads: int, optional

    :return: dictionary of generated files
    :rtype: dict
    """
    results = {}
    adatas = []
    matrix_paths = matrix_paths or []
    barcodes_paths = barcodes_paths or []
    genes_paths = genes_paths or []
    ec_paths = ec_paths or []
    for matrix_path, barcodes_path, genes_ec_path in zip(
            matrix_paths, barcodes_paths, ec_paths
            if not genes_paths or None in genes_paths else genes_paths):
        logger.info('Reading matrix {}'.format(matrix_path))
        adatas.append(
            import_tcc_matrix_as_anndata(
                matrix_path,
                barcodes_path,
                genes_ec_path,
                txnames_path,
                threads=threads
            ) if tcc else
            import_matrix_as_anndata(matrix_path, barcodes_path, genes_ec_path)
        )
    logger.info('Combining matrices')
    adata = sum_anndatas(*adatas) if nucleus else overlay_anndatas(*adatas)
    if loom:
        loom_path = os.path.join(counts_dir, '{}.loom'.format(ADATA_PREFIX))
        logger.info('Writing matrices to loom {}'.format(loom_path))
        adata.write_loom(loom_path)
        results.update({'loom': loom_path})
    if h5ad:
        h5ad_path = os.path.join(counts_dir, '{}.h5ad'.format(ADATA_PREFIX))
        logger.info('Writing matrices to h5ad {}'.format(h5ad_path))
        adata.write(h5ad_path)
        results.update({'h5ad': h5ad_path})
    return results


def filter_with_bustools(
    bus_path,
    ecmap_path,
    txnames_path,
    t2g_path,
    whitelist_path,
    filtered_bus_path,
    counts_prefix=None,
    tcc=False,
    temp_dir='tmp',
    threads=8,
    memory='4G',
    count=True,
    loom=False,
    h5ad=False
):
    """Generate filtered count matrices with bustools.

    :param bus_path: path to sorted, corrected, sorted BUS file
    :type bus_path: str
    :param ecmap_path: path to matrix ec file
    :type ecmap_path: str
    :param txnames_path: path to list of transcripts
    :type txnames_path: str
    :param t2g_path: path to transcript-to-gene mapping
    :type t2g_path: str
    :param whitelist_path: path to filter whitelist to generate
    :type whitelist_path: str
    :param filtered_bus_path: path to filtered BUS file to generate
    :type filtered_bus_path: str
    :param counts_prefix: prefix of count matrix, defaults to `None`
    :type counts_prefix: str, optional
    :param tcc: whether to generate a TCC matrix instead of a gene count matrix,
                defaults to `False`
    :type tcc: bool, optional
    :param temp_dir: path to temporary directory, defaults to `tmp`
    :type temp_dir: str, optional
    :param threads: number of threads to use, defaults to `8`
    :type threads: int, optional
    :param memory: amount of memory to use, defaults to `4G`
    :type memory: str, optional
    :param loom: whether to convert the final count matrix into a loom file,
                 defaults to `False`
    :type loom: bool, optional
    :param h5ad: whether to convert the final count matrix into a h5ad file,
                 defaults to `False`
    :type h5ad: bool, optional

    :return: dictionary of generated files
    :rtype: dict
    """
    logger.info('Filtering with bustools')
    results = {}
    whitelist_result = bustools_whitelist(bus_path, whitelist_path)
    results.update(whitelist_result)
    capture_result = bustools_correct(
        bus_path,
        os.path.join(temp_dir, os.path.basename(filtered_bus_path)),
        whitelist_result['whitelist'],
    )
    sort_result = bustools_sort(
        capture_result['bus'],
        filtered_bus_path,
        temp_dir=temp_dir,
        threads=threads,
        memory=memory,
    )
    results.update({'bus_scs': sort_result['bus']})

    if count:
        counts_dir = os.path.dirname(counts_prefix)
        make_directory(counts_dir)
        count_result = bustools_count(
            sort_result['bus'],
            counts_prefix,
            t2g_path,
            ecmap_path,
            txnames_path,
            tcc=tcc,
        )
        results.update(count_result)

        if loom or h5ad:
            results.update(
                convert_matrix(
                    counts_dir,
                    count_result['mtx'],
                    count_result['barcodes'],
                    genes_path=count_result.get('genes'),
                    ec_path=count_result.get('ec'),
                    txnames_path=txnames_path,
                    loom=loom,
                    h5ad=h5ad,
                    tcc=tcc
                )
            )

    return results


def stream_fastqs(fastqs, temp_dir='tmp'):
    """Given a list of fastqs (that may be local or remote paths), stream any
    remote files. Internally, calls utils.

    :param fastqs: list of (remote or local) fastq paths
    :type fastqs: list
    :param temp_dir: temporary directory
    :type temp_dir: str

    :return: all remote paths substituted with a local path
    :rtype: list
    """
    return [
        stream_file(fastq, os.path.join(temp_dir, os.path.basename(fastq)))
        if urlparse(fastq).scheme in ('http', 'https', 'ftp', 'ftps') else fastq
        for fastq in fastqs
    ]


def copy_or_create_whitelist(technology, bus_path, out_dir):
    """Copies a pre-packaged whitelist if it is provided. Otherwise, runs
    `bustools whitelist` to generate a whitelist.

    :param technology: single-cell technology used
    :type technology: str
    :param bus_path: path to BUS file generate the whitelist from
    :type bus_path: str
    :param out_dir: path to output directory
    :type out_dir: str

    :return: path to copied or generated whitelist
    :rtype: str
    """
    if whitelist_provided(technology):
        logger.info(
            'Copying pre-packaged {} whitelist to {}'.format(
                technology.upper(), out_dir
            )
        )
        return copy_whitelist(technology, out_dir)
    else:
        return bustools_whitelist(
            bus_path, os.path.join(out_dir, WHITELIST_FILENAME)
        )['whitelist']


def count(
    index_path,
    t2g_path,
    technology,
    out_dir,
    fastqs,
    whitelist_path=None,
    tcc=False,
    filter=None,
    temp_dir='tmp',
    threads=8,
    memory='4G',
    overwrite=False,
    loom=False,
    h5ad=False,
):
    """Generates count matrices for single-cell RNA seq.

    :param index_path: path to kallisto index
    :type index_path: str
    :param t2g_path: path to transcript-to-gene mapping
    :type t2g_path: str
    :param technology: single-cell technology used
    :type technology: str
    :param out_dir: path to output directory
    :type out_dir: str
    :param fastqs: list of FASTQ file paths
    :type fastqs: list
    :param whitelist_path: path to whitelist, defaults to `None`
    :type whitelist_path: str, optional
    :param tcc: whether to generate a TCC matrix instead of a gene count matrix,
                defaults to `False`
    :type tcc: bool, optional
    :param filter: filter to use to generate a filtered count matrix,
                   defaults to `None`
    :type filter: str, optional
    :param temp_dir: path to temporary directory, defaults to `tmp`
    :type temp_dir: str, optional
    :param threads: number of threads to use, defaults to `8`
    :type threads: int, optional
    :param memory: amount of memory to use, defaults to `4G`
    :type memory: str, optional
    :param overwrite: overwrite an existing index file, defaults to `False`
    :type overwrite: bool, optional
    :param loom: whether to convert the final count matrix into a loom file,
                 defaults to `False`
    :type loom: bool, optional
    :param h5ad: whether to convert the final count matrix into a h5ad file,
                 defaults to `False`
    :type h5ad: bool, optional

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    results = {}

    make_directory(out_dir)
    unfiltered_results = results.setdefault('unfiltered', {})

    bus_result = {
        'bus': os.path.join(out_dir, BUS_FILENAME),
        'ecmap': os.path.join(out_dir, ECMAP_FILENAME),
        'txnames': os.path.join(out_dir, TXNAMES_FILENAME),
    }
    if any(not os.path.exists(path)
           for name, path in bus_result.items()) or overwrite:
        # Pipe any remote files.
        fastqs = stream_fastqs(fastqs, temp_dir=temp_dir)
        bus_result = kallisto_bus(
            fastqs, index_path, technology, out_dir, threads=threads
        )
    else:
        logger.info(
            'Skipping kallisto bus because output files already exist. Use the --overwrite flag to overwrite.'
        )
    unfiltered_results.update(bus_result)

    sort_result = bustools_sort(
        bus_result['bus'],
        os.path.join(temp_dir, BUS_S_FILENAME),
        temp_dir=temp_dir,
        threads=threads,
        memory=memory
    )
    if not whitelist_path:
        logger.info('Whitelist not provided')
        whitelist_path = copy_or_create_whitelist(
            technology, sort_result['bus'], out_dir
        )
        unfiltered_results.update({'whitelist': whitelist_path})

    inspect_result = bustools_inspect(
        sort_result['bus'], os.path.join(out_dir, INSPECT_FILENAME),
        whitelist_path, bus_result['ecmap']
    )
    unfiltered_results.update(inspect_result)
    correct_result = bustools_correct(
        sort_result['bus'], os.path.join(temp_dir, BUS_SC_FILENAME),
        whitelist_path
    )
    sort2_result = bustools_sort(
        correct_result['bus'],
        os.path.join(out_dir, BUS_UNFILTERED_FILENAME),
        temp_dir=temp_dir,
        threads=threads,
        memory=memory
    )
    unfiltered_results.update({'bus_scs': sort2_result['bus']})

    counts_dir = os.path.join(out_dir, UNFILTERED_COUNTS_DIR)
    make_directory(counts_dir)
    counts_prefix = os.path.join(
        counts_dir, TCC_PREFIX if tcc else COUNTS_PREFIX
    )
    count_result = bustools_count(
        sort2_result['bus'],
        counts_prefix,
        t2g_path,
        bus_result['ecmap'],
        bus_result['txnames'],
        tcc=tcc,
    )
    unfiltered_results.update(count_result)

    # Convert outputs.
    if loom or h5ad:
        unfiltered_results.update(
            convert_matrix(
                counts_dir,
                count_result['mtx'],
                count_result['barcodes'],
                genes_path=count_result.get('genes'),
                ec_path=count_result.get('ec'),
                txnames_path=bus_result['txnames'],
                loom=loom,
                h5ad=h5ad,
                tcc=tcc,
            )
        )

    if filter == 'bustools':
        filtered_counts_prefix = os.path.join(
            out_dir, FILTERED_COUNTS_DIR, TCC_PREFIX if tcc else COUNTS_PREFIX
        )
        filtered_whitelist_path = os.path.join(
            out_dir, FILTER_WHITELIST_FILENAME
        )
        filtered_bus_path = os.path.join(out_dir, BUS_FILTERED_FILENAME)
        results['filtered'] = filter_with_bustools(
            sort2_result['bus'],
            bus_result['ecmap'],
            bus_result['txnames'],
            t2g_path,
            filtered_whitelist_path,
            filtered_bus_path,
            counts_prefix=filtered_counts_prefix,
            tcc=tcc,
            temp_dir=temp_dir,
            threads=threads,
            memory=memory,
            loom=loom,
            h5ad=h5ad
        )

    return results


def count_velocity(
    index_path,
    t2g_path,
    cdna_t2c_path,
    intron_t2c_path,
    technology,
    out_dir,
    fastqs,
    whitelist_path=None,
    tcc=False,
    filter=None,
    temp_dir='tmp',
    threads=8,
    memory='4G',
    overwrite=False,
    loom=False,
    h5ad=False,
    nucleus=False,
):
    """Generates RNA velocity matrices for single-cell RNA seq.

    :param index_path: path to kallisto index
    :type index_path: str
    :param t2g_path: path to transcript-to-gene mapping
    :type t2g_path: str
    :param cdna_t2c_path: path to cDNA transcripts-to-capture file
    :type cdna_t2c_path: str
    :param intron_t2c_path: path to intron transcripts-to-capture file
    :type intron_t2c_path: str
    :param technology: single-cell technology used
    :type technology: str
    :param out_dir: path to output directory
    :type out_dir: str
    :param fastqs: list of FASTQ file paths
    :type fastqs: list
    :param whitelist_path: path to whitelist, defaults to `None`
    :type whitelist_path: str, optional
    :param tcc: whether to generate a TCC matrix instead of a gene count matrix,
                defaults to `False`
    :type tcc: bool, optional
    :param filter: filter to use to generate a filtered count matrix,
                   defaults to `None`
    :type filter: str, optional
    :param temp_dir: path to temporary directory, defaults to `tmp`
    :type temp_dir: str, optional
    :param threads: number of threads to use, defaults to `8`
    :type threads: int, optional
    :param memory: amount of memory to use, defaults to `4G`
    :type memory: str, optional
    :param overwrite: overwrite an existing index file, defaults to `False`
    :type overwrite: bool, optional
    :param loom: whether to convert the final count matrix into a loom file,
                 defaults to `False`
    :type loom: bool, optional
    :param h5ad: whether to convert the final count matrix into a h5ad file,
                 defaults to `False`
    :type h5ad: bool, optional
    :param nucleus: whether this is a single-nucleus experiment. if `True`, the
                    spliced and unspliced count matrices will be summed,
                    defaults to `False`
    :type nucleus: bool, optional

    :return: dictionary containing path to generated index
    :rtype: dict
    """
    results = {}

    make_directory(out_dir)
    unfiltered_results = results.setdefault('unfiltered', {})

    bus_result = {
        'bus': os.path.join(out_dir, BUS_FILENAME),
        'ecmap': os.path.join(out_dir, ECMAP_FILENAME),
        'txnames': os.path.join(out_dir, TXNAMES_FILENAME),
    }
    if any(not os.path.exists(path)
           for name, path in bus_result.items()) or overwrite:
        fastqs = stream_fastqs(fastqs, temp_dir=temp_dir)
        bus_result = kallisto_bus(
            fastqs, index_path, technology, out_dir, threads=threads
        )
    else:
        logger.info(
            'Skipping kallisto bus because output files already exist. Use the --overwrite flag to overwrite.'
        )
    unfiltered_results.update(bus_result)

    sort_result = bustools_sort(
        bus_result['bus'],
        os.path.join(temp_dir, BUS_S_FILENAME),
        temp_dir=temp_dir,
        threads=threads,
        memory=memory
    )
    if not whitelist_path:
        logger.info('Whitelist not provided')
        whitelist_path = copy_or_create_whitelist(
            technology, sort_result['bus'], out_dir
        )
        unfiltered_results.update({'whitelist': whitelist_path})

    inspect_result = bustools_inspect(
        sort_result['bus'], os.path.join(out_dir, INSPECT_FILENAME),
        whitelist_path, bus_result['ecmap']
    )
    unfiltered_results.update(inspect_result)
    correct_result = bustools_correct(
        sort_result['bus'], os.path.join(temp_dir, BUS_SC_FILENAME),
        whitelist_path
    )
    sort2_result = bustools_sort(
        correct_result['bus'],
        os.path.join(out_dir, BUS_UNFILTERED_FILENAME),
        temp_dir=temp_dir,
        threads=threads,
        memory=memory
    )
    unfiltered_results.update({'bus_scs': sort2_result['bus']})

    prefixes = [BUS_CDNA_PREFIX, BUS_INTRON_PREFIX]
    # The prefix and t2cs are swapped because we call bustools capture with
    # the --complement flag.
    prefix_to_t2c = {
        BUS_CDNA_PREFIX: intron_t2c_path,
        BUS_INTRON_PREFIX: cdna_t2c_path,
    }
    counts_dir = os.path.join(out_dir, UNFILTERED_COUNTS_DIR)
    make_directory(counts_dir)
    for prefix, t2c_path in prefix_to_t2c.items():
        capture_result = bustools_capture(
            sort2_result['bus'],
            os.path.join(temp_dir, '{}.bus'.format(prefix)), t2c_path,
            bus_result['ecmap'], bus_result['txnames']
        )
        sort_result = bustools_sort(
            capture_result['bus'],
            os.path.join(out_dir, '{}{}'.format(prefix, BUS_UNFILTERED_SUFFIX)),
            temp_dir=temp_dir,
            threads=threads,
            memory=memory
        )

        if prefix not in unfiltered_results:
            unfiltered_results[prefix] = {}
        unfiltered_results[prefix].update(sort_result)

        counts_prefix = os.path.join(counts_dir, prefix)
        count_result = bustools_count(
            sort_result['bus'],
            counts_prefix,
            t2g_path,
            bus_result['ecmap'],
            bus_result['txnames'],
            tcc=tcc,
        )
        unfiltered_results[prefix].update(count_result)

    if loom or h5ad:
        unfiltered_results.update(
            convert_matrices(
                counts_dir,
                [unfiltered_results[prefix]['mtx'] for prefix in prefixes],
                [unfiltered_results[prefix]['barcodes'] for prefix in prefixes],
                genes_paths=[
                    unfiltered_results[prefix].get('genes')
                    for prefix in prefixes
                ],
                ec_paths=[
                    unfiltered_results[prefix].get('ec') for prefix in prefixes
                ],
                txnames_path=bus_result['txnames'],
                loom=loom,
                h5ad=h5ad,
                tcc=tcc,
                nucleus=nucleus
            )
        )

    if filter:
        filtered_results = results.setdefault('filtered', {})
        if filter == 'bustools':
            filtered_results.update(
                filter_with_bustools(
                    sort2_result['bus'],
                    bus_result['ecmap'],
                    bus_result['txnames'],
                    t2g_path,
                    os.path.join(out_dir, FILTER_WHITELIST_FILENAME),
                    os.path.join(out_dir, BUS_FILTERED_FILENAME),
                    count=False
                )
            )

            for prefix, t2c_path in prefix_to_t2c.items():
                filtered_capture_result = bustools_capture(
                    filtered_results['bus_scs'],
                    os.path.join(temp_dir, '{}.bus'.format(prefix)), t2c_path,
                    bus_result['ecmap'], bus_result['txnames']
                )
                filtered_sort_result = bustools_sort(
                    filtered_capture_result['bus'],
                    os.path.join(
                        out_dir, '{}{}'.format(prefix, BUS_FILTERED_SUFFIX)
                    ),
                    temp_dir=temp_dir,
                    threads=threads,
                    memory=memory
                )

                filtered_results.setdefault(prefix,
                                            {}).update(filtered_sort_result)

                filtered_counts_dir = os.path.join(out_dir, FILTERED_COUNTS_DIR)
                make_directory(filtered_counts_dir)
                filtered_counts_prefix = os.path.join(
                    filtered_counts_dir, prefix
                )
                count_result = bustools_count(
                    filtered_sort_result['bus'],
                    filtered_counts_prefix,
                    t2g_path,
                    bus_result['ecmap'],
                    bus_result['txnames'],
                    tcc=tcc,
                )
                filtered_results[prefix].update(count_result)

        if loom or h5ad:
            filtered_results.update(
                convert_matrices(
                    filtered_counts_dir,
                    [filtered_results[prefix]['mtx'] for prefix in prefixes],
                    [
                        filtered_results[prefix]['barcodes']
                        for prefix in prefixes
                    ],
                    genes_paths=[
                        filtered_results[prefix].get('genes')
                        for prefix in prefixes
                    ],
                    ec_paths=[
                        filtered_results[prefix].get('ec')
                        for prefix in prefixes
                    ],
                    txnames_path=bus_result['txnames'],
                    loom=loom,
                    h5ad=h5ad,
                    tcc=tcc,
                    nucleus=nucleus,
                )
            )

    return results