{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":84969,"databundleVersionId":10033515,"sourceType":"competition"}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Overview\r\nThis notebook provides an exploratory data analysis of the CZII Cryo-ET Particle Detection Challenge dataset. It focuses on:\r\n- Loading and visualizing tomogram data\r\n- Analyzing particle distributions\r\n- Visualizing different particle types:\r\n  - apo-ferritin (easy)\r\n  - beta-amylase (impossible, not scored)\r\n  - beta-galactosidase (hard)\r\n  - ribosome (easy)\r\n  - thyroglobulin (hard)\r\n  - virus-like-particle (easy)","metadata":{}},{"cell_type":"markdown","source":"## Table of Contents\n1. Loading Required Libraries\n2. Reading Tomogram Data\n3. Loading Particle Coordinates\n4. Visualizing Particles in Tomogram Slices\n5. Statistical Analysis of Particle Distribution","metadata":{}},{"cell_type":"code","source":"!pip install -q zarr\n!pip install -q ome-zarr\n!pip install -q copick","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:16.015593Z","iopub.execute_input":"2024-11-07T01:49:16.016003Z","iopub.status.idle":"2024-11-07T01:49:52.407997Z","shell.execute_reply.started":"2024-11-07T01:49:16.015965Z","shell.execute_reply":"2024-11-07T01:49:52.406523Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Loading Required Libraries","metadata":{}},{"cell_type":"code","source":"# First, let's import the necessary libraries\nimport zarr\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom pathlib import Path","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:52.410756Z","iopub.execute_input":"2024-11-07T01:49:52.411214Z","iopub.status.idle":"2024-11-07T01:49:52.527293Z","shell.execute_reply.started":"2024-11-07T01:49:52.411169Z","shell.execute_reply":"2024-11-07T01:49:52.526007Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Reading Tomogram Data","metadata":{}},{"cell_type":"code","source":"# Define the path to the zarr file\nzarr_path = Path('/kaggle/input/czii-cryo-et-object-identification/train/static/ExperimentRuns/TS_5_4/VoxelSpacing10.000/denoised.zarr')\n\n# Open the zarr array\nzarr_store = zarr.open(str(zarr_path))\n\n# Print basic information about the zarr store\nprint(\"Zarr store structure:\")\nprint(zarr_store.tree())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:52.528988Z","iopub.execute_input":"2024-11-07T01:49:52.530083Z","iopub.status.idle":"2024-11-07T01:49:52.611319Z","shell.execute_reply.started":"2024-11-07T01:49:52.530020Z","shell.execute_reply":"2024-11-07T01:49:52.610080Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Get the highest resolution data (scale 0)\ntomogram = zarr_store['0'][:]\n\nprint(f\"Tomogram shape: {tomogram.shape}\")\nprint(f\"Data type: {tomogram.dtype}\")\nprint(f\"Min value: {tomogram.min()}\")\nprint(f\"Max value: {tomogram.max()}\")\nprint(f\"Mean value: {tomogram.mean()}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:52.612717Z","iopub.execute_input":"2024-11-07T01:49:52.613078Z","iopub.status.idle":"2024-11-07T01:49:53.337967Z","shell.execute_reply.started":"2024-11-07T01:49:52.613042Z","shell.execute_reply":"2024-11-07T01:49:53.335892Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Reading and printing the structure of the apo-ferritin JSON\nimport json\nfrom pathlib import Path\n\n# Read the JSON file for apo-ferritin\njson_path = Path('/kaggle/input/czii-cryo-et-object-identification/train/overlay/ExperimentRuns/TS_5_4/Picks/apo-ferritin.json')\n\nwith open(json_path, 'r') as f:\n    data = json.load(f)\n    \n# Examine the structure\nprint(\"Keys in the JSON file:\", data.keys())\nprint(\"\\nFirst few points:\")\nprint(json.dumps(data['points'][:2], indent=2))  # Print first 2 points for examination","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:53.342104Z","iopub.execute_input":"2024-11-07T01:49:53.342942Z","iopub.status.idle":"2024-11-07T01:49:53.355967Z","shell.execute_reply.started":"2024-11-07T01:49:53.342873Z","shell.execute_reply":"2024-11-07T01:49:53.352707Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Loading Particle Coordinates","metadata":{}},{"cell_type":"code","source":"import json\nimport numpy as np\nfrom pathlib import Path\n\ndef load_apo_ferritin_coordinates(experiment_name='TS_5_4'):\n    \"\"\"Load apo-ferritin coordinates from JSON file.\"\"\"\n    json_path = Path('/kaggle/input/czii-cryo-et-object-identification/train/overlay/ExperimentRuns') / experiment_name / 'Picks/apo-ferritin.json'\n    \n    try:\n        with open(json_path, 'r') as f:\n            data = json.load(f)\n            \n        # Extract coordinates from the points array\n        coords = []\n        for point in data['points']:\n            coords.append([\n                point['location']['z'],\n                point['location']['y'],\n                point['location']['x']\n            ])\n        \n        coords = np.array(coords)\n        print(f\"Loaded {len(coords)} apo-ferritin coordinates\")\n        \n        # Print some basic statistics\n        print(\"\\nCoordinate ranges:\")\n        print(f\"Z range: {coords[:, 0].min():.1f} to {coords[:, 0].max():.1f}\")\n        print(f\"Y range: {coords[:, 1].min():.1f} to {coords[:, 1].max():.1f}\")\n        print(f\"X range: {coords[:, 2].min():.1f} to {coords[:, 2].max():.1f}\")\n        \n        return coords\n        \n    except Exception as e:\n        print(f\"Error reading coordinates: {e}\")\n        return np.array([])\n\n# Load the coordinates\napo_ferritin_coords = load_apo_ferritin_coordinates()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:53.358041Z","iopub.execute_input":"2024-11-07T01:49:53.358749Z","iopub.status.idle":"2024-11-07T01:49:53.377365Z","shell.execute_reply.started":"2024-11-07T01:49:53.358688Z","shell.execute_reply":"2024-11-07T01:49:53.375790Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Tomogram shape:\", tomogram.shape)\nprint(\"\\nCoordinate ranges before scaling:\")\nprint(f\"Z range: {apo_ferritin_coords[:, 0].min():.1f} to {apo_ferritin_coords[:, 0].max():.1f}\")\nprint(f\"Y range: {apo_ferritin_coords[:, 1].min():.1f} to {apo_ferritin_coords[:, 1].max():.1f}\")\nprint(f\"X range: {apo_ferritin_coords[:, 2].min():.1f} to {apo_ferritin_coords[:, 2].max():.1f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:53.379466Z","iopub.execute_input":"2024-11-07T01:49:53.380117Z","iopub.status.idle":"2024-11-07T01:49:53.389748Z","shell.execute_reply.started":"2024-11-07T01:49:53.380055Z","shell.execute_reply":"2024-11-07T01:49:53.388564Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Let's create a function to scale the coordinates\ndef scale_coordinates(coords, tomogram_shape):\n    \"\"\"Scale coordinates to match tomogram dimensions.\"\"\"\n    scaled_coords = coords.copy()\n    \n    # Scale factors for each dimension\n    scale_z = tomogram_shape[0] / coords[:, 0].max()\n    scale_y = tomogram_shape[1] / coords[:, 1].max()\n    scale_x = tomogram_shape[2] / coords[:, 2].max()\n    \n    # Apply scaling\n    scaled_coords[:, 0] = coords[:, 0] * scale_z\n    scaled_coords[:, 1] = coords[:, 1] * scale_y\n    scaled_coords[:, 2] = coords[:, 2] * scale_x\n    \n    return scaled_coords\n\n# Scale the coordinates\nscaled_coords = scale_coordinates(apo_ferritin_coords, tomogram.shape)\n\nprint(\"\\nCoordinate ranges after scaling:\")\nprint(f\"Z range: {scaled_coords[:, 0].min():.1f} to {scaled_coords[:, 0].max():.1f}\")\nprint(f\"Y range: {scaled_coords[:, 1].min():.1f} to {scaled_coords[:, 1].max():.1f}\")\nprint(f\"X range: {scaled_coords[:, 2].min():.1f} to {scaled_coords[:, 2].max():.1f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:49:53.391337Z","iopub.execute_input":"2024-11-07T01:49:53.391784Z","iopub.status.idle":"2024-11-07T01:49:53.405537Z","shell.execute_reply.started":"2024-11-07T01:49:53.391743Z","shell.execute_reply":"2024-11-07T01:49:53.404305Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Visualizing Particles in Tomogram Slices","metadata":{}},{"cell_type":"code","source":"from mpl_toolkits.axes_grid1 import ImageGrid\n\n# Updated visualization function\ndef visualize_apo_ferritin(tomogram, coords, n_slices=3, slice_thickness=10):\n    \"\"\"\n    Visualize apo-ferritin particles in tomogram slices.\n    \"\"\"\n    fig = plt.figure(figsize=(20, 10))\n    grid = ImageGrid(fig, 111,\n                    nrows_ncols=(1, n_slices),\n                    axes_pad=0.3,\n                    share_all=True,\n                    cbar_location=\"right\",\n                    cbar_mode=\"single\",\n                    cbar_size=\"5%\",\n                    cbar_pad=0.1)\n    \n    # Normalize tomogram data\n    vmin, vmax = np.percentile(tomogram, (1, 99))\n    normalized_tomogram = np.clip((tomogram - vmin) / (vmax - vmin), 0, 1)\n    \n    # Calculate evenly spaced z-positions\n    z_positions = np.linspace(0, tomogram.shape[0]-1, n_slices, dtype=int)\n    \n    # Plot each slice\n    for idx, ax in enumerate(grid):\n        z = z_positions[idx]\n        \n        # Show tomogram slice\n        im = ax.imshow(normalized_tomogram[z, :, :], cmap='gray', vmin=0, vmax=1)\n        \n        # Find particles near this slice\n        mask = np.abs(coords[:, 0] - z) < slice_thickness\n        if np.any(mask):\n            ax.scatter(coords[mask, 2], coords[mask, 1],\n                      color='red', marker='o', s=100, \n                      facecolors='none', linewidth=2,\n                      label='apo-ferritin')\n        \n        ax.set_title(f'Slice Z={z}\\n({np.sum(mask)} particles visible)')\n        ax.grid(False)\n        \n        # Set the axes limits to match the tomogram dimensions\n        ax.set_xlim(0, tomogram.shape[2])\n        ax.set_ylim(tomogram.shape[1], 0)  # Inverted y-axis to match image coordinates\n    \n    # Add colorbar and title\n    grid.cbar_axes[0].colorbar(im)\n    \n    plt.suptitle('Apo-ferritin Particles in Tomogram Slices\\n' + \n                 f'Showing particles within ±{slice_thickness} units of each slice',\n                 fontsize=16, y=1.05)\n    \n    # Add legend to the first subplot\n    grid[0].legend(bbox_to_anchor=(1.5, 1.0))\n    \n    plt.show()\n    \n    # return fig","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:50:47.323384Z","iopub.execute_input":"2024-11-07T01:50:47.323856Z","iopub.status.idle":"2024-11-07T01:50:47.351782Z","shell.execute_reply.started":"2024-11-07T01:50:47.323813Z","shell.execute_reply":"2024-11-07T01:50:47.350689Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Create the visualization with scaled coordinates\nvisualize_apo_ferritin(tomogram, scaled_coords)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:50:49.336322Z","iopub.execute_input":"2024-11-07T01:50:49.336774Z","iopub.status.idle":"2024-11-07T01:50:52.130043Z","shell.execute_reply.started":"2024-11-07T01:50:49.336732Z","shell.execute_reply":"2024-11-07T01:50:52.128300Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Visualizing 6 Particles in Tomogram Slices with difficulty levels","metadata":{}},{"cell_type":"code","source":"# Define particle types with their properties\nPARTICLE_TYPES = {\n    'apo-ferritin': {'color': 'red', 'marker': 'o', 'difficulty': 'easy'},\n    'beta-amylase': {'color': 'gray', 'marker': 's', 'difficulty': 'impossible'},\n    'beta-galactosidase': {'color': 'blue', 'marker': '^', 'difficulty': 'hard'},\n    'ribosome': {'color': 'green', 'marker': 'D', 'difficulty': 'easy'},\n    'thyroglobulin': {'color': 'purple', 'marker': 'p', 'difficulty': 'hard'},\n    'virus-like-particle': {'color': 'orange', 'marker': '*', 'difficulty': 'easy'}\n}\n\ndef load_all_particle_coordinates(experiment_name='TS_5_4'):\n    \"\"\"Load coordinates for all particle types.\"\"\"\n    base_path = Path('/kaggle/input/czii-cryo-et-object-identification/train/overlay/ExperimentRuns')\n    particle_coords = {}\n    \n    for particle_type in PARTICLE_TYPES.keys():\n        json_path = base_path / experiment_name / 'Picks' / f'{particle_type}.json'\n        try:\n            with open(json_path, 'r') as f:\n                data = json.load(f)\n                coords = []\n                for point in data['points']:\n                    coords.append([\n                        point['location']['z'],\n                        point['location']['y'],\n                        point['location']['x']\n                    ])\n                particle_coords[particle_type] = np.array(coords)\n                print(f\"Loaded {len(coords)} {particle_type} coordinates\")\n        except Exception as e:\n            print(f\"Error reading {particle_type} coordinates: {e}\")\n            particle_coords[particle_type] = np.array([])\n    \n    return particle_coords","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:52:37.410530Z","iopub.execute_input":"2024-11-07T01:52:37.410975Z","iopub.status.idle":"2024-11-07T01:52:37.420951Z","shell.execute_reply.started":"2024-11-07T01:52:37.410932Z","shell.execute_reply":"2024-11-07T01:52:37.419785Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Updated particle types with high-contrast colors\nPARTICLE_TYPES = {\n    'apo-ferritin': {'color': '#FF3333', 'marker': 'o', 'difficulty': 'easy'},          # Bright red\n    'beta-amylase': {'color': '#FFFFFF', 'marker': 's', 'difficulty': 'impossible'},     # White\n    'beta-galactosidase': {'color': '#33FFFF', 'marker': '^', 'difficulty': 'hard'},    # Cyan\n    'ribosome': {'color': '#33FF33', 'marker': 'D', 'difficulty': 'easy'},              # Bright green\n    'thyroglobulin': {'color': '#FF33FF', 'marker': 'p', 'difficulty': 'hard'},         # Magenta\n    'virus-like-particle': {'color': '#FFFF33', 'marker': '*', 'difficulty': 'easy'}     # Yellow\n}\n\ndef visualize_all_particles(tomogram, particle_coords, n_slices=3, slice_thickness=20):  # Increased slice thickness\n    \"\"\"\n    Visualize all particle types in tomogram slices with an overlay legend.\n    \"\"\"\n    fig = plt.figure(figsize=(20, 10))\n    grid = ImageGrid(fig, 111,\n                    nrows_ncols=(1, n_slices),\n                    axes_pad=0.3,\n                    share_all=True,\n                    cbar_location=\"right\",\n                    cbar_mode=\"single\",\n                    cbar_size=\"5%\",\n                    cbar_pad=0.1)\n    \n    # Normalize tomogram data\n    vmin, vmax = np.percentile(tomogram, (1, 99))\n    normalized_tomogram = np.clip((tomogram - vmin) / (vmax - vmin), 0, 1)\n    \n    # Find z-positions with maximum particle density\n    all_z_coords = []\n    for coords in particle_coords.values():\n        if len(coords) > 0:\n            all_z_coords.extend(coords[:, 0])\n    \n    if all_z_coords:\n        z_coords = np.array(all_z_coords)\n        z_density = np.histogram(z_coords, bins=50)[0]\n        highest_density_indices = np.argsort(z_density)[-n_slices:]\n        z_positions = np.linspace(z_coords.min(), z_coords.max(), 51)[highest_density_indices]\n    else:\n        z_positions = np.linspace(0, tomogram.shape[0]-1, n_slices, dtype=int)\n    \n    # Plot each slice\n    for idx, ax in enumerate(grid):\n        z = int(z_positions[idx])\n        \n        # Show tomogram slice\n        im = ax.imshow(normalized_tomogram[z, :, :], cmap='gray', vmin=0, vmax=1)\n        \n        # Plot each particle type\n        particles_in_slice = 0\n        particle_counts = {}\n        \n        for particle_type, coords in particle_coords.items():\n            if len(coords) > 0:\n                # Find particles near this slice\n                mask = np.abs(coords[:, 0] - z) < slice_thickness\n                if np.any(mask):\n                    style = PARTICLE_TYPES[particle_type]\n                    ax.scatter(coords[mask, 2], coords[mask, 1],\n                             color=style['color'], marker=style['marker'],\n                             s=100, facecolors='none', linewidth=2,\n                             label=f\"{particle_type}\\n({style['difficulty']})\")\n                    count = np.sum(mask)\n                    particles_in_slice += count\n                    particle_counts[particle_type] = count\n        \n        # Create detailed title showing counts for each particle type\n        title_parts = [f'Slice Z={z}']\n        if particle_counts:\n            for ptype, count in particle_counts.items():\n                if count > 0:\n                    title_parts.append(f'{ptype}: {count}')\n        title = '\\n'.join(title_parts)\n        ax.set_title(title, fontsize=8)\n        \n        ax.grid(False)\n        \n        # Set the axes limits to match the tomogram dimensions\n        ax.set_xlim(0, tomogram.shape[2])\n        ax.set_ylim(tomogram.shape[1], 0)  # Inverted y-axis to match image coordinates\n        \n        # Add legend with semi-transparent background for better visibility\n        if idx == 0:  # Only add legend to first subplot\n            handles, labels = ax.get_legend_handles_labels()\n            legend = ax.legend(handles, labels,\n                             bbox_to_anchor=(0.02, 0.98), \n                             loc='upper left',\n                             borderaxespad=0.,\n                             framealpha=0.8,\n                             facecolor='black',\n                             edgecolor='white',\n                             labelcolor='white',\n                             fontsize=8)\n            \n            for handle in handles:\n                handle.set_edgecolor('black')\n                handle.set_linewidth(1.5)\n    \n    # Add colorbar and title\n    grid.cbar_axes[0].colorbar(im)\n    \n    plt.suptitle('All Particle Types in Tomogram Slices\\n' + \n                 f'Showing particles within ±{slice_thickness} units of each slice',\n                 fontsize=16, y=1.05)\n    \n    plt.show()\n    \n    # Print overall particle statistics\n    print(\"\\nOverall Particle Statistics:\")\n    print(\"-\" * 50)\n    for particle_type, coords in particle_coords.items():\n        if len(coords) > 0:\n            print(f\"\\n{particle_type} ({PARTICLE_TYPES[particle_type]['difficulty']}):\")\n            print(f\"Total particles: {len(coords)}\")\n            print(f\"Z range: {coords[:, 0].min():.1f} to {coords[:, 0].max():.1f}\")\n    \n    # return fig","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:50:56.778208Z","iopub.execute_input":"2024-11-07T01:50:56.778651Z","iopub.status.idle":"2024-11-07T01:50:56.801536Z","shell.execute_reply.started":"2024-11-07T01:50:56.778608Z","shell.execute_reply":"2024-11-07T01:50:56.800179Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load all particle coordinates\nall_particle_coords = load_all_particle_coordinates()\n\n# Scale coordinates for each particle type\nscaled_particle_coords = {\n    particle_type: scale_coordinates(coords, tomogram.shape)\n    for particle_type, coords in all_particle_coords.items()\n}\n\n# Create visualization with all particle types\nvisualize_all_particles(tomogram, scaled_particle_coords)\n\n# Print statistics for each particle type\nprint(\"\\nParticle Statistics:\")\nprint(\"-\" * 50)\nfor particle_type, coords in scaled_particle_coords.items():\n    if len(coords) > 0:\n        print(f\"\\n{particle_type} ({PARTICLE_TYPES[particle_type]['difficulty']}):\")\n        print(f\"Number of particles: {len(coords)}\")\n        print(f\"Z range: {coords[:, 0].min():.1f} to {coords[:, 0].max():.1f}\")\n        print(f\"Y range: {coords[:, 1].min():.1f} to {coords[:, 1].max():.1f}\")\n        print(f\"X range: {coords[:, 2].min():.1f} to {coords[:, 2].max():.1f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T01:52:40.659998Z","iopub.execute_input":"2024-11-07T01:52:40.660461Z","iopub.status.idle":"2024-11-07T01:52:43.598168Z","shell.execute_reply.started":"2024-11-07T01:52:40.660413Z","shell.execute_reply":"2024-11-07T01:52:43.596576Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}