Metadata-Version: 2.4
Name: bamurai
Version: 0.14.0
Summary: A package for splitting reads in BAM/FASTQ files into smaller fragments.
Author: Shian Su
License: Apache-2.0
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: pysam
Requires-Dist: numpy
Requires-Dist: pandas
Requires-Dist: tqdm
Provides-Extra: dev
Requires-Dist: pytest; extra == "dev"
Requires-Dist: black; extra == "dev"
Requires-Dist: mypy; extra == "dev"
Requires-Dist: coverage; extra == "dev"
Requires-Dist: pytest-cov; extra == "dev"
Dynamic: license-file

# Bamurai

A Python toolkit for manipulating BAM, CRAM and FASTQ files, designed to split reads into smaller fragments, extract statistics, validate files, and manage multi-sample data.

**For BAM/SAM/CRAM files, Bamurai counts only primary alignments**, so each sequenced read is counted once. Secondary and supplementary alignments are other mappings of the same read, not additional reads.

What happens to secondary and supplementary records depends on the command:

- **Commands that output reads, rows or statistics** (`split`, `divide`, `chunk`, `stats`, `extract_tag`) ignore them. Including them would duplicate reads and inflate counts.
- **Commands that write a subset of an alignment file** (`add_tags`, `assign_samples`, `extract_sample`, `split_samples`, `split_by_tag`, `extract_regions_by_group`) keep them, so the output still works with tools that rely on them, such as quantifiers that resolve multi-mapping reads from secondary alignments. Each record is handled by its own barcode or tag exactly as a primary read would be: matched, tagged, or sent to the unmatched/undetermined output by the same rules. The one exception is `extract_sample`, which drops secondary and supplementary records with no barcode tag and warns with a count. Summaries count primary reads only. Pass `--primary-only` to write no secondary or supplementary records.

## Description

Bamurai is a command-line tool for working with long-read BAM/CRAM/FASTQ files: splitting reads into smaller fragments, summarising and validating files, and splitting, extracting or tagging alignments by cell barcode metadata.

These are the current features of Bamurai:

1. Splitting reads in a file to a target length (`split`)
2. Splitting reads in a file to a target number of pieces per read (`divide`)
3. Getting statistics from a BAM/CRAM or FASTQ(.gz) file (`stats`)
4. Basic validation of BAM, CRAM and FASTQ(.gz) files (`validate`)
5. Splitting a file into size-based chunks (`chunk`)
6. Splitting or extracting alignments by donor ID using a barcode mapping
   (`split_samples`, `extract_sample`, `assign_samples`)
7. Extracting hashtag oligo (HTO) information from 10x FASTQ pairs (`get_hto`)
8. Tagging alignments with arbitrary values from a TSV keyed by read ID or barcode (`add_tags`)
9. Extracting a BAM tag's value per read into a TSV (`extract_tag`)
10. Splitting a BAM file into one file per distinct value of a tag (`split_by_tag`)
11. Extracting alignments in genomic regions and splitting them by barcode metadata
    such as donor or cell type (`extract_regions_by_group`)

The `split` command splits reads into a target length, each read will be split into fragments as close to the target length as possible. Reads shorter than the target length will not be split.

The `divide` command splits reads into a target number of pieces (default 2), each read will be split into the number of pieces specified. A further minimum length can be specified to ensure that reads are not split if the resultant fragments are less than the minimum length (default 100).

The `stats` command will output the following information by default:
```
Statistics for input.bam:
  Total reads: 8160
  Average read length: 30638
  Throughput (Gb): 0.25
  N50: 82547
```

It can be used with the `--tsv` argument to output the statistics in a tab-separated format for computational analysis. Note that `throughput` is reported in raw bases here, rather than the gigabases of the default output.
```bash
file_name       total_reads     avg_read_len    throughput      n50
input.bam      8160    30638   250006998       82547
```

