STAR BAM sort memory: --outBAMsortingBinsN is not the lever, read order is

[lhqing]

STAR BAM coordinate-sort memory: what the source actually proves

Sources. All code citations are against alexdobin/STAR at commit
b1edc1208d91a53bf40ebae8669f71d50b994851, which is both master HEAD (checked 2026-08-03) and the
2.7.11b release tag — i.e. exactly the version you measured, and still the newest STAR that exists.
Manual citations are doc/STARmanual.pdf from that tree. Issue citations are GitHub issues on that
repo, restricted to comments authored by alexdobin.

Permalink prefix: https://github.com/alexdobin/STAR/blob/b1edc1208d91a53bf40ebae8669f71d50b994851/

---

⚠ REVISION NOTICE (round 2)

A real Log.out from a sorted STARsolo run was obtained after this report was first written. It
confirms the root cause (non-random read order) but corrects the mechanism sketch in two
places and adds an important caveat to the total/(nBins-1) formula. See
§ Round 2 at the end — read it
before acting on §"Answer up front" items 3, 4 and 6. Everything in §Evidence 1-7 stands unchanged.

Answer up front

1. Your hypothesis is wrong about the mechanism. Bin boundaries are NOT equal genomic length.
They are equal-count quantiles of observed alignment start positions — STAR sorts a sample of
start positions and cuts at equally-spaced ranks. Confidence: certain (BAMoutput.cpp:138-142,
quoted below). A coverage pileup spread over ~10 kb is therefore not, by itself, able to defeat
binning: the quantile cuts will subdivide the pileup itself.

