Skip to navigation
Tool ReferenceTools

minimap2

View as Markdown

Run a GPU-accelerated minimap2.

This tool aligns long read sequences against a large reference database using accelerated seeding, chaining, and KSW2 to convert FASTQ to BAM/CRAM.

Refer to the minimap2 Reference section for a detailed listing of all available options.

Quick Start

# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun minimap2 \
    --ref /workdir/${REFERENCE_FILE} \
    --in-fq /workdir/${INPUT_FASTQ} \
    --out-bam /outputdir/${OUTPUT_BAM}

Compatible CPU-based minimap2, GATK4 Commands

The commands below are the minimap2-v2.31 and GATK4 counterpart of the Clara Parabricks command above. The output from these commands will be identical to the output from the above command. Refer to the Output Comparison page for comparing the results. You may need to increase the Java heap size based on your dataset, or decrease the number of —MAX_RECORDS_IN_RAM.

# Run minimap2 and pipe the output to create a sorted BAM.
$ minimap2 -ax map-pbmm2 \
    <INPUT_DIR>/${REFERENCE_FILE} \
    <INPUT_DIR>/${INPUT_FASTQ} | \
  gatk SortSam \
    --java-options -Xmx30g \
    --MAX_RECORDS_IN_RAM 5000000 \
    -I /dev/stdin \
    -O cpu.bam \
    --SORT_ORDER coordinate

Note that three changes must be made to the baseline minimap2 code in order to match the results exactly:

Firstly, a new preset must be made in options.c in the mm_set_opt function that tries to replicate the preset of pbmm2 by setting these parameters as a new preset named “map-pbmm2”:

io->k = 19;
io->w = 19;
io->batch_size = 0x7fffffffffffffffL; // always build a uni-part index
mo->flag |= MM_F_CIGAR;
mo->flag |= MM_F_SOFTCLIP;
mo->flag |= MM_F_LONG_CIGAR;
mo->flag |= MM_F_EQX;
mo->flag |= MM_F_NO_PRINT_2ND; // Allow secondaries with enforced mapping, but disable per default!
mo->zdrop = 400;
mo->zdrop_inv = 50;
mo->a = 1;
mo->b = 4;
mo->q = 6;
mo->q2 = 26;
mo->e = 2;
mo->e2 = 1;
mo->bw = 2000;
mo->max_gap = 10000;
mo->occ_dist = 500;
mo->min_mid_occ = 50;
mo->max_mid_occ = 500;
mo->min_dp_max = 500;

Secondly, a fix must be made to the baseline KSW2 code to round the loop fission start and end points by changing them to st and en respectively. If the start point (st0) is a number below 16, but greater than 0, its scoring values will not be initialized correctly, but will still be used later when computing the actual alignment. This can be fixed by rounding the start and end points to multiples of 16.

To make this fix, change the following code in ksw2_extd2_sse.c:

// loop fission: set scores first
if (!(flag & KSW_EZ_GENERIC_SC)) {
    for (t = st0; t <= en0; t += 16) {
        __m128i sq, st, tmp, mask;
        sq = _mm_loadu_si128((__m128i*)&sf[t]);
        st = _mm_loadu_si128((__m128i*)&qrr[t]);
        mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_));
        tmp = _mm_cmpeq_epi8(sq, st);
    #ifdef __SSE4_1__
                    tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp);
                    tmp = _mm_blendv_epi8(tmp,     sc_N_,   mask);
    #else
                    tmp = _mm_or_si128(_mm_andnot_si128(tmp,  sc_mis_), _mm_and_si128(tmp,  sc_mch_));
                    tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp),     _mm_and_si128(mask, sc_N_));
    #endif
                    _mm_storeu_si128((__m128i*)((int8_t*)s + t), tmp);
                }
            } else {
                for (t = st0; t <= en0; ++t)
                    ((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]];
            }

Fixed version that uses lf_start and lf_en:

