rmarkello / abagen

A toolbox for working with Allen Human Brain Atlas microarray expression data
https://abagen.readthedocs.io
BSD 3-Clause "New" or "Revised" License
93 stars 41 forks source link

Donor-specific atlases and tolerance=0 fails for two donors #190

Closed rmarkello closed 3 years ago

rmarkello commented 3 years ago

The issue

When using donor-specific atlases and providing tolerance=0 to the abagen.get_expression_data() workflow, two of the donors ('9861' and '10021') have no matched samples. This is wrong and is due entirely to this line:

https://github.com/rmarkello/abagen/blob/ed4ce09d7e628b82a5f52d4c03f3fa61250a87e0/abagen/matching.py#L316

This line is important for when we're using a group-level atlas and trying to determine whether a sample (using its MNI coordinates) falls directly within a voxel of the provided atlas; unfortunately, the x-offset for the two aforementioned donors is a decimal (0.5), so the call to np.floor changes the coordinates such that they no longer align with the donor-specific atlases.

To reproduce:

>>> import abagen
>>> atlases = abagen.fetch_desikan_killiany(native=True)['image']
>>> expression = abagen.get_expression_data(atlases, tolerance=0, return_donors=True)
>>> np.all(np.isnan(expression[0]))
True
>>> np.all(np.isnan(expression[1]))
True

Proposed solution

The call to np.floor is a bit too basic. It assumes that all the provided atlases will always be on a regular grid (i.e., 1mm isotropic) and that the offsets will always be integer numbers—and while this may generally be the case, it is impossible to say that it will always be true.

The np.floor line should be removed entirely, and, perhaps somewhere in samples_.update_coords() where we first handle the coordinates, we should update the MNI coordinates such that they align with the provided atlas grid. This will require maintaining some information about the image grid for the relevant atlas(es) (or at the very least voxel sizes).

(This does raise another point about whether we should be rounding the MNI coordinates down to the nearest atlas voxel at all. I maintain that this is reasonable to do, but it is a design choice that could be argued against.)