2. --outBAMsortingBinsN genuinely does reduce peak sort memory — it is the primary knob.
--limitBAMsortRAM is compared against the largest single bin, not the total. So the sort
requirement is total / (nBins-1) when binning works. Dobin states the same formula himself
(issue #457): "the max RAM required for sorting can be approximated as the size of unzipped fastq
files ... divided by the --outBAMsortingBinsN value." Confidence: certain for the code,
high for the practical formula.

3. Memory is therefore NOT linear in total records — it is linear in the largest bin. It is
bounded, and the bound is under your control, provided the quantile sample is representative.

4. The reason binning did nothing for you is almost certainly NOT the rDNA pileup. It is almost
certainly that your FASTQ has non-random read order. The quantile boundaries are estimated from
the first ~50 MB of BAM output — i.e. the first few hundred thousand reads of the file and
nothing else. If those reads are in genome-coordinate order (the universal signature of a FASTQ
regenerated from a CellRanger/coordinate-sorted BAM via bamtofastq), every boundary lands inside
one tiny genomic window, and all the rest of the genome falls past the last boundary into a single
bin. Then max-bin ≈ total, and no value of --outBAMsortingBinsN can help — which is exactly what
you measured. Dobin has diagnosed this same failure four separate times (#289, #1136, #1762, #2070),
and in #2070 the user's requirement dropped from 400 GB to 60 GB purely by shuffling the FASTQ.
Confidence: high (see the arithmetic in "Does it explain our measurement?" — your numbers are
~74x above what working bins would produce, which is the signature of total collapse, not of a 44%
pileup).

5. There is no spill path. None. At all. The sort is strictly in-RAM per bin: one new char[binS+1]
for the whole bin plus one new uint[binN*3], qsort, write, free. Refusal-or-bad_alloc are the
only outcomes. Dobin, 2018: "robust sorting is high on my TODO list" — still open in 2025.

6. For your 2.44-billion-record sample: yes, it is sortable in sane RAM — but the BAM sort is not
your binding constraint. With working bins at --outBAMsortingBinsN 200 the sort needs ~2 GB per
bin. The real floor is STARsolo's readInfo array, 16 bytes × 2.23e9 input reads = ~36 GB, which
is invisible to --limitBAMsortRAM and is held live during the sort. Plus rGeneUMI at ~24 GB.
Budget ~70-90 GB for that sample, dominated by Solo, not by sorting. There is also a hard architectural
ceiling at 2^32 = 4.295 billion reads (read index is packed into the top 32 bits of a 64-bit word);
you are at 52% of it.

---

Evidence

1. How bin boundaries are computed — quantiles of a sample, not genomic length

source/BAMoutput.cpp:118-144 (BAMoutput::coordBins) — the decisive function.
permalink

void BAMoutput::coordBins() {//define genomic starts for bins
    nBins=P.outBAMcoordNbins;//this is the true number of bins
    ...
    if (P.outBAMsortingBinStart[0]!=0) {//it's set to 0 only after the bin sizes are determined
        //extract coordinates and sort
        uint *startPos = new uint [binTotalN[0]+1];//array of aligns start positions
        for (uint ib=0,ia=0;ia<binTotalN[0];ia++) {
            uint32 *bamIn32=(uint32*) (binStart[0]+ib);
            startPos[ia]  =( ((uint) bamIn32[1]) << 32) | ( (uint)bamIn32[2] );
            ib+=bamIn32[0]+sizeof(uint32)+sizeof(uint);
        };
        qsort((void*) startPos, binTotalN[0], sizeof(uint), funCompareUint1);

        //determine genomic starts of the bins
        P.inOut->logMain << "BAM sorting: "<<binTotalN[0]<< " mapped reads\n";
        P.inOut->logMain << "BAM sorting bins genomic start loci:\n";

        P.outBAMsortingBinStart[0]=0;
        for (uint32 ib=1; ib<(nBins-1); ib++) {
            P.outBAMsortingBinStart[ib]=startPos[binTotalN[0]/(nBins-1)*ib];
            P.inOut->logMain << ib <<"\t"<< (P.outBAMsortingBinStart[ib]>>32) << "\t" << ((P.outBAMsortingBinStart[ib]<<32)>>32) <<endl;
            //how to deal with equal boundaries???
        };

Read startPos[binTotalN[0]/(nBins-1)*ib]: startPos is the sorted array of alignment start
positions, and the index is an equally spaced rank. That is a quantile cut. Nothing here
references chrStart, nGenome, or chromosome lengths. Your "equal genomic length" hypothesis is
refuted by this line. The design intent is explicitly equal counts per bin.

The author's own inline comment //how to deal with equal boundaries??? is an acknowledged
unhandled case — see §"tie floor" below.

Note also the encoding: alignG = (chr_index << 32) | position, a single 64-bit genome-global
coordinate, so ordering across chromosomes is well defined.

The sample that defines the boundaries is a fixed ~50 MB and does not grow with the dataset.
source/BAMoutput.cpp:9-33:

    nBins=P.outBAMcoordNbins;
    binSize=P.chunkOutBAMsizeBytes/nBins;
    bamArraySize=binSize*nBins;
    ...
    binSize1=binStart[nBins-1]-binStart[0];
    nBins=1;//start with one bin to estimate genomic bin sizes

nBins starts at 1, so every alignment goes to bin 0 until that buffer (binSize1, ≈
chunkOutBAMsizeBytes*(nBins-1)/nBins) fills; then coordBins() runs
(BAMoutput.cpp:97-105) and re-bins retroactively. chunkOutBAMsizeBytes = limitIObufferSize[1]
(Parameters.cpp:1163), default 50000000 (parametersDefault:232). At ~160 B/record that is
~300,000 alignment records — for a 2 M-record run or a 2.44 G-record run, identically.

P.outBAMsortingBinStart is global and computed once, by whichever thread fills first, under a
mutex, guarded by outBAMsortingBinStart[0]!=0 (Parameters.cpp:655 initialises it to 1). So the
boundaries come from one thread's first ~50 MB, i.e. effectively the head of the FASTQ.

Dobin confirms this verbatim in issue #289 (2017-07-07):

"my suspicion is that your FASTQ files were derived from coordinate-sorted BAM files. This would
explain why STAR needs so much RAM for sorting - it decides on the genomic bin size for sorting
based on the first 100000 reads."

Bin assignment — source/BAMoutput.cpp:86-93:

        bamIn32=(uint32*) bamIn;
        alignG=( ((uint) bamIn32[1]) << 32 ) | ( (uint)bamIn32[2] );
        if (bamIn32[1] == ((uint32) -1) ) {//unmapped
            iBin=P.outBAMcoordNbins-1;
        } else if (nBins>1) {//bin starts have already been determined
            iBin=binarySearch1a <uint64> (alignG, P.outBAMsortingBinStart, (int32) (nBins-1));
        };

Mapped reads occupy bins 0 .. nBins-2 (that is nBins-1 bins); the last bin is reserved for
unmapped reads. This is why every memory formula below divides by nBins-1, not nBins.

The tie floor. source/serviceFuns.cpp:239-263, binarySearch1a "returns the last X element
that is <= x", and explicitly walks forward over duplicates:

    while (i1<N-1 && x==X[i1+1])
        ++i1; //go forward to check for equals

Consequence, proven by code: if a single genomic start coordinate holds more than 1/(nBins-1)
of the sample, several consecutive quantile boundaries take the same value; the bins between them
are empty (skipped by binS==0 continue, bamSortByCoordinate.cpp:59) and every record at that
coordinate lands in one bin. So:

max_bin_records >= (number of alignment records sharing the single most common start coordinate)

That is a hard floor no --outBAMsortingBinsN can go below. It is the mechanism behind
TyberiusPrime's report in issue #413 (2025-06-16): a targeted-panel scRNA-seq where "STAR wants
142 gigs of RAM for sorting, no matter that I set --outBAMSortingBinsN to 6400."

2. What --limitBAMsortRAM is compared against — the LARGEST BIN

source/bamSortByCoordinate.cpp:14-33.
permalink

        //check max size needed for sorting
        uint maxMem=0;
        for (uint32 ibin=0; ibin<nBins-1; ibin++) {//check all bins
            uint binS=0;
            for (int it=0; it<P.runThreadN; it++) {//collect sizes from threads
                binS += RAchunk[it]->chunkOutBAMcoord->binTotalBytes[ibin]+24*RAchunk[it]->chunkOutBAMcoord->binTotalN[ibin];
            };
            if (binS>maxMem) maxMem=binS;
        };
        ...
        P.inOut->logMain << "Max memory needed for sorting = "<<maxMem<<endl;
        if (maxMem>P.limitBAMsortRAM) {
            ostringstream errOut;
            errOut <<"EXITING because of fatal ERROR: not enough memory for BAM sorting: \n";
            errOut <<"SOLUTION: re-run STAR with at least --limitBAMsortRAM " <<maxMem+1000000000;

if (binS>maxMem) maxMem=binS; — this is a max over bins, summed across mapping threads within
each bin. It is not the total and not a running allocation. ibin<nBins-1 excludes the
unmapped bin.

Per-record cost: binTotalBytes accumulates bamSize + sizeof(uint) (BAMoutput.cpp:113), i.e. the
BAM record plus an 8-byte read index; the check adds 24 more bytes/record for the sort key triple.
So memory/record = BAM_record_bytes + 32, and this matches the actual allocation exactly
(BAMbinSortByCoordinate.cpp:11-12: new char[binS+1] + new uint[binN*3]).

3. What "Max memory needed for sorting" means

It is maxMem = the largest single bin, as printed at bamSortByCoordinate.cpp:28. Not the
total. The number in the error message is that same value + 1 GB of slack
(bamSortByCoordinate.cpp:32).

4. Is there any spill / adaptive / second-pass path? — No.

source/BAMbinSortByCoordinate.cpp:7-52 is the entire sort:

    //allocate arrays
    char *bamIn=new char[binS+1];
    uint *startPos=new uint[binN*3];
    ...
    //sort
    qsort((void*) startPos, binN, sizeof(uint)*3, funCompareArrays<uint,3>);

Whole bin into RAM, single qsort, stream out via bgzf, delete[]. No chunked merge, no temp
spill, no re-binning, no retry with more bins. The only "adaptivity" anywhere is the one-shot
coordBins() re-bin of the first buffer.

--limitBAMsortRAM has a second, distinct role as a concurrency budget
(bamSortByCoordinate.cpp:64-77):

                    uint newMem=binS+binN*24;
                    bool boolWait=true;
                    while (boolWait) {
                        #pragma omp critical
                        if (totalMem+newMem < P.limitBAMsortRAM) {
                            boolWait=false;
                            totalMem+=newMem;
                        };
                        sleep(0.1);
                    };

So actual peak RSS during sorting ≈ min(limitBAMsortRAM, Σ concurrently-sorted bins), whereas
the refusal threshold is max-bin. Raising --limitBAMsortRAM past what is needed therefore lets
STAR consume more. (Two incidental defects: sleep(0.1) truncates to sleep(0) — a busy spin at
100% CPU; and if maxMem == limitBAMsortRAM exactly, the > check passes but the strict < gate
never opens — an infinite hang. Neither should bite you if you leave headroom.)

The unmapped bin is not a memory risk: BAMbinSortUnmapped.cpp:19-78 is a streaming k-way
merge over 2*nThreads file handles with one BAMoutput_oneAlignMaxBytes buffer each — O(1) in
record count. --outSAMunmapped Within costs nothing here.

5. What the manual says

doc/STARmanual.pdf §17 (Output: SAM/BAM), option list — verbatim and complete:

--outBAMsortingBinsN
default: 50
int: >0: number of genome bins for coordinate-sorting

The manual states no memory relationship at all. (It also says "genome bins", which is the
misleading phrase that probably seeded your equal-genomic-length hypothesis — the code does quantiles.)

The only other manual text is under --outSAMtype ... SortedByCoordinate:

"sorted by coordinate. This option will allocate extra memory for sorting which can be specified by
--limitBAMsortRAM."

and --limitBAMsortRAM:

"default: 0 / int>=0: maximum available RAM (bytes) for sorting BAM. If =0, it will be set to the
genome index size. 0 value can only be used with --genomeLoad NoSharedMemory option."

The memory claim lives in CHANGES.md:371 (STAR 2.6.0a, 2018/04/23), where the option was
introduced:

"--outBAMsortingBinsN option to control the number of sorting bins. Increasing this number
reduces the amount of RAM required for sorting."

6. Known issues and what Dobin recommends

The most quantitative statement, issue #457, alexdobin, 2020-03-21:

"If the reads in your file are randomized, which is normally the case for raw files that come out
of the sequencer, then the max RAM required for sorting can be approximated as the size of
unzipped fastq files (sum of two mates for PE) divided by the --outBAMsortingBinsN value.
However, sometimes the reads do not get distributed uniformly into the bins, so to be safe, I would
double --outBAMsortingBinsN."

That formula is exactly total / nBins and is the author confirming binning divides the requirement.

Issue #2070 ("STARsolo BAM sorting required too much RAM", user needed 450 GB). alexdobin,
2024-02-23:

"Were these FASTQ files preprocessed in any way? If not, then it means there is a locus in the
genome that contains most of the reads, which will require large RAM for sorting."

alexdobin, 2024-03-01 — the decisive diagnosis, and a precise description of the code path above:

"It appears the reads are ordered for some reason. The first ~70k reads all map on chr1 at
~1,000,000 base, so STAR makes a lot of sorting bins around that locus, but most of the reads end
up in the last bin outside that locus, and so take a lot of RAM and time to sort."

Outcome, reported by the user 2024-03-04: FASTQ had been converted from a BAM; after
fastq-shuffle.pl, RAM dropped from 400 GB+ to 60 GB. Note this user was already running
--outBAMsortingBinsN 400 and it had not helped — same symptom as yours.

Issue #289, alexdobin, 2017-07-10 — the diagnostic recipe:

"The way to check whether the reads are sorted is to map the first ~100,000 reads
(--readMapNumber 100000) with default parameters, and check which chromosomes are present in the
SAM output, i.e. cat Aligned.out.sam | grep -v "^@" | cut -f3 | sort | uniq -c. If one
chromosome is dominant, the reads are likely to be sorted and STAR sorting should not be used."

Issue #1762, alexdobin, 2023-02-21 (STARsolo, FASTQ from a CellRanger BAM):

"The original BAM was sorted, which caused trouble for STAR BAM sorting. You have a couple of
options: 1). Generate Unsorted BAM with STAR, then sort it with samtools. 2). Shuffle the original
BAM file before generating BAMs with samtools collate."

Issue #1136, alexdobin, 2021-02-11: same diagnosis ("FASTQs recreated from sorted BAM files"),
same two remedies.

Issue #870, alexdobin, 2020-04-10 — the pileup variant:

"You can try to increase --outBAMsortingBinsN from the default 50 to 100 or even 200, though it
will require a large number of open files. It seems like there a few loci that are very highly
expressed, which STAR sorting cannot deal nicely with... So it may be a safer option to use the
samtools sort."
(The reporter later confirmed: "I was able to solve this using --outBAMsortingBinsN.")

Issue #1792, alexdobin, 2023-03-17 — directly on your exact situation (STARsolo CB/UB forcing the
internal sort):

"Unfortunately, at the moment STAR can output CB/UB tags only into sorted BAM files. You can try to
reduce RAM usage by increasing --outBAMsortingBinsN to 100 or 200 from the default 50."

Issue #413, alexdobin, 2018-05-04 — no spill path, and none planned then:

"robust sorting is high on my TODO list, it should be doable with some extra space and time for some
marginal cases. Presently, I would recommend scaling up --outBAMsortingBinsN parameter from the
default value of 50, which should solve most of the problems with sorting at the cost of increased
number (not size) of temp files."

That issue is still open as of the 2025-06-16 comment. Nothing shipped.

Issue #942, alexdobin, 2020-06-16 and #1923, 2023-10-13: the general fallback is
--outSAMtype BAM Unsorted + samtools sort — which for you is blocked by CB/UB.

7. Interaction with --outBAMsortingThreadN and thread count

source/Parameters.cpp:648-653:

                if (outBAMsortingThreadN==0) {
                    outBAMsortingThreadNactual=min(6, runThreadN);
                } else {
                    outBAMsortingThreadNactual=outBAMsortingThreadN;
                };
                outBAMcoordNbins=max((uint32)outBAMsortingThreadNactual*3,outBAMsortingBinsN);

Three real interactions:

• A floor on bins. nBins >= 3 * sortingThreads. --outBAMsortingThreadN 20 silently forces

  =60 bins. It can only raise the bin count, never lower it, so it never hurts max-bin.

• Peak RSS scales with sorting threads, max-bin does not. The refusal threshold is max-bin
  (thread-count independent), but actual concurrent allocation is up to
  outBAMsortingThreadNactual bins at once, gated only by limitBAMsortRAM
  (bamSortByCoordinate.cpp:64-77). More sorting threads = higher real RSS, same refusal
  threshold. If you are RSS-constrained rather than refusal-constrained, lower
  --outBAMsortingThreadN.
• File descriptors, not memory, are the cost of many bins. BAMoutput.cpp:27 opens one
  ofstream per bin per mapping thread in the constructor. The mapping-stage buffer
  (bamArraySize) is fixed at chunkOutBAMsizeBytes regardless of nBins — bins cost zero extra
  mapping RAM. Dobin, issue STAR invents introns longer than any that exist: no map module sets alignIntronMax, and the default is mammal-calibrated #457, 2018-08-03:

  "the number of temp files used for sorting scales as (Nthreads+1)*Nbins, and default Linux
  ulimit -n is usually set at 1024, so with 12 threads, you hit the limit somewhere around 80."
  His fix: ulimit -n 10000 before STAR. Without it, high bin counts segfault in
  outBAMwriteHeader (issue STAR invents introns longer than any that exist: no map module sets alignIntronMax, and the default is mammal-calibrated #457 traceback).

---

Does it explain our measurement?

Your pileup hypothesis: refuted as stated, and the arithmetic says it is not the cause either

Two independent reasons.

(a) The mechanism is wrong. Boundaries are count-quantiles (BAMoutput.cpp:139), not genomic
length. A 44%-of-reads pileup spread over ~10 kb attracts ~44% of the boundaries too. With 50 bins,
~21 of the 48 cuts land inside the rDNA array, spaced roughly
10 kb / 21 ≈ 480 bp apart. The pileup gets subdivided, exactly as designed.

(b) The magnitude is wrong. Only an exact-coordinate tie collapses bins (§1, tie floor), and
the floor is the multiplicity of the single most common start position. At your 2 M-record test with
44% rDNA, that is ~880 k records spread over ~10 kb ≈ ~90 records per coordinate, i.e. a floor of
~14 KB. Even 100x concentration gives ~1.4 MB. That is three orders of magnitude below your measured
394 MB.

What your numbers actually say

Take your data at face value, assuming you read "Max memory needed for sorting" from Log.out (or
subtracted the 1 GB slack from the error message) — either way, that number is max-bin:

records	measured max-bin	bytes/record	what working 50-bin quantile binning predicts
1,999,909	394 MB	197	394 MB / 49 ≈ 8.0 MB
9,844,534	1,590 MB	161	1590 MB / 49 ≈ 32 MB

Your measured max-bin is ~49x larger than the prediction — i.e. it equals the TOTAL. Max-bin ≈
total is the arithmetic signature of one bin holding essentially every mapped record. And
bytes/record = 161-197 is precisely the expected BAM_record + 32 for a single un-split pool — it
is a total-bytes rate, not a max-bin rate.

Total collapse plus insensitivity to nBins (50 ≈ 200 ≈ 1000) is the exact fingerprint Dobin
describes in #2070: boundaries all crammed into one small window seen at the head of the file,
everything else past the last boundary in one terminal bin. Adding boundaries just subdivides the
already-empty head window more finely — which is why 1000 bins was marginally worse, not better
(the sample buffer binSize1 = bamArraySize*(nBins-1)/nBins grows slightly with nBins, and the
terminal bin still absorbs everything).

Most likely root cause: your FASTQs are in genome-coordinate order. For 10x data this is
overwhelmingly common — public 10x submissions are frequently deposited as CellRanger BAMs and
regenerated with bamtofastq/cellranger bamtofastq, which preserves the coordinate-sorted order.
This is exactly issues #289, #1136, #1762 and #2070.

This is good news for you: it is a property of the file, not of the biology, and it is fixable by
shuffling — the #2070 user went 400 GB -> 60 GB.

Caveat I must flag. I am inferring root cause from the ratio of your numbers to what the code
predicts. I have not seen your Log.out. There is a free, zero-compute check that settles it
outright — see "cheapest decisive experiment" below.

---

What it means for a 2.44-billion-record sample

Using your own 160 B/record (which is BAM_record + 32, consistent with the code):

If binning works (random read order): max-bin ≈ 2.44e9 * 160 / (nBins-1)

--outBAMsortingBinsN	max-bin ("Max memory needed for sorting")
50 (default)	~8.0 GB
100	~3.9 GB
200	~2.0 GB
500	~0.8 GB

Your 390 GB naive extrapolation is not the linear-in-records law — it is precisely the
bins-fully-collapsed number. With working bins the sort is a non-issue at any sample size.

But the BAM sort is not your binding constraint. STARsolo allocates, and holds live across the
sort:

• SoloFeature_countCBgeneUMI.cpp:16 — readInfo.resize(nReadsInput, ...) where
  nReadsInput = g_statsAll.readN+1 (SoloFeature_sumThreads.cpp:11), i.e. total input reads.
  readInfoStruct is {uint64 cb; uint32 umi;} (SoloCommon.h:7-10) = 16 B after padding.
  → 2.23e9 x 16 B ≈ 36 GB.
• SoloFeature_countCBgeneUMI.cpp:21 — rGeneUMI = new uint32[rguStride*nReadsMapped], with
  rguStride=3 whenever readInfo is on (ParametersSolo.cpp:436 sets readIndexYes = readInfoYes)
  → ~24 GB at ~2e9 mapped reads.

Neither is counted by --limitBAMsortRAM, and both must be resident during sorting, because
SoloFeature::addBAMtags reads readInfo[iread] per record while writing each bin
(SoloFeature_addBAMtags.cpp:13-17; called from BAMbinSortByCoordinate.cpp:70-71). Ordering is
fixed in STAR.cpp: genome freed at :247, soloMain.processAndOutput() at :256, then
bamSortByCoordinate(...) at :272.

Practical recipe for the big sample:

--outSAMtype BAM SortedByCoordinate
--outBAMsortingBinsN 200          # ~2 GB max-bin; Dobin's #457 advice is "double it to be safe"
--outBAMsortingThreadN 4          # caps concurrent bins; also forces nBins >= 12
--limitBAMsortRAM 20000000000     # 20 GB: refusal headroom AND the concurrency budget

plus, in the same shell before STAR:

ulimit -n 8192                    # runThreadN * 200 descriptors; #457 segfaults without this

Total node RAM to request: ~90-110 GB (≈36 GB readInfo + ~24 GB rGeneUMI + ~20 GB sort budget

• genome/index + slack). The genome index is freed before sorting (STAR.cpp:247) so it does not
  stack with the sort, but it does stack with mapping.

Two hard ceilings to be aware of at this scale:

1. 2^32 reads. The read index is packed into the top 32 bits of a 64-bit word:
   outBAMcoord->coordOneAlign(..., (iReadAll<<32) | (iTr<<8) | ...)
   (ReadAlign_outputAlignments.cpp:200, also :228, :247) and recovered as iread >> 32
   (SoloFeature_addBAMtags.cpp:8-9). Above 4,294,967,295 input reads the CB/UB tags silently
   attach to the wrong reads. At 2.23e9 you have ~1.9x headroom — fine, but do not concatenate two
   such samples into one STAR run.
2. No escape hatch in 2.7.11b. ParametersSolo.cpp:405-411 hard-refuses:

   if ( (pP->outSAMattrPresent.CB || pP->outSAMattrPresent.UB) && type!=SoloTypes::CB_samTagOut) {
       samAttrYes=true;
       if (!pP->outBAMcoord) {
           errOut << "EXITING because of fatal PARAMETERS error: CB and/or UB attributes in --outSAMattributes can only be output in the sorted BAM file.\n";

   Your premise is correct and there is no --soloAddTagsToUnsorted in any released STAR (I grepped
   the whole tree; 2.7.11b is both the latest release and master HEAD). The one partial exemption is
   --soloType CB_samTagOut, which emits CR/CB into unsorted BAM but does no UMI counting and
   explicitly forbids UB (ParametersSolo.cpp:413-415) — so it is not a substitute if you need
   UB or the count matrices.

---

What remains uncertain, and the cheapest decisive experiment

Uncertain:

1. Whether your 394 MB / 1590 MB figures are Log.out's Max memory needed for sorting (max-bin) or
   something else (peak RSS, /usr/bin/time). My whole root-cause inference rests on them being
   max-bin. If they are peak RSS they mean something quite different and the analysis needs redoing.
2. Whether your FASTQs are coordinate-ordered. Strongly implied by the arithmetic; not directly
   observed by me.
3. Whether your rDNA array additionally has a single-coordinate tie large enough to matter. The
   code proves such a floor exists; only your data can say how big it is.

Cheapest decisive experiment — free, no new run, settles #1 and #2 immediately:

grep the Log.out you already have from the 9.8 M-record run:

grep -A60 "BAM sorting bins genomic start loci" Log.out
grep "BAM sorting:" Log.out            # sample size used to set boundaries
grep "Max memory needed for sorting" Log.out

STAR prints every bin boundary as index<TAB>chr_index<TAB>position (BAMoutput.cpp:140). Read
the list:

• All boundaries inside one small genomic window / one chromosome -> read-order pathology
  confirmed (Dobin #2070). Fix by shuffling; binning will then work and the monster sample is
  routine. This is my prediction.
• Many boundaries with identical chr+position -> exact-coordinate tie collapse. Binning
  cannot help below that floor; use samtools sort externally, or accept the floor.
• Boundaries spread sensibly across all six chromosomes -> binning worked, and then
  Max memory needed for sorting should have been ≈ total/49. If it was not, my model of the code is
  wrong and I would want to see the file.

Also worth grepping while you are there:
grep "Allocated and initialized readInfo array" Log.out — confirms nReadsInput and lets you
verify the 16 B/read figure against real RSS.

Second-cheapest — Dobin's own test, ~2 minutes (issue #289), if Log.out is gone:

STAR --genomeDir <ce11> --readFilesIn <R2> <R1> --readMapNumber 100000 \
     --outSAMtype SAM --outFileNamePrefix ordercheck.
grep -v '^@' ordercheck.Aligned.out.sam | cut -f3 | sort | uniq -c | sort -rn | head

Dominance by one chromosome (or, better, cut -f3,4 showing monotonically increasing positions)
= reads are coordinate-ordered.

If you want the controlled bin-count discrimination anyway — the experiment that separates
"pileup" from "read order" without touching the monster sample: take one small real sample and run it
three ways at --outBAMsortingBinsN 50 and 500: (a) as-is, (b) after shuf-pairing / samtools collate the reads, (c) after removing the rDNA locus from the reference. Compare
Max memory needed for sorting.

• (b) collapses the number and (a) does not -> read order. (My prediction.)
• (c) collapses it and (b) does not -> pileup / exact-coordinate ties.
• 500 bins helps in (b) roughly 10x versus 50 bins -> binning is healthy and the 2.44 G sample is
  safe with --outBAMsortingBinsN 200.

That is a ~10-minute experiment on a 2 M-record sample and it answers the operational question
completely.

---

---

Round 2: what the real Log.out proves, and what it revises

New observation. STAR 2.7.11b, ce11, GSE208154 / SAMN29720279 lane L001, 12,993,522 reads,
--runThreadN 24, default --outBAMsortingBinsN 50.
Max memory needed for sorting = 2801566182 (2.80 GB). All 48 logged boundaries on chrIdx 0
(chr I, 15,072,434 bp in ce11), in two clusters:

boundaries  1-37 : chr I  3,885,713 - 4,500,516     (615 kb window)
boundaries 38-48 : chr I 15,062,683 - 15,069,919    (7 kb window = tandem rDNA array)

The read-order theory is CONFIRMED, not refuted — the read-name argument was a false negative

The counter-argument was: read names are Illumina tile coordinates
(K00125:218:HCL2KBBXY:1:2106:15818:42829), therefore the file is in instrument order, therefore not
bamtofastq-derived.

That inference does not hold. Read names carry no information about record order. bamtofastq,
samtools fastq and SRA's own BAM→FASTQ conversion all copy the original read-name string verbatim.
A FASTQ regenerated from a coordinate-sorted BAM has original Illumina tile names and
genome-coordinate record order. The two are independent.

Direct external confirmation. ENA Portal API, result=read_run for SAMN29720279:

run_accession  sample_accession  submitted_format  submitted_ftp                                      fastq_ftp
SRR20172067    SAMN29720279      BAM;BAI           .../SRR20172067/possorted_genome_bam_TC2_d15_1.bam.1   (empty)

The submitter uploaded possorted_genome_bam...bam — CellRanger's position-sorted BAM, named
as such — and fastq_ftp is empty, i.e. no FASTQ was ever submitted. Every FASTQ for this run is
a conversion of a coordinate-sorted BAM. This is exactly Dobin's #289 / #1136 / #1762 / #2070
scenario, confirmed from the archive rather than inferred.

(Incidentally 1:2106:... is a surface-2 tile; a genuinely instrument-ordered file starts at tile
1101. Weak evidence on its own, but it points the same way.)

Q1. What exactly is the boundary sample? — one thread's first ~50 MB, as 2-3 disjoint read windows

REVISES round 1, which said "the first ~50 MB of BAM output ... effectively the head of the
FASTQ". More precisely, proven from code:

• One BAMoutput per mapping thread — ReadAlignChunk.cpp:51:
  chunkOutBAMcoord = new BAMoutput (iChunk, P.outBAMsortTmpDir, P);
• Each thread starts with one bin — BAMoutput.cpp:33: nBins=1;//start with one bin to estimate genomic bin sizes. All its alignments accumulate in bin 0's buffer.
• Buffer capacity is binSize1, independent of nBins — BAMoutput.cpp:11-13, 32:
  binSize=P.chunkOutBAMsizeBytes/nBins; bamArraySize=binSize*nBins; ... binSize1=binStart[nBins-1]-binStart[0];
  which is (chunkOutBAMsizeBytes/nBins)*(nBins-1) ≈ chunkOutBAMsizeBytes = limitIObufferSize[1]
  = 50,000,000 for every nBins >= 2. At ~140 B/record in-buffer that is ~350,000 records.
• The first thread to fill it decides for the entire run — BAMoutput.cpp:122-137:

      if (P.runThreadN>1) pthread_mutex_lock(&g_threadChunks.mutexBAMsortBins);
      if (P.outBAMsortingBinStart[0]!=0) {//it's set to 0 only after the bin sizes are determined

  with Parameters.cpp:655 initialising outBAMsortingBinStart[0]=1. The winner computes, sets
  [0]=0, and the other 23 threads contribute nothing — they only re-bin their own buffers
  against the winner's boundaries.
• coordFlush() is per-thread-per-run, not per-chunk — ReadAlignChunk_processChunks.cpp:255-259,
  outside the while (!noReadsLeft) chunk loop:

      };//cycle over input chunks
      if (P.outFilterBySJoutStage!=1 && RA->iRead>0) {
          if (P.outBAMunsorted) chunkOutBAMunsorted->unsortedFlush();
          if (P.outBAMcoord) chunkOutBAMcoord->coordFlush();

  So the only alternative sampling path is a thread that produced <50 MB in the whole run.
• Which reads land in that buffer — chunk intake is serialised, so threads take contiguous
  chunks in file order (ReadAlignChunk_processChunks.cpp:18,22):

      if (P.runThreadN>1) pthread_mutex_lock(&g_threadChunks.mutexInRead);
      ...
      while (chunkInSizeBytesTotal[0] < P.chunkInSizeBytes && ...

  chunkInSizeBytes ≈ limitIObufferSize[0]/readNends = 30 MB/2 = 15 MB per mate
  (Parameters.cpp:1161-1162) ≈ ~60,000 reads per chunk. Filling 50 MB of BAM takes ~2-3 chunks,
  and with 24 threads round-robining, those chunks sit at scattered offsets within roughly the first
  ~72 chunks (~4 M reads) of the file.

Net: the sample is 2-3 disjoint windows of ~60,000 consecutive reads, taken from the head third of
the file, by one thread. For an instrument-ordered file each window is a genome-wide random sample
and everything works. For a coordinate-ordered file each window is one contiguous genomic
interval — and you get exactly as many genomic clusters as the winning thread consumed chunks.

Observed: 2 clusters. Predicted: 2-3. That is the fingerprint. Round 1 predicted "one tiny
window", which was wrong in detail — the round-robin chunk hand-off produces several windows, not one.

Q2. Why only chr I? — the sample is not representative, and the cuts are working correctly

The cuts are doing exactly what we think. The decisive property is that a boundary value is
literally an element of the sample (BAMoutput.cpp:139):

            P.outBAMsortingBinStart[ib]=startPos[binTotalN[0]/(nBins-1)*ib];

startPos holds only sampled alignG values. So the boundary set can never contain a coordinate
that no sampled alignment occupies.

Read the observed boundaries as sample quantiles — boundary ib sits at rank ib/49:

rank	genomic span	share of sample	share of library
0 - 2.0%	< chr I:3,885,713	2.0%	—
2.0 - 75.5%	chr I:3,885,713 - 4,500,516 (615 kb)	73.5%	—
75.5 - 77.6%	chr I:4.50 - 15.06 Mb (10.5 Mb)	2.0%	—
77.6 - 98.0%	chr I:15,062,683 - 15,069,919 (7 kb, rDNA)	20.4%	—
98.0 - 100%	> chr I:15,069,919 — all of chr II/III/IV/V/X/MtDNA	≤2.0%	56%

Boundary 48 is at sample rank 48/49 = 97.96% and sits at chr I:15,069,919, i.e. 2.5 kb from the
end of chr I. So ≥98% of the sampled alignments are on chr I, while chr I is 44% of the library.
The sample over-represents chr I by ~2.2x and, critically, has ≤2% mass on the 85 Mb / 56% of records
that live beyond it.

That is exactly the two-contiguous-windows signature: 73.5% of the sample in one 615 kb interval and
20.4% in the 7 kb rDNA — proportions that are impossible for a genome-wide sample and natural for two
contiguous slices of a coordinate-sorted file.

Answer: the sample is not representative of the file. The equal-rank cuts are not at fault.

Q3. The 11 boundaries inside 7 kb are a SYMPTOM, and are in fact correct behaviour

Two things follow from the code.

1. No tie-collapse is occurring. Round 1 §"tie floor" warned that duplicate boundary values cause
   binarySearch1a to dump everything into one bin (serviceFuns.cpp:261,
   while (i1<N-1 && x==X[i1+1]) ++i1;). The 11 rDNA boundaries are at 11 distinct positions
   spanning 15,062,683-15,069,919. No duplicates ⇒ no collapse ⇒ those 11 bins are real and each
   takes a genuine slice of the rDNA. That failure mode is not what is happening here.
2. Dense boundaries around a pileup are the design working. A region holding a large share of
   records deserves many bins — that is the whole point of equal-count quantiles. If the sample were
   representative, 11 bins on the rDNA would be correct and desirable.

The damage is not the cluster; it is the absence of any boundary beyond chr I:15,069,919. That
single gap forces 56% of the library into one bin. 0.56 × total ≈ 2.80 GB observed ⇒ implied total
sort footprint ≈ 5.0 GB, i.e. ~160 B/record over ~31 M records — internally consistent.

Q4. Remedy — and a REVISION to the total/(nBins-1) formula

Why more bins gave no benefit on the earlier runs — now provable

Every boundary is an element of the sample, and binarySearch1a returns "the last X element that is
<= x". Therefore, if max(startPos) = chrI:15,069,919, then for any nBins:

every alignment with alignG > chrI:15,069,919 — i.e. all of chr II..MtDNA — returns index
nBins-2, the terminal mapped bin.

More bins can only subdivide the coordinate range the sample already covers. They cannot create a
cut where the sample has no mass. Raising nBins also cannot enlarge the sample, because
binSize1 ≈ chunkOutBAMsizeBytes regardless of nBins (Q1). That is the mechanical proof behind
"200 bins identical, 1000 slightly worse".

REVISION to §"Answer up front" item 2 and the 2.44 G table

Round 1 stated max-bin ≈ total/(nBins-1). That formula silently assumes a representative sample.
The correct statement, with K = nBins-1 mapped bins, s = share of the sample lying in a genomic
region, L = share of the library in that same region:

max-bin ≈ total × max over regions of ( L / (s·K) ) — a region gets s·K boundaries but must
absorb L of the data.

When the sample is representative, s = L and this collapses to total/K — round 1's formula, and
Dobin's #457 formula. When it is not, the worst region dominates. Here, for the tail beyond chr I,
s ≈ 0.02 and L = 0.56:

--outBAMsortingBinsN	boundaries landing beyond chr I (s·K)	predicted max-bin
50 (observed)	0.98 → 1 bin holds all 56%	2.80 GB ✓ (matches the log)
200	~4	~0.70 GB
1000	~20	~0.14 GB

So more bins probably WOULD help this particular run — roughly 4x at 200 and up to ~20x at 1000 —
but the ceiling is set by s, which is a property of read order, not of nBins. Note the earlier
"200/1000 bins don't help" measurements were taken on the other, smaller runs (~2 M and ~9.8 M
records), not on this one; if those had s ≈ 0 the flat response is expected. This is a testable
prediction and the cheapest next experiment (below).

The actual lever, in order

1. Restore random read order. This is the fix. The input is a converted possorted_genome_bam;
   un-sort it and s becomes L, the formula collapses to total/K, and default 50 bins is ample.
   • Best: samtools collate -u -O possorted_genome_bam.bam | <bamtofastq> — operate on the BAM you
     already have, no re-download. Dobin's own recommendation, #1762: "Shuffle the original BAM file
     before generating BAMs with samtools collate."
   • Or shuffle the FASTQ pairs (the #2070 user used fastq-shuffle.pl; 400 GB → 60 GB).
   • samtools collate is name-grouping, which is not a uniform shuffle, but it fully destroys the
     coordinate correlation that causes this — which is all that is required.
2. Raise --outBAMsortingBinsN to 200-1000 as a cheap partial mitigation if shuffling is
   impossible. Bounded benefit as tabulated; requires ulimit -n >= runThreadN * nBins (STAR invents introns longer than any that exist: no map module sets alignIntronMax, and the default is mammal-calibrated #457).
3. Enlarge the boundary sample with the second value of --limitIObufferSize (that is
   chunkOutBAMsizeBytes, Parameters.cpp:1163), e.g. --limitIObufferSize 30000000 500000000.
   A 10x sample reads 10x further into the file. Costs runThreadN × that value — at 24 threads
   that is 12 GB of buffers, and on a coordinate-sorted file it still only reaches ~10x deeper rather
   than genome-wide. Poor value; listed for completeness.
4. samtools sort externally — still blocked by CB/UB (ParametersSolo.cpp:405-411).

Is this unavoidable for any skewed-expression RNA-seq library? — No, and that matters

No. Expression skew alone does not trigger this. The failure requires the boundary sample to
miss part of the genome. With instrument-ordered reads, ~350,000 records sampled from anywhere in
the file is a genome-wide sample; equal-rank cuts then track the expression skew correctly, and a
pileup simply receives proportionally many bins (precisely what the 11 correct rDNA boundaries in Q3
demonstrate). Deep sequencing makes the sample more reliable, not less, because the sample size is
constant while the estimate is a quantile.

The only skew-driven failure that survives a random read order is the exact-coordinate tie floor
(round 1 §1): max-bin >= records sharing the single most common start coordinate. That is a real but
much smaller effect, and it is what TyberiusPrime hit in #413 with a targeted panel.

However, there is a systemic finding worth carrying into seqforge. The trigger — a dataset whose
reads were re-derived from a coordinate-sorted BAM — is common in public 10x data, because
submitters routinely upload possorted_genome_bam.bam rather than FASTQ. It is also detectable
without downloading a byte: ENA's filereport exposes submitted_format and submitted_ftp
directly, as demonstrated above. For a compiler processing ~10⁴ public datasets, a check of the form
"submitted_format contains BAM, or submitted_ftp basename matches *sorted*bam*" flags the affected
datasets up front and lets the recipe insert a samtools collate step. That is a far more actionable
result than a per-sample RAM tuning exercise.

Revised uncertainty, and the cheapest next experiment

Now settled: boundaries come from one thread's fixed ~50 MB buffer as 2-3 contiguous read windows
(Q1, code); the cuts are correct and the sample is not (Q2, code + log arithmetic); the rDNA cluster
is a symptom and no tie-collapse is occurring (Q3, distinct boundary values); more bins cannot create
cuts outside the sampled range (Q4, code); the input derives from a position-sorted BAM (ENA
submitted_format).

Still open: the exact value of s (sample mass beyond chr I). The boundary list only bounds it at
≤2%; it cannot distinguish "2% on chr II+" from "2% still on chr I". s is what decides whether
remedy 2 buys 4-20x or nothing.

Cheapest decisive experiment — one number from the run you already have:

grep "BAM sorting:" Log.out          #  "BAM sorting: <N> mapped reads"  = binTotalN[0], the sample size

That is printed at BAMoutput.cpp:134, immediately above the boundary list, and gives the sample size
directly. If it is ~350,000, the model above is confirmed end to end. Please capture this — it costs
nothing and it is the one number missing.

Cheapest experiment that discriminates the remaining question (~15 min, this same 13 M-read lane):

Re-run at --outBAMsortingBinsN 1000 (with ulimit -n 32768), change nothing else, and read
Max memory needed for sorting:

• ≈0.1-0.2 GB → s ≈ 2%, binning helps ~20x, remedy 2 is viable as a stopgap.
• still ≈2.8 GB → s ≈ 0: the sample contains zero alignments beyond chr I, binning is
  structurally incapable of helping, and shuffling is the only fix.

Either outcome is decisive, and neither requires the 2.23 G-read sample.

Confirmatory run (~20 min), worth doing once to validate the recipe for the whole corpus:
samtools collate the source BAM → regenerate FASTQ → re-run the same lane at default 50 bins.
Prediction: boundaries spread across all six chromosomes plus MtDNA, and
Max memory needed for sorting ≈ 5.0 GB / 49 ≈ 0.10 GB (a ~27x reduction). If that holds, the
2.44 G-record sample needs ≈ 2.44e9 × 160 / 49 ≈ 8 GB for sorting at default settings — and the
binding constraint reverts to STARsolo's readInfo/rGeneUMI arrays (~60 GB), exactly as §"What it
means for a 2.44-billion-record sample" concluded.

---

Round 3 — the prediction, measured

Round 2 predicted more bins would help this run substantially (~20x at 1000 bins). Measured on the
same lane, that is wrong. Same input, same 24 threads, only --outBAMsortingBinsN changed:

bins	boundary sample	max memory needed	vs 50 bins
50 (default)	142,760 reads	2,801,566,182 (2.80 GB)	—
1000	205,513 reads	2,007,006,631 (2.01 GB)	1.40x

A 20x increase in bins bought 1.4x, not 20x. And the reason is visible in the log: 998 of the
1000 boundaries are still on chrIdx 0 (chromosome I) — the same confinement as at 50 bins. Nothing
past the last chr-I boundary gets a cut, so the terminal bin still swallows chromosomes II/III/IV/V/X/
MtDNA, i.e. 56% of the library.

The small gain is fully explained by the sample growing (142,760 -> 205,513 reads), because buffer
capacity binSize1 = (chunkOutBAMsizeBytes/nBins)*(nBins-1) rises slightly as nBins grows. A bigger
sample subdivides chr I a little more finely. It does not, and cannot, place a boundary where the
sample has no mass.

Conclusion, measured rather than modelled: --outBAMsortingBinsN is not the lever for this data.
It is a real and correctly-implemented knob — it bounds max-bin when the boundary sample is
representative — but here the sample is confined to one chromosome by read order, and no bin count
repairs that. The fix has to act on the input: shuffle (samtools collate), or avoid the condition by
detecting it at resolve time.

What this means for seqforge

The condition is detectable for free, before any download. ENA's filereport exposes what the
submitter actually uploaded; verified independently for this run:

$ curl -s "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=SRR20172067&result=read_run\
&fields=run_accession,submitted_format,submitted_ftp,fastq_ftp&format=tsv"
run_accession  submitted_format  submitted_ftp                          fastq_ftp
SRR20172067    BAM;BAI           possorted_genome_bam_TC2_d15_1.bam.1   (empty)

submitted_format contains BAM, fastq_ftp is empty, and the file is CellRanger's
position-sorted BAM — so every FASTQ for this run is a conversion of a coordinate-sorted file.
seqforge io resolve already queries ENA, so this is one extra field on a call that already happens.

Sizing resources.mem_gb per sample treats the symptom; flagging genome-ordered input at resolve time
and collating treats the cause, and is cheaper for a corpus of ~10^4 datasets.

[lhqing]

Resolution — what we decided, and what we did NOT do

Closing this out. The investigation above is correct and worth keeping, but the decision went a
different way than its final section implies, so read that section as analysis rather than
plan.

Decided:

1. Keep CB/UB in the CRAM (option A, as shipped in starsolo: the retained CRAM carries the barcode, and the counts become CellRanger-comparable #200). No switch to raw CR/UR tags,
   even though those avoid the sort entirely.
2. Keep Velocyto in soloFeatures.
3. Escalate memory on retry instead of sizing every job for the worst sample — starsolo: escalate memory on retry, so a job that does not fit fails loudly instead of being OOM-killed #205.
4. A job that still does not fit fails loudly, and that is an accepted outcome.

NOT adopted, despite being where the analysis pointed:

• Shuffling genome-ordered input (samtools collate), and detecting it at resolve time from
  ENA's submitted_format. This remains the correct fix for the root cause — the boundary-sample
  pathology documented above — and the detection is free. It was judged more machinery than the
  problem currently justifies. If sort memory becomes a recurring operational cost rather than a
  single-sample annoyance, this is the thing to build, and the analysis is already done.
• Deriving --limitBAMsortRAM from read count. More correct than the 3/4-of-request heuristic,
  but it needs read_count plumbed from the archive record into compose, and retry escalation
  covers the same ground for less.

One correction to the round-2/round-3 analysis

The report attributes the per-read readInfo array to the CB/UB request. That is only one of
two independent triggers — ParametersSolo.cpp sets readInfoYes from Velocyto as well:

if (featureYes[Velocyto] || featureYes[VelocytoSimple])  readInfoYes[Gene] = true;
if (samAttrYes)                                          readInfoYes[samAttrFeature] = true;

Our soloFeatures has always included Velocyto, so seqforge has been paying 16 B x every input
read since well before #198. CB/UB consume that array; they did not cause it. This matters for
anyone revisiting the CR/UR idea: raw tags avoid the sort, not readInfo.

For comparison: what scRecounter does

Same problem, different constraints (config/process.config):

withLabel:process_high { cpus = 8; memory = { 72.GB * task.attempt }; maxRetries = 3 }
errorStrategy { task.attempt <= maxRetries ? 'retry' : 'ignore' }

They sidestep the sort completely with --outSAMtype None — no BAM, no bins, --limitBAMsortRAM
irrelevant. They still pay readInfo (they run Velocyto too), and handle it with 72 GB escalating to
216 GB, then 'ignore' — the sample is silently dropped from the corpus.

Two things worth taking: 8 CPUs to 72 GB is the right shape for this workload (memory-heavy,
not core-heavy), and retry escalation is the mechanism we adopted in #205. The 'ignore' is not
available to us — silently dropping a dataset is defensible for a count-matrix corpus, not for a
compiler whose contract is an exit code.