20919343111

Committed 12 Jan 2026 12:27PM UTC coverage: 70.772% (+0.2%) from 70.6%

Build # 20919343111

Build Type

Pull #160

github

Committed by

web-flow

Commit Message

Merge f9b616da0 into dd823d605

Pull Request Pull Request #160: Optimize align_and_fix runtime with sambamba and skip_realign option

Coverage Stats

16 of 18 new or added lines in 1 file covered. (88.89%)

4649 of 6569 relevant lines covered (70.77%)

0.71 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

77.32

/read_utils.py

#!/usr/bin/env python3
"""
Utilities for working with sequence reads, such as converting formats and
fixing mate pairs.
"""

__author__ = "irwin@broadinstitute.org, dpark@broadinstitute.org"
__commands__ = []

import argparse
import logging
import math
import os
import sqlite3
import subprocess
import tempfile
import time
import shutil
import sys
import concurrent.futures
from contextlib import contextmanager
import functools

from Bio import SeqIO
import pysam

import util.cmd
import util.file
import util.misc
from util.file import mkstempfname
import tools.bbmap
import tools.bwa
import tools.cdhit
import tools.picard
import tools.samtools
import tools.minimap2
import tools.mvicuna
import tools.prinseq
import tools.novoalign
import tools.gatk
import tools.sambamba

log = logging.getLogger(__name__)


# ===============================
# ***  index_fasta_samtools   ***
# ===============================


def parser_index_fasta_samtools(parser=argparse.ArgumentParser()):
    parser.add_argument('inFasta', help='Reference genome, FASTA format.')
    util.cmd.common_args(parser, (('loglevel', None), ('version', None)))
    util.cmd.attach_main(parser, main_index_fasta_samtools)
    return parser


def main_index_fasta_samtools(args):
    '''Index a reference genome for Samtools.'''
    tools.samtools.SamtoolsTool().faidx(args.inFasta, overwrite=True)
    return 0


__commands__.append(('index_fasta_samtools', parser_index_fasta_samtools))

# =============================
# ***  index_fasta_picard   ***
# =============================


def parser_index_fasta_picard(parser=argparse.ArgumentParser()):
    parser.add_argument('inFasta', help='Input reference genome, FASTA format.')
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.CreateSequenceDictionaryTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[],
        nargs='*',
        help='Optional arguments to Picard\'s CreateSequenceDictionary, OPTIONNAME=value ...'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_index_fasta_picard)
    return parser


def main_index_fasta_picard(args):
    '''Create an index file for a reference genome suitable for Picard/GATK.'''
    tools.picard.CreateSequenceDictionaryTool().execute(
        args.inFasta, overwrite=True,
        picardOptions=args.picardOptions,
        JVMmemory=args.JVMmemory
    )
    return 0


__commands__.append(('index_fasta_picard', parser_index_fasta_picard))

# =============================
# ***  mkdup_picard   ***
# =============================


def parser_mkdup_picard(parser=argparse.ArgumentParser()):
    parser.add_argument('inBams', help='Input reads, BAM format.', nargs='+')
    parser.add_argument('outBam', help='Output reads, BAM format.')
    parser.add_argument('--outMetrics', help='Output metrics file. Default is to dump to a temp file.', default=None)
    parser.add_argument(
        "--remove",
        help="Instead of marking duplicates, remove them entirely (default: %(default)s)",
        default=False,
        action="store_true",
        dest="remove"
    )
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.MarkDuplicatesTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[],
        nargs='*',
        help='Optional arguments to Picard\'s MarkDuplicates, OPTIONNAME=value ...'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_mkdup_picard)
    return parser


def main_mkdup_picard(args):
    '''Mark or remove duplicate reads from BAM file.'''
    opts = list(args.picardOptions)
    if args.remove:
        opts = ['REMOVE_DUPLICATES=true'] + opts
    tools.picard.MarkDuplicatesTool().execute(
        args.inBams, args.outBam,
        args.outMetrics, picardOptions=opts,
        JVMmemory=args.JVMmemory
    )
    return 0


__commands__.append(('mkdup_picard', parser_mkdup_picard))

# =============================
# ***  revert_bam_picard   ***
# =============================

def parser_revert_sam_common(parser=argparse.ArgumentParser()):
    parser.add_argument('--clearTags', dest='clear_tags', default=False, action='store_true', 
                        help='When supplying an aligned input file, clear the per-read attribute tags')
    parser.add_argument("--tagsToClear", type=str, nargs='+', dest="tags_to_clear", default=["XT", "X0", "X1", "XA", 
                        "AM", "SM", "BQ", "CT", "XN", "OC", "OP"], 
                        help='A space-separated list of tags to remove from all reads in the input bam file (default: %(default)s)')
    parser.add_argument(
        '--doNotSanitize',
        dest="do_not_sanitize",
        action='store_true',
        help="""When being reverted, picard's SANITIZE=true
                is set unless --doNotSanitize is given. 
                Sanitization is a destructive operation that removes reads so the bam file is consistent.
                From the picard documentation:
                    'Reads discarded include (but are not limited to) paired reads with missing mates, duplicated records, 
                    records with mismatches in length of bases and qualities.'
                For more information see:
                    https://broadinstitute.github.io/picard/command-line-overview.html#RevertSam
                """
    )
    try:
        parser.add_argument(
            '--JVMmemory',
            default=tools.picard.RevertSamTool.jvmMemDefault,
            help='JVM virtual memory size (default: %(default)s)'
        )
    except argparse.ArgumentError:
        # if --JVMmemory is already present in the parser, continue on...
        pass
    return parser

def parser_revert_bam_picard(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format.')
    parser.add_argument('outBam', help='Output reads, BAM format.')
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.RevertSamTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[],
        nargs='*',
        help='Optional arguments to Picard\'s RevertSam, OPTIONNAME=value ...'
    )
    parser = parser_revert_sam_common(parser)
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_revert_bam_picard, split_args=True)
    return parser


def main_revert_bam_picard(inBam, outBam, clear_tags=False, tags_to_clear=None, picardOptions=None, do_not_sanitize=False, JVMmemory=None):
    '''Revert BAM to raw reads'''
    picardOptions = picardOptions or []
    tags_to_clear = tags_to_clear or []

    if clear_tags:
        for tag in tags_to_clear:
            picardOptions.append("ATTRIBUTE_TO_CLEAR={}".format(tag))
            
            if not do_not_sanitize and not any(opt.startswith("SANITIZE") for opt in picardOptions):
                picardOptions.append('SANITIZE=true')

    tools.picard.RevertSamTool().execute(inBam, outBam, picardOptions=picardOptions, JVMmemory=JVMmemory)
    return 0
__commands__.append(('revert_bam_picard', parser_revert_bam_picard))

@contextmanager
def revert_bam_if_aligned(inBam, revert_bam=None, clear_tags=True, tags_to_clear=None, picardOptions=None, JVMmemory="4g", sanitize=False):
    '''
        revert the bam file if aligned. 
        If a revertBam file path is specified, it is used, otherwise a temp file is created.
    '''

    revertBamOut = revert_bam if revert_bam else util.file.mkstempfname('.bam')
    picardOptions = picardOptions or []
    tags_to_clear = tags_to_clear or []

    outBam = inBam
    with pysam.AlignmentFile(inBam, 'rb', check_sq=False) as bam:
        # if it looks like the bam is aligned, revert it
        if 'SQ' in bam.header and len(bam.header['SQ'])>0:
            if clear_tags:
                for tag in tags_to_clear:
                    picardOptions.append("ATTRIBUTE_TO_CLEAR={}".format(tag))

            if not any(opt.startswith("SORT_ORDER") for opt in picardOptions):
                picardOptions.append('SORT_ORDER=queryname')

            if sanitize and not any(opt.startswith("SANITIZE") for opt in picardOptions):
                picardOptions.append('SANITIZE=true')

            tools.picard.RevertSamTool().execute(
                inBam, revertBamOut, picardOptions=picardOptions 
            )
            outBam = revertBamOut
        else:
            # if we don't need to produce a revertBam file
            # but the user has specified one anyway
            # simply touch the output
            if revert_bam:
                log.warning("An output was specified for 'revertBam', but the input is unaligned, so RevertSam was not needed. Touching the output.")
                util.file.touch(revertBamOut)
                # TODO: error out? run RevertSam anyway?

    yield outBam

    if revert_bam == None:
        os.unlink(revertBamOut)

# =========================
# ***  generic picard   ***
# =========================


def parser_picard(parser=argparse.ArgumentParser()):
    parser.add_argument('command', help='picard command')
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.PicardTools.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[], nargs='*',
        help='Optional arguments to Picard, OPTIONNAME=value ...'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_picard)
    return parser


def main_picard(args):
    '''Generic Picard runner.'''
    tools.picard.PicardTools().execute(args.command, picardOptions=args.picardOptions, JVMmemory=args.JVMmemory)
    return 0


__commands__.append(('picard', parser_picard))

# ===================
# ***  sort_bam   ***
# ===================


def parser_sort_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input bam file.')
    parser.add_argument('outBam', help='Output bam file, sorted.')
    parser.add_argument(
        'sortOrder',
        help='How to sort the reads. [default: %(default)s]',
        choices=tools.picard.SortSamTool.valid_sort_orders,
        default=tools.picard.SortSamTool.default_sort_order
    )
    parser.add_argument(
        "--index",
        help="Index outBam (default: %(default)s)",
        default=False,
        action="store_true",
        dest="index"
    )
    parser.add_argument(
        "--md5",
        help="MD5 checksum outBam (default: %(default)s)",
        default=False,
        action="store_true",
        dest="md5"
    )
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.SortSamTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[],
        nargs='*',
        help='Optional arguments to Picard\'s SortSam, OPTIONNAME=value ...'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_sort_bam)
    return parser