The `validate` command will check the integrity of a BAM, CRAM or FASTQ(.gz) file and output the following information if the file is valid:
```bash
input.bam is a valid BAM file with 8160 records.
```

Unlike the read-processing commands, `validate` inspects every record in the file, so the count it reports includes any secondary and supplementary alignments.

If the file is invalid, `validate` prints the first problem it finds and exits with a non-zero status, as it does for a file in an unsupported format, so it can gate later steps in a pipeline.

## Installation

To install the released version of Bamurai from PyPI

```bash
pip install bamurai
```

To install the latest version of Bamurai from GitHub

```bash
pip install git+https://github.com/Shians/Bamurai.git
```

## Usage

To get help on the command-line interface and list available commands
```bash
bamurai --help
```

To get help on a specific command
```bash
bamurai <command> --help
```

A command that fails on bad input, such as a missing file or an unknown TSV
column, prints `bamurai: error: <message>` and exits with status 1. To see the
full Python traceback, for example when reporting a bug, put `--debug` before
the command:
```bash
bamurai --debug stats input.bam
```

### CRAM files

Every command that takes a BAM also takes a CRAM. CRAM stores reads as
differences from a reference genome, so reading one needs the reference FASTA
it was written against. Pass it with `--reference`:
```bash
bamurai stats input.cram --reference genome.fa
```

When reading, `--reference` can be left out if htslib can find the reference
itself, through the path in the CRAM header or the `REF_PATH` and `REF_CACHE`
environment variables. Writing CRAM always needs `--reference`.

Commands that take several inputs accept any mix of BAM and CRAM, but every
input must have the same reference sequences, as their reads are merged into
shared output files. One `--reference` covers all inputs. Bamurai checks each
input before writing anything, and stops with an error if a CRAM cannot be
read or if an input's sequence names or lengths differ from the reference.

Commands that write alignments choose the output format as follows:

| Command | Output format |
| --- | --- |
| `add_tags`, `extract_sample`, `assign_samples` | From the `--output` extension (`.cram`, `.sam`, otherwise BAM) |
| `split_by_tag` | Same as the input, or `--output-format bam\|cram\|sam` |
| `split_samples`, `extract_regions_by_group` | BAM, or `--output-format cram` |

### Splitting reads to target size

To split a file into 10,000 bp reads
```bash
bamurai split input.bam --len-target 10000 --output output.fastq
```

An output name ending in `.gz` is written gzipped. Without `--output`, reads
go to stdout:
```bash
bamurai split input.bam --len-target 10000 --output output.fastq.gz
bamurai split input.bam --len-target 10000 | gzip > output.fastq.gz
```

### Dividing reads into a target number of pieces

To divide reads into 2 pieces
```bash
bamurai divide input.bam --num-fragments 2 --output output.fastq
```

To divide reads into 2 pieces unless resultant fragments are less than 1000 bp
```bash
bamurai divide input.bam --num-fragments 2 --min-length 1000 --output output.fastq
```

### Getting statistics from a BAM/CRAM or FASTQ file

To get stats from a BAM file
```bash
bamurai stats input.bam
```

To get stats from a FASTQ file or Gzipped FASTQ file
```bash
bamurai stats input.fastq
bamurai stats input.fastq.gz
```

### Validating BAM, CRAM or FASTQ files

To validate a BAM file
```bash
bamurai validate input.bam
```

As the SAM spec allows, a record may have no sequence or no qualities (`*`),
as aligners such as minimap2 write secondary alignments. Such a file is valid.
If a primary read lacks either, `validate` gives a warning with a count,
because some commands can't use it: `split`, `divide` and `chunk` stop with
an error at a read without qualities, since FASTQ needs them, and every
command that reads sequences stops at a read without a sequence. `stats`
counts reads without qualities as normal.

### Splitting a file into chunks

To split a file into chunks of at least 1 GB each
```bash
bamurai chunk input.bam --size 1G
```

Sizes are given in human-readable notation (`1G`, `100M`, `1000K`). Each chunk is at least the requested size, and may be slightly larger because reads are never split across chunks.

