Repository: rmarkello/abagen

Issue:
The hippocampal formation: cortex or subcortex?
## The issue

**tl;dr** The Allen Institute ontology classifies hippocampus as part of cortex, not subcortex, which could cause problems for matching some microarray samples to ROIs.

---

When users provide a file or dataframe to the `atlas_info` parameter in `abagen.get_expression_data()` they are required to specify a broad structural class for each region in their `atlas` (in a column labelled 'structure' in the file/dataframe). The current options for this structural class include:

1. cortex,
2. subcortex,
3. cerebellum,
4. brainstem,
5. white matter, and
6. other (i.e., ventricles and such)

We match these designations with the information from the Allen ontology such that samples that don't fall directly within a region in the `atlas` aren't incorrectly assigned to regions in `atlas` across hemispheric / structural boundaries. 

That is, if one of the samples from the Allen Institute is labelled as having come from the **left hemisphere subcortex** we make sure to only assign it to a region in the user-specified `atlas` labelled as belonging to the **left hemisphere subcortex**. This impacts only a minority of samples (i.e., we don't currently check whether this is the case for those samples having coordinates directly within a region in the atlas), but a significant minority, nonetheless.

While matching these designations seems like a reasonable approach in most cases, the one point of contention that a general user might have is that the Allen Institute ontology classifies the hippocampal formation (including the subiculum, dentate gyrus, and CA1-4) as part of "cortex" rather than "subcortex".  Specifically, their ontology specifies:

```
brain 
└─ gray matter
   └─ telencephalon
      └─ cerebral cortex
         └─ limbic lobe
            └─ hippocampal formation
```

Thus, if a researcher provides an atlas where they label all their hippocampal ROIs as "subcortex" they're liable to get vastly different results than if they label all their hippocampal ROIs as "cortex."

While I have it on good authority that the hippocampus is often considered part of "allocortex," I'm hesitant to add this as a permissible structural class to `abagen` since it seems quite a bit more specific than the current (rather broad) structural designations listed above (1-6).

## Proposed solution

I genuinely don't know! It would be _great_ to allow either specification for the hippocampus (i.e., "cortex" or "subcortex"), but the current framework for getting these structural classes from the Allen ontology doesn't allow for this hedging. I can think about how to modify it for this one instance in particular, but in the interim it would be great to come up with alternatives.

One option that might be worthwhile is to simply allow users to specify either (or both) of the expected 'hemisphere' and 'structure' information in `atlas_info` and just use whatever is available. Then, users who have hippocampal ROIs can refrain from specifying the 'structure' for their ROIs and we'll do our best to ensure samples simply don't cross hemispheric boundaries. This isn't necessarily ideal because there's the possibility that samples will get incorrectly assigned across e.g., cortical/subcortical boundaries for regions that _aren't_ the hippocampus (but we might still consider this option outside of the current problem!).

Alternatively (and perhaps most immediately appealing), we can add a warning on the documentation about this designation and inform users to specify that their hippocampal ROIs are part of "cortex" (not "subcortex") when they provide `atlas_info`.

Discussion:

Patch:
diff --git a/abagen/samples.py b/abagen/samples.py
index 38bda53..2341b5b 100644
--- a/abagen/samples.py
+++ b/abagen/samples.py
@@ -22,7 +22,8 @@ ONTOLOGY = Recoder(
      ('9131', 'pons', 'brainstem'),
      ('9512', 'myelencephalon', 'brainstem'),
      ('9218', 'white matter', 'white matter'),
-     ('9352', 'sulci & spaces', 'other')),
+     ('9352', 'sulci & spaces', 'other'),
+     ('4219', 'hippocampal formation', 'subcortex')),
     fields=('id', 'name', 'structure')
 )
 
@@ -122,11 +123,12 @@ def _get_struct(structure_path):
         Structure, or None if unable to identify a corresponding structure
     """
 
-    structure_path = set(structure_path.strip('/').split('/'))
-    ids = list(set(ONTOLOGY.value_set('id')) & structure_path)
+    structures = set(structure_path.strip('/').split('/'))
+    ids = sorted(set(ONTOLOGY.value_set('id')) & structures,
+                 key=lambda x: structure_path.index(x))
 
     try:
-        return ONTOLOGY.structure[ids[0]]
+        return ONTOLOGY.structure[ids[-1]]
     except IndexError:
         return
 
diff --git a/docs/user_guide/expression.rst b/docs/user_guide/expression.rst
index cea8279..ee85161 100644
--- a/docs/user_guide/expression.rst
+++ b/docs/user_guide/expression.rst
@@ -82,13 +82,13 @@ expression data normalized and aggregated across donors:
     >>> print(expression)
     gene_symbol      A1BG  A1BG-AS1       A2M  ...       ZYX     ZZEF1      ZZZ3
     label                                      ...
-    1            0.533064  0.792167  0.447196  ...  0.521263  0.224733  0.396009
-    2            0.421408  0.716915  0.460235  ...  0.419497  0.117964  0.558332
+    1            0.533876  0.764106  0.437582  ...  0.515566  0.211259  0.391993
+    2            0.426346  0.690736  0.438414  ...  0.418789  0.106776  0.546173
     3                 NaN       NaN       NaN  ...       NaN       NaN       NaN
     ...               ...       ...       ...  ...       ...       ...       ...
-    81           0.387847  0.045987  0.473161  ...  0.092883  0.596697  0.410640
-    82           0.796091  0.134167  0.372935  ...  0.140178  0.279248  0.716110
-    83           0.356067  0.022527  0.802499  ...  0.279927  0.886430  0.470536
+    81           0.385569  0.050425  0.409866  ...  0.087828  0.607935  0.367567
+    82           0.822947  0.154579  0.371812  ...  0.140344  0.288985  0.708077
+    83           0.356487  0.027261  0.804702  ...  0.282122  0.882226  0.458286
     <BLANKLINE>
     [83 rows x 15656 columns]
 


Test patch:
diff --git a/abagen/tests/test_samples.py b/abagen/tests/test_samples.py
index 044dfcb..d23508a 100644
--- a/abagen/tests/test_samples.py
+++ b/abagen/tests/test_samples.py
@@ -145,7 +145,7 @@ def test_drop_mismatch_samples(mm_annotation, ontology):
     expected = pd.DataFrame(dict(hemisphere=['L', 'R', np.nan],
                                  mni_x=[-10, 30, 0],
                                  structure_acronym=['S', 'Cl', 'CC'],
-                                 structure=['cortex', 'subcortex', 'other'],
+                                 structure=['subcortex', 'subcortex', 'other'],
                                  structure_id=[4251, 4323, 9422],
                                  structure_name=['subiculum, left',
                                                  'claustrum, right',