Replace the blocking blastp -remote call with NCBI's BLAST URL API and a bounded timeout

[mrubash1]

Stacked on #101 (which is stacked on #100). The two commits new to this PR are the last two; review from Replace blastp -remote with NCBI's URL API onward.

The problem

run_blast shells out to blastp -remote, which blocks indefinitely polling NCBI's shared public queue. Observed during this work: 0.31 s of CPU over 51 minutes, no open socket, and a 0-byte output file. Reproduced identically on macOS arm64 and on x86-64 Linux, so it is not a local or an architecture problem.

The cause is queue congestion, not the pipeline's arguments. NCBI's own time-to-completion estimate (RTOE) came back at 7838 s (131 minutes) for a single 375-residue actin query. -word_size, -outfmt, -max_target_seqs, and shell=True were each ruled out by probing RTOE against a control: the same parameters submitted six times consecutively returned RTOE=20 every time, versus 8655 for those same parameters minutes earlier. RTOE tracks global server load, not the query.

There is no way to bound how long blastp -remote will block. That is what this PR changes.

The change

run_blast now drives NCBI's BLAST URL API with requests (already a dependency of envs/blast.yml), keeping the nr database and the existing retry loop:

• Submit — CMD=Put with PROGRAM=blastp, DATABASE=nr, mapping max_target_seqs→HITLIST_SIZE, word_size→WORD_SIZE, evalue→EXPECT, plus tool and email so NCBI can identify the client.
• Poll — CMD=Get&FORMAT_OBJECT=SearchInfo, first after the returned RTOE, then no more often than every 60 seconds, as NCBI asks.
• Fetch — CMD=Get&FORMAT_TYPE=JSON2_S, converted to the pipeline's existing 17-column -outfmt 6 layout.

run_blast still returns a CompletedProcess whose stderr is prefixed Error: on failure, so blast_call_failed() and the word-size backoff added in #101 are unchanged.

Why JSON2_S rather than Tabular

FORMAT_TYPE=Tabular would need no conversion, but NCBI does not actually deliver it. Fetching a completed search (RID 7TJ2P2DR014) returns a 62-byte body containing only the status block and no hits, with or without FORMAT_OBJECT=Alignment. The same RID returns real results in every other format:

Tabular  :      62 bytes (status block only)
Text     : 157,935 bytes
XML2_S   : 511,934 bytes
JSON2_S  : 422,458 bytes

JSON2_S also carries the subject accession and the taxonomy id directly, so staxids is populated rather than dropped.

Timeout

blast_timeout_seconds (default 1800, i.e. 30 minutes) bounds the total time spent on one protein across all retry attempts, not each attempt: a single deadline is established before the retry loop and each attempt receives only what remains. A retry cannot reset the budget. On expiry the search fails with a message reporting how long it waited and what NCBI estimated.

The timeout bounds the waiting, not the requests: a poll already in flight when it expires runs to completion, so a call can overrun by up to one HTTP_TIMEOUT_SECONDS (60s). An earlier revision of this description claimed the timeout was honoured exactly, which was wrong.

Database

blast_database config parameter, defaulting to nr — what the pipeline has always searched, so no existing run changes. It is a config parameter rather than a constant because which database is searched decides which homologs are found, so it belongs with the other parameters that determine a map's contents.

The documented alternative is nr_cluster_seq (ClusteredNR), not one of the RefSeq-only databases: it is built by clustering nr, so it does not drop the GenBank CDS translations, PDB, PIR and PRF entries that refseq_protein excludes outright. Measured against NCBI on 2026-08-13, nr and refseq_protein received identical time estimates and refseq_protein finished no faster, so the smaller RefSeq databases do not appear to buy the speed they are usually chosen for.

Failure handling

Three failure modes that this PR originally got wrong, each now covered by a test that fails without its fix:

• A search whose estimate exceeded the budget was never polled at all. The deadline was re-checked between waiting and polling, so with the default 1800s timeout and a measured RTOE of 9168s the code slept the full 30 minutes, sent zero requests, and reported a timeout — for a search NCBI may well have already finished. The status is now always checked after waiting.
• A poll that failed to reach NCBI ended the search. One refused connection discarded the request id and everything already spent waiting, and the retry submitted a fresh search to the back of the queue. Poll failures are now retried, bounded by MAX_CONSECUTIVE_POLL_FAILURES and by the timeout. The search is queued on NCBI's side, so failing to ask about it says nothing about the search.
• An unusable results response was written as an empty results file. A report whose BlastOutput2 held no search, or a body that parsed as JSON but was not an object, produced an empty file and rc=0 — resurfacing downstream as the misleading "no hits were returned" that this PR exists to remove. Both are now reported as failures.

Other

• blast_email config parameter, overridable with PROTEINCARTOGRAPHY_BLAST_EMAIL. The default is a project address, not an individual's.
• resources: ncbi_remote=1 on the run_blast rule. Inert unless snakemake is given a budget for it (--resources ncbi_remote=1), which serializes BLAST so that a run with many input proteins does not submit concurrent searches to a shared queue. Documented in the README rather than forced on.

