# HCP Suite

## Description

This repository provides

- a set of scripts to perform common operations and analyses on the Human Connectome Project's neuroimaging data with a focus on task-based fMRI
- a Python-based framework for connectome-based analysis (CPM) using [Ray](https://ray.io) to parallelize CPM which has become the main focus of this project

## Installation

Other than cloning this repository, you need to have bash installed (which is most likely the case if you use Linux, *BSD or even MacOS).
For the Python code, the arguably easiest and cleanest way is to set up a Python virtual environment and install the dependencies there:

``` shell
$ python3 -m venv ./hcp-suite-venv # Setup the virtual environment
$ source ./hcp-suite-venv/bin/activate # Activate the virtual environment
$ pip install pandas pingouin networkx nilearn nibabel ray # Install dependencies within the virtual environment
$ pip install ipython jupyterlab # Install interactive Python shells to run hcp-suite code in
```

## Usage: CPM tutorial

The following tutorial uses the gambling task as an example, the variable to be predicted is BMI.

### Overview

Python code in this repository is written to be used in an interactive Python shell. Either Jupyter Lab or Jupyer Notebook is recommended as plots and images are conveniently displayed in-line but any Python shell (e.g. iPython) should work.

Generally speaking, the procedure is as follows:

1. Downloading data of subjects to be included in the analysis
2. Parcellation of CIFTI files with `task.sh prepare` (at the code's current state, the HCP's folder structure is assumed)
3. Extraction of fMRI time series based on EV files provided by the HCP with `get_ev_timeseries()` from `hcpsuite.py`
4. Computing of correlation matrices via `compute_correlations()` from `hcpsuite.py`
5. Performing of CPM, the functions are provided by `cpm_ray.py`
6. Establishing of statistical significance by way of permutation testing (also `cpm_ray.py`)

The following code snippets are to be run in an interactive Python shell (e.g. the "Console" of Jupyter Lab):

#### Downloading data

We will use Amazon Web Services to download HCP data. Set up access to HCP data via Amazon Web Services by following their [documentation](https://wiki.humanconnectome.org/docs/How%20To%20Connect%20to%20Connectome%20Data%20via%20AWS.html). You should be provided with an AWS Access Key ID and Secret Access Key we are going to put into Python variables (in quotation marks) for easy access:

``` python
aws_access_key_id="Replace with your access key ID"
aws_secret_access_key="Replace with your secret access key"

# Also, specify a location you want to download files to. We will use this variable repeatedly.
data_dir = "./HCP_1200"
```

We now need a list of subject IDs to be included in the analysis. Save them in a simple text file with one line per subject ID. For the rest of this tutorial, we will use the gambling task as an example (the list of IDs will be called ```gambling_subjects.ids```). As preprocessing for resting-state data differs, we will provide resting-state specific instructions when needed (the list of IDs will be called ```rest_subjects.ids```).

``` python
import sys

sys.path.append("../hcp-suite/lib") # This assumes the hcp-suite directory is in the current working directory's parent directory. You can also use absolute paths.
from hcpsuite import *
from cpm_ray import *

subj_ids = load_list_from_file("./gambling_subjects.ids") # Change the path to gambling_subjects.ids if it resides elsewhere

# Check the number of subjects
len(subj_ids)
```

Next, we will use the function download_hcp_files() to download all needed files from the HCP's AWS storage. This function takes the following arguments in this order: List of subject IDs, download location, task name, AWS access key ID, and AWS secret access key. It will return a list of missing subjects we capture in the variable missing_subjs:

``` python
missing_subjs = download_hcp_files(subj_ids, data_dir, "gambling", aws_access_key_id, aws_secret_access_key)
```

The previous command will print a list of missing subjects (if any) since not all data from all subjects is available on AWS, unfortunately. You can either obtain data from other sources or remove the missing subjects from our list of subject IDs:


``` python
for missing_subj in missing_subjs:
    subj_ids.remove(missing_subj)
```

#### Parcellation

Parcellation involves combining voxels of the "raw" CIFTI files into parcels as specified by the combined cortical, subcortical, and cerebellar parcellation we dubbed "RenTianGlasser" after the authors of the individual parcellations.

``` python
parcellate_ciftis(subj_ids=subj_ids,
                  parcellation_fname="./hcp-suite/data/parcellations/RenTianGlasser.dlabel.nii",
                  task="gambling",
                  data_dir=data_dir) # If you have installed Connectome Workbench in a non-standard way,
                                     # so that 'wb_command' is not in your $PATH, add
                                     # 'wb_command="/replace/with/path/to/wb_command"'
```

#### Extraction of time series


``` python
ev_data_dict = get_ev_timeseries(subj_ids, ["win.txt"], task="gambling",
                                 runs=('LR', 'RL'),
                                 parcellation="RenTianGlasser",
                                 data_dir=data_dir) # If you have installed Connectome Workbench in a non-standard way,
                                                    # so that 'wb_command' is not in your $PATH, add
                                                    # 'wb_command="/replace/with/path/to/wb_command"'

# Now, we save the extraced time series as text files in a directory of our choice
# (in this case: ./GAMBLING_win)
save_data_dict(ev_data_dict, path="./GAMBLING_win")
```

#### Computation of correlation matrices

We continue in our Python shell to compute correlation matrices:

``` python
# We load the saved time series from the previous step
cd GAMBLING_win # save_data_dict() writes file names into file "ts_files" but without paths,
                # thus the easiest way is to change into the directory containing the files
time_series, time_series_files = load_time_series("./ts_files") # ... and read them from there

correlation_measure, correlation_matrices = compute_correlations(time_series, kind='tangent')
# Tangent in our experience provides the best results, but there are alternatives:
# https://nilearn.github.io/dev/modules/generated/nilearn.connectome.ConnectivityMeasure.html

# We then save the matrices into a single file in the NIfTI format for downstream processing
save_matrix(cifti_dim_to_nifti_dim(correlation_matrices), "./GAMBLING_win-tangent.nii.gz")
```

#### CPM

For actual CPM, we need to install and run [Ray](https://ray.io) (run this e.g. in your Python virtual environment as described in [Installation](#installation)):

``` shell
$ pip install ray
$ ray start --head --num-cpus=16 # Run this on your main node.
                                 # Processes take up a lot of RAM, be careful not to use too many CPUs
```

Optionally add more Ray nodes to form a cluster, see the [Ray documentation](https://docs.ray.io/en/latest/cluster/vms/user-guides/launching-clusters/on-premises.html) for details.

##### Merging behavioral/biometric data in Python

For our analysis, we need values from both the unrestricted and [restricted HCP data](https://www.humanconnectome.org/study/hcp-young-adult/document/restricted-data-usage), which are available as separate CSV files. For easier handling, we merge them into a single CSV file:

``` python
import pandas as pd
unrestricted = pd.read_csv("/path/to/unrestricted.csv")
restricted = pd.read_csv("/path/to/restricted.csv")
merged = pd.merge(restricted, unrestricted, on="Subject")
merged.to_csv("./merged.csv") # Save the merged DataFrame as a CSV file in the current directory
```

##### Loading and preparing data in Python

``` python
behav = 'BMI' # Which variable do we want to predict?
covars = ['Age_in_Yrs', 'Gender' ] # What variables do we want to correct for?

behav_data = get_behav_data("./path/to/merged.csv", ids) # Loading of behavioral and biometrical data as a Pandas DataFrame from a CSV file
# We need to binarize gender for our analysis
behav_data["Gender"] = behav_data["Gender"].replace("F", 0)
behav_data["Gender"] = behav_data["Gender"].replace("M", 1)

functional_data = convert_matrices_to_dataframe(nifti_dim_to_cifti_dim(get_nimg_data("./path/to/GAMBLING_win-tangent.nii.gz")), subj_ids) # Loading of correlation matrices as Pandas DataFrame
```

##### Starting the Ray handler

ray_handler() is a Python class through which data management and the starting and coordination of Ray Actors (i.e. the processes working in parallel) is being handled.

``` python
n_folds = 128 # In this example, we use 128 folds, which is a good starting point
n_perm = 1000 # How many permutations are we planning to do later on?

ray_handler = RayHandler(
    functional_data.copy(),
    behav_data.copy(),
    behav,
    covars,
    address="auto", # We assume that the main Ray process runs on the same host
    n_perm=n_perm,
) # You can safely ignore the PerformanceWarning messages

ray_handler.add_kfold_indices(n_folds, clean=True) # By setting "clean" to True, we remove twins from the fold so they don't predict each other
ray_handler.upload_data() # Functional and behavioral data are uploaded into the shared storage to which Ray Actors have access
```

##### Starting the analysis

First we define the jobs for the actors; a job is a Python list object consisting of the following items: "job type", "fold number", "permutation number". The permutation number for the actual, i.e. unpermutated, data is "-1".

``` python
job_list = [["fselection_and_prediction", fold, perm] for perm in [-1] for fold in range(n_folds)] # This is the job list without permutations

# If we wanted to also compute the permutations (which takes a very long time), the job list can be created as follows:
#job_list = [["fselection_and_prediction", fold, perm] for perm in [-1, *range(n_perm)] for fold in range(n_folds)]

ray_handler.start_actors(job_list) # Start computing
```

##### Monitoring progress and retrieving results

``` python
n_done, n_remaining, n_held = ray_handler.status() # Prints a status report (see screenshot)

results = ray_handler.get_results(n=100) # Retrieving a huge number of results (e.g. when performing permutation analysis)
                                         # and especially from distributed Ray actors can take a long time. Specifying the
                                         # number of results (e.g. n=100) to be retrieved from the results store at once
                                         # allows for a display of progress

# Optional: Save results into a file (NumPy format)
np.save("./GAMBLING_win-results.npy", results)

# This file can be used to restore results
results = np.load("./GAMBLING_win-results.npy", allow_pickle=True).item()
```

You might consider fetching results and saving them periodically with a simple loop:

``` python
sleep_time = 1200 # Sleep for 20 minutes and then rerun loop
n_remaining = 1 # Set to something > 0 to get the loop started
results_path = get_safe_path("./GAMBLING_win-results", ".npy")
while n_remaining > 0: # Run until no jobs to be fetched are remaining
    n_done, n_remaining, n_held = ray_handler.status()
    if n_held > 0:
        results = ray_handler.get_results(n=100)
        # BEWARE: the file in results_path will be overwritten without asking
        # but we have used get_safe_path for risk mitigation
        print("\nSaving results to {}...".format(results_path), end='', flush=True)
        np.save(results_path, results)
        print(" done.")
    else:
        print("\nNo results to fetch and save.")
    print("Sleeping for {} seconds...".format(sleep_time))
    sleep(sleep_time) # Will sleep for sleep_time seconds
```

Screenshot of ray_handler.status() showing 1773 actors computing simultaneously:
<img src=./data/img/screenshot-ray_handler.status.png alt="Screenshot of ray_handler.status() showing 1773 actors computing simultaneously" width="50%" height="50%">

#### Check for completion

Rarely single jobs or actors die before completion. Having run your analyses, you can check your results for completion as follows and rerun analyses as needed (the function ```check_for_completion``` will advise you on how do this):

``` python
incomplete_jobs = check_for_completion(results)
```

#### Presenting results

To plot observed values against values as predicted by the GLM, use plot_predictions():

``` python
perm = -1 # Permutation -1 selects unpermutated results
g = plot_predictions(results['prediction_results'][perm], tail="glm", color="green")
```

Example output:

<img src=./data/img/example_prediction_results_plot.png alt="Example output of plot_predictions()" width="25%" height="25%">

---

To plot permutation results, use plot_permutation_results(), which will take either a list of paths to saved results files (f.i. when you want to combine multiple permutation runs) or the prediction results dictionary, which we will use in the following example:

``` python
# plot_permutation_results() will automatically remove incomplete permutations
# and return/print basic descriptive statistics (minimum and maximum r value, total
# number of permutations, and the minimum p value to be achieved with these permutations

plot, min_value, max_value, count, min_p = plot_permutation_results(results['prediction_results'])
```

Example output:

<img src=./data/img/screenshot-plot_permutation_results.png alt="Example output of plot_permutation_results()" width="25%" height="25%">

For more presentation and plotting examples including static and interactive plots of predictive networks, see ```save.py``` in folder ```utils```.

#### Overlap of predictive networks

cpm_ray.py has a function ```get_overlap()``` to create overlap networks of different predictive networks. For example, if you have two different CPM results as generated by ```results = ray_handler.get_results()``` and saved in two files ```results_A.npy``` and ```results_B.npy```, you can do the following:

``` python
# Load results into Python
results_A = np.load("/path/to/results_A.npy", allow_pickle=True).item()
results_B = np.load("/path/to/results_B.npy", allow_pickle=True).item()

# Load coordinates of parcels for plotting
coords = np.loadtxt("/path/to/hcp-suite/data/parcellations/RenTianGlasser.coords")

# Use overlap() function, which takes as its main input a list of results, in this case
# [results_A, results_B]. You can specify a result as the odd one out (e.g. odd_one_out=0
# for results_A), which means that result_A's tails will be switched, so that the positive
# tail result A will be overlapped with the negative tail of result B and vice versa.
#
# overlap() returns a bunch of dictionaries with the tail (positive and negative predictive
# networks)  as a primary key and usually several subkeys. These are more or less
# self-explanatory.

overlap, degrees_sorted, top_n_edges, plot_dict, m_cons = get_overlap([results_A, results_B], odd_one_out=None, coords=coords, plot_top_nodes=True, top_n_edges=50)

```

## Task.sh usage

``` shell
Usage: task.sh <action> [action-specific arguments] [common options]
  Actions:
    prepare
      -s <arg>: specify file with subject IDs or a space-separated list of subjects
      -t <arg>: specify a task
      -p <arg>: specify a parcellation
      You might want to consider using -e to prevent exiting when one subject fails.
    design
      -s <arg>: specify file with subject IDs or a space-separated list of subjects
      -o <arg>: specify an output file (e.g. design.mat)
      -V <arg>: specify variables (e.g. "-V Age_in_Yrs -V Gender -V BMI"; repeatable)
      -f <arg>: specify CSV file containing the specified variables (repeatable; e.g.
                "-f restricted.csv -f unrestricted.csv")
    eb
      -s <arg>: specify file with subject IDs or a space-separated list of subjects
      -o <arg>: specify an output file (e.g. eb.csv)
      -f <arg>: specify CSV file containing the specified variables (repeatable; e.g.
                "-f restricted.csv -f unrestricted.csv")
    analyse
      -s <arg>: specify file with subject IDs or a space-separated list of subjects
      -t <arg>: specify a task
      -c <arg>: specify a COPE number
      -p <arg>: specify a parcellation
      -d <arg>: specify target directory
      [-n <arg>]: number of permutations (default: 5000)
      [-o <arg>]: additional PALM arguments
      Specify design, contrast and EB via additional PALM arguments or place them in
      <target directory>/../ as design_<task>.mat, design.con and eb_<task>.csv

  Common options:
    [-D <arg>]: specify the directory containing HCP subject data (overrides HCP variable in ../lib/common.sh
    [-l <log file>]: specify log file instead of default /tmp/${0}-${RANDOM}.log
    [-e]: turn off error handling (do not exit on error)
```
