age_conn.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import os.path as op\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import nibabel as nib\n",
    "from nimare.transforms import p_to_z\n",
    "\n",
    "import nilearn\n",
    "from nilearn import plotting\n",
    "from nilearn import datasets\n",
    "import nilearn.reporting\n",
    "from nilearn.image import threshold_img\n",
    "\n",
    "from nilearn.maskers import NiftiSpheresMasker\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = \"./dset\"\n",
    "age_dir = op.join(data_dir, \"age-effect5-21old/habenula\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_theme(color_codes=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get peaks\n",
    "brik_fn = op.join(age_dir, \"sub-group_task-rest_desc-agehabenula_briks+tlrc.BRIK\")\n",
    "table_fn = op.join(age_dir, \"sub-group_task-rest_desc-Agehabenula_table.txt\")\n",
    "nii_1s_fn = op.join(age_dir, \"sub-group_task-rest_desc-ageDiff_zmap.nii.gz\")\n",
    "cluster_fn = op.join(age_dir, \"clustsim\", \"clustsim_out.NN3_bisided.1D\")\n",
    "\n",
    "column_names = [\".10000\", \".05000\", \".02000\", \".01000\"]\n",
    "cluster_df = pd.read_table(\n",
    "    cluster_fn, skiprows=8, delim_whitespace=True, names=column_names\n",
    ")\n",
    "cluster_df = cluster_df.reset_index()\n",
    "cluster_df.rename(columns={\"index\": \"pthr\"}, inplace=True)\n",
    "print(cluster_df)\n",
    "\n",
    "cmap = \"cool\"\n",
    "brik_idx = [11]\n",
    "nii_fns = [nii_1s_fn]\n",
    "tests = [\"1s\"]\n",
    "alpha = \".10000\"\n",
    "pthrs = [0.0050]\n",
    "\n",
    "cohen_thresh = 0\n",
    "\n",
    "data_df = pd.read_csv(table_fn, sep=\"\\t\")\n",
    "n_sub = data_df.groupby(\"group\").size().sum()\n",
    "n_sub_1, n_sub_2 = data_df.groupby(\"group\").size().values\n",
    "\n",
    "for brik_i, nii_fn, test, pthr in zip(brik_idx, nii_fns, tests, pthrs):\n",
    "    convert = f\"3dAFNItoNIFTI \\\n",
    "        -prefix {nii_fn} \\\n",
    "        {brik_fn}[{brik_i}]\"\n",
    "    os.system(convert)\n",
    "\n",
    "    nii_img = nib.load(nii_fn)\n",
    "    z_thresh = p_to_z(pthr)\n",
    "    info = nii_img.get_fdata()\n",
    "\n",
    "    clust_ext = cluster_df.loc[cluster_df[\"pthr\"] == pthr, alpha].values[0]\n",
    "    nii_thr_img = threshold_img(nii_img, z_thresh, cluster_threshold=clust_ext)\n",
    "    print(clust_ext, pthr, z_thresh)\n",
    "\n",
    "    clusters = nilearn.reporting.get_clusters_table(\n",
    "        nii_thr_img, z_thresh, two_sided=True\n",
    "    )\n",
    "    print(clusters)  # coordinates are same as affine of input (MNI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import datasets\n",
    "from nilearn.image import load_img\n",
    "from nilearn.input_data import NiftiLabelsMasker\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Load the Juelich atlas\n",
    "'''juelich_atlas = datasets.fetch_atlas_juelich(\"maxprob-thr0-1mm\")\n",
    "atlas_filename = juelich_atlas.maps\n",
    "atlas_labels = juelich_atlas.labels'''\n",
    "\n",
    "# Load the Harvard-Oxford atlas\n",
    "harvard_oxford_atlas = datasets.fetch_atlas_harvard_oxford(\"sub-maxprob-thr25-1mm\")\n",
    "atlas_filename = harvard_oxford_atlas.maps\n",
    "atlas_labels = harvard_oxford_atlas.labels\n",
    "\n",
    "# Load the atlas image\n",
    "atlas_img = load_img(atlas_filename)\n",
    "\n",
    "\n",
    "# Function to find the region for a given coordinate\n",
    "def find_region(coord):\n",
    "    # Transform MNI coordinates to voxel indices\n",
    "    voxel_indices = np.round(\n",
    "        nib.affines.apply_affine(np.linalg.inv(atlas_img.affine), coord)\n",
    "    ).astype(int)\n",
    "    # Get the label index for the voxel\n",
    "    label_idx = atlas_img.get_fdata()[tuple(voxel_indices)]\n",
    "    region_name = (\n",
    "        atlas_labels[int(label_idx)]\n",
    "        if int(label_idx) < len(atlas_labels)\n",
    "        else \"Unknown\"\n",
    "    )\n",
    "    return region_name\n",
    "\n",
    "\n",
    "# Extract the coordinates from the clusters table and map to regions\n",
    "def map_clusters_to_regions(clusters_df):\n",
    "    coordinates = clusters_df[[\"X\", \"Y\", \"Z\"]].values\n",
    "    regions = [find_region(coord) for coord in coordinates]\n",
    "    clusters_df[\"Region\"] = regions\n",
    "    return clusters_df\n",
    "\n",
    "\n",
    "# Iterate over your clusters and map to regions\n",
    "for brik_i, nii_fn, test, pthr in zip(brik_idx, nii_fns, tests, pthrs):\n",
    "    convert = f\"3dAFNItoNIFTI -prefix {nii_fn} {brik_fn}[{brik_i}]\"\n",
    "    os.system(convert)\n",
    "\n",
    "    nii_img = nib.load(nii_fn)\n",
    "    z_thresh = p_to_z(pthr)\n",
    "    info = nii_img.get_fdata()\n",
    "\n",
    "    clust_ext = cluster_df.loc[cluster_df[\"pthr\"] == pthr, alpha].values[0]\n",
    "    nii_thr_img = threshold_img(nii_img, z_thresh, cluster_threshold=clust_ext)\n",
    "    print(clust_ext, pthr, z_thresh)\n",
    "\n",
    "    clusters = nilearn.reporting.get_clusters_table(nii_thr_img, stat_threshold=z_thresh, two_sided=True)\n",
    "    clusters_with_regions = map_clusters_to_regions(clusters)\n",
    "    print(clusters_with_regions)\n",
    "\n",
    "    output_filename = op.join(age_dir, f\"clusters_with_subregions_{test}_pthr_{pthr}.csv\")\n",
    "    clusters_with_regions.to_csv(output_filename, index=False)\n",
    "    print(f\"Saved clusters with regions to {output_filename}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "import pandas as pd\n",
    "import nibabel as nib\n",
    "import numpy as np\n",
    "from nilearn.image import threshold_img\n",
    "from nilearn.reporting import get_clusters_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Function to execute the AFNI whereami command and parse its output\n",
    "def find_region(coord):\n",
    "    # Format the coordinates properly for the whereami command\n",
    "    coord_str = f\"{coord[0]} {coord[1]} {coord[2]}\"\n",
    "    # Command to query AFNI's whereami\n",
    "    whereami_cmd = f\"whereami [{coord_str}] -atlas MNI_Glasser_HCP_v1.0 -classic\"\n",
    "    try:\n",
    "        result = subprocess.run(\n",
    "            whereami_cmd,\n",
    "            shell=True,\n",
    "            check=True,\n",
    "            stdout=subprocess.PIPE,\n",
    "            stderr=subprocess.PIPE,\n",
    "            text=True,\n",
    "        )\n",
    "        output = result.stdout\n",
    "        # Print the full output for debugging purposes\n",
    "        print(f\"whereami output for {coord}:\")\n",
    "        print(output)\n",
    "        # Parse the output to extract the region name\n",
    "        for line in output.split(\"\\n\"):\n",
    "            if \"Focus point (LPI)\" in line:\n",
    "                parts = line.split()\n",
    "                if len(parts) > 3:\n",
    "                    return parts[3]  # Extract the region name\n",
    "        return \"Unknown\"\n",
    "    except subprocess.CalledProcessError as e:\n",
    "        print(\n",
    "            f\"Error running whereami command: {e}\\nOutput: {e.output}\\nError: {e.stderr}\"\n",
    "        )\n",
    "        return \"Unknown\"\n",
    "\n",
    "\n",
    "# Extract the coordinates from the clusters table and map to regions\n",
    "def map_clusters_to_regions(clusters_df):\n",
    "    coordinates = clusters_df[[\"X\", \"Y\", \"Z\"]].values\n",
    "    regions = [find_region(coord) for coord in coordinates]\n",
    "    clusters_df[\"Region\"] = regions\n",
    "    return clusters_df\n",
    "\n",
    "\n",
    "# Iterate over your clusters and map to regions\n",
    "for brik_i, nii_fn, test, pthr in zip(brik_idx, nii_fns, tests, pthrs):\n",
    "    convert = f\"3dAFNItoNIFTI -prefix {nii_fn} {brik_fn}[{brik_i}]\"\n",
    "    os.system(convert)\n",
    "\n",
    "    nii_img = nib.load(nii_fn)\n",
    "    z_thresh = p_to_z(pthr)\n",
    "    info = nii_img.get_fdata()\n",
    "\n",
    "    clust_ext = cluster_df.loc[cluster_df[\"pthr\"] == pthr, alpha].values[0]\n",
    "    nii_thr_img = threshold_img(nii_img, z_thresh, cluster_threshold=clust_ext)\n",
    "\n",
    "    clusters = get_clusters_table(nii_thr_img, stat_threshold=z_thresh, two_sided=True)\n",
    "    clusters_with_regions = map_clusters_to_regions(clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter the DataFrame to include only rows where 'Cluster Size (mm3)' is empty\n",
    "subclusters = clusters[clusters[\"Cluster Size (mm3)\"] == \"\"]\n",
    "\n",
    "# Extract the main cluster ID from the 'Cluster ID' column for grouping\n",
    "subclusters[\"Main Cluster ID\"] = subclusters[\"Cluster ID\"].str.extract(\n",
    "    r\"(\\d+)\", expand=False\n",
    ")\n",
    "\n",
    "# Group by the main cluster ID and extract coordinates\n",
    "grouped_coordinates = subclusters.groupby(\"Main Cluster ID\").apply(\n",
    "    lambda df: df[[\"X\", \"Y\", \"Z\"]].values.tolist()\n",
    ")\n",
    "\n",
    "# Convert the grouped coordinates to a dictionary for easier access\n",
    "coordinates_dict = grouped_coordinates.to_dict()\n",
    "\n",
    "# Print the coordinates dictionary\n",
    "for cluster_id, coords in coordinates_dict.items():\n",
    "    print(f\"Cluster {cluster_id}: {coords}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a list of all directories in the base directory\n",
    "dset = os.listdir(data_dir)\n",
    "\n",
    "# Filter out only subject directories (assuming they start with 'sub-')\n",
    "sub_dir = [d for d in dset if d.startswith(\"sub-\")]\n",
    "\n",
    "# Loop over each subject directory\n",
    "# on hpc will use partipant.tsv\n",
    "for subject in sub_dir:\n",
    "    subject_path = os.path.join(data_dir, subject, \"func\")\n",
    "\n",
    "    if os.path.isdir(subject_path):\n",
    "        print(f\"Processing subject: {subject}\")\n",
    "        subject_brik = op.join(\n",
    "            subject_path,\n",
    "            f\"{subject}_task-rest_space-MNI152NLin2009cAsym_res-2_desc-norm_bucketREML+tlrc.BRIK\",\n",
    "        )\n",
    "        subject_nii = op.join(\n",
    "            subject_path,\n",
    "            f\"{subject}_task-rest_space-MNI152NLin2009cAsym_res-2_desc-norm_bucketREML+tlrc.nii.gz\",\n",
    "        )\n",
    "\n",
    "        convert = f\"3dAFNItoNIFTI \\\n",
    "        -prefix {subject_nii} \\\n",
    "        {subject_brik}[{1}]\"\n",
    "        os.system(convert)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = []\n",
    "\n",
    "for cluster_id, coords in coordinates_dict.items():\n",
    "    print(f\"Cluster {cluster_id}: {coords}\")\n",
    "    for subject in sub_dir:\n",
    "        subject_path = os.path.join(data_dir, subject, \"func\")\n",
    "\n",
    "        if os.path.isdir(subject_path):\n",
    "            print(f\"Processing subject: {subject}\")\n",
    "\n",
    "            subject_nii = nib.load(op.join(\n",
    "                subject_path,\n",
    "                f\"{subject}_task-rest_space-MNI152NLin2009cAsym_res-2_desc-norm_bucketREML+tlrc.nii.gz\",\n",
    "            ))\n",
    "\n",
    "            mask_img = datasets.load_mni152_brain_mask(resolution=1)\n",
    "\n",
    "            for coord in coords:\n",
    "                print(coord)\n",
    "                masker = NiftiSpheresMasker([coord], radius = 5, mask_img=mask_img)\n",
    "                masker = masker.fit_transform(subject_nii)\n",
    "                print(masker)\n",
    "\n",
    "                results.append(\n",
    "                    {\"Coordinate\": coord, \"Subject\": subject, \"Correlation\": masker}\n",
    "                )\n",
    "\n",
    "corr_df = pd.DataFrame(results)\n",
    "print(corr_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cluster_id, coords in coordinates_dict.items():\n",
    "    print(f\"Cluster {cluster_id}: {coords}\")\n",
    "\n",
    "coords_all = []\n",
    "for coord_list in coordinates_dict.values():\n",
    "    coords_all.extend(coord_list)\n",
    "\n",
    "# Create an empty image\n",
    "brain_mask = datasets.load_mni152_template(resolution=1)\n",
    "spheres_img = np.zeros(brain_mask.shape)\n",
    "\n",
    "# Fill the spheres in the empty image\n",
    "affine = brain_mask.affine\n",
    "for coord in coords_all:\n",
    "    # Convert from MNI coordinates to voxel indices\n",
    "    x, y, z = nib.affines.apply_affine(np.linalg.inv(affine), coord).astype(int)\n",
    "    radius = 5  # radius in voxels\n",
    "\n",
    "    # Create a sphere\n",
    "    for i in range(-radius, radius + 1):\n",
    "        for j in range(-radius, radius + 1):\n",
    "            for k in range(-radius, radius + 1):\n",
    "                if i**2 + j**2 + k**2 <= radius**2:\n",
    "                    if (\n",
    "                        0 <= x + i < spheres_img.shape[0]\n",
    "                        and 0 <= y + j < spheres_img.shape[1]\n",
    "                        and 0 <= z + k < spheres_img.shape[2]\n",
    "                    ):\n",
    "                        spheres_img[x + i, y + j, z + k] = 1\n",
    "\n",
    "# Convert the numpy array to a NIfTI image\n",
    "spheres_nifti = nib.Nifti1Image(spheres_img, affine)\n",
    "\n",
    "# Plot the spheres on a glass brain\n",
    "plotting.plot_glass_brain(\n",
    "    spheres_nifti, display_mode=\"ortho\", colorbar=True, title=f\"{radius}mm Spheres for Age Effect\", cmap='cool'\n",
    ")\n",
    "\n",
    "# Show the plot\n",
    "plotting.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cluster_id, coords in coordinates_dict.items():\n",
    "    print(f\"Cluster {cluster_id}: {coords}\")\n",
    "\n",
    "coords_all = []\n",
    "for coord_list in coordinates_dict.values():\n",
    "    coords_all.extend(coord_list)\n",
    "\n",
    "# Create an empty image\n",
    "brain_mask = datasets.load_mni152_template(resolution=1)\n",
    "spheres_img = np.zeros(brain_mask.shape)\n",
    "\n",
    "# Fill the spheres in the empty image\n",
    "affine = brain_mask.affine\n",
    "for coord in coords_all:\n",
    "    # Convert from MNI coordinates to voxel indices\n",
    "    x, y, z = nib.affines.apply_affine(np.linalg.inv(affine), coord).astype(int)\n",
    "    radius = 5  # radius in voxels\n",
    "\n",
    "    # Create a sphere\n",
    "    for i in range(-radius, radius + 1):\n",
    "        for j in range(-radius, radius + 1):\n",
    "            for k in range(-radius, radius + 1):\n",
    "                if i**2 + j**2 + k**2 <= radius**2:\n",
    "                    if (\n",
    "                        0 <= x + i < spheres_img.shape[0]\n",
    "                        and 0 <= y + j < spheres_img.shape[1]\n",
    "                        and 0 <= z + k < spheres_img.shape[2]\n",
    "                    ):\n",
    "                        spheres_img[x + i, y + j, z + k] = 1\n",
    "\n",
    "# Convert the numpy array to a NIfTI image\n",
    "spheres_nifti = nib.Nifti1Image(spheres_img, affine)\n",
    "\n",
    "# Plot the spheres on a glass brain\n",
    "plotting.plot_glass_brain(\n",
    "    spheres_nifti, display_mode=\"ortho\", colorbar=True, title=f\"{radius}mm Spheres for Age Effect\", cmap='cool'\n",
    ")\n",
    "\n",
    "# Show the plot\n",
    "plotting.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coordinates_dict = {\n",
    "    \"1\": [[25.5, -60.5, -6.5], [23.5, -68.5, -2.5], [27.5, -56.5, 7.5]],\n",
    "    \"2\": [[-30.5, -54.5, -16.5], [-32.5, -56.5, -6.5], [-20.5, -62.5, -4.5]],\n",
    "}\n",
    "\n",
    "corr_df = pd.read_csv(op.join(age_dir, \"age-correlation.csv\"))\n",
    "corr_df[\"Coordinate\"] = corr_df[\"Coordinate\"].apply(eval)\n",
    "\n",
    "for cluster_id, coords in coordinates_dict.items():\n",
    "    for coord in coords:\n",
    "        coord_df = corr_df[corr_df[\"Coordinate\"].apply(lambda x: x == coord)]\n",
    "        print(f\"Condensed DataFrame for Coordinate {coord}:\")\n",
    "        print(coord_df)\n",
    "\n",
    "        g = sns.lmplot(\n",
    "            data=coord_df,\n",
    "            x=\"Age\",\n",
    "            y=\"Correlation\",\n",
    "            hue=\"Group\",\n",
    "            markers = [\"+\", \"d\"],\n",
    "            height=7,\n",
    "        )\n",
    "\n",
    "        plt.ylim(-0.6, 0.6)\n",
    "        \n",
    "        # Add title to the plot\n",
    "        plt.title(f\"Age Effect: {radius}mm Sphere for MNI {coord}\")\n",
    "\n",
    "        # Show the plot\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "def r2(x, y):\n",
    "    return stats.pearsonr(x, y)[0] ** 2\n",
    "\n",
    "\n",
    "coordinates_dict = {\n",
    "    \"1\": [[25.5, -60.5, -6.5], [23.5, -68.5, -2.5], [27.5, -56.5, 7.5]],\n",
    "    \"2\": [[-30.5, -54.5, -16.5], [-32.5, -56.5, -6.5], [-20.5, -62.5, -4.5]],\n",
    "}\n",
    "\n",
    "# Assuming corr_df is your DataFrame loaded from \"age-correlation.csv\"\n",
    "corr_df = pd.read_csv(op.join(age_dir, \"age-correlation.csv\"))\n",
    "\n",
    "# Group by 'Subject' and 'Cluster', and calculate the mean of 'Correlation'\n",
    "corr_df_mean = corr_df.groupby([\"Subject\", \"Cluster\"], as_index=False)[\n",
    "    \"Correlation\"\n",
    "].mean()\n",
    "corr_df_mean.rename(columns={\"Correlation\": \"Correlation (Fisher's Z)\"}, inplace=True)\n",
    "\n",
    "# Merge averaged correlations back into the original DataFrame and replace 'Correlation' column\n",
    "corr_df = pd.merge(corr_df, corr_df_mean, on=[\"Subject\", \"Cluster\"], how=\"left\")\n",
    "corr_df[\"Correlation\"] = corr_df[\n",
    "    \"Correlation (Fisher's Z)\"\n",
    "]  # Replace original 'Correlation' with 'Avg_Correlation'\n",
    "corr_df.drop(\n",
    "    columns=[\"Correlation\", \"Coordinate\"], inplace=True\n",
    ")  # Drop the temporary column\n",
    "\n",
    "corr_df.drop_duplicates(subset=[\"Subject\", \"Cluster\"], keep=\"first\", inplace=True)\n",
    "\n",
    "corr_df[\"Group\"] = corr_df[\"Group\"].replace({\"asd\": \"ASD\", \"td\": \"NT\"})\n",
    "\n",
    "clusters = [1, 2]\n",
    "custom_palette = {\"ASD\": \"#A1AF58\", \"NT\": \"#73AFF1\"}  # Custom palette with colors for each group\n",
    "custom_palette2 = {\"ASD\": \"#A1AF58\", \"NT\": \"#73AFF1\"}  # Map Group names to colors\n",
    "\n",
    "regions = [\"Right Parahippocampal Gyrus\", \"Left Parahippocampal Gyrus\"]\n",
    "\n",
    "# Print the number of subjects and average age for each group\n",
    "group_counts = corr_df[\"Group\"].value_counts()\n",
    "group_means = corr_df.groupby(\"Group\")[\"Age\"].mean()\n",
    "overall_mean_age = corr_df[\"Age\"].mean()\n",
    "\n",
    "print(f\"Number of ASD subjects: {group_counts['ASD']}\")\n",
    "print(f\"Number of NT subjects: {group_counts['NT']}\")\n",
    "print(f\"Average age of ASD subjects: {group_means['ASD']:.2f}\")\n",
    "print(f\"Average age of NT subjects: {group_means['NT']:.2f}\")\n",
    "\n",
    "\n",
    "for cluster, region in zip(clusters, regions):\n",
    "    cluster_df = corr_df[corr_df[\"Cluster\"] == cluster].copy()\n",
    "\n",
    "    # Print the number of subjects and average age for each group\n",
    "    group_counts = cluster_df[\"Group\"].value_counts()\n",
    "    group_means = cluster_df.groupby(\"Group\")[\"Age\"].mean()\n",
    "    overall_mean_age = cluster_df[\"Age\"].mean()\n",
    "    overall_std_age = cluster_df[\"Age\"].std()\n",
    "\n",
    "    print(f\"Number of ASD subjects: {group_counts['ASD']}\")\n",
    "    print(f\"Number of NT subjects: {group_counts['NT']}\")\n",
    "    print(f\"Average age of ASD subjects: {group_means['ASD']:.2f}\")\n",
    "    print(f\"Average age of NT subjects: {group_means['NT']:.2f}\")\n",
    "    print(\n",
    "        f\"Average age of all subjects: {overall_mean_age:.2f} (SD: {overall_std_age:.2f})\"\n",
    "    )\n",
    "\n",
    "    # Create a joint plot\n",
    "    g = sns.jointplot(\n",
    "        data=cluster_df,\n",
    "        x=\"Age\",\n",
    "        y=\"Correlation (Fisher's Z)\",\n",
    "        hue=\"Group\",\n",
    "        palette=custom_palette,\n",
    "        kind=\"scatter\",\n",
    "        height=7,\n",
    "    )\n",
    "    # Move the legend to the bottom left\n",
    "    g.ax_joint.legend(title=\"Group\", loc=\"lower left\", handletextpad=0.5, fontsize=10)\n",
    "\n",
    "    g.ax_joint.collections[0].set_sizes([15])\n",
    "\n",
    "    # Adding regression lines for each group\n",
    "    for group_name, group_data in cluster_df.groupby(\"Group\"):\n",
    "        sns.regplot(\n",
    "            x=\"Age\",\n",
    "            y=\"Correlation (Fisher's Z)\",\n",
    "            data=group_data,\n",
    "            scatter=False,\n",
    "            ax=g.ax_joint,\n",
    "            color=custom_palette2[group_name],  # Match color with the group\n",
    "            label=f\"Group {group_name}\",\n",
    "        )\n",
    "\n",
    "    # Set limits for y-axis\n",
    "    g.ax_joint.set_ylim(-0.2, 0.2)\n",
    "    g.ax_joint.set_xlim(5, 21)\n",
    "\n",
    "    # Add a title to the plot\n",
    "    g.ax_joint.set_xlabel(f\"Age in Years\")\n",
    "    g.fig.suptitle(f\"{region}\", y=1.02)\n",
    "\n",
    "    # Display the plot\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}