def main_sort_bam(args):
    '''Sort BAM file'''
    opts = list(args.picardOptions)
    if args.index:
        opts = ['CREATE_INDEX=true'] + opts
    if args.md5:
        opts = ['CREATE_MD5_FILE=true'] + opts
    tools.picard.SortSamTool().execute(
        args.inBam, args.outBam, args.sortOrder,
        picardOptions=opts, JVMmemory=args.JVMmemory
    )
    return 0


__commands__.append(('sort_bam', parser_sort_bam))

# ====================
# ***  downsample_bams  ***
# ====================

def parser_downsample_bams(parser=argparse.ArgumentParser()):
    parser.add_argument('in_bams', help='Input bam files.', nargs='+')
    parser.add_argument('--outPath', dest="out_path", type=str, help="""Output path. If not provided, 
                        downsampled bam files will be written to the same paths as each 
                        source bam file""")
    parser.add_argument('--readCount', dest="specified_read_count", type=int, help='The number of reads to downsample to.')
    group = parser.add_mutually_exclusive_group()
    group.add_argument('--deduplicateBefore', action='store_true', dest="deduplicate_before", help='de-duplicate reads before downsampling.')
    group.add_argument('--deduplicateAfter', action='store_true', dest="deduplicate_after", help='de-duplicate reads after downsampling.')
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.DownsampleSamTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[],
        nargs='*',
        help='Optional arguments to Picard\'s DownsampleSam, OPTIONNAME=value ...'
    )
    util.cmd.common_args(parser, (('threads', None), ('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_downsample_bams, split_args=True)
    return parser


def main_downsample_bams(in_bams, out_path, specified_read_count=None, deduplicate_before=False, deduplicate_after=False, picardOptions=None, threads=None, JVMmemory=None):
    '''Downsample multiple bam files to the smallest read count in common, or to the specified count.'''
    if picardOptions is None:
        picardOptions = []

    opts = list(picardOptions) + []

    def get_read_counts(bams):
        # get read counts for bam files provided
        read_counts = {}
        with concurrent.futures.ProcessPoolExecutor(max_workers=util.misc.sanitize_thread_count(threads)) as executor:
            samtools = tools.samtools.SamtoolsTool()
            for x, result in zip(bams, executor.map(samtools.count, bams)):
                read_counts[x] = int(result)

        return read_counts

    def get_downsample_target_count(bams, readcount_requested):
        read_counts = get_read_counts(bams)
        downsample_target = 0
        if readcount_requested is not None:
            downsample_target = readcount_requested
            if downsample_target > min(read_counts.values()):
                raise ValueError("The smallest file has %s reads, which is less than the downsample target specified, %s. Please reduce the target count or omit it to use the read count of the smallest input file." % (min(read_counts.values()), downsample_target))
        else:
            downsample_target = min(read_counts.values())
        return downsample_target

    def downsample_bams(data_pairs, downsample_target, threads=None, JVMmemory=None):
        workers = min(util.misc.sanitize_thread_count(threads), len(data_pairs))
        JVMmemory = JVMmemory if JVMmemory else tools.picard.DownsampleSamTool.jvmMemDefault
        jvm_worker_memory_mb = str(int(util.misc.convert_size_str(JVMmemory,"m")[:-1])//workers)+"m"
        with concurrent.futures.ProcessPoolExecutor(max_workers=workers) as executor:
            downsamplesam = tools.picard.DownsampleSamTool()
            future_to_file = {executor.submit(downsamplesam.downsample_to_approx_count, *(list(fp)+[downsample_target]), JVMmemory=jvm_worker_memory_mb): fp[0] for fp in data_pairs}
            for future in concurrent.futures.as_completed(future_to_file):
                f = future_to_file[future]
                try:
                    data = future.result()
                except Exception as exc:
                    raise

    def dedup_bams(data_pairs, threads=None):
        workers = min(util.misc.sanitize_thread_count(threads), len(data_pairs))
        with concurrent.futures.ProcessPoolExecutor(max_workers=workers) as executor:
            future_to_file = {executor.submit(rmdup_mvicuna_bam, *fp): fp[0] for fp in data_pairs}
            for future in concurrent.futures.as_completed(future_to_file):
                f = future_to_file[future]
                try:
                    data = future.result()
                except Exception as exc:
                    raise

    data_pairs = []
    if deduplicate_before:
        with util.file.tempfnames(suffixes=[ '{}.dedup.bam'.format(os.path.splitext(os.path.basename(x))[0]) for x in in_bams]) as deduped_bams:
            data_pairs = list(zip(in_bams, deduped_bams))
            dedup_bams(data_pairs, threads=threads)
            downsample_target = get_downsample_target_count(deduped_bams, specified_read_count)
            
            if out_path:
                util.file.mkdir_p(out_path)
                data_pairs = list(zip(deduped_bams, [os.path.join(out_path, os.path.splitext(os.path.basename(x))[0]+".dedup.downsampled-{}.bam".format(downsample_target)) for x in in_bams]))
            else:
                data_pairs = list(zip(deduped_bams, [os.path.splitext(x)[0]+".dedup.downsampled-{}.bam".format(downsample_target) for x in in_bams]))
            downsample_bams(data_pairs, downsample_target, threads=threads, JVMmemory=JVMmemory)
    else:
        downsample_target = get_downsample_target_count(in_bams, specified_read_count)
        if deduplicate_after:
            log.warning("--deduplicateAfter has been specified. Read count of output files is not guaranteed.")
            with util.file.tempfnames(suffixes=[ '{}.downsampled-{}.bam'.format(os.path.splitext(os.path.basename(x))[0], downsample_target) for x in in_bams]) as downsampled_tmp_bams:
                data_pairs = list(zip(in_bams, downsampled_tmp_bams))
                downsample_bams(data_pairs, downsample_target, threads=threads, JVMmemory=JVMmemory)
                
                if out_path:
                    util.file.mkdir_p(out_path)
                    data_pairs = list(zip(downsampled_tmp_bams, [os.path.join(out_path, os.path.splitext(os.path.basename(x))[0]+".downsampled-{}.dedup.bam").format(downsample_target) for x in in_bams]))
                else:
                    data_pairs = list(zip(downsampled_tmp_bams, [os.path.splitext(x)[0]+".downsampled-{}.dedup.bam".format(downsample_target) for x in in_bams]))
                dedup_bams(data_pairs, threads=threads)
        else:
            if out_path:
                util.file.mkdir_p(out_path)
                data_pairs = list(zip(in_bams, [os.path.join(out_path, os.path.splitext(os.path.basename(x))[0]+".downsampled-{}.bam".format(downsample_target)) for x in in_bams]))
            else:
                data_pairs = list(zip(in_bams, [os.path.splitext(x)[0]+".downsampled-{}.bam".format(downsample_target) for x in in_bams]))
            downsample_bams(data_pairs, downsample_target, threads=threads, JVMmemory=JVMmemory)
    return 0

__commands__.append(('downsample_bams', parser_downsample_bams))

# ====================
# ***  merge_bams  ***
# ====================


def parser_merge_bams(parser=argparse.ArgumentParser()):
    parser.add_argument('inBams', help='Input bam files.', nargs='+')
    parser.add_argument('outBam', help='Output bam file.')
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.MergeSamFilesTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[],
        nargs='*',
        help='Optional arguments to Picard\'s MergeSamFiles, OPTIONNAME=value ...'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_merge_bams)
    return parser


def main_merge_bams(args):
    '''Merge multiple BAMs into one'''
    opts = list(args.picardOptions) + ['USE_THREADING=true']
    tools.picard.MergeSamFilesTool().execute(args.inBams, args.outBam, picardOptions=opts, JVMmemory=args.JVMmemory)
    return 0


__commands__.append(('merge_bams', parser_merge_bams))

# ====================
# ***  filter_bam  ***
# ====================


def parser_filter_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input bam file.')
    parser.add_argument('readList', help='Input file of read IDs.')
    parser.add_argument('outBam', help='Output bam file.')
    parser.add_argument(
        "--exclude",
        help="""If specified, readList is a list of reads to remove from input.
            Default behavior is to treat readList as an inclusion list (all unnamed
            reads are removed).""",
        default=False,
        action="store_true",
        dest="exclude"
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_filter_bam)
    return parser


def main_filter_bam(args):
    '''Filter BAM file by read name'''
    with util.file.tmp_dir(suffix='_filter_bam') as tmpdir:
        db_path = os.path.join(tmpdir, 'read_ids.db')
        with ReadIdStore(db_path) as store:
            store.add_from_readlist(args.readList)
            store.filter_bam_by_ids(args.inBam, args.outBam, include=not args.exclude)
    return 0


__commands__.append(('filter_bam', parser_filter_bam))



# =======================
# ***  fastq_to_bam   ***
# =======================


def fastq_to_bam(
    inFastq1,
    inFastq2,
    outBam,
    sampleName=None,
    header=None,
    JVMmemory=None,  # Kept for backwards compatibility, ignored (was for Picard)
    picardOptions=None,  # Kept for backwards compatibility, parsed for RG tags
    threads=None
):
    ''' Convert a pair of fastq paired-end read files and optional text header
        to a single bam file.

        Uses samtools import for multi-threaded FASTQ to BAM conversion.
        The JVMmemory parameter is kept for backwards compatibility but ignored.
        The picardOptions parameter is parsed for Picard-style RG tag options.
    '''
    picardOptions = picardOptions or []

    if header:
        fastqToSamOut = mkstempfname('.bam')
    else:
        fastqToSamOut = outBam
    if sampleName is None:
        sampleName = 'Dummy'    # Will get overwritten by rehead command

    # Parse picardOptions to extract RG tag values
    # Supports Picard-style options like LIBRARY_NAME=Alexandria
    opts_dict = {}
    for opt in picardOptions:
        if '=' in opt:
            key, value = opt.split('=', 1)
            opts_dict[key.upper()] = value

    # Map Picard option names to import_fastq parameters
    library_name = opts_dict.get('LIBRARY_NAME')
    platform = opts_dict.get('PLATFORM', 'illumina')
    platform_unit = opts_dict.get('PLATFORM_UNIT')
    sequencing_center = opts_dict.get('SEQUENCING_CENTER')
    run_date = opts_dict.get('RUN_DATE')
    read_group_id = opts_dict.get('READ_GROUP_NAME')

    tools.samtools.SamtoolsTool().import_fastq(
        inFastq1, inFastq2, fastqToSamOut,
        sample_name=sampleName,
        library_name=library_name,
        platform=platform,
        read_group_id=read_group_id,
        platform_unit=platform_unit,
        sequencing_center=sequencing_center,
        run_date=run_date,
        threads=threads
    )

    if header:
        tools.samtools.SamtoolsTool().reheader(fastqToSamOut, header, outBam)

    return 0


def parser_fastq_to_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inFastq1', help='Input fastq file; 1st end of paired-end reads.')
    parser.add_argument('inFastq2', help='Input fastq file; 2nd end of paired-end reads.')
    parser.add_argument('outBam', help='Output bam file.')
    headerGroup = parser.add_mutually_exclusive_group(required=True)
    headerGroup.add_argument('--sampleName', help='Sample name to insert into the read group header.')
    headerGroup.add_argument('--header', help='Optional text file containing header.')
    parser.add_argument(
        '--JVMmemory',
        default=None,
        help='Deprecated: kept for backwards compatibility, ignored. (Was for Picard)'
    )
    parser.add_argument(
        '--picardOptions',
        default=[],
        nargs='*',
        help='''Read group options in Picard format (OPTIONNAME=value).
                Supported: LIBRARY_NAME, PLATFORM, PLATFORM_UNIT, SEQUENCING_CENTER,
                RUN_DATE, READ_GROUP_NAME. Note that header-related options will be
                overwritten by HEADER if present.'''
    )
    parser.add_argument(
        '--threads',
        type=int,
        default=None,
        help='Number of threads for BAM compression (default: auto)'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, fastq_to_bam, split_args=True)
    return parser


__commands__.append(('fastq_to_bam', parser_fastq_to_bam))


def join_paired_fastq(
        inFastqs,
        output,
        outFormat
):
    ''' Join paired fastq reads into single reads with Ns
    '''
    inFastqs = list(inFastqs)
    if output == '-':
        output = sys.stdout
    SeqIO.write(util.file.join_paired_fastq(inFastqs, output_format=outFormat), output, outFormat)
    return 0


def parser_join_paired_fastq(parser=argparse.ArgumentParser()):
    parser.add_argument('output', help='Output file.')
    parser.add_argument('inFastqs', nargs='+', help='Input fastq file (2 if paired, 1 if interleaved)')
    parser.add_argument('--outFormat', default='fastq', help='Output file format.')
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, join_paired_fastq, split_args=True)
    return parser

__commands__.append(('join_paired_fastq', parser_join_paired_fastq))


# ======================
# ***  split_reads   ***
# ======================
defaultIndexLen = 2
defaultMaxReads = 1000
defaultFormat = 'fastq'


def split_bam(inBam, outBams):
    '''Split BAM file equally into several output BAM files. '''
    samtools = tools.samtools.SamtoolsTool()
    picard = tools.picard.PicardTools()

    # get totalReadCount and maxReads
    # maxReads = totalReadCount / num files, but round up to the nearest
    # even number in order to keep read pairs together (assuming the input
    # is sorted in query order and has no unmated reads, which can be
    # accomplished by Picard RevertSam with SANITIZE=true)
    totalReadCount = samtools.count(inBam)
    maxReads = int(math.ceil(float(totalReadCount) / len(outBams) / 2) * 2)
    log.info("splitting %d reads into %d files of %d reads each", totalReadCount, len(outBams), maxReads)

    # load BAM header into memory
    header = samtools.getHeader(inBam)
    if 'SO:queryname' not in header[0]:
        raise Exception('Input BAM file must be sorted in queryame order')

    # dump to bigsam
    bigsam = mkstempfname('.sam')
    samtools.view(['-@', '3'], inBam, bigsam)

    # split bigsam into little ones
    with util.file.open_or_gzopen(bigsam, 'rt') as inf:
        for outBam in outBams:
            log.info("preparing file " + outBam)
            tmp_sam_reads = mkstempfname('.sam')
            with open(tmp_sam_reads, 'wt') as outf:
                for row in header:
                    outf.write('\t'.join(row) + '\n')
                for _ in range(maxReads):
                    line = inf.readline()
                    if not line:
                        break
                    outf.write(line)
                if outBam == outBams[-1]:
                    for line in inf:
                        outf.write(line)
            picard.execute(
                "SamFormatConverter", [
                    'INPUT=' + tmp_sam_reads, 'OUTPUT=' + outBam, 'VERBOSITY=WARNING'
                ],
                JVMmemory='512m'
            )
            os.unlink(tmp_sam_reads)
    os.unlink(bigsam)


def parser_split_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input BAM file.')
    parser.add_argument('outBams', nargs='+', help='Output BAM files')
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, split_bam, split_args=True)
    return parser


__commands__.append(('split_bam', parser_split_bam))

# =======================
# ***  reheader_bam   ***
# =======================


def parser_reheader_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format.')
    parser.add_argument('rgMap', help='Tabular file containing three columns: field, old, new.')
    parser.add_argument('outBam', help='Output reads, BAM format.')
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_reheader_bam)
    return parser