Output field fidelity

Converting the real JSON response reproduces the real blastp -remote row for a hit that both outputs share — 16 of 17 fields identical:

field                                                              match
qseqid, sseqid, pident, length, mismatch, gapopen                  OK
qstart, qend, sstart, send, evalue                                 OK
sacc, saccver, sgi, staxids, scomnames                             OK
bitscore                788 vs 787                                 differs

bitscore differs by less than 1 because the two searches ran against different nr snapshots, and the bit score depends on the effective search space. The test asserts this rather than hiding it.

scomnames is written as N/A. NCBI's JSON provides sciname (a scientific name), and the column is defined as the common name, so populating it would mislabel the data. sgi is written as 0, matching what the CLI writes for retired GIs. Only sacc is read downstream (extract_blast_hits.py:50).

What is and is not verified

Verified against live NCBI, end to end in a single run — on macOS arm64 and, independently, on x86-64 Linux (Modal):

Attempt 1/1 to run a blast search (using a word size of 5, with 14400s of the 14400s budget remaining)
NCBI queued the blast search as '7VPVAHA6014' and estimated that it will take 20s
RUN_BLAST_EXIT=0                                          # 84 seconds, macOS arm64

NCBI queued the blast search as '7VRWRVHH014' and estimated that it will take 31s
The blast search '7VRWRVHH014' is READY after 91s         # 93 seconds, x86-64 Linux

• Submit, RID/RTOE parsing, polling, fetch, and conversion, all in one invocation
• 10 hits written, every row with all 17 fields; staxids populated (89673, 61221, 9764, ...)
• The repo's own extract_blast_hits.py run against that output returns 10 accessions, several of which (RLW01512, XP_007129366, NP_001009784) match the accessions in the repo's committed id-mapping fixture
• The timeout path, verified separately: a search whose queue estimate was 7838 s failed at the configured deadline with a clear message rather than blocking
• Gapped alignments: a 500-hit search produced 5 rows with gapopen > 0 and 444 rows with mismatches, so the gap-open counting and the mismatch calculation are exercised against real NCBI output. Note the repo's committed blastresults.tsv artifact contains no gapped rows at all, so this is broader coverage than that fixture provides.

Not verified against live NCBI:

• Non-zero e-value formatting. Every hit in every real response obtained here had evalue=0, and the committed CLI artifact likewise contains only 0.0, so there is no real-world reference to compare against. format_evalue implements blastp's formatting thresholds and is unit-tested, but only against synthetic values. Note that evalue is not read downstream; only sacc is (extract_blast_hits.py:50), and that is verified above.

Test results

Run on macOS arm64 and independently on x86-64 Linux (Modal), both with CI's pinned tool versions:

74 passed, 3 warnings in 124.34s              # make test
26 of 26 steps (100%) done, 1 passed          # make smoke-test
ruff check .                                  # clean
ruff format --check .                         # 43 files left unchanged
snakefmt --check .                            # All 1 file(s) would be left unchanged

The 3 warnings are pre-existing deprecation notices from snakemake and stopit, not from this PR.

Unit tests mock every HTTP request and fake the clock, so multi-hour timeouts are exercised instantly and without sleeping. A trimmed 9.2 KB excerpt of a real NCBI response is committed as a fixture so the parser is tested against genuine output rather than only hand-written mocks; hits and descriptions were dropped, never edited.

Note for reviewers

This does not make BLAST faster. NCBI's queue is the bottleneck and remains so. What changes is that the pipeline waits a bounded, configurable amount of time and then fails with a clear message, instead of blocking indefinitely.

NCBI's congestion is intermittent, not permanent: the RTOE estimate for the same query was 7838 s on one day and 20 s the next. blastp -remote would therefore often work fine. The point of this change is the days when it does not, where it currently blocks for hours with no bound.

Alternative back ends were benchmarked (EBI's REST service returned 1000 hits against UniProtKB in ~8.5 minutes while NCBI returned nothing in any configuration) but deliberately not adopted: EBI has no nr, caps at 1000 hits against a configured max_blast_hits: 3000, and returns UniProt rather than RefSeq accessions. Changing which database is searched changes the results, which is a scientific decision rather than an engineering one.

[mrubash1]

Closing in favour of #107.

This PR replaced the BLAST transport wholesale. #103 has since merged its own www/QBlast implementation to main, so the choice is no longer between two proposals — it would now mean deleting working code. Keeping that transport and porting this PR's improvements onto it is the better trade.

#107 carries the substance across: the bounded total timeout, the budget spanning attempts, the consecutive-poll-failure cap, response validation, config parameters instead of environment variables, and a testable heartbeat. It also fixes a defect found while reading the merged implementation — _www_get treated any unrecognised response as READY, so a transient NCBI error page ended the search and was reported as a success with no hits.

What is deliberately not carried across is the JSON2_S parsing and the ~1200 lines of tests built around it. Those are validated row-for-row against real blastp CLI output and are stronger than the hand-rolled XML path, but swapping the parser is the wholesale replacement being avoided here. Worth revisiting only if the XML path is found to produce a wrong row.