diff --git a/_freeze/fcn/execute-results/html.json b/_freeze/fcn/execute-results/html.json
index 01263dd..e2b6b15 100644
--- a/_freeze/fcn/execute-results/html.json
+++ b/_freeze/fcn/execute-results/html.json
@@ -1,8 +1,8 @@
{
- "hash": "1ea02656b247e6be6413bed43f4e2641",
+ "hash": "8e142f10111df79d0c061abcbc192f41",
"result": {
"engine": "knitr",
- "markdown": "# Functional Connectivity {#sec-functional-connectivity}\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(arrow)\nlibrary(dplyr)\nlibrary(fs)\nlibrary(ggplot2)\nlibrary(mgc)\nlibrary(purrr)\nlibrary(readr)\nlibrary(stringr)\nlibrary(tidyr)\n```\n:::\n\n\nFunctional connectivity is a measurement of how activity in different regions of the brain relates to one another over time. It is calculated by reducing the voxel-wise BOLD timeseries into a smaller set of regional timecourses, and then calculating a measure of correlation between those regional timecourses. \n\nTo define these regions, researchers rely on \"brain parcellations\" or atlases. Some atlases divide the brain into discrete, non-overlapping anatomical regions. Other methods, like [Dictionary Learning for Functional Regions (DiFuMo)](https://nilearn.github.io/modules/generated/nilearn.regions.fetch_atlas_difumo.html), do not create a strict parcellation. Instead, DiFuMo extracts a set of \"soft\", overlapping functional networks based on massive amounts of fMRI data. When using these different approaches, it is important to remember that you are correlating activity between fundamentally different kinds of brain divisions.\n\nThis kit reviews the A2CPS functional connectivity derivatives that have been created by pre-defined parcellations and soft atlases. For functional connectivity derived from subject-specific independent components, see @sec-gift.\n\n## Starting Project\n\n### Locate data\n\nOn TACC, the data are stored underneath the releases. For example, data release v{{< var release.version >}} is underneath\n\n```bash\n{{< var release.root >}}/{{< var release.dir >}}\n```\n\nPaths in this kit are written relative to that release folder, so the neuroimaging data are underneath `mris`.\n\n\nThe functional connectivity derivatives are underneath `mris/derivatives/fcn`\n\n```bash\n$ ls mris/derivatives/fcn/\ncleaned confounds confounds.json connectivity connectivity.json timeseries timeseries.json\n```\n\nThe functional connectivity data is in a tabular format in the folder `connectivity`. The A2CPS dataset includes functional connectivity from several atlases, estimated using multiple methods. For details on the atlases, see the data dictionary `connectivity.json`.\n\nTo enable flexibility, the timeseries parcellations are also provided in the folder `timeseries`, with data dictionary `timeseries.json`. These enable analyses such as dynamic connectivity or calculation of alternative measures of connectivity (e.g., partial correlation).\n\nThe tabular data comprise parquet files that have been partitioned in a [hive style](https://hive.apache.org/). That is, subfolder names contain column information (e.g., participant label). \n\n```bash\n$ tree connectivity | head -n 20\nconnectivity\n├── sub=10003\n│ └── ses=V1\n│ ├── task=cuff\n│ │ └── run=1\n│ │ ├── atlas=difumo_dimension-1024_resolution-2mm\n│ │ │ ├── estimator=empirical\n│ │ │ │ └── part-0.parquet\n│ │ │ └── estimator=leodit_wolf\n│ │ │ └── part-0.parquet\n│ │ ├── atlas=difumo_dimension-64_resolution-2mm\n│ │ │ ├── estimator=empirical\n│ │ │ │ └── part-0.parquet\n│ │ │ └── estimator=leodit_wolf\n│ │ │ └── part-0.parquet\n│ │ ├── atlas=dmn\n│ │ │ ├── estimator=empirical\n│ │ │ │ └── part-0.parquet\n│ │ │ └── estimator=leodit_wolf\n│ │ │ └── part-0.parquet\n```\n\nThe timeseries were extracted from the NIfTI files in `cleaned`, which are the outputs of fMRIPrep after band-pass filtering (0.01 - 0.1Hz), quadratic detrending, and nuisance regression (@sadil_acute_2024). The nuisance regressors are stored in the folder `confounds`.\n\n### Extract data\n\nMost users will start with the functional connectivity data. For example, here we grab the connectivity associated with the Default Mode Network (DMN) atlas (the nodal coordinates were derived from @baliki_corticostriatal_2012). In this case, we'll restrict results to connectivities generated with the [empirical estimator](https://scikit-learn.org/stable/modules/generated/sklearn.covariance.EmpiricalCovariance.html#sklearn.covariance.EmpiricalCovariance).\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndmn <- open_dataset(\"data/pre-surgery/mris/derivatives/fcn/connectivity\") |>\n filter(atlas == \"dmn\") |>\n filter(estimator == \"empirical\") |>\n arrange(sub, ses, task, run, source, target) |> # for reproducibility\n select(sub, source, target, connectivity) |>\n collect()\nhead(dmn)\n```\n\n::: {.cell-output-display}\n\n\n| sub| source| target| connectivity|\n|-----:|------:|------:|------------:|\n| 10003| 1| 2| 0.1972440|\n| 10003| 1| 3| 0.0254901|\n| 10003| 1| 4| 0.1471590|\n| 10003| 2| 3| 0.0555853|\n| 10003| 2| 4| 0.4430639|\n| 10003| 3| 4| -0.0632595|\n\n
\n:::\n:::\n\n\nNotice that the `source` and `target` fields are simply integer indices for this atlas. As specified in `connectivity.json`, information about these regions is available in one of the A2CPS GitHub repos. That table can be read directly from a URL.\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndmn_labels <- read_csv(\n \"https://raw.githubusercontent.com/a2cps/functional_connectivity/3aa91a6c10d14dcc7d1fe9890e7a6db95d2aad8b/src/functional_connectivity/data/baliki.csv\"\n)\ndmn_labels\n```\n\n::: {.cell-output-display}\n\n\n| region|label | x| y| z|\n|------:|:-------|---:|---:|--:|\n| 1|mPFC | 2| 52| -2|\n| 2|rNAC | 10| 12| -8|\n| 3|rInsula | 40| -6| -2|\n| 4|S1/M1 | -32| -34| 66|\n\n
\n:::\n:::\n\n\nAfter reading in the labels, they can be merged with the functional connectivity results. \n\n\n::: {.cell}\n\n```{.r .cell-code}\ndmn_labeled <- dmn |>\n left_join(dmn_labels, by = join_by(source == region)) |>\n select(-source, -x, -y, -z) |>\n rename(source = label) |>\n left_join(dmn_labels, by = join_by(target == region)) |>\n select(-target, -x, -y, -z) |>\n rename(target = label)\n\nhead(dmn_labeled)\n```\n\n::: {.cell-output-display}\n\n\n| sub| connectivity|source |target |\n|-----:|------------:|:-------|:-------|\n| 10003| 0.1972440|mPFC |rNAC |\n| 10003| 0.0254901|mPFC |rInsula |\n| 10003| 0.1471590|mPFC |S1/M1 |\n| 10003| 0.0555853|rNAC |rInsula |\n| 10003| 0.4430639|rNAC |S1/M1 |\n| 10003| -0.0632595|rInsula |S1/M1 |\n\n
\n:::\n:::\n\n\n### Example Analysis: Discriminability\n\nIn this section, we show how the connectivity results could be used to calculate discriminability\\ [@bridgeford_eliminating_2021], which is a multivariate measure of replicability like the intra-class correlation coefficient or fingerprinting. In the context of connectomes, it answers the question: \"Is the functional connectivity profile of a participant distinct enough to identify them among a group of people?\" Discriminability ranges from 0 - 1, with 0.5 indicating something like an equal chance that the two scans from the same participants are as similar as the two scans from different participants, and 1 indicating that the two scans from the same participant are always more similar than two scans from differing participants.\n\nWe are going to assess whether a person's connectivity matrix is consistent across runs, even across runs of different types (e.g., CUFF1 vs CUFF2, REST1 vs CUFF1).\n\nFirst, we need a list of participants that have all four scan types.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsubs_with_all_runs <- open_dataset(\n \"data/pre-surgery/mris/derivatives/fcn/connectivity\"\n) |>\n distinct(sub, task, run) |>\n count(sub) |>\n filter(n == 4) |>\n select(sub) |>\n collect()\n```\n:::\n\n\nAs usual, we should also restrict analyses to only those runs that are not \"red\".\n\n\n::: {.cell}\n\n```{.r .cell-code}\nred_fmri <- read_tsv(\n dir_ls(\"data/pre-surgery/mris/bids\", glob = \"*scans.tsv\", recurse = TRUE),\n na = \"n/a\"\n) |>\n filter(str_detect(filename, \"task\")) |>\n filter(rating == \"red\") |>\n mutate(\n sub = str_extract(filename, \"(?<=sub-)[[:digit:]]{5}\") |> as.integer()\n ) |>\n distinct(sub)\nhead(red_fmri)\n```\n\n::: {.cell-output-display}\n\n\n| sub|\n|-----:|\n| 10058|\n| 10074|\n| 10077|\n| 10086|\n| 10097|\n| 10103|\n\n
\n:::\n:::\n\n\nOf the list of participants with all four functional scans, filter out the participants for which any of the scans were red.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsubs_with_all_runs_ok <- subs_with_all_runs |>\n anti_join(red_fmri, by = join_by(sub))\n```\n:::\n\n\nUse this to filter the connectivity results. We'll select just one of the smaller DiFuMo atlases\\ [@dadi_fine_2020]. As before, we'll stick with the empirical estimator.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nfcn <- open_dataset(\"data/pre-surgery/mris/derivatives/fcn/connectivity\") |>\n filter(atlas == \"difumo_dimension-64_resolution-2mm\") |>\n filter(estimator == \"empirical\") |>\n semi_join(subs_with_all_runs_ok, by = join_by(sub)) |>\n arrange(sub, task, run, source, target) |> # for reproducibility\n mutate(scan = str_c(task, run)) |>\n select(-atlas, -ses, -task, -run, -estimator) |>\n collect()\nhead(fcn)\n```\n\n::: {.cell-output-display}\n\n\n| source| target| connectivity| sub|scan |\n|------:|------:|------------:|-----:|:-----|\n| 1| 2| 0.0832526| 10010|cuff1 |\n| 1| 3| 0.6190895| 10010|cuff1 |\n| 1| 4| 0.4174389| 10010|cuff1 |\n| 1| 5| 0.1506552| 10010|cuff1 |\n| 1| 6| -0.0530966| 10010|cuff1 |\n| 1| 7| 0.6338826| 10010|cuff1 |\n\n
\n:::\n:::\n\n\nNext, define some helper functions to break up the different parts of the analysis pipeline.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nget_scan_combinations <- function(\n scans = c(\"rest1\", \"rest2\", \"cuff1\", \"cuff2\"),\n .col1 = scan1,\n .col2 = scan2\n) {\n combn(scans, 2) |>\n t() |>\n as_tibble() |>\n rename({{ .col1 }} := V1, {{ .col2 }} := V2)\n}\n\njoin_fcn_to_combinations <- function(.data, fcn) {\n fcn_nested <- group_nest(fcn, scan)\n .data |>\n left_join(fcn_nested, by = join_by(scan1 == scan)) |>\n left_join(fcn_nested, by = join_by(scan2 == scan)) |>\n mutate(\n data = map2(\n data.x,\n data.y,\n \\(x, y) left_join(x, y, by = join_by(source, target, sub))\n )\n ) |>\n select(-starts_with(\"data.\"))\n}\n\nget_discr <- function(.data) {\n d <- .data |>\n mutate(feature = interaction(source, target)) |>\n select(-source, -target) |>\n pivot_longer(starts_with(\"connectivity\")) |>\n pivot_wider(names_from = \"feature\")\n\n discr.stat(as.matrix(select(d, -sub, -name)), as.matrix(select(d, sub)))$discr\n}\n```\n:::\n\n\nApply these helper functions to the functional connectivity data, calculating discriminability.\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndiscriminability <- get_scan_combinations() |>\n join_fcn_to_combinations(fcn) |>\n mutate(discr = map_dbl(data, get_discr)) |>\n select(-data)\n```\n:::\n\n\nTo review the results, plot the data. When plotting, let's color the points based on whether the two scans are of the same type.\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndiscriminability |>\n mutate(\n same_type = (str_detect(scan1, \"rest\") & str_detect(scan2, \"rest\")) |\n (str_detect(scan1, \"cuff\") & str_detect(scan2, \"cuff\")),\n scans = interaction(scan1, scan2)\n ) |>\n ggplot(aes(x = scans, y = discr, color = same_type)) +\n geom_point() +\n coord_flip() +\n ylim(0, 1)\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nOverall, discriminability is around 0.8, with at most minor differences between different pairings of runs.\n\n## Considerations While Working on the Project\n\n### Variability Across Scanners\n\nMany MRI biomarkers exhibit variability across the scanners, which may confound some analyses. For an up-to-date assessment of the issue and overview of current thinking, please see @sec-mri-harmonization.\n\n\n### Data Quality\n\nAs with any MRI derivative, all pipeline derivatives have been included. This means that products were included regardless of their quality, and so some products may have been generated from images that are known to have poor quality---rated \"red\", or incomparable. For details on the ratings and how to exclude them, see @sec-rawdata-mri-qc-joining. Additionally, extensive QC has not yet been performed on the derivatives themselves, and so there may be cases where pipelines produced atypical outputs. For an overview of planned checks, see [Confluence](https://a2cps.atlassian.net/wiki/spaces/DOC/pages/262471688/Image+Processing+Checks).\n\n\n### Citations\n\nIn publications or presentations including data from A2CPS, please include the following statement as attribution:\n\n> Data were provided (in part) by the A2CPS Consortium funded by the National Institutes of Health (NIH) Common Fund, which is managed by the Office of the Director (OD)/Office of Strategic Coordination (OSC). Consortium components and their associated funding sources include Clinical Coordinating Center (U24NS112873), Data Integration and Resource Center (U54DA049110), Omics Data Generation Centers (U54DA049116, U54DA049115, U54DA049113), Multi-site Clinical Center 1 (MCC1) (UM1NS112874), and Multi-site Clinical Center 2 (MCC2) (UM1NS118922).\n\n::: {.callout-note}\nThe following published papers should be cited when referring to A2CPS Protocol and Biomarkers: @sluka_predicting_2023 @berardi_multi_2022\n:::\n\n\nWhen using neuroimaging derivatives, please also cite @sadil_acute_2024.\n",
+ "markdown": "# Functional Connectivity {#sec-functional-connectivity}\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(arrow)\nlibrary(dplyr)\nlibrary(fs)\nlibrary(readr)\nlibrary(ggplot2)\nlibrary(mgc)\nlibrary(purrr)\nlibrary(readr)\nlibrary(stringr)\nlibrary(tidyr)\n```\n:::\n\n\nFunctional connectivity is a measurement of how activity in different regions of the brain relates to one another over time. It is calculated by reducing the voxel-wise BOLD timeseries into a smaller set of regional timecourses, and then calculating a measure of correlation between those regional timecourses.\n\nTo define these regions, researchers rely on \"brain parcellations\" or atlases. Some atlases divide the brain into discrete, non-overlapping anatomical regions. Other methods, like [Dictionary Learning for Functional Regions (DiFuMo)](https://nilearn.github.io/modules/generated/nilearn.regions.fetch_atlas_difumo.html), do not create a strict parcellation. Instead, DiFuMo extracts a set of \"soft\", overlapping functional networks based on massive amounts of fMRI data. When using these different approaches, it is important to remember that you are correlating activity between fundamentally different kinds of brain divisions.\n\nThis kit reviews the A2CPS functional connectivity derivatives that have been created by pre-defined parcellations and soft atlases. For functional connectivity derived from subject-specific independent components, see @sec-gift.\n\n## Starting Project\n\n### Locate data\n\nOn TACC, the data are stored underneath the releases. For example, data release v{{< var release.version >}} is underneath\n\n```bash\n{{< var release.root >}}/{{< var release.dir >}}\n```\n\nPaths in this kit are written relative to that release folder, so the neuroimaging data are underneath `mris`.\n\n\nThe functional connectivity derivatives are underneath `mris/derivatives/fcn`\n\n```bash\n$ ls mris/derivatives/fcn/\ncleaned confounds confounds.json connectivity connectivity.json timeseries timeseries.json\n```\n\nThe functional connectivity data is in a tabular format in the folder `connectivity`. The A2CPS dataset includes functional connectivity from several atlases, estimated using multiple methods. For details on the atlases, see the data dictionary `connectivity.json`.\n\nTo enable flexibility, the timeseries parcellations are also provided in the folder `timeseries`, with data dictionary `timeseries.json`. These enable analyses such as dynamic connectivity or calculation of alternative measures of connectivity (e.g., partial correlation).\n\nThe tabular data comprise parquet files that have been partitioned in a [hive style](https://hive.apache.org/). That is, subfolder names contain column information (e.g., participant label).\n\n```bash\n$ tree connectivity | head -n 20\nconnectivity\n├── sub=10003\n│ └── ses=V1\n│ ├── task=cuff\n│ │ └── run=1\n│ │ ├── atlas=difumo_dimension-1024_resolution-2mm\n│ │ │ ├── estimator=empirical\n│ │ │ │ └── part-0.parquet\n│ │ │ └── estimator=leodit_wolf\n│ │ │ └── part-0.parquet\n│ │ ├── atlas=difumo_dimension-64_resolution-2mm\n│ │ │ ├── estimator=empirical\n│ │ │ │ └── part-0.parquet\n│ │ │ └── estimator=leodit_wolf\n│ │ │ └── part-0.parquet\n│ │ ├── atlas=dmn\n│ │ │ ├── estimator=empirical\n│ │ │ │ └── part-0.parquet\n│ │ │ └── estimator=leodit_wolf\n│ │ │ └── part-0.parquet\n```\n\nThe timeseries were extracted from the NIfTI files in `cleaned`, which are the outputs of fMRIPrep after band-pass filtering (0.01 - 0.1Hz), quadratic detrending, and nuisance regression (@sadil_acute_2024). The nuisance regressors are stored in the folder `confounds`.\n\n### Extract data\n\nMost users will start with the functional connectivity data. For example, here we grab the connectivity associated with the Default Mode Network (DMN) atlas (the nodal coordinates were derived from @baliki_corticostriatal_2012). In this case, we'll restrict results to connectivities generated with the [empirical estimator](https://scikit-learn.org/stable/modules/generated/sklearn.covariance.EmpiricalCovariance.html#sklearn.covariance.EmpiricalCovariance).\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndmn <- open_dataset(\"data/pre-surgery/mris/derivatives/fcn/connectivity\") |>\n filter(atlas == \"dmn\") |>\n filter(estimator == \"empirical\") |>\n arrange(sub, ses, task, run, source, target) |> # for reproducibility\n select(sub, source, target, connectivity) |>\n collect()\nhead(dmn)\n```\n\n::: {.cell-output-display}\n\n\n| sub| source| target| connectivity|\n|-----:|------:|------:|------------:|\n| 10003| 1| 2| 0.1972440|\n| 10003| 1| 3| 0.0254901|\n| 10003| 1| 4| 0.1471590|\n| 10003| 2| 3| 0.0555853|\n| 10003| 2| 4| 0.4430639|\n| 10003| 3| 4| -0.0632595|\n\n
\n:::\n:::\n\n\nNotice that the `source` and `target` fields are simply integer indices for this atlas. As specified in `connectivity.json`, information about these regions is available in one of the A2CPS GitHub repos. That table can be read directly from a URL.\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndmn_labels <- read_csv(\n \"https://raw.githubusercontent.com/a2cps/functional_connectivity/3aa91a6c10d14dcc7d1fe9890e7a6db95d2aad8b/src/functional_connectivity/data/baliki.csv\"\n)\ndmn_labels\n```\n\n::: {.cell-output-display}\n\n\n| region|label | x| y| z|\n|------:|:-------|---:|---:|--:|\n| 1|mPFC | 2| 52| -2|\n| 2|rNAC | 10| 12| -8|\n| 3|rInsula | 40| -6| -2|\n| 4|S1/M1 | -32| -34| 66|\n\n
\n:::\n:::\n\n\nAfter reading in the labels, they can be merged with the functional connectivity results.\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndmn_labeled <- dmn |>\n left_join(dmn_labels, by = join_by(source == region)) |>\n select(-source, -x, -y, -z) |>\n rename(source = label) |>\n left_join(dmn_labels, by = join_by(target == region)) |>\n select(-target, -x, -y, -z) |>\n rename(target = label)\n\nhead(dmn_labeled)\n```\n\n::: {.cell-output-display}\n\n\n| sub| connectivity|source |target |\n|-----:|------------:|:-------|:-------|\n| 10003| 0.1972440|mPFC |rNAC |\n| 10003| 0.0254901|mPFC |rInsula |\n| 10003| 0.1471590|mPFC |S1/M1 |\n| 10003| 0.0555853|rNAC |rInsula |\n| 10003| 0.4430639|rNAC |S1/M1 |\n| 10003| -0.0632595|rInsula |S1/M1 |\n\n
\n:::\n:::\n\n\n### Example Analysis: Discriminability\n\nIn this section, we show how the connectivity results could be used to calculate discriminability\\ [@bridgeford_eliminating_2021], which is a multivariate measure of replicability like the intra-class correlation coefficient or fingerprinting. In the context of connectomes, it answers the question: \"Is the functional connectivity profile of a participant distinct enough to identify them among a group of people?\" Discriminability ranges from 0 - 1, with 0.5 indicating something like an equal chance that the two scans from the same participants are as similar as the two scans from different participants, and 1 indicating that the two scans from the same participant are always more similar than two scans from differing participants.\n\nWe are going to assess whether a person's connectivity matrix is consistent across runs, even across runs of different types (e.g., CUFF1 vs CUFF2, REST1 vs CUFF1).\n\nFirst, we need a list of participants that have all four scan types.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsubs_with_all_runs <- open_dataset(\n \"data/pre-surgery/mris/derivatives/fcn/connectivity\"\n) |>\n distinct(sub, task, run) |>\n count(sub) |>\n filter(n == 4) |>\n select(sub) |>\n collect()\n```\n:::\n\n\nAs usual, we should also restrict analyses to only those runs that are not \"red\".\n\n\n::: {.cell}\n\n```{.r .cell-code}\nred_fmri <- read_tsv(\n dir_ls(\"data/pre-surgery/mris/bids\", glob = \"*scans.tsv\", recurse = TRUE),\n na = \"n/a\"\n) |>\n filter(str_detect(filename, \"task\")) |>\n filter(rating == \"red\") |>\n mutate(\n sub = str_extract(filename, \"(?<=sub-)[[:digit:]]{5}\") |> as.integer()\n ) |>\n distinct(sub)\nhead(red_fmri)\n```\n\n::: {.cell-output-display}\n\n\n| sub|\n|-----:|\n| 10058|\n| 10074|\n| 10077|\n| 10086|\n| 10097|\n| 10103|\n\n
\n:::\n:::\n\n\nOf the list of participants with all four functional scans, filter out the participants for which any of the scans were red.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsubs_with_all_runs_ok <- subs_with_all_runs |>\n anti_join(red_fmri, by = join_by(sub))\n```\n:::\n\n\nUse this to filter the connectivity results. We'll select just one of the smaller DiFuMo atlases\\ [@dadi_fine_2020]. As before, we'll stick with the empirical estimator.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nfcn <- open_dataset(\"data/pre-surgery/mris/derivatives/fcn/connectivity\") |>\n filter(atlas == \"difumo_dimension-64_resolution-2mm\") |>\n filter(estimator == \"empirical\") |>\n semi_join(subs_with_all_runs_ok, by = join_by(sub)) |>\n arrange(sub, task, run, source, target) |> # for reproducibility\n mutate(scan = str_c(task, run)) |>\n select(-atlas, -ses, -task, -run, -estimator) |>\n collect()\nhead(fcn)\n```\n\n::: {.cell-output-display}\n\n\n| source| target| connectivity| sub|scan |\n|------:|------:|------------:|-----:|:-----|\n| 1| 2| 0.0832526| 10010|cuff1 |\n| 1| 3| 0.6190895| 10010|cuff1 |\n| 1| 4| 0.4174389| 10010|cuff1 |\n| 1| 5| 0.1506552| 10010|cuff1 |\n| 1| 6| -0.0530966| 10010|cuff1 |\n| 1| 7| 0.6338826| 10010|cuff1 |\n\n
\n:::\n:::\n\n\nNext, define some helper functions to break up the different parts of the analysis pipeline.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nget_scan_combinations <- function(\n scans = c(\"rest1\", \"rest2\", \"cuff1\", \"cuff2\"),\n .col1 = scan1,\n .col2 = scan2\n) {\n combn(scans, 2) |>\n t() |>\n as_tibble() |>\n rename({{ .col1 }} := V1, {{ .col2 }} := V2)\n}\n\njoin_fcn_to_combinations <- function(.data, fcn) {\n fcn_nested <- group_nest(fcn, scan)\n .data |>\n left_join(fcn_nested, by = join_by(scan1 == scan)) |>\n left_join(fcn_nested, by = join_by(scan2 == scan)) |>\n mutate(\n data = map2(\n data.x,\n data.y,\n \\(x, y) left_join(x, y, by = join_by(source, target, sub))\n )\n ) |>\n select(-starts_with(\"data.\"))\n}\n\nget_discr <- function(.data) {\n d <- .data |>\n mutate(feature = interaction(source, target)) |>\n select(-source, -target) |>\n pivot_longer(starts_with(\"connectivity\")) |>\n pivot_wider(names_from = \"feature\")\n\n discr.stat(as.matrix(select(d, -sub, -name)), as.matrix(select(d, sub)))$discr\n}\n```\n:::\n\n\nApply these helper functions to the functional connectivity data, calculating discriminability.\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndiscriminability <- get_scan_combinations() |>\n join_fcn_to_combinations(fcn) |>\n mutate(discr = map_dbl(data, get_discr)) |>\n select(-data)\n```\n:::\n\n\nTo review the results, plot the data. When plotting, let's color the points based on whether the two scans are of the same type.\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndiscriminability |>\n mutate(\n same_type = (str_detect(scan1, \"rest\") & str_detect(scan2, \"rest\")) |\n (str_detect(scan1, \"cuff\") & str_detect(scan2, \"cuff\")),\n scans = interaction(scan1, scan2)\n ) |>\n ggplot(aes(x = scans, y = discr, color = same_type)) +\n geom_point() +\n coord_flip() +\n ylim(0, 1)\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\nOverall, discriminability is around 0.8, with at most minor differences between different pairings of runs.\n\n## Considerations While Working on the Project\n\n### Mislabeled Schaefer Atlas\n\nTwo of the pre-computed connectivity matrices are based on [Schaefer's atlas](https://github.com/ThomasYeoLab/CBIG/tree/v0.39.0-Zhang2026_EIAD/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal). Unfortunately, the labels are stored incorrectly, as described in [this issue](https://github.com/a2cps/biomarker-extractor/issues/16). The fix is relatively straightforward. Here, we download the file referenced in that issue, which contains a mapping between the current and correct labels.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlabel_mapping <- arrow::read_tsv_arrow(\n \"https://github.com/user-attachments/files/22352624/labels_mapping.tsv\",\n schema = schema(current = int64(), correct = int64()), # specifying to make the join below easier\n skip = 1,\n as_data_frame = FALSE\n)\nhead(label_mapping)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nTable\n6 rows x 2 columns\n$current \n$correct \n```\n\n\n:::\n:::\n\n\nWe can use those labels to remap the existing source and target columns, as in the code chunk below. Note that the chunk below does this for the `schaefer_nrois-400_resolution-2_networks-7` atlas, but it will also need to be done for the `schaefer_nrois-400_resolution-2_networks-17` atlas.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nopen_dataset(\n \"data/pre-surgery/mris/derivatives/fcn/connectivity\"\n) |>\n filter(atlas == \"schaefer_nrois-400_resolution-2_networks-7\") |>\n left_join(label_mapping, by = join_by(source == current)) |>\n left_join(\n label_mapping,\n by = join_by(target == current),\n suffix = c(\".source\", \".target\")\n ) |>\n select(-source, -target) |>\n rename(source = correct.source, target = correct.target) |>\n head() |>\n collect()\n```\n\n::: {.cell-output-display}\n\n\n| connectivity| sub|ses |task | run|atlas |estimator | source| target|\n|------------:|-----:|:---|:----|---:|:------------------------------------------|:-----------|------:|------:|\n| 0.6768990| 10003|V1 |cuff | 1|schaefer_nrois-400_resolution-2_networks-7 |ledoit_wolf | 1| 2|\n| 0.6125918| 10003|V1 |cuff | 1|schaefer_nrois-400_resolution-2_networks-7 |ledoit_wolf | 1| 3|\n| 0.6882472| 10003|V1 |cuff | 1|schaefer_nrois-400_resolution-2_networks-7 |ledoit_wolf | 1| 4|\n| 0.5958827| 10003|V1 |cuff | 1|schaefer_nrois-400_resolution-2_networks-7 |ledoit_wolf | 1| 5|\n| 0.4730263| 10003|V1 |cuff | 1|schaefer_nrois-400_resolution-2_networks-7 |ledoit_wolf | 1| 6|\n| 0.1177192| 10003|V1 |cuff | 1|schaefer_nrois-400_resolution-2_networks-7 |ledoit_wolf | 1| 7|\n\n
\n:::\n:::\n\n\n### Variability Across Scanners\n\nMany MRI biomarkers exhibit variability across the scanners, which may confound some analyses. For an up-to-date assessment of the issue and overview of current thinking, please see @sec-mri-harmonization.\n\n\n### Data Quality\n\nAs with any MRI derivative, all pipeline derivatives have been included. This means that products were included regardless of their quality, and so some products may have been generated from images that are known to have poor quality---rated \"red\", or incomparable. For details on the ratings and how to exclude them, see @sec-rawdata-mri-qc-joining. Additionally, extensive QC has not yet been performed on the derivatives themselves, and so there may be cases where pipelines produced atypical outputs. For an overview of planned checks, see [Confluence](https://a2cps.atlassian.net/wiki/spaces/DOC/pages/262471688/Image+Processing+Checks).\n\n\n### Citations\n\nIn publications or presentations including data from A2CPS, please include the following statement as attribution:\n\n> Data were provided (in part) by the A2CPS Consortium funded by the National Institutes of Health (NIH) Common Fund, which is managed by the Office of the Director (OD)/Office of Strategic Coordination (OSC). Consortium components and their associated funding sources include Clinical Coordinating Center (U24NS112873), Data Integration and Resource Center (U54DA049110), Omics Data Generation Centers (U54DA049116, U54DA049115, U54DA049113), Multi-site Clinical Center 1 (MCC1) (UM1NS112874), and Multi-site Clinical Center 2 (MCC2) (UM1NS118922).\n\n::: {.callout-note}\nThe following published papers should be cited when referring to A2CPS Protocol and Biomarkers: @sluka_predicting_2023 @berardi_multi_2022\n:::\n\n\nWhen using neuroimaging derivatives, please also cite @sadil_acute_2024.\n",
"supporting": [],
"filters": [
"rmarkdown/pagebreak.lua"
diff --git a/fcn.qmd b/fcn.qmd
index 4c0f5da..705cc1d 100644
--- a/fcn.qmd
+++ b/fcn.qmd
@@ -6,6 +6,7 @@
library(arrow)
library(dplyr)
library(fs)
+library(readr)
library(ggplot2)
library(mgc)
library(purrr)
@@ -14,7 +15,7 @@ library(stringr)
library(tidyr)
```
-Functional connectivity is a measurement of how activity in different regions of the brain relates to one another over time. It is calculated by reducing the voxel-wise BOLD timeseries into a smaller set of regional timecourses, and then calculating a measure of correlation between those regional timecourses.
+Functional connectivity is a measurement of how activity in different regions of the brain relates to one another over time. It is calculated by reducing the voxel-wise BOLD timeseries into a smaller set of regional timecourses, and then calculating a measure of correlation between those regional timecourses.
To define these regions, researchers rely on "brain parcellations" or atlases. Some atlases divide the brain into discrete, non-overlapping anatomical regions. Other methods, like [Dictionary Learning for Functional Regions (DiFuMo)](https://nilearn.github.io/modules/generated/nilearn.regions.fetch_atlas_difumo.html), do not create a strict parcellation. Instead, DiFuMo extracts a set of "soft", overlapping functional networks based on massive amounts of fMRI data. When using these different approaches, it is important to remember that you are correlating activity between fundamentally different kinds of brain divisions.
@@ -37,7 +38,7 @@ The functional connectivity data is in a tabular format in the folder `connectiv
To enable flexibility, the timeseries parcellations are also provided in the folder `timeseries`, with data dictionary `timeseries.json`. These enable analyses such as dynamic connectivity or calculation of alternative measures of connectivity (e.g., partial correlation).
-The tabular data comprise parquet files that have been partitioned in a [hive style](https://hive.apache.org/). That is, subfolder names contain column information (e.g., participant label).
+The tabular data comprise parquet files that have been partitioned in a [hive style](https://hive.apache.org/). That is, subfolder names contain column information (e.g., participant label).
```bash
$ tree connectivity | head -n 20
@@ -92,7 +93,7 @@ dmn_labels <- read_csv(
dmn_labels
```
-After reading in the labels, they can be merged with the functional connectivity results.
+After reading in the labels, they can be merged with the functional connectivity results.
```{r}
#| label: dmn-labeled
@@ -246,6 +247,41 @@ Overall, discriminability is around 0.8, with at most minor differences between
## Considerations While Working on the Project
+### Mislabeled Schaefer Atlas
+
+Two of the pre-computed connectivity matrices are based on [Schaefer's atlas](https://github.com/ThomasYeoLab/CBIG/tree/v0.39.0-Zhang2026_EIAD/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal). Unfortunately, the labels are stored incorrectly, as described in [this issue](https://github.com/a2cps/biomarker-extractor/issues/16). The fix is relatively straightforward. Here, we download the file referenced in that issue, which contains a mapping between the current and correct labels.
+
+```{r}
+#| label: read-remap-labels
+
+label_mapping <- arrow::read_tsv_arrow(
+ "https://github.com/user-attachments/files/22352624/labels_mapping.tsv",
+ schema = schema(current = int64(), correct = int64()), # specifying to make the join below easier
+ skip = 1,
+ as_data_frame = FALSE
+)
+head(label_mapping)
+```
+
+We can use those labels to remap the existing source and target columns, as in the code chunk below. Note that the chunk below does this for the `schaefer_nrois-400_resolution-2_networks-7` atlas, but it will also need to be done for the `schaefer_nrois-400_resolution-2_networks-17` atlas.
+
+```{r}
+open_dataset(
+ "data/pre-surgery/mris/derivatives/fcn/connectivity"
+) |>
+ filter(atlas == "schaefer_nrois-400_resolution-2_networks-7") |>
+ left_join(label_mapping, by = join_by(source == current)) |>
+ left_join(
+ label_mapping,
+ by = join_by(target == current),
+ suffix = c(".source", ".target")
+ ) |>
+ select(-source, -target) |>
+ rename(source = correct.source, target = correct.target) |>
+ head() |>
+ collect()
+```
+
### Variability Across Scanners
{{< include _snippets/mri-scanner-variability.qmd >}}