def main_reheader_bam(args):
    ''' Copy a BAM file (inBam to outBam) while renaming elements of the BAM header.
        The mapping file specifies which (key, old value, new value) mappings. For
        example:
            LB        lib1        lib_one
            SM        sample1        Sample_1
            SM        sample2        Sample_2
            SM        sample3        Sample_3
            CN        broad        BI
    '''
    # read mapping file
    mapper = dict((a + ':' + b, a + ':' + c) for a, b, c in util.file.read_tabfile(args.rgMap))
    # read and convert bam header
    header_file = mkstempfname('.sam')
    with open(header_file, 'wt') as outf:
        for row in tools.samtools.SamtoolsTool().getHeader(args.inBam):
            if row[0] == '@RG':
                row = [mapper.get(x, x) for x in row]
            outf.write('\t'.join(row) + '\n')
    # write new bam with new header
    tools.samtools.SamtoolsTool().reheader(args.inBam, header_file, args.outBam)
    os.unlink(header_file)
    return 0


__commands__.append(('reheader_bam', parser_reheader_bam))


def parser_reheader_bams(parser=argparse.ArgumentParser()):
    parser.add_argument('rgMap', help='Tabular file containing three columns: field, old, new.')
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_reheader_bams)
    return parser


def main_reheader_bams(args):
    ''' Copy BAM files while renaming elements of the BAM header.
        The mapping file specifies which (key, old value, new value) mappings. For
        example:
            LB  lib1  lib_one
            SM  sample1 Sample_1
            SM  sample2 Sample_2
            SM  sample3 Sample_3
            CN  broad   BI
            FN  in1.bam out1.bam
            FN  in2.bam out2.bam
    '''
    # read mapping file
    mapper = dict((a + ':' + b, a + ':' + c) for a, b, c in util.file.read_tabfile(args.rgMap) if a != 'FN')
    files = list((b, c) for a, b, c in util.file.read_tabfile(args.rgMap) if a == 'FN')
    header_file = mkstempfname('.sam')
    # read and convert bam headers
    for inBam, outBam in files:
        if os.path.isfile(inBam):
            with open(header_file, 'wt') as outf:
                for row in tools.samtools.SamtoolsTool().getHeader(inBam):
                    if row[0] == '@RG':
                        row = [mapper.get(x, x) for x in row]
                    outf.write('\t'.join(row) + '\n')
            # write new bam with new header
            tools.samtools.SamtoolsTool().reheader(inBam, header_file, outBam)
    os.unlink(header_file)
    return 0


__commands__.append(('reheader_bams', parser_reheader_bams))

# ============================
# ***  dup_remove_mvicuna  ***
# ============================


