piscem is a next-generation index + mapper tool consisting of a compacted de Bruijn graph constructor (with tiling information), and a very efficient index and read mapper. It is a thin command-line wrapper around piscem-rs (indexing and mapping) and cf1-rs (cDBG construction); both are pure Rust, so building piscem needs no C or C++ toolchain.

Notes

• If you are primarily interested in simply using piscem, you can obtain pre-compiled binaries from the GitHub releases page for linux x86-64, linux ARM, OSX x86-64, and OSX ARM. Likewise, piscem can be installed using bioconda. The instructions below are primarily for those who need to build piscem from source.

• Please ensure that the user file handle limit is set to 2048. This may already be set (and should be fine already on OSX), but you can accomplish this by executing:

$ ulimit -n 2048

before running piscem.

Building

git clone https://github.com/COMBINE-lab/piscem.git
cd piscem
cargo build --release

The checked-in .cargo/config.toml selects a portable baseline per target (x86-64-v3 + AVX2 on x86-64, neoverse-n1 on linux ARM, apple-a14 on OSX ARM). It deliberately does not use target-cpu=native, because source packagers such as bioconda honour this file, and a binary built for the build machine's exact CPU crashes with SIGILL elsewhere. If you are building only for the machine you are on and want the last few percent, override it locally:

RUSTFLAGS="-C target-cpu=native" cargo build --release

Usage

Fast, space-efficient k-mer-based read mapper for bulk RNA-seq, scRNA-seq, and scATAC-seq

Usage: piscem [OPTIONS] <COMMAND>

Commands:
  build        Index a reference sequence
  map-sc       map reads for single-cell processing
  map-bulk     map reads for bulk processing
  map-sc-atac  map reads for scAtac processing
  help         Print this message or the help of the given subcommand(s)

Options:
  -q, --quiet    be quiet (no effect yet for cDBG building phase of indexing)
  -h, --help     Print help
  -V, --version  Print version

Each subcommand is summarised below. Run piscem <subcommand> --help for the authoritative and complete option list — the options are numerous enough that reproducing them here only invites drift.

build

Indexes one or more reference sequences, building a piscem index over them. References are supplied through exactly one of --ref-seqs (a , separated list of FASTA files), --ref-lists (files that each list FASTA paths), or --ref-dirs (directories whose FASTA files are all indexed, non-recursively).

The -k parameter is the k-mer size for the underlying colored compacted de Bruijn graph (must be odd and <= 31), and -m is the minimizer size used to build the sshash data structure (must be < k). Passing --decoy-paths adds poison k-mer information from the given decoy sequences to the index.

Two directories are involved, and they are not the same:

• --work-dir (default ./workdir.noindex) holds the temporary files for cDBG construction.
• --tmp-dir holds SSHash's external minimizer-sort scratch, a later phase. Unset, SSHash uses sshash_tmp in the current directory.

--ram-limit-gib caps the RAM that external minimizer sort may use (default 8 GiB); a smaller value spills to disk sooner. Both are useful on HPC systems where /tmp is small or the job runs under a cgroup memory limit. The index build also detects available memory automatically, walking the cgroup chain so a container limit is respected rather than the host's total.

Note If you are running on an Apple silicon machine, it is highly recommended that you set -t to be less than or equal to the number of high performance cores that you have (rather than the total number of cores including efficiency cores), as using efficiency cores in the piscem build step has been observed to severely degrade performance.

map-sc

Maps single-cell sequencing reads against a piscem index and produces a RAD format output file that can be processed by alevin-fry.

You can provide multiple files to -1 and -2 as a , separated list; it is important to provide them in the same order to both flags. The --geometry flag specifies the geometry of the UMIs and cell barcodes for the reads; see geometry below.

map-bulk

Maps bulk sequencing reads against a piscem index. The tool performs non-spliced alignment, and is therefore applicable to e.g. metagenomic reads against a set of metagenomes, DNA-seq alignment against one or more references, or RNA-seq alignment against a transcriptome (but not a genome). It produces a bulk RAD format output file that can be processed by piscem-infer to estimate the abundances of all references in the index given the mapped reads.

