{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Getting started with pyAFQ - GroupAFQ\n\nThere are two ways to :doc:`use pyAFQ </tutorials/index>`: through the\ncommand line interface, and by writing Python code. This tutorial will walk you\nthrough the basics of the latter, using pyAFQ's Python Application Programming\nInterface (API).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os.path as op\n\nimport matplotlib.pyplot as plt\nimport nibabel as nib\nimport plotly\nimport pandas as pd\n\nfrom AFQ.api.group import GroupAFQ\nimport AFQ.data.fetch as afd\nimport AFQ.viz.altair as ava\nimport AFQ.definitions.image as afm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Example data\npyAFQ can be called using GroupAFQ to handle data\norganized in a BIDS compliant directory.\nIf this is not the case, refer to the Participant AFQ example.\nTo get users started with this tutorial, we will download some example\ndata and organize it in a BIDS compliant way (for more details on how\nBIDS is used in pyAFQ, refer to :doc:`plot_006_bids_layout`).\n\nThe following call downloads a a single subject's data from the Healthy Brain\nNetwork Processed Open Diffusion Derivatives dataset (HBN-POD2) [1]_, [2]_\nand organizes it in BIDS in the user's home directory under::\n\n  ``~/AFQ_data/HBN/``\n\nThe data is also placed in a derivatives directory, signifying that it has\nalready undergone the required preprocessing necessary for pyAFQ to run.\n\nThe clear_previous_afq is used to remove any previous runs of the afq object\nstored in the `~/AFQ_data/HBN/` BIDS directory. Set it to None if\nyou want to use the results of previous runs.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bids_path = afd.fetch_hbn_preproc(\n    [\"NDARAA948VFH\"],\n    clear_previous_afq=\"all\")[1]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Set tractography parameters (optional)\nWe make create a `tracking_params` variable, which we will pass to the\nGroupAFQ object which specifies that we want 100,000 seeds randomly\ndistributed in the white matter. We only do this to make this example faster\nand consume less space; normally, we use more seeds.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tracking_params = dict(n_seeds=1e5,\n                       random_seeds=True,\n                       rng_seed=2025,\n                       trx=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Define PVE images (optional)\nTo improve segmentation and tractography results, we can provide\npartial volume estimate (PVE) images for the cerebrospinal fluid (CSF),\ngray matter (GM), and white matter (WM). Here, we define these images\nusing the AFQ.definitions.image.PVEImages class, which takes as input\nthree AFQ.definitions.image.ImageFile objects, one for each tissue type.\nOne can also provide a single PVE image with all three tissue types\nusing the AFQ.definitions.image.PVEImage class. Finally, by default,\nif no PVE images are provided, pyAFQ will use SynthSeg2 to compute\nthese images.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pve = afm.PVEImages(\n    afm.ImageFile(\n        suffix=\"probseg\", filters={\"scope\": \"qsiprep\", \"label\": \"CSF\"}),\n    afm.ImageFile(\n        suffix=\"probseg\", filters={\"scope\": \"qsiprep\", \"label\": \"GM\"}),\n    afm.ImageFile(\n        suffix=\"probseg\", filters={\"scope\": \"qsiprep\", \"label\": \"WM\"}))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Brain Mask Definition (optional)\n\nBy default, pyAFQ will compute a brain mask from the T1. However,\nthis requires onnxruntime to be installed. If you do not have onnxruntime\ninstalled, or if you want to use a different brain mask, you can specify\nit here.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "brain_mask_definition = afm.ImageFile(\n    suffix=\"mask\", filters={\"desc\": \"brain\", \"scope\": \"qsiprep\"})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initialize a GroupAFQ object:\n\nCreates a GroupAFQ object, that encapsulates tractometry. This object can be\nused to manage the entire :doc:`/explanations/index`, including:\n\n- Tractography\n- Registration\n- Segmentation\n- Cleaning\n- Profiling\n- Visualization\n\nThis will also create an output folder for the corresponding AFQ derivatives\nin the AFQ data directory: ``AFQ_data/HBN/derivatives/afq/``\n\nTo initialize this object we will pass in the path location to our BIDS\ncompliant data, the name of the preprocessing pipeline we want to use, \nthe name of the t1 preprocessing pipeline we want to use (in this case,\nits the same, qsiprep [3]), the participant labels we want to process\n(in this case, just a single subject), the PVE images we defined above, and\nthe tracking parameters we defined above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "myafq = GroupAFQ(\n    bids_path=op.join(afd.afq_home, 'HBN'),\n    dwi_preproc_pipeline='qsiprep',\n    t1_preproc_pipeline='qsiprep',\n    participant_labels=['NDARAA948VFH'],\n    pve=pve,\n    brain_mask_definition=brain_mask_definition,\n    tracking_params=tracking_params)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Calculating DKI FA (Diffusion Kurtosis Imaging Fractional Anisotropy)\nThe GroupAFQ object has a method called `export`, which allows the user\nto calculate various derived quantities from the data.\n\nFor example, FA can be computed using the DKI model, by explicitly\ncalling `myafq.export(\"dki_fa\")`. This triggers the computation of DKI\nparameters for all subjects in the dataset, and stores the results in\nthe AFQ derivatives directory. In addition, it calculates the FA\nfrom these parameters and stores it in a different file in the same\ndirectory.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The AFQ API computes quantities lazily. This means that DKI parameters\n   are not computed until they are required. This means that the first\n   line below is the one that requires time.</p></div>\n\nThe result of the call to `export` is a dictionary, with the subject\nIDs as keys, and the filenames of the corresponding files as values.\nThis means that to extract the filename corresponding to the FA of the first\nsubject, we can do:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "FA_fname = myafq.export(\"dki_fa\", collapse=False)[\"NDARAA948VFH\"][\"HBNsiteRU\"]\n\n# We will then use `nibabel` to load the deriviative file and retrieve the\n# data array.\n\nFA_img = nib.load(FA_fname)\nFA = FA_img.get_fdata()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualize the result with Matplotlib\nAt this point `FA` is an array, and we can use standard Python tools to\nvisualize it or perform additional computations with it.\n\nIn this case we are going to take an axial slice halfway through the\nFA data array and plot using a sequential color map.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The data array is structured as a xyz coordinate system.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1)\nax.matshow(FA[:, :, FA.shape[-1] // 2], cmap='viridis')\nax.axis(\"off\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Recognizing the bundles and calculating tract profiles:\nTypically, users of pyAFQ are interested in calculating not only an overall\nmap of the FA, but also the major white matter pathways (or bundles) and\ntract profiles of tissue properties along their length. To trigger the\npyAFQ pipeline that calculates the profiles, users can call the\n`export('profiles')` method:\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>Running the code below triggers the full pipeline of operations\n   leading to the computation of the tract profiles. Therefore, it\n   takes a little while to run (about 40 minutes, typically).</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "myafq.export('profiles')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualizing the bundles and calculating tract profiles:\nThe pyAFQ API provides several ways to visualize bundles and profiles.\n\nFirst, we will run a function that exports an html file that contains\nan interactive visualization of the bundles that are segmented.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>By default we resample a 100 points within a bundle, however to reduce\n   processing time we will only resample 50 points.</p></div>\n\nOnce it is done running, it should pop a browser window open and let you\ninteract with the bundles.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>You can hide or show a bundle by clicking the legend, or select a\n   single bundle by double clicking the legend. The interactive\n   visualization will also all you to pan, zoom, and rotate.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bundle_html = myafq.export(\"all_bundles_figure\", collapse=False)\nplotly.io.show(bundle_html[\"NDARAA948VFH\"][\"HBNsiteRU\"][0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also visualize the tract profiles in all of the bundles. These\nplots show both FA (left) and MD (right) laid out anatomically.\nTo make this plot, it is required that you install with\n`pip install pyAFQ[plot]` so that you have the necessary dependencies.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig_files = myafq.export(\"tract_profile_plots\", collapse=False)[\n    \"NDARAA948VFH\"][\"HBNsiteRU\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: {{ fig_files[0] }}\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can even use altair to visualize the tract profiles in all\nof the bundles. We provide a more customizable interface for visualizing\nthe tract profiles using altair.\nAgain, to make this plot, it is required that you install with\n`pip install pyAFQ[plot]` so that you have the necessary dependencies.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "profiles_df = myafq.combine_profiles()\naltair_df = ava.combined_profiles_df_to_altair_df(\n    profiles_df,\n    tissue_properties=['dki_fa', 'dki_md'])\naltair_chart = ava.altair_df_to_chart(altair_df)\naltair_chart.display()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Exporting citations\nFinally, we can export the citations for the some of methods used in this\nanalysis. These are not guaranteed to be comprehensive, but they\nshould be a good starting point.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "myafq.export(\"citations\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can check the number of streamlines per bundle, to make sure\nevery bundle is found with a reasonable amount of streamlines.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bundle_counts = pd.read_csv(\n    myafq.export(\"sl_counts\", collapse=False)[\n        \"NDARAA948VFH\"][\"HBNsiteRU\"], index_col=[0])\nfor ind in bundle_counts.index:\n    if ind == \"Total Recognized\":\n        threshold = 2500\n    else:\n        threshold = 20\n    if bundle_counts[\"n_streamlines\"][ind] < threshold:\n        raise ValueError((\n            \"Small number of streamlines found \"\n            f\"for bundle(s):\\n{bundle_counts}\"))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Alternative way to initialize GroupAFQ when using QSIPrep data\nAs a final note, if you are using QSIPrep preprocessed data,\nyou can also initialize the GroupAFQ object using the\n`from_qsiprep` class method. This method will automatically set\nthe appropriate BIDS filters to find the preprocessed DWI data.\nAdditionally, it will find and use the brain masks and PVE images\nthat QSIPrep generates. Outside of BIDS filters, the arguments\nare the same as those used when initializing the GroupAFQ object\ndirectly.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "myafq = GroupAFQ.from_qsiprep(\n    qsi_dir=op.join(afd.afq_home, 'HBN'),\n    participant_labels=['NDARAA948VFH'],\n    tracking_params=tracking_params)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n.. [1] Alexander LM, Escalera J, Ai L, et al. An open resource for\n    transdiagnostic research in pediatric mental health and learning\n    disorders. Sci Data. 2017;4:170181.\n\n.. [2] Richie-Halford A, Cieslak M, Ai L, et al. An analysis-ready and quality\n    controlled resource for pediatric brain white-matter research. Scientific\n    Data. 2022;9(1):1-27.\n\n.. [3] Cieslak M, Cook PA, He X, et al. QSIPrep: an integrative platform for\n    preprocessing and reconstructing diffusion MRI data. Nat Methods.\n    2021;18(7):775-778.\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.13.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}