class ReadIdStore:
    """
    SQLite-backed store for read IDs with O(1) memory footprint.

    Uses SQLite for disk-based storage with automatic deduplication via
    UNIQUE constraint. Supports random sampling for downsampling operations.
    """

    def __init__(self, db_path):
        """
        Initialize the read ID store.

        Args:
            db_path: Path to SQLite database file (typically in tmpdir)
        """
        self.db_path = db_path
        self.conn = sqlite3.connect(db_path)
        # Performance pragmas (safe since db_path is in tmpdir)
        self.conn.execute('PRAGMA synchronous = OFF')
        self.conn.execute('PRAGMA journal_mode = OFF')
        self.conn.execute('PRAGMA cache_size = -64000')  # 64MB cache (negative = KiB)
        self.conn.execute('PRAGMA temp_store = MEMORY')
        self.conn.execute('''
            CREATE TABLE IF NOT EXISTS read_ids (
                id INTEGER PRIMARY KEY AUTOINCREMENT,
                read_id TEXT UNIQUE NOT NULL
            )
        ''')
        self.conn.commit()

    def add_from_fastq(self, fastq_path):
        """
        Extract read IDs from a FASTQ file and add to the store.

        Handles paired-end reads by stripping /1 or /2 suffixes.
        Duplicates are automatically ignored via INSERT OR IGNORE.

        Args:
            fastq_path: Path to FASTQ file (can be gzipped)

        Returns:
            Number of new read IDs added
        """
        cursor = self.conn.cursor()
        count_before = len(self)

        with util.file.open_or_gzopen(fastq_path, 'rt') as inf:
            line_num = 0
            batch = []
            for line in inf:
                if (line_num % 4) == 0:
                    idVal = line.rstrip('\n')[1:]
                    # Remove /1 or /2 suffix if present (paired-end reads)
                    if idVal.endswith('/1') or idVal.endswith('/2'):
                        idVal = idVal[:-2]
                    batch.append((idVal,))

                    # Batch inserts for performance
                    if len(batch) >= 10000:
                        cursor.executemany(
                            'INSERT OR IGNORE INTO read_ids (read_id) VALUES (?)',
                            batch
                        )
                        batch = []
                line_num += 1

            # Insert remaining batch
            if batch:
                cursor.executemany(
                    'INSERT OR IGNORE INTO read_ids (read_id) VALUES (?)',
                    batch
                )

        self.conn.commit()
        return len(self) - count_before

    def add_from_readlist(self, readlist_path):
        """
        Add read IDs from a text file (one ID per line).

        Args:
            readlist_path: Path to text file with one read ID per line.
                           Empty lines and whitespace are stripped.
        """
        with util.file.open_or_gzopen(readlist_path, 'rt') as f:
            self.extend(line.strip() for line in f if line.strip())

    def __len__(self):
        """Return the number of unique read IDs in the store."""
        cursor = self.conn.execute('SELECT COUNT(*) FROM read_ids')
        return cursor.fetchone()[0]

    def __iter__(self):
        """Iterate over all read IDs in insertion order. O(1) memory."""
        cursor = self.conn.execute('SELECT read_id FROM read_ids ORDER BY id')
        for row in cursor:
            yield row[0]

    def __contains__(self, read_id):
        """Check if a read ID exists in the store."""
        cursor = self.conn.execute(
            'SELECT 1 FROM read_ids WHERE read_id = ? LIMIT 1',
            (read_id,)
        )
        return cursor.fetchone() is not None

    def contains_batch(self, read_ids, batch_size=10000):
        """
        Check multiple read IDs for membership in batched queries.

        Args:
            read_ids: List of read ID strings to check
            batch_size: Number of IDs to query at once (default 10000)

        Returns:
            Set of read IDs that exist in the store
        """
        found = set()
        batch = []
        for read_id in read_ids:
            batch.append(read_id)
            if len(batch) >= batch_size:
                placeholders = ','.join('?' * len(batch))
                cursor = self.conn.execute(
                    f'SELECT read_id FROM read_ids WHERE read_id IN ({placeholders})',
                    batch
                )
                found.update(row[0] for row in cursor)
                batch = []
        if batch:
            placeholders = ','.join('?' * len(batch))
            cursor = self.conn.execute(
                f'SELECT read_id FROM read_ids WHERE read_id IN ({placeholders})',
                batch
            )
            found.update(row[0] for row in cursor)
        return found

    def add(self, read_id):
        """
        Add a single read ID to the store.

        Args:
            read_id: The read ID string to add

        Returns:
            True if the read ID was added, False if it already existed
        """
        cursor = self.conn.execute(
            'INSERT OR IGNORE INTO read_ids (read_id) VALUES (?)',
            (read_id,)
        )
        self.conn.commit()
        return cursor.rowcount > 0

    def extend(self, read_ids):
        """
        Add multiple read IDs to the store efficiently.

        Uses executemany for batch insertion. O(1) memory - processes
        the input iterable in chunks.

        Args:
            read_ids: Iterable of read ID strings

        Returns:
            Number of new read IDs added
        """
        count_before = len(self)
        cursor = self.conn.cursor()

        # Process in batches for memory efficiency
        batch = []
        for read_id in read_ids:
            batch.append((read_id,))
            if len(batch) >= 10000:
                cursor.executemany(
                    'INSERT OR IGNORE INTO read_ids (read_id) VALUES (?)',
                    batch
                )
                batch = []

        # Insert remaining batch
        if batch:
            cursor.executemany(
                'INSERT OR IGNORE INTO read_ids (read_id) VALUES (?)',
                batch
            )

        self.conn.commit()
        return len(self) - count_before

    def __delitem__(self, read_id):
        """
        Delete a read ID from the store.

        Args:
            read_id: The read ID string to delete

        Raises:
            KeyError: If the read ID does not exist
        """
        cursor = self.conn.execute(
            'DELETE FROM read_ids WHERE read_id = ?',
            (read_id,)
        )
        self.conn.commit()
        if cursor.rowcount == 0:
            raise KeyError(read_id)

    def discard(self, read_id):
        """
        Remove a read ID if it exists, without raising an error if absent.

        Args:
            read_id: The read ID string to remove
        """
        self.conn.execute(
            'DELETE FROM read_ids WHERE read_id = ?',
            (read_id,)
        )
        self.conn.commit()

    def write_to_file(self, out_path, max_reads=None):
        """
        Write read IDs to a file.

        Args:
            out_path: Output file path
            max_reads: If specified, randomly sample this many read IDs

        Returns:
            Number of read IDs written
        """
        if max_reads is not None and max_reads < len(self):
            # Random sampling using ORDER BY RANDOM()
            query = 'SELECT read_id FROM read_ids ORDER BY RANDOM() LIMIT ?'
            cursor = self.conn.execute(query, (max_reads,))
        else:
            # Return all in insertion order
            query = 'SELECT read_id FROM read_ids ORDER BY id'
            cursor = self.conn.execute(query)

        count = 0
        with open(out_path, 'wt') as outf:
            for row in cursor:
                outf.write(row[0] + '\n')
                count += 1

        return count

    def shrink_to_subsample(self, n):
        """
        Destructively shrink the store to a random subsample of n elements.

        This modifies the store in-place, keeping exactly n randomly selected
        read IDs and deleting all others. If the store has n or fewer elements,
        no changes are made.

        Args:
            n: Number of read IDs to keep

        Note:
            Uses SQLite's ORDER BY RANDOM() for random selection. The operation
            is O(N log N) where N is the current store size due to random sorting.
        """
        current_size = len(self)
        if current_size <= n:
            return  # Nothing to do

        # Delete all rows except n randomly selected ones
        self.conn.execute('''
            DELETE FROM read_ids
            WHERE id NOT IN (
                SELECT id FROM read_ids
                ORDER BY RANDOM()
                LIMIT ?
            )
        ''', (n,))
        self.conn.commit()

    def filter_bam_by_ids(self, inBam, outBam, include=True):
        """
        Filter a BAM file to reads matching (or not matching) IDs in this store.

        This is a faster, O(1) memory replacement for picard.FilterSamReadsTool.
        Headers are fully preserved from input to output.

        Args:
            inBam: Input BAM file path
            outBam: Output BAM file path
            include: If True, keep only reads with IDs in the store (inclusion list).
                     If False, keep only reads with IDs NOT in the store (exclusion list).

        Note:
            For paired-end reads, both mates share the same QNAME so both will be
            kept or excluded together based on whether the QNAME is in the store.
        """
        # Handle empty input BAM
        samtools = tools.samtools.SamtoolsTool()
        if samtools.isEmpty(inBam):
            shutil.copyfile(inBam, outBam)
            return

        # Handle empty read ID store
        if len(self) == 0:
            if include:
                # Include mode with empty list = empty output (keep header only)
                with pysam.AlignmentFile(inBam, 'rb', check_sq=False) as inb:
                    with pysam.AlignmentFile(outBam, 'wb', header=inb.header) as outf:
                        pass  # Write header only, no reads
            else:
                # Exclude mode with empty list = keep all reads
                shutil.copyfile(inBam, outBam)
            return

        # Hybrid approach: samtools subprocess pipes + batched SQLite queries
        start = time.time()
        read_count = 0
        write_count = 0
        lookup_time = 0.0
        batch_size = 50000

        samtools_path = samtools.install_and_get_path()

        # Read: samtools view -h --no-PG input.bam -> SAM text (with header, no PG added)
        # Write: samtools view -b --no-PG -o output.bam - (converts SAM to BAM, no PG added)
        read_proc = subprocess.Popen(
            [samtools_path, 'view', '-h', '--no-PG', inBam],
            stdout=subprocess.PIPE)
        write_proc = subprocess.Popen(
            [samtools_path, 'view', '-b', '-@', '2', '--no-PG', '-o', outBam, '-'],
            stdin=subprocess.PIPE)

        try:
            batch = []  # [(line, qname), ...]
            for line in read_proc.stdout:
                # Header lines (start with @) are passed through unchanged
                if line.startswith(b'@'):
                    write_proc.stdin.write(line)
                    continue

                read_count += 1

                # Parse only QNAME (first field before tab)
                tab_pos = line.find(b'\t')
                qname = line[:tab_pos].decode()
                batch.append((line, qname))

                if len(batch) >= batch_size:
                    t0 = time.time()
                    found = self.contains_batch([q for _, q in batch])
                    lookup_time += time.time() - t0

                    for ln, q in batch:
                        if (include and q in found) or (not include and q not in found):
                            write_proc.stdin.write(ln)
                            write_count += 1
                    batch = []

            # Process remaining batch
            if batch:
                t0 = time.time()
                found = self.contains_batch([q for _, q in batch])
                lookup_time += time.time() - t0

                for ln, q in batch:
                    if (include and q in found) or (not include and q not in found):
                        write_proc.stdin.write(ln)
                        write_count += 1
        finally:
            write_proc.stdin.close()
            write_proc.wait()
            read_proc.wait()

        elapsed = time.time() - start
        log.debug(f"PERF: filter_time={elapsed:.2f}s lookup_time={lookup_time:.2f}s "
                  f"reads={read_count} written={write_count} batch_size={batch_size}")

    def close(self):
        """Close the database connection."""
        self.conn.close()

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        self.close()
        return False


