> For the complete documentation index, see [llms.txt](https://time-series-features.gitbook.io/pyspi/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://time-series-features.gitbook.io/pyspi/installing-and-using-pyspi/usage/walkthrough-tutorials/neuroimaging-fmri-time-series.md).

# Neuroimaging: fMRI Time Series

In this tutorial, we will be applying *pyspi* to data derived from the [UCLA Consortium for Neuropsychiatric Phenomics LA5c Study](https://openneuro.org/datasets/ds000030/versions/00016), which is an open-access data source. We will focus on blood-oxygen level-dependent signalling (BOLD) time-series obtained from functional magnetic resonance imaging (fMRI).

The data used in this tutorial is **four sets of four brain** regions recorded from an anonymous participant. All fMRI time-series in this dataset are comprised of **152 time points** each (i.e., T = 152).

The following cells are intended to be run sequentially, from top to bottom.

## 1. Preparing the environment

We will first load all the modules needed for this demonstration. You may need to `pip install` one or more packages first.

```python
import pandas as pd
import numpy as np
import requests
import dill
from os import chdir, getcwd
from scipy.stats import zscore
import matplotlib.pyplot as plt
from scipy.stats import zscore
import seaborn as sns
from pyspi.calculator import Calculator
from matplotlib import colors, cm
from copy import deepcopy
```

***

## 2. Preparing the fMRI multivariate time-series (MTS) data

Now that we have imported the necessary libraries, let's read in the four sets of fMRI data, each corresponding to a separate region of interest (ROI), into *pandas* DataFrames:

<pre class="language-python"><code class="lang-python"><strong># Store the set names
</strong><strong>sets = ["BOLD_fMRI_TS_set1", "BOLD_fMRI_TS_set2", "BOLD_fMRI_TS_set3", "BOLD_fMRI_TS_set4"]
</strong>base_url = "https://raw.githubusercontent.com/anniegbryant/CNS_2022/main/pyspi_tutorial/tutorial_example_data/"

# Download each of the sets and save locally
for se in sets:
    url = base_url + se + ".csv"
    r = requests.get(url)
    fname = se + ".csv"
    with open(fname, 'wb') as f:
        f.write(r.content)
    f.close()

# Load in the four sets (each 4 ROIs with 152 time points) of raw BOLD fMRI time-series
BOLD_fMRI_TS_set1 = pd.read_csv("BOLD_fMRI_TS_set1.csv",
                   header = None)

BOLD_fMRI_TS_set2 = pd.read_csv("BOLD_fMRI_TS_set2.csv",
                   header = None)

BOLD_fMRI_TS_set3 = pd.read_csv("BOLD_fMRI_TS_set3.csv",
                   header = None)

BOLD_fMRI_TS_set4 = pd.read_csv("BOLD_fMRI_TS_set4.csv",
                   header = None)
# We can arbitrarily label our brain regions ROI1, ROI2, ROI3, ROI4
region_labels = ["ROI1", "ROI2", "ROI3", "ROI4"]
</code></pre>

Now, let's inspect the first set (set1) to see how the data is set up:

```
BOLD_fMRI_TS_set1
```

You should obtain the following dataframe. As per the output, the dataframe is 4 rows by 152 columns. Each row corresponds to one of the four brain regions, and each column corresponds to a measurement at a single time point:

<figure><img src="https://636746764-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FIw3ORxNbDkeyBcdB5svU%2Fuploads%2FHxt0ttVpxnBFWduY02Ok%2FScreenshot%202023-12-16%20at%204.16.23%20pm.png?alt=media&amp;token=5b7c7fd1-7046-449d-8bd3-06521b23934b" alt=""><figcaption></figcaption></figure>

This is the format that *pyspi* expects: **processes** (e.g., brain regions) **as the rows** and **time-points as the columns**. If you have different input data that is transposed, you can easily transpose the pandas DataFrame using the `transpose()` function.

To continue our exploratory data analysis, we can view the raw time series values for the first set of regions:

```python
def plot_data_lines(data, labels, title):
    plt.figure()
    data.transpose().plot(colormap=cm.jet)

    ax = plt.gca()
    ax.legend(labels=labels, loc='center left', bbox_to_anchor=(1, 0.5))
    plt.title(title)
    plt.xlabel('Time (fMRI frame)')
    plt.ylabel('BOLD Signal')
    plt.show()
    
plot_data_lines(BOLD_fMRI_TS_set1, title="Region Set 1", labels = region_labels)
```

<figure><img src="https://636746764-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FIw3ORxNbDkeyBcdB5svU%2Fuploads%2FczIZWWiL0uzCPRMpX9Pf%2FBOLD_fMRI_TS_Plot.svg?alt=media&amp;token=4d92cfa9-98c0-402b-b6be-792bfcc54d3f" alt="" width="461"><figcaption></figcaption></figure>

Alternatively, we can visualise the time-series for this first set of brain regions as a heatmap:

```python
# Plotting the data as a heatmap
def plot_data_heatmap(data, labels, title):
    plt.subplots()
    plt.pcolormesh(data,cmap=sns.color_palette('icefire',as_cmap=True))
    plt.colorbar()
    plt.title(title)
    ticks = [t+0.5 for t in range(len(labels))]
    plt.yticks(ticks=ticks, labels=labels)
    plt.xlabel('Time (fMRI frame)')
    plt.ylabel('Brain Region')
    plt.show()

plot_data_heatmap(BOLD_fMRI_TS_set1, labels = region_labels, title="Region Set 1")
```

You should obtain a heatmap like the following:

<figure><img src="https://636746764-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FIw3ORxNbDkeyBcdB5svU%2Fuploads%2FkTOgMRm5NoeAFVtS1VhA%2FBOLD_fMRI_TS_heatamp.svg?alt=media&amp;token=55a46795-f4f9-4e14-bc96-83bbbea46943" alt="" width="563"><figcaption></figcaption></figure>

You can also use the same two plotting functions -- `plot_data_lines()` and `plot_data_heatmap()` to visualise the raw time-series data for the other three sets of regions.

***

## 3. Run *pyspi* on a single dataset

Our data is now downloaded and ready to process with *pyspi.* There are two main steps to run *pyspi* and extract all of our pairwise statistics:

1. Initialise the `Calculator` object; and
2. Compute the `Calculator` object.

In the first step, you have some flexibility to play around with which SPIs you want to compute from the input time-series data. By default, *pyspi* will compute all available SPIs. However, this can take a long time depending on how many observations you have. If you want to give pyspi a try with a reduced feature set that runs quickly, you can pass `subset='fast'` as an additional parameter when instantiating the `Calculator:`

```python
# First, ensure there is no NAN's inthe dataset
BOLD_fMRI_TS_set1 = np.nan_to_num(BOLD_fMRI_TS_set1)

# Instantiate calculator object
calc_set1 = Calculator(BOLD_fMRI_TS_set1, subset="fast")

# Compute the instantiated calculator object
calc_set1.compute()
```

***OPTIONAL:** If you want to use the pre-computed calc\_set1 object, you can simply download the file: pyspi\_calc\_set1.pkl and **save** it in the same location as your jupyter notebook.*

{% file src="/files/EJPawcfTikhzC68LWNXL" %}
Download the pre-computed calculator object for set 1.
{% endfile %}

To load the .*pkl* file, do the following:

```python
with open('pyspi_calc_set1.pkl', 'rb') as f:
    calc_set1 = dill.load(f)
```

We can then inspect the resulting statistical pairwise interactions (SPIs) using the `calc_set1.table` object:

```python
print(calc_set1.table)
```

This output contains the resulting values from all of the SPIs concatenated together. We can view the list of SPIs that we calculated using the `.keys()` method:

```python
calc_set1.spis.keys()
```

Using their corresponding key/identifier, we can isolate one of the SPIs and visualise the results across the brain regions. For example, let's look at the matrix of pairwise interactions (MPI) for the covariance (identifier `cov_EmpiricalCovariance`) :

```python
def plot_mpi(S,identifier,labels,ax=None):
    """ Plot a given matrix of pairwise interactions, annotating the process labels and identifier
    """
    if ax is None:
        _, ax = plt.subplots()
    plt.sca(ax)

    # Use a diverging cmap if our statistic goes negative (and a sequential cmap otherwise)
    if np.nanmin(S) < 0.:
        maxabsval = max(abs(np.nanmin(S)),abs(np.nanmax(S)))
        norm = colors.Normalize(vmin=-maxabsval, vmax=maxabsval)
        plt.imshow(S,cmap='coolwarm',norm=norm)
    else:
        plt.imshow(S,cmap='Reds',vmin=0)

    plt.xticks(ticks=range(len(labels)),labels=labels,rotation=90)
    plt.yticks(ticks=range(len(labels)),labels=labels)
    plt.xlabel('Brain Region')
    plt.ylabel('Brain Region')
    plt.title(identifier)
    plt.colorbar()


# Plot this dataframe
plot_mpi(S = calc_set1.table["cov_EmpiricalCovariance"], identifier = "cov_EmpiricalCovariance", labels = region_labels)
```

You should obtain the following plot of the MPI:

<figure><img src="https://636746764-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FIw3ORxNbDkeyBcdB5svU%2Fuploads%2FvP4d6ofX4oyTvNCCuSrA%2FMPI_cov.svg?alt=media&amp;token=bd85aa1e-dab1-4b8c-8aa5-c6741b6cbf4c" alt=""><figcaption></figcaption></figure>

***

## 4. Running *pyspi* for multiple datasets

In the previous section, we ran through the *pyspi* pipeline for one dataset. In practice, you may have several datasets that you wish to process with *pyspi*, so here we cover tips for iterative processing.

We begin by creating a dictionary containing our other three datasets to process:

```python
datasets_to_process = {"Set2": BOLD_fMRI_TS_set2,
            "Set3": BOLD_fMRI_TS_set3,
            "Set4": BOLD_fMRI_TS_set4}

# We can initialise one calculator and copy it for each dataset to save time
calc = Calculator(subset="fast")

results = {}

# Iterate over each dataset
for key in datasets_to_process:

    # copy over the top-level calculator
    mycalc = deepcopy(calc)
    
    # Name the calculator after the ROI set
    mycalc.name = key[0]
    
    # Ensure we have no NaNs in our dataset
    dataset = np.nan_to_num(datasets_to_process[key])
    
    # Load the dataset
    mycalc.load_dataset(dataset)
    
    # Compute all pairwise interactions
    mycalc.compute()
    
    # Store our results
    results[key] = mycalc.table
    
```

This combines the SPI tables from each set into a single dictionary. This step takes about 2 minutes to run, so if you want to skip this part, you can load the pre-computed results dictionary as follows:

```python
with open('three_sets_results.pkl','rb') as f:
    results = dill.load(f)
```

We can then concatenate the dictionary of dataframes into one large dataframe as follows:

```
df_all_results = pd.concat(results, axis=0)
df_all_results
```

***

## 5. Downstream Analysis

### Exporting data to R

For users who use *R* for data wrangling and visualisation, we can save the *pyspi* output to a pickle file (`.pkl`) and write a custom function to load this data into *R* with the `reticulate()` package. We can practice with the `calc_set1` results, by first saving `calc_set1.table` to its own `.pkl` file:

```python
with open('pyspi_calc_set1_table.pkl', 'wb') as f:
    dill.dump(calc_set1.table, f)
```

We can then define a separate python script containing a function to extract the SPI data from our `pyspi_calc_table.pkl` file -- found in this repo as `pickle_reader_for_R.py`. Here are the contents:

```python
from pygments import highlight
from pygments.lexers import PythonLexer
from pygments.formatters import HtmlFormatter
import IPython

with open('R_interface/pickle_reader_for_R.py') as f:
    code = f.read()

formatter = HtmlFormatter()
IPython.display.HTML('<style type="text/css">{}</style>{}'.format(
    formatter.get_style_defs('.highlight'),
    highlight(code, PythonLexer(), formatter)))
```


---

# Agent Instructions
This documentation is published with GitBook. GitBook is the documentation platform designed so that both humans and AI agents can read, navigate, and reason over technical content effectively. Learn more at gitbook.com.

## Querying This Documentation
If you need additional information that is not directly available in this page, you can query the documentation by asking a question.

Perform an HTTP GET request on the following URL with the `ask` and `goal` query parameters:

```
GET https://time-series-features.gitbook.io/pyspi/installing-and-using-pyspi/usage/walkthrough-tutorials/neuroimaging-fmri-time-series.md?ask=<question>&goal=<user_goal>
```

`ask` is the immediate question: it should be specific, self-contained, and written in natural language.
`goal` is what the user is ultimately trying to achieve, the reason they need the answer. Sharing it helps GitBook give you a better, more relevant answer. A goal is most helpful when it describes the outcome the user wants rather than restating the question. For example, with `ask=how do I create an API token`, a goal like `automate deployments from our CI pipeline` lets GitBook tailor the answer to that use case.

The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.