pheno_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 nilearn import datasets\n",
    "from nimare.transforms import p_to_z\n",
    "\n",
    "import nilearn\n",
    "from nilearn import plotting\n",
    "import nilearn.reporting\n",
    "from nilearn.image import threshold_img\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = \"./dset\"\n",
    "abide_dir = op.join(data_dir, \"group/habenula\")\n",
    "clust_dir = op.join(abide_dir, \"clustsim\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_theme(color_codes=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Use Nilearn to identify peak clusters from z thresholded group difference correlation maps that will be used for masking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get peaks\n",
    "brik_fn = op.join(\n",
    "    abide_dir, \"sub-group_task-rest_desc-1S2StTesthabenula_briks+tlrc.BRIK\"\n",
    ")\n",
    "table_fn = op.join(abide_dir, \"sub-group_task-rest_desc-1S2StTesthabenula_table.txt\")\n",
    "nii_1s_fn = op.join(abide_dir, \"sub-group_task-rest_desc-2SampletTest_zmap.nii.gz\")\n",
    "cluster_fn = op.join(clust_dir, \"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 = \".02000\"\n",
    "pthrs = [0.0001]\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(nii_thr_img, z_thresh, two_sided=True)\n",
    "    print(clusters)  # coordinates are same as affine of input (MNI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Atlases\n",
    "#### load Harvard and Juliech Atlases "
   ]
  },
  {
   "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": [
    "# 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(\"cort-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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\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(\n",
    "        nii_thr_img, stat_threshold=z_thresh, two_sided=True\n",
    "    )\n",
    "    clusters_with_regions = map_clusters_to_regions(clusters)\n",
    "    print(clusters_with_regions)\n",
    "\n",
    "    output_filename = op.join(\n",
    "        abide_dir, f\"clusters_with_regions_{test}_pthr_{pthr}.csv\"\n",
    "    )\n",
    "    clusters_with_regions.to_csv(output_filename, index=False)\n",
    "    print(f\"Saved clusters with regions to {output_filename}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create dataframes for Statistical Significance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load phenotype correlation data\n",
    "corr_df = pd.read_csv(op.join(abide_dir, \"pheno4clstr-correlation2.csv\"))\n",
    "corr_df = corr_df[corr_df[\"score\"] != -9999.0]  # Remove rows with score == -9999.0\n",
    "\n",
    "\n",
    "# Calculate the mean Correlation for each Subject, cluster, and phenotype\n",
    "mean_corr_df = corr_df.groupby(\n",
    "    [\"Subject\", \"Group\", \"Age\", \"Cluster\", \"phenotype\", \"score\"], as_index=False\n",
    ")[\"Correlation\"].mean()\n",
    "\n",
    "# Pivot the DataFrame to have phenotypes as columns with the mean score values\n",
    "pivot_df = mean_corr_df.pivot(\n",
    "    index=[\"Subject\", \"Group\", \"Age\", \"Cluster\", \"Correlation\"],\n",
    "    columns=\"phenotype\",\n",
    "    values=\"score\",\n",
    ")\n",
    "\n",
    "# Flatten the columns and reset the index to get a clean DataFrame\n",
    "pivot_df.columns.name = None\n",
    "pivot_df.reset_index(inplace=True)\n",
    "\n",
    "pivot_df.rename(\n",
    "    columns={\n",
    "        \"Correlation\": \"RSFC\",\n",
    "        \"ADOS_GOTHAM_SOCAFFECT\": \"Phen1\",\n",
    "        \"ADOS_GOTHAM_RRB\": \"Phen2\",\n",
    "        \"SRS_MOTIVATION\": \"Phen3\",\n",
    "        \"VINELAND_DAILYLVNG_STANDARD\": \"Phen4\",\n",
    "        \"VINELAND_COPING_V_SCALED\": \"Phen5\",\n",
    "        \"VINELAND_RECEPTIVE_V_SCALED\": \"Phen6\",\n",
    "        \"VINELAND_EXPRESSIVE_V_SCALED\": \"Phen7\",\n",
    "        \"VINELAND_WRITTEN_V_SCALED\": \"Phen8\",\n",
    "        \"VINELAND_COMMUNICATION_STANDARD\": \"Phen9\",\n",
    "        \"SRS_COMMUNICATION\": \"Phen10\",\n",
    "    },\n",
    "    inplace=True,\n",
    ")\n",
    "# Load participant.tsv\n",
    "participants_df = pd.read_csv(\n",
    "    op.join(data_dir, \"participants.tsv\"),\n",
    "    sep=\"\\t\",\n",
    "    usecols=[\"participant_id\", \"SITE_ID\", \"SEX\"],\n",
    ")\n",
    "# Rename the columns\n",
    "participants_df.rename(\n",
    "    columns={\"participant_id\": \"Subject\", \"SITE_ID\": \"Site\", \"SEX\": \"Sex\"}, inplace=True\n",
    ")\n",
    "\n",
    "# Merge the two DataFrames on Subject\n",
    "merged_df = pd.merge(pivot_df, participants_df, on=\"Subject\")\n",
    "\n",
    "# Reorder the columns\n",
    "ordered_columns = [\n",
    "    \"Subject\",\n",
    "    \"Cluster\",\n",
    "    \"Group\",\n",
    "    \"Age\",\n",
    "    \"Sex\",\n",
    "    \"Site\",\n",
    "    \"RSFC\",\n",
    "    \"Phen1\",\n",
    "    \"Phen2\",\n",
    "    \"Phen3\",\n",
    "    \"Phen4\",\n",
    "    \"Phen5\",\n",
    "    \"Phen6\",\n",
    "    \"Phen7\",\n",
    "    \"Phen8\",\n",
    "    \"Phen9\",\n",
    "    \"Phen10\",\n",
    "]\n",
    "\n",
    "pheno_df = merged_df[ordered_columns]\n",
    "\n",
    "# Display the reordered DataFrame\n",
    "print(pheno_df)\n",
    "clusters = [1, 2, 3, 4]\n",
    "\n",
    "for cluster in clusters:\n",
    "    cluster_df = pheno_df[\n",
    "        pheno_df[\"Cluster\"] == cluster\n",
    "    ].copy()  # Filter for current cluster\n",
    "\n",
    "    # Drop the cluster column\n",
    "    cluster_df.drop(columns=[\"Cluster\"], inplace=True)\n",
    "\n",
    "    # Drop rows with NaN values\n",
    "    # cluster_df.dropna(axis=0, inplace=True)\n",
    "\n",
    "    # Save the filtered DataFrame to a CSV file\n",
    "    output_path = op.join(abide_dir, f\"cluster-{cluster}_data.csv\")\n",
    "    cluster_df.to_csv(output_path, index=False)\n",
    "\n",
    "    # Print the filtered DataFrame for verification\n",
    "    print(cluster_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Correlations as a function of Phenotypic Score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the r2 function to calculate R^2\n",
    "def r2(x, y):\n",
    "    return stats.pearsonr(x, y)[0] ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename columns of pheno df to what you want to be in your plots\n",
    "\n",
    "pheno_df.rename(\n",
    "    columns={\n",
    "        \"RSFC\": \"Correlation (Fisher's Z)\",\n",
    "        \"Phen1\": \"ADOS_GOTHAM_SOCAFFECT\",\n",
    "        \"Phen2\": \"Restricted and Repetitive Behaviors Total Subscore for Gotham Algorithm of the ADOS\",\n",
    "        \"Phen3\": \"Social Responsivenes Scale Social Motivation Subscore Raw Total\",\n",
    "        \"Phen4\": \"Vineland Adaptive Behavior Scales Daily Living Skills Standard Score\",\n",
    "        \"Phen5\": \"VINELAND_COPING_V_SCALED\",\n",
    "        \"Phen6\": \"Vineland Adaptive Behavior Scales Receptive Language V Scaled Score\",\n",
    "        \"Phen7\": \"VINELAND_EXPRESSIVE_V_SCALED\",\n",
    "        \"Phen8\": \"VINELAND_WRITTEN_V_SCALED\",\n",
    "        \"Phen9\": \"VINELAND_COMMUNICATION_STANDARD\",\n",
    "        \"Phen10\": \"SRS_COMMUNICATION\",\n",
    "    },\n",
    "    inplace=True,\n",
    ")\n",
    "\n",
    "print(pheno_df)\n",
    "\n",
    "# Melt the DataFrame to long format\n",
    "plotcorr_df = pd.melt(\n",
    "    pheno_df,\n",
    "    id_vars=[\n",
    "        \"Subject\",\n",
    "        \"Cluster\",\n",
    "        \"Group\",\n",
    "        \"Age\",\n",
    "        \"Sex\",\n",
    "        \"Site\",\n",
    "        \"Correlation (Fisher's Z)\",\n",
    "    ],  # Assuming 'Subject' is the identifier column\n",
    "    value_vars=[\n",
    "        \"ADOS_GOTHAM_SOCAFFECT\",\n",
    "        \"Restricted and Repetitive Behaviors Total Subscore for Gotham Algorithm of the ADOS\",\n",
    "        \"Social Responsivenes Scale Social Motivation Subscore Raw Total\",\n",
    "        \"Vineland Adaptive Behavior Scales Daily Living Skills Standard Score\",\n",
    "        \"VINELAND_COPING_V_SCALED\",\n",
    "        \"Vineland Adaptive Behavior Scales Receptive Language V Scaled Score\",\n",
    "        \"VINELAND_EXPRESSIVE_V_SCALED\",\n",
    "        \"VINELAND_WRITTEN_V_SCALED\",\n",
    "        \"VINELAND_COMMUNICATION_STANDARD\",\n",
    "        \"SRS_COMMUNICATION\",\n",
    "    ],\n",
    "    var_name=\"Phenotype\",\n",
    "    value_name=\"Score\",\n",
    ")\n",
    "\n",
    "print(plotcorr_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_cluster_correlation(\n",
    "    phenotype_df, phenotype, cluster, region, custom_palette, custom_palette2\n",
    "):\n",
    "    # Calculate x-axis limits based on score column\n",
    "    xlim_min = phenotype_df[\"Score\"].min() - 1  # Minus 1 from the lowest score\n",
    "    xlim_max = phenotype_df[\"Score\"].max() + 1  # Plus 1 to the highest score\n",
    "\n",
    "    # Create a joint plot\n",
    "    g = sns.jointplot(\n",
    "        data=phenotype_df,\n",
    "        x=\"Score\",\n",
    "        y=\"Correlation (Fisher's Z)\",\n",
    "        hue=\"Group\",\n",
    "        palette=custom_palette,\n",
    "        kind=\"scatter\",\n",
    "        height=7,\n",
    "    )\n",
    "\n",
    "    # Adding regression lines for each group\n",
    "    for i, (group_name, group_data) in enumerate(phenotype_df.groupby(\"Group\")):\n",
    "        sns.regplot(\n",
    "            x=\"Score\",\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",
    "        )\n",
    "\n",
    "        # Calculate R^2 for the current group's regression line\n",
    "        r_squared = r2(group_data[\"Score\"], group_data[\"Correlation (Fisher's Z)\"])\n",
    "\n",
    "        # Annotate the plot with the R^2 value, with box around it\n",
    "        # Position annotations vertically\n",
    "        y_offset = i * 0.06  # Adjust vertical spacing between annotations\n",
    "        g.ax_joint.annotate(\n",
    "            f\"$R^2$ = {r_squared:.4f}\",\n",
    "            xy=(0.1, 0.9 - y_offset),  # Adjust the position of the annotation (x, y)\n",
    "            xycoords=\"axes fraction\",\n",
    "            ha=\"left\",\n",
    "            va=\"top\",  # Horizontal and vertical alignment\n",
    "            fontsize=10,\n",
    "            color=custom_palette2[group_name],  # Match color with the group\n",
    "            bbox=dict(\n",
    "                boxstyle=\"round,pad=0.3\",\n",
    "                edgecolor=custom_palette2[group_name],\n",
    "                facecolor=\"white\",\n",
    "                alpha=0.9,\n",
    "            ),  # Box style and colors\n",
    "        )\n",
    "\n",
    "    # Set limits for y-axis\n",
    "    g.ax_joint.set_ylim(-0.3, 0.3)\n",
    "    g.ax_joint.set_xlim(xlim_min, xlim_max)\n",
    "\n",
    "    # Add a title to the plot\n",
    "    g.fig.suptitle(f\"{phenotype}: {region}\", y=1.02)\n",
    "\n",
    "    # Display the plot\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = [1, 2, 3 ,4]  # Example clusters\n",
    "regions = [\n",
    "    \"Right Middle Temporal Gyrus\",\n",
    "    \"Left Superior Temporal Gyrus\",\n",
    "    \"Left Precentral Gyrus\",\n",
    "    \"Left Angular Gyrus\",\n",
    "]  # Example regions\n",
    "custom_palette = [\"#A1AF58\", \"#73AFF1\"]  # Custom palette with colors for each group\n",
    "custom_palette2 = {\"asd\": \"#A1AF58\", \"td\": \"#73AFF1\"}  # Map Group names to colors\n",
    "\n",
    "for phenotype in plotcorr_df[\"Phenotype\"].unique():\n",
    "    # Filter plotcorr_df for the current phenotype\n",
    "    phenotype_df = plotcorr_df[plotcorr_df[\"Phenotype\"] == phenotype].copy()\n",
    "    # print(phenotype_df)\n",
    "\n",
    "    for cluster, region in zip(clusters, regions):\n",
    "        cluster_df = phenotype_df[\n",
    "            phenotype_df[\"Cluster\"] == cluster\n",
    "        ].copy()  # Filter for current cluster\n",
    "\n",
    "        #  Drop rows with NaN in 'Score' column\n",
    "        cluster_df = cluster_df.dropna(subset=[\"Score\"])\n",
    "        print(cluster_df)\n",
    "\n",
    "        plot_cluster_correlation(\n",
    "            cluster_df, phenotype, cluster, region, custom_palette, custom_palette2\n",
    "        )"
   ]
  }
 ],
 "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
}