fw-heudiconv.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using ``heudiconv``-style heuristics to get accurate BIDS \n",
    "\n",
    "BIDS curation is not a straightforward process and can be very tricky when non-standard names are used at the scanner.\n",
    "\n",
    "The python package ``heudiconv`` approaches this problem by having users write a \"heuristic\" python script that assigns scans to BIDS names based on properties in their dicom headers. When a user runs ``heudiconv`` they supply the path to dicoms, the subject ID, session ID and the path to a python script that implements their heuristic. The software then groups the dicoms based on studyUID and other shared properties in their headers, uses the user's script to determine what their names should be in BIDS, then runs dcm2niix and renames the files.\n",
    "\n",
    "The first and last parts of this process are already fully automated in flywheel. Therefore, one possible solution for getting accurate BIDS curation would be to build a system that can query a session's acquisitions and adjust their metadata according to the user-provided heuristic.\n",
    "\n",
    "## Heuristic files\n",
    "\n",
    "At the end of the day, the heuristic file loops over acquisitions and sends their metadata through some user-provided logic. The scans themselves are provided as part of a ``SeqInfo`` object. The scans in the SeqInfo expose attributes to the user that can be used to classify those scans. Here is the complete list:\n",
    "\n",
    "```python\n",
    "seqinfo_fields = [\n",
    "    'total_files_till_now',  # 0\n",
    "    'example_dcm_file',      # 1\n",
    "    'series_id',             # 2\n",
    "    'dcm_dir_name',          # 3\n",
    "    'unspecified2',          # 4\n",
    "    'unspecified3',          # 5\n",
    "    'dim1', 'dim2', 'dim3', 'dim4', # 6, 7, 8, 9\n",
    "    'TR', 'TE',              # 10, 11\n",
    "    'protocol_name',         # 12\n",
    "    'is_motion_corrected',   # 13\n",
    "    'is_derived',            # 14\n",
    "    'patient_id',            # 15\n",
    "    'study_description',     # 16\n",
    "    'referring_physician_name', # 17\n",
    "    'series_description',    # 18\n",
    "    'sequence_name',         # 19\n",
    "    'image_type',            # 20\n",
    "    'accession_number',      # 21\n",
    "    'patient_age',           # 22\n",
    "    'patient_sex',           # 23\n",
    "    'date',                  # 24\n",
    "    'series_uid'             # 25\n",
    "]\n",
    "```\n",
    "\n",
    "### The ``infotodict`` function and its output\n",
    "\n",
    "This function has to be implemented in the user's heuristic script. It takes a ``SeqInfo`` object as its only argument and returns a dictionary that maps a formatting tuple to lists of scan ids. The formatting tuple is defined by the user in the script, based on a boilerplate that heudiconv provides in its template. For reference, here is the function that creates the keys used in the ``infotodict`` output:\n",
    "\n",
    "```python\n",
    "def create_key(template, outtype=('nii.gz',), annotation_classes=None):\n",
    "    if template is None or not template:\n",
    "        raise ValueError('Template must be a valid format string')\n",
    "    return template, outtype, annotation_classes\n",
    "```\n",
    "\n",
    "The return value is hashable and can be used as a key. The user creates one key per scan type at the beginning of their implementation of ``infotodict``. Here is an example from the GRMPY_822831 project in the BBL group:\n",
    "\n",
    "```python\n",
    "def infotodict(seqinfo):\n",
    "    \"\"\"Heuristic evaluator for determining which runs belong where\n",
    "\n",
    "    allowed template fields - follow python string module:\n",
    "\n",
    "    item: index within category\n",
    "    subject: participant id\n",
    "    seqitem: run number during scanning\n",
    "    subindex: sub index within group\n",
    "    \"\"\"\n",
    "\n",
    "    last_run = len(seqinfo)\n",
    "    \n",
    "    # Create Keys\n",
    "    t1w = create_key(\n",
    "       'sub-{subject}/{session}/anat/sub-{subject}_{session}_T1w')\n",
    "    t2w = create_key(\n",
    "       'sub-{subject}/{session}/anat/sub-{subject}_{session}_T2w')\n",
    "    dwi = create_key(\n",
    "       'sub-{subject}/{session}/dwi/sub-{subject}_{session}_acq-multiband_dwi')\n",
    "    \n",
    "    # Field maps\n",
    "    b0_phase = create_key(\n",
    "       'sub-{subject}/{session}/fmap/sub-{subject}_{session}_phasediff')\n",
    "    b0_mag = create_key(\n",
    "       'sub-{subject}/{session}/fmap/sub-{subject}_{session}_magnitude')\n",
    "    pe_rev = create_key(\n",
    "        'sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-multiband_dir-j_epi')\n",
    "    # fmri scans\n",
    "    rest_mb = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-multiband_bold')\n",
    "    rest_sb = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-singleband_bold')\n",
    "    fracback = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-fracback_acq-singleband_bold')\n",
    "    face = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-face_acq-singleband_bold')\n",
    "\n",
    "    # ASL scans\n",
    "    asl = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_asl')\n",
    "    asl_dicomref = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_acq-ref_asl')\n",
    "    m0 = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_m0')\n",
    "    mean_perf = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_mean-perfusion')\n",
    "```\n",
    "\n",
    "\n",
    "The object that gets returned can then be defined as \n",
    "\n",
    "```python\n",
    "    info = {t1w:[], t2w:[], dwi:[], b0_phase:[], \n",
    "            b0_mag:[], pe_rev:[], rest_mb:[], rest_sb:[],\n",
    "            fracback:[], asl_dicomref:[], face:[], asl:[],\n",
    "            m0:[], mean_perf:[]}\n",
    "\n",
    "```\n",
    "\n",
    "Finally, the ``SeqInfo`` object is looped over and the scan_ids are added to the lists in the ``info`` dictionary. Here is part of an example loop:\n",
    "\n",
    "```python\n",
    "\n",
    "    def get_latest_series(key, s):\n",
    "        if len(info[key]) == 0:\n",
    "            info[key].append(s.series_id)\n",
    "        else:\n",
    "            if info[key] < s.series_id:\n",
    "                info[key] = [s.series_id]\n",
    "\n",
    "    for s in seqinfo:\n",
    "        protocol = s.protocol_name.lower()\n",
    "        if \"mprage\" in protocol:\n",
    "            get_latest_series(t1w,s)\n",
    "        elif \"t2_sag\" in protocol:\n",
    "            get_latest_series(t2w,s)\n",
    "        elif \"b0map\" in protocol and \"M\" in s.image_type:\n",
    "            info[b0_mag].append(s.series_id)\n",
    "        elif \"b0map\" in protocol and \"P\" in s.image_type:\n",
    "            info[b0_phase].append(s.series_id)\n",
    "        elif \"topup_ref\" in protocol:\n",
    "            get_latest_series(pe_rev, s)\n",
    "        elif \"dti_multishell\" in protocol and not s.is_derived:\n",
    "            get_latest_series(dwi, s)\n",
    "                elif s.series_description.endswith(\"_ASL\"):\n",
    "            get_latest_series(asl, s)\n",
    "        elif protocol.startswith(\"asl_dicomref\"):\n",
    "            get_latest_series(asl_dicomref, s)\n",
    "        elif s.series_description.endswith(\"_M0\"):\n",
    "            get_latest_series(m0, s)\n",
    "        elif s.series_description.endswith(\"_MeanPerf\"):\n",
    "            get_latest_series(mean_perf, s)\n",
    "\n",
    "        elif \"fracback\" in protocol:\n",
    "            get_latest_series(fracback, s)\n",
    "        elif \"face\" in protocol:\n",
    "            get_latest_series(face,s)\n",
    "        elif \"rest\" in protocol:\n",
    "            if \"MB\" in s.image_type:\n",
    "                get_latest_series(rest_mb,s)\n",
    "            else:\n",
    "                get_latest_series(rest_mb,s)\n",
    "    return info\n",
    "```\n",
    "\n",
    "Note that the ``get_latest_series`` function is used to address the situation where the a scan is interrupted and the interrupted scan is then followed by a complete version of the scan. It is common to have both the aborted and complete scan in a session, so this will always use the latest version of the scan. It is also possible to append multiple series_ids to the list and their ordering is formatted with their position in the list (like ``run-{item:02d}`` in the format string). \n",
    "\n",
    "It is then straightforward to take the ``info`` dictionary and update the bids fields in flywheel. \n",
    "\n",
    "\n",
    "## Goals\n",
    "\n",
    "If possible, it would be a major achievement to write a gear or python package that can take a heuristic file that works with ``heudiconv`` and have it work without any changes in flywheel. Since it won't have to do dicom conversion or actual file renaming, the only actual work to do is to:\n",
    "\n",
    "  1. Write an implementation of the ``SeqInfo`` class that is built from flywheel scan info\n",
    "  2. Write an implementation of the file renaming that updates scan metadata on flywheel\n",
    "  \n",
    "## Reference materials\n",
    "\n",
    "Here is the full ``grmpy_heuristic.py`` file:\n",
    "\n",
    "```python\n",
    "import os\n",
    "\n",
    "\n",
    "def create_key(template, outtype=('nii.gz',), annotation_classes=None):\n",
    "    if template is None or not template:\n",
    "        raise ValueError('Template must be a valid format string')\n",
    "    return template, outtype, annotation_classes\n",
    "\n",
    "\n",
    "def infotodict(seqinfo):\n",
    "    \"\"\"Heuristic evaluator for determining which runs belong where\n",
    "\n",
    "    allowed template fields - follow python string module:\n",
    "\n",
    "    item: index within category\n",
    "    subject: participant id\n",
    "    seqitem: run number during scanning\n",
    "    subindex: sub index within group\n",
    "    \"\"\"\n",
    "\n",
    "    last_run = len(seqinfo)\n",
    "    \n",
    "    # Create Keys\n",
    "    t1w = create_key(\n",
    "       'sub-{subject}/{session}/anat/sub-{subject}_{session}_T1w')\n",
    "    t2w = create_key(\n",
    "       'sub-{subject}/{session}/anat/sub-{subject}_{session}_T2w')\n",
    "    dwi = create_key(\n",
    "       'sub-{subject}/{session}/dwi/sub-{subject}_{session}_acq-multiband_dwi')\n",
    "    \n",
    "    # Field maps\n",
    "    b0_phase = create_key(\n",
    "       'sub-{subject}/{session}/fmap/sub-{subject}_{session}_phasediff')\n",
    "    b0_mag = create_key(\n",
    "       'sub-{subject}/{session}/fmap/sub-{subject}_{session}_magnitude')\n",
    "    pe_rev = create_key(\n",
    "        'sub-{subject}/{session}/fmap/sub-{subject}_{session}_acq-multiband_dir-j_epi')\n",
    "\n",
    "    # fmri scans\n",
    "    rest_mb = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-multiband_bold')\n",
    "    rest_sb = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-rest_acq-singleband_bold')\n",
    "    fracback = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-fracback_acq-singleband_bold')\n",
    "    face = create_key(\n",
    "       'sub-{subject}/{session}/func/sub-{subject}_{session}_task-face_acq-singleband_bold')\n",
    "\n",
    "    # ASL scans\n",
    "    asl = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_asl')\n",
    "    asl_dicomref = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_acq-ref_asl')\n",
    "    m0 = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_m0')\n",
    "    mean_perf = create_key(\n",
    "       'sub-{subject}/{session}/asl/sub-{subject}_{session}_mean-perfusion')\n",
    "\n",
    "    info = {t1w:[], t2w:[], dwi:[], b0_phase:[], \n",
    "            b0_mag:[], pe_rev:[], rest_mb:[], rest_sb:[],\n",
    "            fracback:[], asl_dicomref:[], face:[], asl:[],\n",
    "            m0:[], mean_perf:[]}\n",
    "\n",
    "    def get_latest_series(key, s):\n",
    "        if len(info[key]) == 0:\n",
    "            info[key].append(s.series_id)\n",
    "        else:\n",
    "            if info[key] < s.series_id:\n",
    "                info[key] = [s.series_id]\n",
    "\n",
    "    for s in seqinfo:\n",
    "        protocol = s.protocol_name.lower()\n",
    "        if \"mprage\" in protocol:\n",
    "            get_latest_series(t1w,s)\n",
    "        elif \"t2_sag\" in protocol:\n",
    "            get_latest_series(t2w,s)\n",
    "        elif \"b0map\" in protocol and \"M\" in s.image_type:\n",
    "            info[b0_mag].append(s.series_id)\n",
    "        elif \"b0map\" in protocol and \"P\" in s.image_type:\n",
    "            info[b0_phase].append(s.series_id)\n",
    "        elif \"topup_ref\" in protocol:\n",
    "            get_latest_series(pe_rev, s)\n",
    "        elif \"dti_multishell\" in protocol and not s.is_derived:\n",
    "            get_latest_series(dwi, s)\n",
    "\n",
    "        elif s.series_description.endswith(\"_ASL\"):\n",
    "            get_latest_series(asl, s)\n",
    "        elif protocol.startswith(\"asl_dicomref\"):\n",
    "            get_latest_series(asl_dicomref, s)\n",
    "        elif s.series_description.endswith(\"_M0\"):\n",
    "            get_latest_series(m0, s)\n",
    "        elif s.series_description.endswith(\"_MeanPerf\"):\n",
    "            get_latest_series(mean_perf, s)\n",
    "\n",
    "        elif \"fracback\" in protocol:\n",
    "            get_latest_series(fracback, s)\n",
    "        elif \"face\" in protocol:\n",
    "            get_latest_series(face,s)\n",
    "        elif \"rest\" in protocol:\n",
    "            if \"MB\" in s.image_type:\n",
    "                get_latest_series(rest_mb,s)\n",
    "            else:\n",
    "                get_latest_series(rest_mb,s)\n",
    "    return info\n",
    "```\n",
    "\n",
    "And here is an example of the fields that go into the ``SeqInfo`"
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}