Allen Common Coordinate Framework (CCF)¶

import datajoint as dj
import numpy as np
import pandas as pd
from pathlib import Path
import urllib.request

schema = dj.Schema('tutorial_allen_ccf')

DATA_DIR = Path('./data')
DATA_DIR.mkdir(exist_ok=True)

ONTOLOGY_URL = (
    "http://api.brain-map.org/api/v2/data/query.csv?"
    "criteria=model::Structure,rma::criteria,[ontology_id$eq1],"
    "rma::options[order$eq%27structures.graph_order%27][num_rows$eqall]"
)
ONTOLOGY_FILE = DATA_DIR / 'allen_structure_graph.csv'

if not ONTOLOGY_FILE.exists():
    print("Downloading Allen brain structure ontology...")
    urllib.request.urlretrieve(ONTOLOGY_URL, ONTOLOGY_FILE)
    print(f"Downloaded to {ONTOLOGY_FILE}")
else:
    print(f"Using cached {ONTOLOGY_FILE}")

ontology = pd.read_csv(ONTOLOGY_FILE)
print(f"Loaded {len(ontology)} brain regions")
ontology.head()

dj.Diagram(schema)

# Insert CCF metadata
CCF.insert1(
    {
        'ccf_id': 1,
        'ccf_version': 'CCFv3',
        'ccf_resolution': 25.0,
        'ccf_description': 'Allen Mouse CCF v3 (25µm resolution)'
    },
    skip_duplicates=True
)

CCF()

# Populate brain regions
BrainRegion.populate(display_progress=True)

# View sample regions
BrainRegion() & 'region_id < 100'

# Populate parent relationships
RegionParent.populate(display_progress=True)

# Find primary visual cortex
BrainRegion & {'acronym': 'VISp'}

# Get VISp region ID
visp = (BrainRegion & {'acronym': 'VISp'}).fetch1()
visp_id = visp['region_id']

# Find all children (direct descendants)
children = BrainRegion * (RegionParent & f'parent_id = {visp_id}' & f'region_id != {visp_id}')
children.proj('acronym', 'region_name')

def get_ancestors(region_acronym, ccf_id=1):
    """Get the path from a region to the root."""
    region = (BrainRegion & {'acronym': region_acronym, 'ccf_id': ccf_id}).fetch1()
    region_id = region['region_id']
    
    path = [region['acronym']]
    
    while True:
        parent_id = (RegionParent & {'ccf_id': ccf_id, 'region_id': region_id}).fetch1('parent_id')
        if parent_id == region_id:  # reached root
            break
        parent = (BrainRegion & {'ccf_id': ccf_id, 'region_id': parent_id}).fetch1()
        path.append(parent['acronym'])
        region_id = parent_id
    
    return ' → '.join(reversed(path))

# Show path from VISp layer 1 to root
print("Path from VISp1 to root:")
print(get_ancestors('VISp1'))

# Aggregate regions by depth in hierarchy
import matplotlib.pyplot as plt

depths = (RegionParent & 'ccf_id = 1').to_arrays('depth')
unique, counts = np.unique(depths, return_counts=True)

plt.figure(figsize=(10, 4))
plt.bar(unique, counts)
plt.xlabel('Depth in Hierarchy')
plt.ylabel('Number of Regions')
plt.title('Brain Regions by Hierarchy Depth')
plt.xticks(unique)

# Find all visual-related regions
(BrainRegion & "region_name LIKE '%Visual%'").proj('acronym', 'region_name')

@schema
class RecordingSite(dj.Manual):
    definition = """
    # Recording electrode location
    recording_id : int32
    ---
    ap : float32                    # anterior-posterior (µm from bregma)
    dv : float32                    # dorsal-ventral (µm from brain surface)
    ml : float32                    # medial-lateral (µm from midline)
    -> BrainRegion                  # assigned brain region (includes ccf_id)
    """

# Insert example recording sites
RecordingSite.insert([
    {'recording_id': 1, 'ccf_id': 1, 'ap': -3500, 'dv': 500, 'ml': 2500, 
     'region_id': (BrainRegion & {'acronym': 'VISp'}).fetch1('region_id')},
    {'recording_id': 2, 'ccf_id': 1, 'ap': -1800, 'dv': 1200, 'ml': 1500,
     'region_id': (BrainRegion & {'acronym': 'CA1'}).fetch1('region_id')},
], skip_duplicates=True)

# View recordings with region info
RecordingSite * BrainRegion.proj('acronym', 'region_name')

# Final schema diagram
dj.Diagram(schema)