def mvicuna_fastqs_to_readlist(inFastq1, inFastq2, readList):
    # Run M-Vicuna on FASTQ files
    outFastq1 = mkstempfname('.1.fastq')
    outFastq2 = mkstempfname('.2.fastq')
    if inFastq2 is None or os.path.getsize(inFastq2) < 10:
        tools.mvicuna.MvicunaTool().rmdup_single(inFastq1, outFastq1)
    else:
        tools.mvicuna.MvicunaTool().rmdup((inFastq1, inFastq2), (outFastq1, outFastq2), None)

    # Make a list of reads to keep
    with open(readList, 'at') as outf:
        for fq in (outFastq1, outFastq2):
            with util.file.open_or_gzopen(fq, 'rt') as inf:
                line_num = 0
                for line in inf:
                    if (line_num % 4) == 0:
                        idVal = line.rstrip('\n')[1:]
                        if idVal.endswith('/1'):
                            outf.write(idVal[:-2] + '\n')
                        # single-end reads do not have /1 /2 mate suffix
                        # so pass through their IDs
                        if not (idVal.endswith('/1') or idVal.endswith('/2')):
                            outf.write(idVal + '\n')
                    line_num += 1
    os.unlink(outFastq1)
    os.unlink(outFastq2)


def rmdup_cdhit_bam(inBam, outBam, max_mismatches=None, jvm_memory=None):
    ''' Remove duplicate reads from BAM file using cd-hit-dup.
    '''
    max_mismatches = max_mismatches or 4
    tmp_dir = tempfile.mkdtemp()

    tools.picard.SplitSamByLibraryTool().execute(inBam, tmp_dir)

    s2fq_tool = tools.picard.SamToFastqTool()
    cdhit = tools.cdhit.CdHit()
    out_bams = []
    for f in os.listdir(tmp_dir):
        out_bam = mkstempfname('.bam')
        out_bams.append(out_bam)
        library_sam = os.path.join(tmp_dir, f)

        in_fastqs = mkstempfname('.1.fastq'), mkstempfname('.2.fastq')

        s2fq_tool.execute(library_sam, in_fastqs[0], in_fastqs[1])
        if not os.path.getsize(in_fastqs[0]) > 0 and not os.path.getsize(in_fastqs[1]) > 0:
            continue

        out_fastqs = mkstempfname('.1.fastq'), mkstempfname('.2.fastq')
        options = {
            '-e': max_mismatches,
        }
        if in_fastqs[1] is not None and os.path.getsize(in_fastqs[1]) > 10:
            options['-i2'] = in_fastqs[1]
            options['-o2'] = out_fastqs[1]

        log.info("executing cd-hit-dup on library " + library_sam)
        # cd-hit-dup cannot operate on piped fastq input because it reads twice
        # Run cd-hit-dup synchronously (not in background) to ensure output files are complete
        # before FastqToSamTool tries to read them
        cdhit.execute('cd-hit-dup', in_fastqs[0], out_fastqs[0], options=options, background=False)

        tools.samtools.SamtoolsTool().import_fastq(
            out_fastqs[0], out_fastqs[1], out_bam,
            sample_name=f,
        )
        for fn in in_fastqs:
            os.unlink(fn)

    with util.file.fifo(name='merged.sam') as merged_bam:
        merge_opts = ['SORT_ORDER=queryname']
        tools.picard.MergeSamFilesTool().execute(out_bams, merged_bam, picardOptions=merge_opts, JVMmemory=jvm_memory, background=True)
        tools.picard.ReplaceSamHeaderTool().execute(merged_bam, inBam, outBam, JVMmemory=jvm_memory)


def parser_rmdup_cdhit_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format.')
    parser.add_argument('outBam', help='Output reads, BAM format.')
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.FilterSamReadsTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)',
        dest='jvm_memory'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, rmdup_cdhit_bam, split_args=True)
    return parser


__commands__.append(('rmdup_cdhit_bam', parser_rmdup_cdhit_bam))

def _merge_fastqs_and_mvicuna(lb, files):
    readList = mkstempfname('.keep_reads.txt')
    log.info("executing M-Vicuna DupRm on library " + lb)

    # create merged FASTQs per library
    infastqs = (mkstempfname('.1.fastq'), mkstempfname('.2.fastq'))
    for d in range(2):
        with open(infastqs[d], 'wt') as outf:
            for fprefix in files:
                fn = '%s_%d.fastq' % (fprefix, d + 1)

                if os.path.isfile(fn):
                    with open(fn, 'rt') as inf:
                        for line in inf:
                            outf.write(line)
                    os.unlink(fn)
                else:
                    log.warning(
                        """no reads found in %s,
                                assuming that's because there's no reads in that read group""", fn
                    )

    # M-Vicuna DupRm to see what we should keep (append IDs to running file)
    if os.path.getsize(infastqs[0]) > 0 or os.path.getsize(infastqs[1]) > 0:
        mvicuna_fastqs_to_readlist(infastqs[0], infastqs[1], readList)
    for fn in infastqs:
        os.unlink(fn)

    return readList

def rmdup_mvicuna_bam(inBam, outBam, threads=None):
    ''' Remove duplicate reads from BAM file using M-Vicuna. The
        primary advantage to this approach over Picard's MarkDuplicates tool
        is that Picard requires that input reads are aligned to a reference,
        and M-Vicuna can operate on unaligned reads.
    '''

    # Convert BAM -> FASTQ pairs per read group and load all read groups
    tempDir = tempfile.mkdtemp()
    tools.picard.SamToFastqTool().per_read_group(inBam, tempDir, picardOptions=['VALIDATION_STRINGENCY=LENIENT'])
    read_groups = [x[1:] for x in tools.samtools.SamtoolsTool().getHeader(inBam) if x[0] == '@RG']
    read_groups = [dict(pair.split(':', 1) for pair in rg) for rg in read_groups]

    # Collect FASTQ pairs for each library
    lb_to_files = {}
    for rg in read_groups:
        lb_to_files.setdefault(rg.get('LB', 'none'), set())
        fname = rg['ID']
        lb_to_files[rg.get('LB', 'none')].add(os.path.join(tempDir, fname))
    log.info("found %d distinct libraries and %d read groups", len(lb_to_files), len(read_groups))

    # Create ReadIdStore and collect read IDs from all libraries
    with util.file.tmp_dir(suffix='_mvicuna_filter') as filter_tmpdir:
        db_path = os.path.join(filter_tmpdir, 'read_ids.db')
        with ReadIdStore(db_path) as store:
            # For each library, merge FASTQs and run M-Vicuna, collecting read IDs
            with concurrent.futures.ProcessPoolExecutor(
                    max_workers=threads or util.misc.available_cpu_count()) as executor:
                futures = [executor.submit(_merge_fastqs_and_mvicuna, lb, files)
                           for lb, files in lb_to_files.items()]
                for future in concurrent.futures.as_completed(futures):
                    log.info("mvicuna finished processing library")
                    try:
                        readList = future.result()
                        # Stream read IDs directly into store (no intermediate concat)
                        with open(readList, 'rt') as f:
                            store.extend(line.strip() for line in f if line.strip())
                        os.unlink(readList)  # Clean up per-library file immediately
                    except Exception as exc:
                        log.error('mvicuna process call generated an exception: %s' % (exc))
                        raise

            # Filter original input BAM against keep-list
            store.filter_bam_by_ids(inBam, outBam, include=True)

    return 0


