MDAnalysis / GridDataFormats

GridDataFormats is a pure Python library to handle data on a regular grid using commonly used file formats in molecular simulations.
https://mdanalysis.org/GridDataFormats
GNU Lesser General Public License v3.0
29 stars 18 forks source link

Impossible reading .ccp4 data #122

Closed plongo-fzj closed 1 year ago

plongo-fzj commented 1 year ago

I'm trying to read a .ccp4 file. This file can be handled by ChimeraX, PyMOL and VMD nevertheless whenever I try to use gridData I get an error. Here's the code I use:

from gridData import Grid
grid = Grid('water_density.ccp4')

And here's the error I get

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 grid = Grid('water_density.ccp4')

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/gridData/core.py:241, in Grid.__init__(self, grid, edges, origin, delta, metadata, interpolation_spline_order, file_format)
    238         filename = str(grid)
    240 if filename is not None:
--> 241     self.load(filename, file_format=file_format)
    242 else:
    243     self._load(grid, edges, metadata, origin, delta)

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/gridData/core.py:548, in Grid.load(self, filename, file_format)
    546     raise IOError(errno.ENOENT, "file not found", filename)
    547 loader = self._get_loader(filename, file_format=file_format)
--> 548 loader(filename)

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/gridData/core.py:559, in Grid._load_mrc(self, filename)
    557 def _load_mrc(self, filename):
    558     """Initializes Grid from a MRC/CCP4 file."""
--> 559     mrcfile = mrc.MRC(filename)
    560     grid, edges = mrcfile.histogramdd()
    561     self._load(grid=grid, edges=edges, metadata=self.metadata)

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/gridData/mrc.py:92, in MRC.__init__(self, filename)
     90 self.filename = filename
     91 if filename is not None:
---> 92     self.read(filename)

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/gridData/mrc.py:98, in MRC.read(self, filename)
     96 if filename is not None:
     97     self.filename = filename
---> 98 with mrcfile.open(filename) as mrc:
     99     if not mrc.is_volume():                           #pragma: no cover
    100         raise ValueError(
    101             "MRC file {} is not a volumetric density.".format(filename))

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/load_functions.py:139, in open(name, mode, permissive, header_only)
    137         elif start[:2] == b'BZ':
    138             NewMrc = Bzip2MrcFile
--> 139 return NewMrc(name, mode=mode, permissive=permissive,
    140               header_only=header_only)

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/mrcfile.py:115, in MrcFile.__init__(self, name, mode, overwrite, permissive, header_only, **kwargs)
    113         self._create_default_attributes()
    114     else:
--> 115         self._read(header_only)
    116 except Exception:
    117     self._close_file()

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/mrcfile.py:131, in MrcFile._read(self, header_only)
    129 """Override _read() to move back to start of file first."""
    130 self._iostream.seek(0)
--> 131 super(MrcFile, self)._read(header_only)

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/mrcinterpreter.py:176, in MrcInterpreter._read(self, header_only)
    174 self._read_extended_header()
    175 if not header_only:
--> 176     self._read_data()

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/mrcfile.py:139, in MrcFile._read_data(self)
    136 header_size = self.header.nbytes + self.header.nsymbt
    137 remaining_file_size = file_size - header_size
--> 139 super(MrcFile, self)._read_data(max_bytes=remaining_file_size)
    141 # Check if the file is the expected size.
    142 if self.data is not None:

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/mrcinterpreter.py:322, in MrcInterpreter._read_data(self, max_bytes)
    301 """Read the data array from the stream.
    302 
    303 This method uses information from the header to set the data array's
   (...)
    319         :data:`True`.
    320 """
    321 try:
--> 322     dtype = utils.data_dtype_from_header(self.header)
    323 except ValueError as err:
    324     if self._permissive:

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/utils.py:66, in data_dtype_from_header(header)
     47 """Return the data dtype indicated by the given header.
     48 
     49 This function calls :func:`dtype_from_mode` to get the basic dtype, and
   (...)
     63         mode.
     64 """
     65 mode = header.mode
---> 66 return dtype_from_mode(mode).newbyteorder(mode.dtype.byteorder)

File ~/miniconda3/envs/standard/lib/python3.11/site-packages/mrcfile/utils.py:172, in dtype_from_mode(mode)
    170     return np.dtype(_mode_to_dtype[mode])
    171 else:
--> 172     raise ValueError("Unrecognised mode '{0}'".format(mode))

ValueError: Unrecognised mode '33554432'

Again I don't believe it is a problem of the file itself since it doesn't create any problem for ChimeraX, VMD and PyMOL.

plongo-fzj commented 1 year ago

probably an error from the ccp4 file used although many visualization tools (ChimeraX, VMD, Pymol) could open it

IAlibay commented 1 year ago

@plongo-fzj sorry we didn't get around to looking at this issue. Could you expand on why you say it might be an error with the file itself?

plongo-fzj commented 1 year ago

I couldn't read it using another python library like gemmi: script:

import gemmi
gemmi.read_ccp4_map(filename)

error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[30], line 1
----> 1 gemmi.read_ccp4_map(filename)

RuntimeError: Unexpected axis value in word 17: 16777216

Unfortunately I couldn't compare my ccp4 with other ccp4 files built from other softwares therefore I cannot say that for sure and this is why I wrote 'probably'. I'll take the occasion to congrats with all developers for this amazing software that is MDAnalysis!

IAlibay commented 1 year ago

Thanks @plongo-fzj - I'm happy to go with gemmi as the sanity check that it's probably a file level issue. However, please do re-open this issue if/when needed!