Reads are given either as paired -1/-2 lists or as a single --reads list.

Note For map-bulk, -o is an output file stem, not a directory: -o some/dir/sample writes some/dir/sample.rad and some/dir/sample.map_info.json. This differs from map-sc and map-sc-atac, where -o names a directory that map.rad is written into.

map-sc-atac

Maps single-cell ATAC-seq reads against a piscem index, producing a RAD file for alevin-fry atac. Reads may be paired (-1/-2) or single (--reads), and the cell barcodes are supplied separately through --barcode.

threading and decompression

For all three mapping subcommands, -t is a single execution-slot budget shared between mapping and gzip decoding — not a mapping-thread count. The effective budget is capped by the parallelism actually available, so a cpuset or cgroup limit is respected; both the requested and effective values are recorded in map_info.json.

--decoder selects how gzip input is decompressed:

value	behaviour
auto (default)	opens seekable gzip through the shared parallel decoder pool and lets a live thread broker adapt the mapping/decode split during the run
serial	disables the parallel pool and gives mapping the full budget
parallel	forces the parallel path wherever the input supports positional reads, still adapting the split
parallel=N	fixes N decode slots per decoder-capable input and disables adaptation

Non-regular inputs — FIFOs, process substitutions such as -1 <(zcat …) — always use the serial path, because the parallel decoder requires positional reads. The downgrade is per file, so one FIFO among the inputs does not cost the regular files their parallel decoder.

--thread-policy <FILE> points at a JSON document overriding thread and decoder policy. Every field is optional and falls back to a measured default, so a file need only name what it changes; an unknown field is a hard error rather than a silent no-op. Currently:

{ "parallel_decode": { "min_threads_per_stream": 8 } }

That value is how many threads must be available per gzip input before the parallel decoder is engaged at all. Below it, the serial decoder is faster, because it decodes inline on threads that also map and so is never idle, while a dedicated decode slot cannot map.

Note that with more than one mapping thread the byte layout of the RAD file varies from run to run, because records are grouped into chunks by whichever worker finishes first. The mapping results are unaffected. Use -t 1 if you need byte-identical output to compare against.

For the full rationale and the measurements behind these defaults, see docs/threading-and-decompression.md in piscem-rs.

geometry

The geometry parameter --geometry|-g can take either a specific geometry name, or a generic specifier string. The current valid names are chromium_v2, chromium_v2_5p, chromium_v3, chromium_v3_5p, and chromium_v4_3p. The custom format is as follows: you must specify the content of read 1 and read 2 in terms of the barcode, UMI, and mappable read sequence. A specification looks like this:

1{b[16]u[12]x:}2{r:}

In particular, this is how one would specify the 10x Chromium v3 geometry using the custom syntax. The format string says that the read pair should be interpreted as read 1 1{...} followed by read 2 2{...}. The syntax inside the {} says how the read should be interpreted. Here b[16]u[12]x: means that the first 16 bases constitute the barcode, the next 12 constitute the UMI, and anything that comes after that (if it exists) until the end of read 1 should be discarded (x). For read 2, we have 2{r:}, meaning that we should interpret read 2, in its full length, as biological sequence.

It is possible to have pieces of geometry repeated, in which case they will be extracted and concatenated together. For example, 1{b[16]u[12]b[4]x:} would mean that we should obtain the barcode by extracting bases 1-16 (1-based indexing) and 29-32 and concatenating them together to obtain the full barcode. A specification that is followed by a specific length (i.e. a number in [] like b[10] or x[4]) is said to be bounded. The specification string can have many bounded pieces, but only one unbounded piece (an unbounded piece is a specifier like r or x, followed by :). Likewise, since the : specifier means to extract this piece until the end of the string, the unbounded specifier must be the last specifier in the description of each read (if it occurs).