def parser_rmdup_mvicuna_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format.')
    parser.add_argument('outBam', help='Output reads, BAM format.')
    util.cmd.common_args(parser, (('threads',None), ('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, rmdup_mvicuna_bam, split_args=True)
    return parser


__commands__.append(('rmdup_mvicuna_bam', parser_rmdup_mvicuna_bam))


def rmdup_bbnorm_bam(inBam, outBam,
                      target=None, k=None, passes=None,
                      memory=None, threads=None,
                      min_input_reads=None, max_output_reads=None):
    """
    Remove duplicate/normalize reads from BAM file using BBNorm.

    Convert BAM to interleaved FASTQ, run bbnorm, extract kept read IDs,
    and filter original BAM to those IDs using pysam (O(1) memory).

    Args:
        inBam: Input BAM file
        outBam: Output BAM file
        target: BBNorm target normalization depth (default: bbnorm default of 100)
        k: Kmer length (default: bbnorm default of 31)
        passes: Number of passes (default: bbnorm default of 2)
        memory: Java memory for bbnorm (e.g., "4g")
        threads: Number of threads for bbnorm
        min_input_reads: Skip processing if input has fewer reads (copy input to output)
        max_output_reads: Randomly downsample keep-list if larger than this
    """
    samtools = tools.samtools.SamtoolsTool()

    # Count input reads
    input_read_count = samtools.count(inBam)
    log.info("Input BAM has %d reads", input_read_count)

    # Skip processing if empty or below min_input_reads threshold
    min_threshold = min_input_reads if min_input_reads is not None else 1
    if input_read_count < min_threshold:
        if min_input_reads is not None:
            log.info("Input read count %d below min_input_reads %d, copying input to output",
                     input_read_count, min_input_reads)
        else:
            log.info("Input BAM is empty, copying to output")
        shutil.copyfile(inBam, outBam)
        return 0

    with util.file.tmp_dir(suffix='_bbnorm') as tmpdir:
        # Convert BAM to interleaved FASTQ using samtools bam2fq
        inFastq = os.path.join(tmpdir, 'input.fastq')
        samtools.bam2fq(inBam, inFastq)  # Single file = interleaved output

        # Run bbnorm
        outFastq = os.path.join(tmpdir, 'output.fastq')
        tools.bbmap.BBMapTool().bbnorm(
            inFastq, outFastq,
            tmpdir=tmpdir, target=target, k=k, passes=passes,
            mindepth=0, threads=threads, memory=memory
        )

        # Extract read IDs from bbnorm output FASTQ using SQLite for O(1) memory
        db_path = os.path.join(tmpdir, 'read_ids.db')

        with ReadIdStore(db_path) as store:
            store.add_from_fastq(outFastq)
            num_ids = len(store)
            log.info("BBNorm retained %d read IDs", num_ids)

            # Downsample if needed (modifies store in-place via SQL)
            if max_output_reads is not None and num_ids > max_output_reads:
                log.info("Downsampling from %d to %d read IDs", num_ids, max_output_reads)
                store.shrink_to_subsample(max_output_reads)

            # Filter original BAM directly using ReadIdStore
            store.filter_bam_by_ids(inBam, outBam, include=True)

    # Count output reads
    output_read_count = samtools.count(outBam)
    log.info("Output BAM has %d reads (%.1f%% of input)",
             output_read_count, 100.0 * output_read_count / max(input_read_count, 1))

    return 0


def parser_rmdup_bbnorm_bam(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format.')
    parser.add_argument('outBam', help='Output reads, BAM format.')
    parser.add_argument(
        '--target',
        type=int,
        default=None,
        help='BBNorm target normalization depth (default: bbnorm default of 100)'
    )
    parser.add_argument(
        '--kmerLength',
        dest='k',
        type=int,
        default=None,
        help='Kmer length for bbnorm (default: bbnorm default of 31)'
    )
    parser.add_argument(
        '--passes',
        type=int,
        default=None,
        help='Number of bbnorm passes (default: bbnorm default of 2)'
    )
    parser.add_argument(
        '--memory',
        default=None,
        help='Java memory for bbnorm (e.g., "4g", "8g")'
    )
    parser.add_argument(
        '--minInputReads',
        dest='min_input_reads',
        type=int,
        default=None,
        help='Skip processing if input has fewer than this many reads'
    )
    parser.add_argument(
        '--maxOutputReads',
        dest='max_output_reads',
        type=int,
        default=None,
        help='Randomly downsample output to at most this many read IDs'
    )
    util.cmd.common_args(parser, (('threads', None), ('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, rmdup_bbnorm_bam, split_args=True)
    return parser


__commands__.append(('rmdup_bbnorm_bam', parser_rmdup_bbnorm_bam))


def parser_rmdup_prinseq_fastq(parser=argparse.ArgumentParser()):
    parser.add_argument('inFastq1', help='Input fastq file; 1st end of paired-end reads.')
    parser.add_argument('inFastq2', help='Input fastq file; 2nd end of paired-end reads.')
    parser.add_argument('outFastq1', help='Output fastq file; 1st end of paired-end reads.')
    parser.add_argument('outFastq2', help='Output fastq file; 2nd end of paired-end reads.')
    parser.add_argument(
        "--includeUnmated",
        help="Include unmated reads in the main output fastq files (default: %(default)s)",
        default=False,
        action="store_true"
    )
    parser.add_argument('--unpairedOutFastq1', default=None, help='File name of output unpaired reads from 1st end of paired-end reads (independent of --includeUnmated)')
    parser.add_argument('--unpairedOutFastq2', default=None, help='File name of output unpaired reads from 2nd end of paired-end reads (independent of --includeUnmated)')
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_rmdup_prinseq_fastq)
    return parser


def main_rmdup_prinseq_fastq(args):
    ''' Run prinseq-lite's duplicate removal operation on paired-end
        reads.  Also removes reads with more than one N.
    '''
    prinseq = tools.prinseq.PrinseqTool()
    prinseq.rmdup_fastq_paired(args.inFastq1, args.inFastq2, args.outFastq1, args.outFastq2, args.includeUnmated, args.unpairedOutFastq1, args.unpairedOutFastq2)
    return 0


__commands__.append(('rmdup_prinseq_fastq', parser_rmdup_prinseq_fastq))


def filter_bam_mapped_only(inBam, outBam):
    ''' Samtools to reduce a BAM file to only reads that are
        aligned (-F 4) with a non-zero mapping quality (-q 1)
        and are not marked as a PCR/optical duplicate (-F 1024).
    '''
    tools.samtools.SamtoolsTool().view(['-b', '-q', '1', '-F', '1028'], inBam, outBam)
    tools.picard.BuildBamIndexTool().execute(outBam)
    return 0


def parser_filter_bam_mapped_only(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input aligned reads, BAM format.')
    parser.add_argument('outBam', help='Output sorted indexed reads, filtered to aligned-only, BAM format.')
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, filter_bam_mapped_only, split_args=True)
    return parser


__commands__.append(('filter_bam_mapped_only', parser_filter_bam_mapped_only))

# ======= Novoalign ========


def parser_novoalign(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format.')
    parser.add_argument('refFasta', help='Reference genome, FASTA format, pre-indexed by Novoindex.')
    parser.add_argument('outBam', help='Output reads, BAM format (aligned).')
    parser.add_argument('--options', default='-r Random', help='Novoalign options (default: %(default)s)')
    parser.add_argument('--min_qual', default=0, help='Filter outBam to minimum mapping quality (default: %(default)s)')
    parser.add_argument(
        '--JVMmemory',
        default=tools.picard.SortSamTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--NOVOALIGN_LICENSE_PATH',
        default=None,
        dest="novoalign_license_path",
        help='A path to the novoalign.lic file. This overrides the NOVOALIGN_LICENSE_PATH environment variable. (default: %(default)s)'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_novoalign)
    return parser


def main_novoalign(args):
    '''Align reads with Novoalign. Sort and index BAM output.'''
    tools.novoalign.NovoalignTool(license_path=args.novoalign_license_path).execute(
        args.inBam,
        args.refFasta,
        args.outBam,
        options=args.options.split(),
        min_qual=args.min_qual,
        JVMmemory=args.JVMmemory
    )
    return 0


__commands__.append(('novoalign', parser_novoalign))


def parser_novoindex(parser=argparse.ArgumentParser()):
    parser.add_argument('refFasta', help='Reference genome, FASTA format.')
    parser.add_argument(
        '--NOVOALIGN_LICENSE_PATH',
        default=None,
        dest="novoalign_license_path",
        help='A path to the novoalign.lic file. This overrides the NOVOALIGN_LICENSE_PATH environment variable. (default: %(default)s)'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None)))
    util.cmd.attach_main(parser, main_novoindex)
    return parser

def main_novoindex(args):
    tools.novoalign.NovoalignTool(license_path=args.novoalign_license_path).index_fasta(args.refFasta)
    return 0


__commands__.append(('novoindex', parser_novoindex))

# ========= GATK ==========


def parser_gatk_ug(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format.')
    parser.add_argument('refFasta', help='Reference genome, FASTA format, pre-indexed by Picard.')
    parser.add_argument(
        'outVcf',
        help='''Output calls in VCF format. If this filename ends with .gz,
        GATK will BGZIP compress the output and produce a Tabix index file as well.'''
    )
    parser.add_argument(
        '--options',
        default='--min_base_quality_score 15 -ploidy 4',
        help='UnifiedGenotyper options (default: %(default)s)'
    )
    parser.add_argument(
        '--JVMmemory',
        default=tools.gatk.GATKTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--GATK_PATH',
        default=None,
        help='A path containing the GATK jar file. This overrides the GATK_ENV environment variable or the GATK conda package. (default: %(default)s)'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_gatk_ug)
    return parser


def main_gatk_ug(args):
    '''Call genotypes using the GATK UnifiedGenotyper.'''
    tools.gatk.GATKTool(path=args.GATK_PATH).ug(
        args.inBam, args.refFasta,
        args.outVcf, options=args.options.split(),
        JVMmemory=args.JVMmemory
    )
    return 0


__commands__.append(('gatk_ug', parser_gatk_ug))


def parser_gatk_realign(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input reads, BAM format, aligned to refFasta.')
    parser.add_argument('refFasta', help='Reference genome, FASTA format, pre-indexed by Picard.')
    parser.add_argument('outBam', help='Realigned reads.')
    parser.add_argument(
        '--JVMmemory',
        default=tools.gatk.GATKTool.jvmMemDefault,
        help='JVM virtual memory size (default: %(default)s)'
    )
    parser.add_argument(
        '--GATK_PATH',
        default=None,
        help='A path containing the GATK jar file. This overrides the GATK_ENV environment variable or the GATK conda package. (default: %(default)s)'
    )
    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_gatk_realign)
    parser.add_argument('--threads', type=int, default=None, help='Number of threads (default: all available cores)')
    return parser


def main_gatk_realign(args):
    '''Local realignment of BAM files with GATK IndelRealigner.'''
    tools.gatk.GATKTool(path=args.GATK_PATH).local_realign(
        args.inBam, args.refFasta,
        args.outBam, JVMmemory=args.JVMmemory,
        threads=args.threads,
    )
    return 0


__commands__.append(('gatk_realign', parser_gatk_realign))

# =========================


def align_and_fix(
    inBam, refFasta,
    outBamAll=None,
    outBamFiltered=None,
    aligner_options='',
    aligner="novoalign",
    bwa_min_score=None,
    novoalign_amplicons_bed=None,
    amplicon_window=4,
    JVMmemory=None,
    threads=None,
    skip_mark_dupes=False,
    skip_realign=False,
    dup_marker='sambamba',
    gatk_path=None,
    novoalign_license_path=None
):
    ''' Take reads, align to reference with Novoalign, minimap2, or BWA-MEM.
        Optionally mark duplicates with Picard or sambamba, optionally realign
        indels with GATK, and optionally filter final file to mapped/non-dupe reads.
    '''
    if not (outBamAll or outBamFiltered):
        log.warning("are you sure you meant to do nothing?")
        return

    assert aligner in ["novoalign", "bwa", "minimap2"]
    samtools = tools.samtools.SamtoolsTool()

    refFastaCopy = mkstempfname('.ref_copy.fasta')
    with util.file.fastas_with_sanitized_ids(refFasta, use_tmp=True) as sanitized_fastas:
        shutil.copyfile(sanitized_fastas[0], refFastaCopy)

    tools.picard.CreateSequenceDictionaryTool().execute(refFastaCopy, overwrite=True, JVMmemory=JVMmemory)
    samtools.faidx(refFastaCopy, overwrite=True)

    if aligner_options is None:
        if aligner=="novoalign":
            aligner_options = '-r Random'
        else:
            aligner_options = '' # use defaults

    if samtools.isEmpty(inBam):
        log.warning("zero reads present in input")

    bam_aligned = mkstempfname('.aligned.bam')
    if aligner=="novoalign":
        if novoalign_amplicons_bed is not None:
            aligner_options += ' --amplicons {} {}'.format(novoalign_amplicons_bed, amplicon_window)

        tools.novoalign.NovoalignTool(license_path=novoalign_license_path).index_fasta(refFastaCopy)

        tools.novoalign.NovoalignTool(license_path=novoalign_license_path).execute(
            inBam, refFastaCopy, bam_aligned,
            options=aligner_options.split(),
            JVMmemory=JVMmemory
        )

    elif aligner=='bwa':
        bwa = tools.bwa.Bwa()
        bwa.index(refFastaCopy)

        opts = aligner_options.split()

        bwa.align_mem_bam(inBam, refFastaCopy, bam_aligned, min_score_to_filter=bwa_min_score, threads=threads, options=opts)

    elif aligner=='minimap2':
        mm2 = tools.minimap2.Minimap2()
        mm2.align_bam(inBam, refFastaCopy, bam_aligned, threads=threads, options=aligner_options.split())

    # Mark duplicates with sambamba (default) or Picard
    if skip_mark_dupes:
        bam_marked = bam_aligned
    else:
        bam_marked = mkstempfname('.mkdup.bam')
        if dup_marker == 'sambamba':
            tools.sambamba.SambambaTool().markdup(
                bam_aligned, bam_marked, threads=threads
            )
        elif dup_marker == 'picard':
            tools.picard.MarkDuplicatesTool().execute(
                [bam_aligned], bam_marked, picardOptions=['CREATE_INDEX=true'],
                JVMmemory=JVMmemory
            )
        else:
            raise ValueError(f"Unknown dup_marker: {dup_marker}. Must be 'sambamba' or 'picard'")
        os.unlink(bam_aligned)

    # Index the marked BAM
    if dup_marker == 'sambamba' and not skip_mark_dupes:
        tools.sambamba.SambambaTool().index(bam_marked, threads=threads)
    else:
        samtools.index(bam_marked)

    # Local realignment with GATK (optional)
    if skip_realign or samtools.isEmpty(bam_marked):
        bam_realigned = bam_marked
    else:
        bam_realigned = mkstempfname('.realigned.bam')
        tools.gatk.GATKTool(path=gatk_path).local_realign(bam_marked, refFastaCopy, bam_realigned, JVMmemory=JVMmemory, threads=threads)
        os.unlink(bam_marked)

    # Generate output files
    if outBamAll:
        shutil.copyfile(bam_realigned, outBamAll)
        if dup_marker == 'sambamba':
            tools.sambamba.SambambaTool().index(outBamAll, threads=threads)
        else:
            tools.picard.BuildBamIndexTool().execute(outBamAll, JVMmemory=JVMmemory)
    if outBamFiltered:
        filtered_any_mapq = mkstempfname('.filtered_any_mapq.bam')
        # filter based on read flags
        samtools.filter_to_proper_primary_mapped_reads(bam_realigned, filtered_any_mapq)
        # remove reads with MAPQ <1
        samtools.view(['-b', '-q', '1'], filtered_any_mapq, outBamFiltered)
        os.unlink(filtered_any_mapq)
        if dup_marker == 'sambamba':
            tools.sambamba.SambambaTool().index(outBamFiltered, threads=threads)
        else:
            tools.picard.BuildBamIndexTool().execute(outBamFiltered, JVMmemory=JVMmemory)
    os.unlink(bam_realigned)


def parser_align_and_fix(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input unaligned reads, BAM format.')
    parser.add_argument('refFasta', help='Reference genome, FASTA format; will be indexed by Picard and Novoalign.')
    parser.add_argument(
        '--outBamAll',
        default=None,
        help='''Aligned, sorted, and indexed reads.  Unmapped and duplicate reads are
                retained. By default, duplicate reads are marked. If "--skipMarkDupes"
                is specified duplicate reads are included in outout without being marked.'''
    )
    parser.add_argument(
        '--outBamFiltered',
        default=None,
        help='''Aligned, sorted, and indexed reads.  Unmapped reads are removed from this file,
                as well as any marked duplicate reads. Note that if "--skipMarkDupes" is provided,
                duplicates will be not be marked and will be included in the output.'''
    )
    parser.add_argument('--aligner_options', default=None, help='aligner options (default for novoalign: "-r Random", bwa: "-T 30"')
    parser.add_argument('--aligner', choices=['novoalign', 'minimap2', 'bwa'], default='novoalign', help='aligner (default: %(default)s)')
    parser.add_argument('--bwa_min_score', type=int, default=None, help='BWA mem on paired reads ignores the -T parameter. Set a value here (e.g. 30) to invoke a custom post-alignment filter (default: no filtration)')
    parser.add_argument('--novoalign_amplicons_bed', default=None, help='Novoalign only: amplicon primer file (BED format) to soft clip')
    parser.add_argument('--amplicon_window', type=int, default=4, help='Novoalign only: amplicon primer window size (default: %(default)s)')
    parser.add_argument('--JVMmemory', default='4g', help='JVM virtual memory size (default: %(default)s)')
    parser.add_argument('--threads', type=int, default=None, help='Number of threads (default: all available cores)')
    parser.add_argument('--skipMarkDupes',
                        help='If specified, duplicate reads will not be marked in the resulting output file.',
                        dest="skip_mark_dupes",
                        action='store_true')
    parser.add_argument('--skipRealign',
                        help='If specified, GATK local realignment will be skipped. '
                             'Recommended for viral genomes where indel realignment provides minimal benefit.',
                        dest="skip_realign",
                        action='store_true')
    parser.add_argument('--dupMarker',
                        choices=['sambamba', 'picard'],
                        default='sambamba',
                        dest="dup_marker",
                        help='Tool to use for marking duplicates. Sambamba is multi-threaded and faster. '
                             'Picard is the legacy option. (default: %(default)s)')
    parser.add_argument(
        '--GATK_PATH',
        default=None,
        dest="gatk_path",
        help='A path containing the GATK jar file. This overrides the GATK_ENV environment variable or the GATK conda package. (default: %(default)s)'
    )
    parser.add_argument(
        '--NOVOALIGN_LICENSE_PATH',
        default=None,
        dest="novoalign_license_path",
        help='A path to the novoalign.lic file. This overrides the NOVOALIGN_LICENSE_PATH environment variable. (default: %(default)s)'
    )

    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, align_and_fix, split_args=True)
    return parser


__commands__.append(('align_and_fix', parser_align_and_fix))

# =========================

def filter_bam_to_proper_primary_mapped_reads(inBam, outBam, doNotRequirePairsToBeProper=False, keepSingletons=False, keepDuplicates=False):
    ''' Take a BAM file and filter to only reads that are properly
        paired and mapped. Optionally reject singletons, and 
        optionally require reads to be properly paired and mapped.

        Output includes reads that are:
            - not flagged as duplicates
            - not secondary or supplementary (split/chimeric reads)
            - For paired-end reads:
            -   marked as proper pair (if require_pairs_to_be_proper=True) OR
                both not unmapped (if require_pairs_to_be_proper=False) OR
                not a member of a pair with a singleton (if reject_singletons=True)
            - For single-end reads:
                mapped
    '''
    samtools = tools.samtools.SamtoolsTool()
    samtools.filter_to_proper_primary_mapped_reads(inBam, 
                                                   outBam,
                                                   require_pairs_to_be_proper=not doNotRequirePairsToBeProper,
                                                   reject_singletons=not keepSingletons,
                                                   reject_duplicates=not keepDuplicates)
    return 0

def parser_filter_bam_to_proper_primary_mapped_reads(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input aligned reads, BAM format.')
    parser.add_argument('outBam', help='Output reads, BAM format.')
    parser.add_argument(
        '--doNotRequirePairsToBeProper',
        help='Do not require reads to be properly paired when filtering (default: %(default)s)',
        action='store_true'
    )
    parser.add_argument(
        '--keepSingletons',
        help='Keep singleton reads when filtering (default: %(default)s)',
        action='store_true'
    )
    parser.add_argument(
        '--keepDuplicates',
        help='When filtering, do not exclude reads due to being flagged as duplicates (default: %(default)s)',
        action='store_true'
    )

    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, filter_bam_to_proper_primary_mapped_reads, split_args=True)
    return parser

# =========================


def minimap2_idxstats(inBam, refFasta, outStats, outReadlist=None, threads=None):
    ''' Align reads to reference with minimap2 and produce idxstats-like counts.

        This uses PAF output format (no SAM/BAM generation) for faster alignment
        and streams the output directly without intermediate files.

        Args:
            inBam: Input reads (BAM format)
            refFasta: Reference genome (FASTA format)
            outStats: Output file in samtools idxstats format
            outReadlist: Optional output file with read IDs that mapped (or None to skip)
            threads: Number of threads for alignment (default: auto-detect)
    '''
    mm2 = tools.minimap2.Minimap2()
    mm2.idxstats(inBam, refFasta, outStats, outReadlist=outReadlist, threads=threads)

    

def parser_minimap2_idxstats(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input unaligned reads, BAM format.')
    parser.add_argument('refFasta', help='Reference genome, FASTA format.')
    parser.add_argument('outStats', help='Output idxstats file (tab-separated: ref_name, ref_length, mapped_count, 0).')
    parser.add_argument('--outReadlist', help='Optional output file listing read IDs that mapped.', default=None)

    util.cmd.common_args(parser, (('threads', None), ('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, minimap2_idxstats, split_args=True)
    return parser

__commands__.append(('minimap2_idxstats', parser_minimap2_idxstats))


def bwamem_idxstats(inBam, refFasta, outBam=None, outStats=None,
        min_score_to_filter=None, aligner_options=None,
        filterReadsAfterAlignment=False, doNotRequirePairsToBeProper=False, keepSingletons=False, keepDuplicates=False):
    ''' Take reads, align to reference with BWA-MEM and perform samtools idxstats.
        Optionally filter reads after alignment, prior to reporting idxstats, to include only those flagged as properly paired.
    '''

    assert outBam or outStats, "Either outBam or outStats must be specified"

    bam_aligned = util.file.mkstempfname('.aligned.bam')
    if outBam is None:
        bam_filtered = mkstempfname('.filtered.bam')
    else:
        bam_filtered = outBam

    samtools = tools.samtools.SamtoolsTool()
    bwa = tools.bwa.Bwa()

    ref_indexed = util.file.mkstempfname('.reference.fasta')
    shutil.copyfile(refFasta, ref_indexed)
    bwa.index(ref_indexed)

    bwa_opts = [] if aligner_options is None else aligner_options.split()
    bwa.mem(inBam, ref_indexed, bam_aligned, options=bwa_opts,
            min_score_to_filter=min_score_to_filter)
    
    if filterReadsAfterAlignment:
        samtools.filter_to_proper_primary_mapped_reads(bam_aligned, 
                                                       bam_filtered, 
                                                       require_pairs_to_be_proper=not doNotRequirePairsToBeProper, 
                                                       reject_singletons=not keepSingletons,
                                                       reject_duplicates=not keepDuplicates)
        os.unlink(bam_aligned)
    else:
        shutil.move(bam_aligned, bam_filtered)

    if outStats is not None:
        samtools.idxstats(bam_filtered, outStats)

    if outBam is None:
        os.unlink(bam_filtered)


def parser_bwamem_idxstats(parser=argparse.ArgumentParser()):
    parser.add_argument('inBam', help='Input unaligned reads, BAM format.')
    parser.add_argument('refFasta', help='Reference genome, FASTA format, pre-indexed by Picard and Novoalign.')
    parser.add_argument('--outBam', help='Output aligned, indexed BAM file', default=None)
    parser.add_argument('--outStats', help='Output idxstats file', default=None)
    parser.add_argument(
        '--minScoreToFilter',
        dest="min_score_to_filter",
        type=int,
        help=("Filter bwa alignments using this value as the minimum allowed "
              "alignment score. Specifically, sum the alignment scores across "
              "all alignments for each query (including reads in a pair, "
              "supplementary and secondary alignments) and then only include, "
              "in the output, queries whose summed alignment score is at least "
              "this value. This is only applied when the aligner is 'bwa'. "
              "The filtering on a summed alignment score is sensible for reads "
              "in a pair and supplementary alignments, but may not be "
              "reasonable if bwa outputs secondary alignments (i.e., if '-a' "
              "is in the aligner options). (default: not set - i.e., do not "
              "filter bwa's output)")
    )
    parser.add_argument(
        '--alignerOptions',
        dest="aligner_options",
        help="bwa options (default: bwa defaults)")
    parser.add_argument(
        '--filterReadsAfterAlignment',
        help=("If specified, reads till be filtered after alignment to include only those flagged as properly paired."
                "This excludes secondary and supplementary alignments."),
        action='store_true'
    )
    parser.add_argument(
        '--doNotRequirePairsToBeProper',
        help='Do not require reads to be properly paired when filtering (default: %(default)s)',
        action='store_true'
    )
    parser.add_argument(
        '--keepSingletons',
        help='Keep singleton reads when filtering (default: %(default)s)',
        action='store_true'
    )
    parser.add_argument(
        '--keepDuplicates',
        help='When filtering, do not exclude reads due to being flagged as duplicates (default: %(default)s)',
        action='store_true'
    )

    util.cmd.common_args(parser, (('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, bwamem_idxstats, split_args=True)
    return parser

__commands__.append(('bwamem_idxstats', parser_bwamem_idxstats))


def parser_extract_tarball(parser=argparse.ArgumentParser()):
    parser.add_argument('tarfile', help='Input tar file. May be "-" for stdin.')
    parser.add_argument('out_dir', help='Output directory')
    parser.add_argument('--compression',
        help='Compression type (default: %(default)s). Auto-detect is incompatible with stdin input unless pipe_hint is specified.',
        choices=('gz', 'bz2', 'lz4', 'zip', 'none', 'auto'),
        default='auto')
    parser.add_argument('--pipe_hint',
        help='If tarfile is stdin, you can provide a file-like URI string for pipe_hint which ends with a common compression file extension if you want to use compression=auto.',
        default=None)
    util.cmd.common_args(parser, (('threads', None), ('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, main_extract_tarball, split_args=True)
    return parser
def main_extract_tarball(*args, **kwargs):
    ''' Extract an input .tar, .tgz, .tar.gz, .tar.bz2, .tar.lz4, or .zip file
        to a given directory (or we will choose one on our own). Emit the
        resulting directory path to stdout.
    '''
    print(util.file.extract_tarball(*args, **kwargs))
__commands__.append(('extract_tarball', parser_extract_tarball))

# =========================

def fasta_read_names(in_fasta, out_read_names):
    """Save the read names of reads in a .fasta file to a text file"""
    with util.file.open_or_gzopen(in_fasta) as in_fasta_f, open(out_read_names, 'wt') as out_read_names_f:
        last_read_name = None
        for line in in_fasta_f:
            if line.startswith('>'):
                read_name = line[1:].strip()
                if read_name.endswith('/1') or read_name.endswith('/2'):
                    read_name = read_name[:-2]
                if read_name != last_read_name:
                    out_read_names_f.write(read_name+'\n')
                last_read_name = read_name


def read_names(in_reads, out_read_names, threads=None):
    """Extract read names from a sequence file"""
    _in_reads = in_reads
    with util.file.tmp_dir(suffix='_read_names.txt') as t_dir:
        if in_reads.endswith('.bam'):
            _in_reads = os.path.join(t_dir, 'reads.fasta')
            tools.samtools.SamtoolsTool().bam2fa(in_reads, _in_reads)
        fasta_read_names(_in_reads, out_read_names)

def parser_read_names(parser=argparse.ArgumentParser()):
    parser.add_argument('in_reads', help='the input reads ([compressed] fasta or bam)')
    parser.add_argument('out_read_names', help='the read names')
    util.cmd.common_args(parser, (('threads', None), ('loglevel', None), ('version', None), ('tmp_dir', None)))
    util.cmd.attach_main(parser, read_names, split_args=True)
    return parser

__commands__.append(('read_names', parser_read_names))


# =========================

def full_parser():
    return util.cmd.make_parser(__commands__, __doc__)


if __name__ == '__main__':
    util.cmd.main_argparse(__commands__, __doc__)

broadinstitute / viral-core / 20919343111

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous