Skip to content

Fetch, cache and display electron density and cryo-EM maps - #1134

Open
aalhossary wants to merge 8 commits into
biojava:masterfrom
aalhossary:aa/electron-density-maps
Open

Fetch, cache and display electron density and cryo-EM maps#1134
aalhossary wants to merge 8 commits into
biojava:masterfrom
aalhossary:aa/electron-density-maps

Conversation

@aalhossary

Copy link
Copy Markdown
Member

Closes #947.

This PR is about a new feature, therefore, I recommend releasing it as BioJava 7.3.0 rather than 6.2.7.
It adds org.biojava.nbio.structure.io.density, which downloads and caches density maps the way LocalPDBDirectory caches coordinate files, and a JmolPanel method that contours the result.

DensityMapCache cache = new DensityMapCache();
DensityMapResult map = cache.getDensityMap(new PdbId("1cbs"), DensityMapKind.TWO_FO_FC);
jmolPanel.loadDensityMap(map);

Why a chain of sources rather than one

edmaps.rcsb.org shut down in October 2024. What RCSB documents in its place —the map coefficients published with the wwPDB validation reports — are structure factors, not a sampled grid, and cannot be displayed without a Fourier transform.
Several other services do serve grids, and they differ enormously in size for the same entry, so rather than picking one, sources are tried in order until one answers:

order
X-ray RCSB density server → PDBe CCP4 → PDBe density server → wwPDB coefficients (off by default)
cryo-EM RCSB density server → PDBe density server → EMDB primary map

Measured sizes for the same data:

entry full archive density server
1cbs (X-ray) 2.1 MB (two CCP4 files) 210 kB at detail 0
EMD-0262 (6hu9) 116 MB 0.5–3.7 MB

wwPDB coefficients are supported for completeness but disabled by default, and DensityFileFormat.isJmolLoadable() lets a viewer refuse them rather than silently drawing nothing.