// loop fission: set scores first
int lf_start = st, lf_en = en;
if (!(flag & KSW_EZ_GENERIC_SC)) {
    for (t = lf_start; t <= lf_en; t += 16) {
        __m128i sq, st, tmp, mask;
        sq = _mm_loadu_si128((__m128i*)&sf[t]);
        st = _mm_loadu_si128((__m128i*)&qrr[t]);
        mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_));
        tmp = _mm_cmpeq_epi8(sq, st);
    #ifdef __SSE4_1__
                    tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp);
                    tmp = _mm_blendv_epi8(tmp,     sc_N_,   mask);
    #else
                    tmp = _mm_or_si128(_mm_andnot_si128(tmp,  sc_mis_), _mm_and_si128(tmp,  sc_mch_));
                    tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp),     _mm_and_si128(mask, sc_N_));
    #endif
                    _mm_storeu_si128((__m128i*)((int8_t*)s + t), tmp);
                }
            } else {
                for (t = lf_start; t <= lf_en; ++t)
                    ((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]];
            }

Thirdly, Parabricks minimap2 uses a stable radix sort for anchors, whereas the default CPU baseline uses an unstable radix sort (‘radix_sort_128x’). This can cause ties between anchors sharing the same reference position to break differently between the two. To obtain bit-for-bit identical results against the GPU output, replace collect_seed_hits and collect_seed_hits_heap in map.c with a stable variant that sorts by (anchor.x, original_index):

typedef struct { uint64_t x; int64_t idx; } mm_stable_key_t;
static int mm_stable_key_cmp(const void *p, const void *q)
{
    const mm_stable_key_t *aa = (const mm_stable_key_t*)p;
    const mm_stable_key_t *bb = (const mm_stable_key_t*)q;
    if (aa->x != bb->x) return aa->x < bb->x ? -1 : 1;
    if (aa->idx != bb->idx) return aa->idx < bb->idx ? -1 : 1;
    return 0;
}
static mm128_t *collect_seed_hits_stable(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
                                int *n_mini_pos, uint64_t **mini_pos)
{
    int i, n_m;
    mm_seed_t *m;
    mm128_t *a, *tmp;
    mm_stable_key_t *keys;
    int64_t j;
    m = mm_collect_matches(km, &n_m, qlen, max_occ, opt->max_max_occ, opt->occ_dist, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
    a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
    for (i = 0, *n_a = 0; i < n_m; ++i) {
        mm_seed_t *q = &m[i];
        const uint64_t *r = q->cr;
        uint32_t k;
        for (k = 0; k < q->n; ++k) {
            int32_t is_self, rpos = (uint32_t)r[k] >> 1;
            mm128_t *p;
            if (skip_seed(opt->flag, r[k], q, qname, qlen, mi, &is_self)) continue;
            p = &a[(*n_a)++];
            if ((r[k]&1) == (q->q_pos&1)) { // forward strand
                p->x = (r[k]&0xffffffff00000000ULL) | rpos;
                p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
            } else if (!(opt->flag & MM_F_QSTRAND)) { // reverse strand and not in the query-strand mode
                p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
                p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
            } else { // reverse strand; query-strand
                int32_t len = mi->seq[r[k]>>32].len;
                p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | (len - (rpos + 1 - q->q_span) - 1); // coordinate only accurate for non-HPC seeds
                p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
            }
            p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
            if (q->is_tandem) p->y |= MM_SEED_TANDEM;
            if (is_self) p->y |= MM_SEED_SELF;
        }
    }
    kfree(km, m);
    keys = (mm_stable_key_t*)kmalloc(km, *n_a * sizeof(mm_stable_key_t));
    tmp  = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
    for (j = 0; j < *n_a; ++j) { keys[j].x = a[j].x; keys[j].idx = j; }
    qsort(keys, *n_a, sizeof(mm_stable_key_t), mm_stable_key_cmp);
    for (j = 0; j < *n_a; ++j) tmp[j] = a[keys[j].idx];
    memcpy(a, tmp, *n_a * sizeof(mm128_t));
    kfree(km, tmp);
    kfree(km, keys);
    return a;
}

Replace both the heap-sort and radix-sort branches in mm_map_frag_core with a single call to collect_seed_hits_stable. The original two-branch code in map.c:

if (opt->flag & MM_F_HEAP_SORT)
    a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv,
                               qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
else
    a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv,
                          qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);

Should be replaced with:

a = collect_seed_hits_stable(b->km, opt, opt->mid_occ, mi, qname, &mv,
                             qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);

Both collect_seed_hits (radix sort) and collect_seed_hits_heap (min-heap merge) produce different tie-breaking orders from each other and from the GPU pipeline, so both must be replaced. collect_seed_hits_stable supersedes both: it builds the anchor array with the same logic as collect_seed_hits but then sorts by the augmented key (anchor.x, original_index) using qsort.