Output files are written as `chunk_1.fastq`, `chunk_2.fastq`, and so on. Use `--prefix` to change the name
```bash
bamurai chunk input.bam --size 100M --prefix sample_a
```

A rerun can produce fewer chunks than an earlier one, which would leave the earlier run's
higher-numbered chunks mixed in with the new ones. For that reason the command fails if any
`<prefix>_<N>.fastq` file already exists. Pass `--overwrite` to replace them: files with the
same names are replaced, and a warning lists any other matching chunk files from the earlier
run, which are left in place. Chunks are moved into place only after the whole input has been
chunked, so a failed run leaves no partial chunks and existing files untouched.

### Working with multi-sample BAM files

Bamurai provides commands for processing BAM or CRAM files with multiple samples based on barcode information.

#### Splitting BAM or FASTQ files by donor ID

To split a BAM/CRAM or FASTQ file into multiple files, one for each donor ID:

```bash
bamurai split_samples --input input.bam --tsv barcode_to_donor.tsv --output-dir donor_bams
bamurai split_samples --input input.fastq.gz --tsv barcode_to_donor.tsv --output-dir donor_fastqs
```

The TSV file should contain a barcode column and a donor ID column, one row per barcode; see [Barcode and donor columns](#barcode-and-donor-columns) below for how those columns are found.

You can process multiple BAM or FASTQ files at once:

```bash
bamurai split_samples --input input1.bam input2.bam --tsv barcode_to_donor.tsv --output-dir donor_bams
bamurai split_samples --input input1.fastq.gz input2.fastq.gz --tsv barcode_to_donor.tsv --output-dir donor_fastqs
```

For BAM input, most of the time goes on compressing the output files. `-t/--threads` (default 1) sets how many CPUs to use in total, however many donors or groups the library is split into. With BAM output, each library is cut into that many parts, which are split in parallel and then joined without being compressed again; the outputs are identical to a single-threaded run, read order included. With CRAM input or output, the threads only decode the input, so they help less. The parts are written to a hidden folder in `--output-dir` and need about as much space as the outputs.

```bash
bamurai split_samples --input input.bam --tsv barcode_to_donor.tsv --output-dir donor_bams --threads 4
```

Every donor in the TSV gets a file, even if it receives no reads, and
reads with no barcode, or one not in the TSV, go to `unmapped.bam` (or
`unmapped.fastq`). For BAM/CRAM input, every alignment is split by its own
barcode, secondary and supplementary ones included; pass `--primary-only` to
write primary alignments only. For FASTQ input, the barcode is read from
tab-separated tag fields in each header line (see
[Which barcode tag is read](#which-barcode-tag-is-read)), and
`--primary-only` has no effect. Inputs must be all BAM/CRAM or all FASTQ.
`--output-dir` must be empty or
new, since files left by a run with a different TSV would mix in with the
new donor files; pass `--overwrite` to write into a non-empty directory
anyway. Files there that this run does not write are reported, not removed.

#### Extracting alignments for a specific donor

To extract all alignments belonging to a specific donor from a BAM file:

```bash
bamurai extract_sample --bam input.bam --tsv barcode_to_donor.tsv --donor-id donor1 --output donor1.bam
```

You can also process multiple BAM files at once, combining all donor-specific alignments into a single output file:

```bash
bamurai extract_sample --bam input1.bam input2.bam input3.bam --tsv barcode_to_donor.tsv --donor-id donor1 --output donor1.bam
```

This writes every alignment whose barcode belongs to the specified donor, secondary and supplementary ones included, to a new BAM file. Secondary and supplementary alignments with no barcode tag are dropped, with a warning giving the count. If the TSV lists no barcodes for the donor, a message is printed and no file is written.

#### Tagging alignments with their donor ID

The `assign_samples` command labels each alignment in a BAM with the donor it came from, using the `RG` (read group) tag. Unlike `split_samples`, which writes one file per donor, this keeps every alignment in a single BAM so that downstream tools can group by read group.

```bash
bamurai assign_samples --bam input.bam --tsv barcode_to_donor.tsv --output assigned.bam
```

This writes `assigned.bam` containing every alignment from the input, where:

- alignments whose cell barcode maps to a donor carry an `RG` tag naming that donor
- alignments with an unknown or absent barcode are written unchanged
- the header gains an `@RG` record per donor, with `ID` and `SM` both set to the donor ID, unless the input already has one with that ID

Secondary and supplementary alignments are assigned by their own barcode,
like primary ones; pass `--primary-only` to write primary alignments only.

A summary of how many reads were assigned to each donor is printed on completion. Note that `--bam` takes a single BAM, unlike `split_samples` and `extract_sample`, which accept several.

#### Barcode and donor columns

`split_samples`, `extract_sample` and `assign_samples` all read the same barcode-to-donor TSV, and all detect its columns the same way. The barcode column is taken from a column named `barcode` or `cell`, and the donor column from one named `donor_id`. If both `barcode` and `cell` are present the file is ambiguous and the command fails, asking you to choose.

Use `--barcode-column` and `--donor-id-column` when your columns are named something else, or to resolve that ambiguity:

```bash
bamurai extract_sample --bam input.bam --tsv mapping.tsv --donor-id donor1 \
  --barcode-column cell_barcode --donor-id-column sample --output donor1.bam
```

#### Which barcode tag is read

An alignment's cell barcode is read from the first of these tags that is present:

| Tag | Meaning |
| --- | --- |
| `CB` | Cell barcode, error-corrected (preferred) |
| `CR` | Cell barcode, raw/uncorrected |
| `XC` | Cell barcode, Drop-seq and early 10x convention (not in the SAM spec) |

`BC` and `RX` are deliberately not consulted: per the SAM specification `BC` identifies a *sample* and `RX` holds the UMI, so neither names the cell that a donor is assigned to. A read carrying only those tags is treated as having no barcode.

`split_samples` reads the same tags from FASTQ input, as tab-separated
fields after the read name in the header line, the form `samtools fastq -T
CB` writes:
```
@read1	CB:Z:AAACCCAAAGGGTTT
```

`split_by_tag` and `extract_regions_by_group` don't use this list: they read
the one tag you name (`--tag`).

### Extracting HTO information

The `get_hto` command pulls hashtag oligo (HTO) information out of a pair of 10x FASTQ files. It assumes R1 holds the cell barcode followed by the UMI, and R2 holds the hashtag sequence.

```bash
bamurai get_hto --r1 sample_R1.fastq.gz --r2 sample_R2.fastq.gz \
  --bc-len 16 --umi-len 12 --output htos.tsv
```

`--bc-len` and `--umi-len` give the barcode and UMI lengths in R1, and are required. The output is a tab-separated file with the columns `read_name`, `cell_barcode`, `umi`, `hto`, `bc_qual`, `umi_qual` and `hto_qual`, where the three `_qual` columns carry the FASTQ quality strings for the preceding sequences.

Where in R2 the hashtag is read from can be adjusted with `--hashtag-len` (how many bases to take, default 15) and `--hashtag-left-buffer` (how many bases to skip first, default 10).

### Tagging reads from a TSV

To tag reads in a BAM file with values from a TSV, keyed by read_id:
```bash
bamurai add_tags --bam input.bam --tsv tags.tsv --by-read-id --output tagged.bam
```

Or keyed by barcode, tagging every read that carries a given barcode:
```bash
bamurai add_tags --bam input.bam --tsv barcode_tags.tsv --by-barcode --output tagged.bam
```

`--by-read-id` and `--by-barcode` are required and mutually exclusive, and
pick which column of the TSV is the key. The TSV must have a header row with
that key column, plus one column per tag to apply, where each tag column
name is the BAM tag to apply (a letter followed by a letter or digit, e.g.
`XX`, `RG`) — every column besides the key column is applied as a tag, and
column names must be unique. Values are strings unless the column header
declares a type with a suffix, as in SAM's `XX:i:5`: `XX:i` for integers
(32-bit), `XX:f` for finite floats (`nan` and `inf` are rejected), `XX:A` for
a single character, and `XX` or `XX:Z` for strings. A value that doesn't fit
its column's type is an error, so a numeric-looking label such as a donor
`007` stays `007` unless you ask for a number. Earlier versions guessed
the type of each cell; TSVs with numeric columns need the suffix to keep
writing numbers. A blank cell for a given tag column just leaves that tag off, and a
key whose tag cells are all blank still counts as matched (its reads are
written with no tags added). Whitespace around header names, keys and values
is stripped, and a UTF-8 byte order mark is ignored.

`--by-read-id` requires a `read_id` column; reads whose read ID is not in
the TSV are unmatched. `--by-barcode` requires a `barcode` or `cell` column;
an alignment's barcode is read from its `CB`, `CR`, or `XC` tag (see
[Which barcode tag is read](#which-barcode-tag-is-read)), and alignments with
no barcode tag, or a barcode not present in the TSV, are unmatched, with a
summary count of each (for primary reads) printed at the end. If the same barcode appears on
multiple TSV rows with different tag values, `add_tags` raises an error
rather than picking one silently. Passing a flag whose required column is
missing from the TSV — e.g. `--by-barcode` on a TSV with only `read_id` —
also raises, rather than falling back to the other mode.

Every alignment is written to `--output`: matched ones with their tags,
unmatched ones unchanged. To keep only some cells, use `extract_sample` or
`split_samples`. Secondary and supplementary alignments are matched exactly
like primary ones (by read ID, or by their own barcode tag), so they get the
same tags. Pass `--primary-only` to write primary alignments only:

```bash
bamurai add_tags --bam input.bam --tsv barcode_tags.tsv --by-barcode \
  --output tagged.bam --primary-only
```

### Extracting a tag to a TSV

To extract the value of a BAM tag for every primary read into a TSV:
```bash
bamurai extract_tag --bam input.bam --tag CB --output cb.tsv
```

The output has a header row `read_id`, `<tag>` and one row per primary read
that carries the tag. By default, reads missing the tag are skipped; pass
`--include-missing` to write them too, with a blank value in the tag column:
```bash
bamurai extract_tag --bam input.bam --tag CB --output cb.tsv --include-missing
```

Because the output uses `read_id` as its key column and the tag name itself
as the value column, it can be fed straight back into `add_tags --by-read-id`
to carry a tag from one BAM file onto reads in another. The value column is
headed with the tag's type unless it is a string (`NH:i`, `AS:f`, `XC:A`), so
integers, floats and characters come back as the same type. A tag must have
the same kind of value (string, character, integer or float) in every read,
and array (`B`) and hex (`H`) tags are an error. `add_tags` applies
it to every alignment of each read, secondary and supplementary ones
included.

### Splitting a BAM file by tag value

To split a BAM/CRAM/SAM file into one output file per distinct value of a
tag:
```bash
bamurai split_by_tag --bam input.bam --tag RG --output-dir output --prefix by_rg
```

This writes `output/by_rg_<value>.bam` for each distinct value the tag takes,
e.g. `by_rg_sample1.bam`, `by_rg_sample2.bam`. Outputs are in the input's
format unless `--output-format` is given. If `--prefix` is omitted it
defaults to `<input_basename>_<tag>`. Secondary and supplementary alignments
are split by their own tag, like primary ones; pass `--primary-only` to write
primary alignments only. By default, alignments without the tag are dropped
rather than written; pass `--unmatched-output` to include them instead, in
`output/by_rg_undetermined.bam`:
```bash
bamurai split_by_tag --bam input.bam --tag RG --output-dir output --prefix by_rg \
  --unmatched-output
```

`split_by_tag` is meant for tags with few values, such as `RG`. It stops
without writing anything if the tag is an array, has more than 200 distinct
values, or has a value that is empty, longer than 100 characters or
`undetermined`. In file names, any character other than letters, digits and
`. _ + -` becomes `_`, so it also stops if two values would get the same file
name, such as `a/b` and `a_b`, or values that differ only in case (the same
file on macOS and Windows). To split by a many-valued tag such
as `CB`, map its values to groups in a TSV and use `split_samples` or
`extract_regions_by_group`.

### Extracting alignments in regions, split by barcode metadata groups

To pull the alignments overlapping some regions out of one or more indexed BAM
or CRAM files and split them by columns of a barcode metadata TSV:
```bash
bamurai extract_regions_by_group --bam sample1.bam sample2.bam --tsv cells.tsv \
  --columns donor,celltype --region chr1:1,000,000-2,000,000 chr6 \
  --output-dir output
```

The TSV needs a `barcode` (or `cell`) column plus the columns named in
`--columns`:
```
barcode	donor	celltype
AAACCCAAAGGGTTT	donor1	T cells
CCCGGGTTTAAACCC	donor2	B cells
```

Each alignment's barcode is taken from the `CB` tag (change with `--tag`) and
the alignment is written to `output/<group>.bam`, where the group joins the barcode's
values in `--columns` with `_`, e.g. `donor1_T_cells.bam`. Spaces in values
become `_`. With `--columns donor` the outputs are `donor1.bam`,
`donor2.bam`, and so on.

Values are read literally, so `NA` or `null` are ordinary values; only blank
cells count as missing. Values in `--columns` may contain only ASCII letters,
digits, `.`, `-`, `_` and spaces. The command stops with an error before
writing anything if:
- a value in `--columns` has any other character, e.g. `CD4+`, `T/cells` or
  `café`;
- a group name starts with `.`, which would make a hidden file, or is too
  long for a file name;
- a barcode appears in more than one row of the TSV;
- two different sets of values would give the same file name, e.g. `T cells`
  and `T_cells`, `Tcell` and `tcell`, or `a_b`+`c` and `a`+`b_c`;
- a group is named `undetermined`, which is reserved for unmatched reads;
- `--output-dir` already exists and is not empty (see below).

Regions are samtools-style strings (1-based, inclusive), a BED file via
`--bed` (0-based, half-open), or both. A region that starts past the end of
its contig is skipped with a warning, which usually means coordinates from a
different genome build. Alignments from every input file and every region go
into the same group files, and an alignment overlapping several regions is
written once. Each alignment is selected by its own position, so a primary
alignment in the regions can be written without a secondary or
supplementary alignment of the same read that lies outside them, and the
other way round. Outputs are coordinate-sorted and indexed. By default,
alignments whose barcode is missing or not in the TSV (or has a blank value in any of
`--columns`) are dropped; pass `--unmatched-output` to write them to
`output/undetermined.bam`. The summary reports how many primary reads were dropped,
and a warning is shown if no read in the regions matched the TSV, which
usually means the barcodes are in a different tag than `--tag`.

Only groups that receive reads get a file, so files left by an earlier run
would mix with the new results. For that reason `--output-dir` must be new
or empty. Pass `--overwrite` to write into a non-empty directory anyway:
files with the same names are replaced, other files are left in place, and
a warning lists any BAM or CRAM files this run did not write. Outputs never
replace an input file, even through a link or a name differing only in case. Outputs
are moved into place only after every group has been sorted and indexed, so
a run that fails part-way leaves the output directory as it was.

Input must be indexed BAM or CRAM; SAM files are rejected, so convert them
with `samtools view -b` and index the result first. With `--output-format
cram`, outputs are `<group>.cram` with a `.crai` index. Each group keeps an
output file open while the inputs are read, so the command raises the
process's open-file limit when there are many groups. If the system limit
is too low, it stops and suggests a `ulimit -n` value.