Design notes

  • A source that has nothing for an entry is skipped; any other transport failure aborts the chain. A network outage must never be reported as "this entry has no density". When every source is exhausted, NoDensityMapException carries the reason from each one, so the UI can explain rather than just fail.
  • A density server response contains both the 2Fo-Fc and Fo-Fc blocks, so the two kinds share one cache entry instead of downloading the identical file twice.
  • Cryo-EM entries are resolved through EMDB's search API, with RCSB as a fallback. That lookup also yields the depositors' recommended contour level, which is how an EM map should be contoured; it is attached to the result whichever source supplied the voxels. The experimental method is never inferred from resolution, which BioJava parses incorrectly for some cryo-EM entries (Resolution lost in cryo-EM #1000).
  • Ccp4Header checks for the MAP stamp at byte 208, so a server answering with an error page and HTTP 200 produces a clean cache miss rather than a corrupt entry.
  • A .meta sidecar fully describes each cached result, so LOCAL_ONLY is served without opening a connection.

Three Jmol behaviours established by experiment

Each of these changed the implementation, and none is documented:

  1. Option order in isosurface is load-bearing. With mesh nofill before the file name, Jmol accepts the command, reports no error, and draws nothing at all.
  2. A negative sigma does not contour at a negative level. Jmol reserves negative sigma internally, so sigma -3.0 silently contours at the default level — the intended red lobe of a difference map came out identical to the blue one. Difference maps are drawn as a single signed surface instead.
  3. Selecting the Fo-Fc block of a cached BinaryCIF file needs the marker in the file name. BCifDensityReader picks the block by testing whether the file name contains &diff=1, which normally arrives in the URL query string. A cached local file has no query string; appending the marker to the file URL fails, and the #diff=1 form the reader's own to-do comment suggests is ignored. Embedding it in the name works, so the cache exposes the difference map under a companion name, hard-linked to the same bytes.

Measured proof of (3): the plain file contours at cutoff 0.356 over a −1.31 to 3.78 range; the companion name contours at 0.374 over −0.69 to 0.85. Different data, so the marker genuinely selects the other block.

Testing

43 unit tests, entirely offline — the chain tests use stub providers, so "404 falls through" and "anything else aborts" are pinned exactly. 6 integration tests against the live services costing a couple of megabytes in total; the cryo-EM path is covered without downloading the 116 MB map by setting the size limit so the guard fires after resolution but before transfer.

Jmol version

No change: BioJava stays on net.sourceforge.jmol:jmol:14.31.10, which is the newest on Maven Central (October 2020). Every behaviour relied on here — MrcBinaryReader, BCifDensityReader, the &diff=1 selection, gzip sniffing — is identical in 14.31.10 and current Jmol, so downstream projects overriding the dependency with a newer build are unaffected.

The download-validation helpers added in 7.0.0 (biojava#979, biojava#980) had several gaps
that only surface once a caller passes a real hash URL or downloads from a
server that can 404. All of them are fixed here, and hash verification is
implemented rather than stubbed.

Correctness fixes:

* createValidationFiles(URL, ...) passed a literal Hash.UNKNOWN to its
  URLConnection overload instead of the caller's argument, so it could never
  write a hash file and threw IllegalArgumentException for any caller that
  supplied a hashURL.
* Neither downloadFile nor createValidationFiles checked the HTTP status, so a
  404 error page was written into the cache as though it were the requested
  file. Because the .size sidecar was then taken from that same error response,
  validateFile subsequently declared it valid. A new HttpStatusException lets
  callers tell "the resource is not there" apart from a transport failure,
  which matters for anything that tries several mirrors in turn.
* downloadFile used FileChannel.transferFrom(rbc, 0, Long.MAX_VALUE), which is
  not guaranteed to drain a socket-backed channel and could silently truncate a
  download. Replaced with Files.copy, which loops to end of stream.
* downloadFile leaked its temporary file on every failure path.
* validateFile threw NullPointerException for a file with no parent directory,
  and again if listFiles() returned null; an empty .size file raised an
  unchecked NoSuchElementException that escaped the surrounding catch.
* validateFile checked only the first hash sidecar it found, ignoring the rest.

New functionality:

* validateFile now really verifies MD5, SHA-1 and SHA-256 instead of throwing
  UnsupportedOperationException. Sidecars are written as bare lowercase hex and
  parsed tolerantly, so a file downloaded verbatim from a server in coreutils
  or BSD layout is also understood. A sidecar that cannot be parsed is skipped
  with a warning rather than failing an otherwise good download.
* ETagPolicy lets an ETag that is a bare hex digest be recorded as a checksum
  without a second request. files.wwpdb.org and files.rcsb.org return the
  content MD5 as the ETag, so every download from the wwPDB archive now gets a
  real integrity check for free. The <mtime>-<size> ETags used by the EBI
  servers contain a dash and can never be misread as a digest; a test pins that.
* downloadFileWithValidation downloads and validates over a single connection.
  The previous pattern opened one connection to read Content-Length and another
  to fetch the bytes, so the recorded size described a different response than
  the one written; if the resource changed in between, the cache entry was left
  permanently failing validation. Content is digested while streaming to a temp
  file and only moved into place once length and checksum check out.

LocalPDBDirectory uses the new single-connection download, and its two-character
directory hash is promoted to the reusable getMiddleHash(String). That hash
deliberately counts from the end of the identifier so that both spellings of an
entry land in the same bucket: 1cbs and pdb_00001cbs both give "cb", where
counting from the start would file the extended form under "db".

Defaults change slightly: the existing four-argument createValidationFiles
overloads now use ETagPolicy.USE_IF_HEX_DIGEST, so callers start recording
checksums where the server offers one. Targeted at 7.3.0.
Exercises the new status checking against the real wwPDB archive: downloadFile
must throw HttpStatusException rather than writing the error page to the
destination, and createValidationFiles must not record a .size for it. Without
the second half the cached error page would pass validateFile, since its
recorded size would match the error body exactly.
It is a general 'has the server got a newer copy' helper with nothing
PDB-specific about it, and other caching code outside this package needs the
same Last-Modified comparison. Protected access only reached subclasses and
the io package itself.
Adds org.biojava.nbio.structure.io.density, which downloads and caches density
maps the same way LocalPDBDirectory caches coordinate files, and hands back a
File that a viewer can contour. Closes biojava#947.

Several services publish density for the PDB and they differ enormously in size
for the same entry, so rather than picking one, sources are tried in order until
one answers. The order is smallest-adequate-first, because the smallest form is
usually perfectly good to look at:

  X-ray:  RCSB density server -> PDBe CCP4 -> PDBe density server
          -> wwPDB map coefficients (disabled by default)
  cryo-EM: RCSB density server -> PDBe density server -> EMDB primary map

For 1cbs a density server slice is about a tenth the size of the equivalent pair
of CCP4 files. For the map behind 6hu9 it is 3.7 MB against a 106 MB primary
map. A size limit, 256 MiB by default, is checked against the size EMDB itself
reports before any of the body is transferred, and exceeding it is not an error:
the chain simply falls back to a smaller representation.

wwPDB map coefficients are supported for completeness, since they are the route
RCSB documents now that edmaps.rcsb.org has shut down, but they are structure
factors rather than a sampled grid and cannot be displayed without a Fourier
transform. They are therefore disabled by default, and DensityFileFormat carries
an isJmolLoadable() flag so that a viewer can refuse them rather than silently
drawing nothing.

Notes on the design:

* A source that has nothing for an entry is skipped and the next is tried, but
  any other transport failure aborts the chain. A network outage must never be
  reported as "this entry has no density". When every source is exhausted,
  NoDensityMapException carries the reason from each one, so a caller can say
  why rather than just that it failed.
* A density server response contains both the 2Fo-Fc and the Fo-Fc blocks, so
  the two kinds share one cache entry instead of downloading the identical file
  twice. Which block to read is a display-time decision.
* Cryo-EM entries are found through their EMDB identifier, looked up from EMDB's
  search API with RCSB as a fallback. That lookup also yields the contour level
  the depositors recommend, which is how an EM map should be contoured; it is
  attached to the result whichever source supplied the voxels. The experimental
  method is never inferred from resolution, which BioJava parses incorrectly for
  some cryo-EM entries (biojava#1000).
* Ccp4Header checks for the MAP stamp at byte 208, so a server that answers with
  an error page and HTTP 200 produces a clean cache miss rather than a corrupt
  cache entry.
* Cached results are fully described by a .meta sidecar, so LOCAL_ONLY requests
  are served without opening a connection.

DemoFetchElectronDensity exercises all three outcomes: an X-ray entry, a cryo-EM
entry resolved through EMDB, and an entry deposited without structure factors.
Adds loadDensityMap and clearDensityMaps to JmolPanel, the single point every
viewer in the module already goes through, so a map fetched by DensityMapCache
can be contoured with one call. Surfaces are given stable ids so they can be
addressed or removed individually, and the state that "Reset Display" restores
is re-saved afterwards; otherwise the button would silently discard the map the
user had just asked for.

Three things had to be established by experiment against Jmol 14.31.10 rather
than assumed, and each changed the implementation:

* Option order in the isosurface command is load-bearing. With "mesh nofill"
  placed before the file name, Jmol accepts the command, reports no error, and
  draws nothing whatsoever. It has to follow the file name.

* A negative sigma does not contour at a negative level. Jmol reserves negative
  sigma for its own internal signalling, so "sigma -3.0" silently contours at
  the default level instead: the intended red lobe of a difference map came out
  identical to the blue one. Difference maps are therefore drawn as a single
  signed surface, which also matches Jmol's own shortcut for them.

* A density server response carries both the 2FO-FC and FO-FC blocks, and Jmol
  chooses between them by testing whether the file NAME contains "&diff=1". That
  marker normally arrives in the URL query string, which a cached local file does
  not have. Appending it to the file URL does not work, and neither does the
  "#diff=1" form the reader's own comment suggests: both were measured and both
  returned the 2FO-FC block or nothing. Embedding the marker in the file name
  does work, so the cache exposes the difference map under a companion name,
  hard-linked to the same bytes where the filesystem allows it.

Verified by contouring the real cached files headlessly: the 2Fo-Fc map at 1
sigma gives cutoff 0.356 over a -1.31 to 3.78 range, and the difference map at 3
sigma gives 0.374 over -0.69 to 0.85 - a different data block, which is what
proves the marker works rather than merely being accepted.

A map that cannot be contoured without a Fourier transform is rejected with an
explanatory exception rather than producing an empty surface.
Puts the feature in reach from the alignment viewer's View menu. The fetch runs
on a SwingWorker: even the smallest source is a few hundred kilobytes and a
full-resolution map can be far larger, so fetching on the event dispatch thread
would freeze the window for the duration.

Requests are built with allowNonRenderableFormats(false), which keeps the map
coefficient source out of the chain automatically rather than relying on the
viewer to notice it cannot draw the result.

When nothing is available the dialog explains why rather than listing HTTP
codes: for an entry whose every source returns 404 the likely reason is that no
structure factors were deposited and there is no associated EMDB map, which is
worth saying plainly. The per-source detail is still shown underneath.

AbstractAlignmentJmol gains getFrame() and setStatus() so a listener in the
neighbouring package can own its dialogs and report progress; both were
previously reachable only as protected fields.

DemoShowElectronDensity displays 1CBS with both maps clipped around the bound
retinoic acid.
Forty-three unit tests run with no network at all. The chain tests use stub
providers, so the behaviour that matters can be pinned exactly: a 404 falls
through to the next source, a too-large map falls through as well, and anything
else aborts. That last one is the point of the design - reporting a dropped
connection as "this entry has no density" would be worse than failing.

The offline set also pins the two-character directory rule against both
spellings of an entry (1cbs and pdb_00001cbs must land in "cb", not "db"), that
every source and kind combination maps to a distinct file, and that a LOCAL_ONLY
request is served entirely from disk with every server pointed at a dead port.

Six integration tests exercise the real services for a couple of megabytes in
total. The cryo-EM path is covered without downloading the 116 MB map: the EMDB
entry is resolved, the author contour level checked against its known value, and
the size guard then declines the full map before any of its body is transferred.
The coefficient test corrupts a downloaded file afterwards to confirm the
ETag-derived MD5 actually catches it rather than merely being recorded.

Writing the header check turned up a real bug: isCcp4 rejected any file shorter
than a CCP4 header, but a gzipped map compresses to a small fraction of the
header it contains, so small EMDB maps would have been rejected as invalid. The
length shortcut is gone; reading decides it.
@aalhossary

Copy link
Copy Markdown
Member Author

Worth discussing: how small is small enough?

This PR deliberately prefers the smallest adequate representation, and I would like a second opinion on where the line should sit.

The density servers return a downsampled grid, and the saving is not marginal. For 1cbs, detail 0 is 210 kB against 2.1 MB for the two full CCP4 files — and one request instead of two, since a single response carries both map kinds. For the cryo-EM map behind 6hu9 it is 0.5–3.7 MB against 116 MB. The current defaults are detail 3 for RCSB and detail 6 for PDBe, matching what each service's own viewer asks for.

Open questions:

  1. Is a downsampled grid adequate for what people actually use BioJava density for?
    For looking at a ligand pocket in a viewer I am confident it is. For anything quantitative — real-space correlation, occupancy refinement, automated density fit scoring — it may well not be, and someone doing that would want the full grid and would not necessarily notice they had been handed a coarser one. The result object does report its source, so the information is available, but "available if you check" is weaker than a sensible default.

  2. Should the default detail level be lower?
    Detail 0 for 1cbs is a tenth the size of detail 3 and still perfectly readable around a ligand. I kept the services' own defaults rather than choosing for them, but I am not attached to that.

  3. Should the size ceiling be the primary mechanism instead of the ordering?
    Right now the 256 MiB ceiling rarely fires: EMD-0262 is 116 MB and would sail under it. What actually avoids the 116 MB download is putting the density servers first. An alternative design would always prefer the full map and fall back only when it is genuinely huge, which is more predictable but much more expensive by default.

  4. Should a quantitative caller get a different default than a viewer?
    There is already allowNonRenderableFormats distinguishing the two audiences; a preferFullResolution flag would be a small addition if the distinction is real.

Happy to change any of this — the ordering is one list and a couple of setters.

The <table summary="..."> form is obsolete: the summary attribute was removed
in HTML5, and javadoc has generated HTML5 since JDK 15, so doclint rejects it
whenever it is switched on. The build sets -Xdoclint:none so this never broke
CI, but it would surface in the release profile and the attribute does nothing
for accessibility any more.

A definition list suits a list of placeholders and their meanings better than a
two-column table in any case.
@aalhossary

Copy link
Copy Markdown
Member Author

Reviewer note: This PR (when approved) should be merged after #1133; because the other PR is included, a prerequisite, and a dependency of this PR.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Fetch, cache, load, and view electron density maps

1 participant