The original index as a tiebreaker makes qsort behave as a stable sort: anchors at the same reference position are ordered by their emission order (query position), which matches the GPU pipeline’s CUB segmented radix sort output.

This change affects only reads that have multiple anchors mapping to the same reference position. In practice, about 0.003% (short reads) - 0.2% (long reads) of reads differ between the original unstable sort and the stable variant.

Options

TypeNameRequired?Description
I/O--ref REFYesPath to the reference file.
I/O--index INDEXNoPath to a minimizer index file generated by vanilla minimap2 to reduce indexing time.
I/O--in-fq [IN_FQ ...]NoPath to query sequence file(s) in fastq or fastq.gz format. This option can be used multiple times. Interleaved FASTQ is not supported. Provide one file for single-ended input or two files for pair-ended input.
I/O--in-bam IN_BAMNoPath to the input BAM/CRAM file.
I/O-j JUMP_BED, --jump-bed JUMP_BEDNoPath to a BED file (.bed) used in the splice:sr preset for jumping junctions in the annotation.
I/O--junc-bed JUNC_BEDNoPath to a BED file (.bed) used in adding bonus scores if an aligned junction matches a junction in the annotation.
I/O--knownSites KNOWNSITESNoPath to a database of known polymorphic sites used to exclude regions around known polymorphisms from analysis. The file can be in vcf or vcf.gz format. This option can be used multiple times.
I/O--interval-file INTERVAL_FILENoPath to an interval file in one of these formats: Picard-style (.interval_list or .picard), GATK-style (.list or .intervals), or BED file (.bed). This option can be used multiple times.
I/O--out-recal-file OUT_RECAL_FILENoPath of a report file after Base Quality Score Recalibration.
I/O--out-bam OUT_BAMYesPath of a BAM/CRAM file after sorting.
I/O--out-duplicate-metrics OUT_DUPLICATE_METRICSNoPath of duplicate metrics file after marking duplicates.
I/O--out-qc-metrics-dir OUT_QC_METRICS_DIRNoPath of the directory where QC metrics will be generated.
Tool--preset PRESETNoWhich preset to apply. Possible values are {map-pbmm2,map-hifi,map-ont,lr:hq,splice,splice:hq,splice:sr}. 'map-pbmm2' is a customized preset that uses pbmm2's default values for PacBio HiFi/CCS genomic reads. 'map-hifi' is minimap2's default preset for PacBio HiFi/CCS genomic reads. 'map-ont' is for Oxford Nanopore genomic reads. 'lr:hq' is for Nanopore Q20 genomic reads. 'splice' is for spliced long reads (strand unknown). 'splice:hq' is for Final PacBio Iso-seq or traditional cDNA. 'splice:sr' is for short RNA-seq reads. (default: map-pbmm2)
Tool--pbmm2NoInclude additional processing to match the format and accuracy of pbmm2. Not compatible with map-ont `--preset` value.
Tool--pbmm2-unmappedNoInclude unmapped records in output of pbmm2. Must be used concurrently with `--pbmm2`. Not compatible with map-ont `--preset` value.
Tool-k MINIMIZER_KMER_LEN, --minimizer-kmer-len MINIMIZER_KMER_LENNoMinimizer k-mer length.
Tool-ub, --both-strandsNoForce minimap2 to consider both strands when finding canonical splicing sites GT-AG.
Tool-uf, --forward-transcript-strandNoForce minimap2 to consider the forward transcript strand only when finding canonical splicing sites GT-AG.
Tool--mdNoOutput the MD tag.
Tool--eqxNoWrite =/X CIGAR operators.
Tool-y, --copy-commentNoAppend FASTQ comment to BAM output via auxiliary tag.
Tool-L INTERVAL, --interval INTERVALNoInterval within which to call bqsr from the input reads. All intervals will have a padding of 100 to get read records, and overlapping intervals will be combined. Interval files should be passed using the --interval-file option. This option can be used multiple times (e.g. "-L chr1 -L chr2:10000 -L chr3:20000+ -L chr4:10000-20000").
Tool-ip INTERVAL_PADDING, --interval-padding INTERVAL_PADDINGNoAmount of padding (in base pairs) to add to each interval you are including.
Tool--standalone-bqsrNoRun standalone BQSR after generating sorted BAM. This option requires both --knownSites and --out-recal-file input parameters.
Tool--read-group-sm READ_GROUP_SMNoSM tag for read groups in this run.
Tool--read-group-lb READ_GROUP_LBNoLB tag for read groups in this run.
Tool--read-group-pl READ_GROUP_PLNoPL tag for read groups in this run.
Tool--read-group-id-prefix READ_GROUP_ID_PREFIXNoPrefix for the ID and PU tags for read groups in this run. This prefix will be used for all pairs of FASTQ files in this run. The ID and PU tags will consist of this prefix and an identifier, that will be unique for a pair of FASTQ files.
Performance--num-threads NUM_THREADSNoNumber of processing threads. (default: 20)
Performance--nstreams NSTREAMSNoNumber of streams to use per GPU. Default: auto (default: auto)
Performance--gpuwriteNoUse one GPU to accelerate writing final BAM/CRAM.
Performance--gpuwrite-deflate-algo GPUWRITE_DEFLATE_ALGONoChoose the nvCOMP DEFLATE algorithm to use with --gpuwrite. Note these options do not correspond to CPU DEFLATE options. Valid options are 1, 2, and 4. Option 1 is fastest, while options 2 and 4 have progressively lower throughput but higher compression ratios. The default value is 1 when the user does not provide an input (i.e., None).
Performance--gpusortNoUse GPUs to accelerate sorting.
Performance--use-gdsNoUse GPUDirect Storage (GDS) to enable a direct data path for direct memory access (DMA) transfers between GPU memory and storage. Must be used concurrently with `--gpuwrite`. Please refer to Parabricks Documentation > Best Performance for information on how to set up and use GPUDirect Storage.
Performance--max-queue-chunks MAX_QUEUE_CHUNKSNoMax number of chunks to allow in the pre-alignment processing stage. Increasing this value may result in faster processing, but it will use more host memory. It is not recommended to use this argument for splice presets.
Performance--max-queue-reads MAX_QUEUE_READSNoMax number of reads to allow in the alignment processing stage. Increasing this value may result in faster processing, but it will use more host memory.
Performance--low-memoryNoUse low memory mode.
Performance--chunk-size CHUNK_SIZENoMax number of reads in a processing chunk. Increasing this value may result in faster processing, but it will use more host memory. Default: auto. (default: auto)
Performance--seeding-chaining-on-gpu SEEDING_CHAINING_ON_GPU, --chaining-on-gpu SEEDING_CHAINING_ON_GPUNoControl GPU seeding+chaining: 'on' always uses GPU, 'off' always uses CPU, 'auto' enables GPU for splice:sr preset and CPU otherwise. Default: auto. (--chaining-on-gpu is a deprecated alias for this option.) (default: auto)
Performance--estimated-read-length ESTIMATED_READ_LENGTHNoEstimated average read length used for pre-allocate device memory.
Performance--estimated-anchor-per-read ESTIMATED_ANCHOR_PER_READNoEstimated average number of anchors per read used to pre-allocate device and host memory for GPU chaining.
Performance--stage2-large-chaining-size STAGE2_LARGE_CHAINING_SIZENoReads with lengths greater than this threshold will be chained on CPU.
Runtime--verboseNoEnable verbose output.
Runtime--x3NoShow full command line arguments.
Runtime--dry-runNoCheck if the current command is valid.
Runtime--logfile LOGFILENoPath to the log file. If not specified, messages will only be written to the standard error output.
Runtime--tmp-dir TMP_DIRNoFull path to the directory where temporary files will be stored. (default: .)
Runtime--with-petagene-dir WITH_PETAGENE_DIRNoFull path to the PetaGene installation directory. By default, this should have been installed at /opt/petagene. Use of this option also requires that the PetaLink library has been preloaded by setting the LD_PRELOAD environment variable. Optionally set the PETASUITE_REFPATH and PGCLOUD_CREDPATH environment variables that are used for data and credentials. Optionally set the PetaLinkMode environment variable that is used to further configure PetaLink, notably setting it to "+write" to enable outputting compressed BAM and .fastq files.
Runtime--keep-tmpNoDo not delete the directory storing temporary files after completion.
Runtime--no-seccomp-overrideNoDo not override seccomp options for docker.
Runtime--versionNoView compatible software versions.
Runtime--preserve-file-symlinksNoOverride default behavior to keep file symlinks intact and *not* resolve the symlink.
Runtime--help-jsonNoPrint out help message in JSON format.
Runtime--num-gpus NUM_GPUSNoNumber of GPUs to use for a run. (default: 1)