diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 00000000..5bf4860b --- /dev/null +++ b/.editorconfig @@ -0,0 +1,26 @@ +root = true + +[*] +charset = utf-8 +insert_final_newline = true +trim_trailing_whitespace = true + +[*.md] +indent_size = 2 +indent_style = space +max_line_length = 100 # Please keep this in sync with bin/lesson_check.py! +trim_trailing_whitespace = false # keep trailing spaces in markdown - 2+ spaces are translated to a hard break (
) + +[*.r] +max_line_length = 80 + +[*.py] +indent_size = 4 +indent_style = space +max_line_length = 79 + +[*.sh] +end_of_line = lf + +[Makefile] +indent_style = tab diff --git a/.github/workflows/README.md b/.github/workflows/README.md new file mode 100755 index 00000000..d6edf88d --- /dev/null +++ b/.github/workflows/README.md @@ -0,0 +1,198 @@ +# Carpentries Workflows + +This directory contains workflows to be used for Lessons using the {sandpaper} +lesson infrastructure. Two of these workflows require R (`sandpaper-main.yaml` +and `pr-recieve.yaml`) and the rest are bots to handle pull request management. + +These workflows will likely change as {sandpaper} evolves, so it is important to +keep them up-to-date. To do this in your lesson you can do the following in your +R console: + +```r +# Install/Update sandpaper +options(repos = c(carpentries = "https://carpentries.r-universe.dev/", + CRAN = "https://cloud.r-project.org")) +install.packages("sandpaper") + +# update the workflows in your lesson +library("sandpaper") +update_github_workflows() +``` + +Inside this folder, you will find a file called `sandpaper-version.txt`, which +will contain a version number for sandpaper. This will be used in the future to +alert you if a workflow update is needed. + +What follows are the descriptions of the workflow files: + +## Deployment + +### 01 Build and Deploy (sandpaper-main.yaml) + +This is the main driver that will only act on the main branch of the repository. +This workflow does the following: + + 1. checks out the lesson + 2. provisions the following resources + - R + - pandoc + - lesson infrastructure (stored in a cache) + - lesson dependencies if needed (stored in a cache) + 3. builds the lesson via `sandpaper:::ci_deploy()` + +#### Caching + +This workflow has two caches; one cache is for the lesson infrastructure and +the other is for the the lesson dependencies if the lesson contains rendered +content. These caches are invalidated by new versions of the infrastructure and +the `renv.lock` file, respectively. If there is a problem with the cache, +manual invaliation is necessary. You will need maintain access to the repository +and you can either go to the actions tab and [click on the caches button to find +and invalidate the failing cache](https://github.blog/changelog/2022-10-20-manage-caches-in-your-actions-workflows-from-web-interface/) +or by setting the `CACHE_VERSION` secret to the current date (which will +invalidate all of the caches). + +## Updates + +### Setup Information + +These workflows run on a schedule and at the maintainer's request. Because they +create pull requests that update workflows/require the downstream actions to run, +they need a special repository/organization secret token called +`SANDPAPER_WORKFLOW` and it must have the `public_repo` and `workflow` scope. + +This can be an individual user token, OR it can be a trusted bot account. If you +have a repository in one of the official Carpentries accounts, then you do not +need to worry about this token being present because the Carpentries Core Team +will take care of supplying this token. + +If you want to use your personal account: you can go to + +to create a token. Once you have created your token, you should copy it to your +clipboard and then go to your repository's settings > secrets > actions and +create or edit the `SANDPAPER_WORKFLOW` secret, pasting in the generated token. + +If you do not specify your token correctly, the runs will not fail and they will +give you instructions to provide the token for your repository. + +### 02 Maintain: Update Workflow Files (update-workflow.yaml) + +The {sandpaper} repository was designed to do as much as possible to separate +the tools from the content. For local builds, this is absolutely true, but +there is a minor issue when it comes to workflow files: they must live inside +the repository. + +This workflow ensures that the workflow files are up-to-date. The way it work is +to download the update-workflows.sh script from GitHub and run it. The script +will do the following: + +1. check the recorded version of sandpaper against the current version on github +2. update the files if there is a difference in versions + +After the files are updated, if there are any changes, they are pushed to a +branch called `update/workflows` and a pull request is created. Maintainers are +encouraged to review the changes and accept the pull request if the outputs +are okay. + +This update is run ~~weekly or~~ on demand. + +### 03 Maintain: Update Package Cache (update-cache.yaml) + +For lessons that have generated content, we use {renv} to ensure that the output +is stable. This is controlled by a single lockfile which documents the packages +needed for the lesson and the version numbers. This workflow is skipped in +lessons that do not have generated content. + +Because the lessons need to remain current with the package ecosystem, it's a +good idea to make sure these packages can be updated periodically. The +update cache workflow will do this by checking for updates, applying them in a +branch called `updates/packages` and creating a pull request with _only the +lockfile changed_. + +From here, the markdown documents will be rebuilt and you can inspect what has +changed based on how the packages have updated. + +## Pull Request and Review Management + +Because our lessons execute code, pull requests are a secruity risk for any +lesson and thus have security measures associted with them. **Do not merge any +pull requests that do not pass checks and do not have bots commented on them.** + +This series of workflows all go together and are described in the following +diagram and the below sections: + +![Graph representation of a pull request](https://carpentries.github.io/sandpaper/articles/img/pr-flow.dot.svg) + +### Pre Flight Pull Request Validation (pr-preflight.yaml) + +This workflow runs every time a pull request is created and its purpose is to +validate that the pull request is okay to run. This means the following things: + +1. The pull request does not contain modified workflow files +2. If the pull request contains modified workflow files, it does not contain + modified content files (such as a situation where @carpentries-bot will + make an automated pull request) +3. The pull request does not contain an invalid commit hash (e.g. from a fork + that was made before a lesson was transitioned from styles to use the + workbench). + +Once the checks are finished, a comment is issued to the pull request, which +will allow maintainers to determine if it is safe to run the +"Receive Pull Request" workflow from new contributors. + +### Recieve Pull Request (pr-recieve.yaml) + +**Note of caution:** This workflow runs arbitrary code by anyone who creates a +pull request. GitHub has safeguarded the token used in this workflow to have no +priviledges in the repository, but we have taken precautions to protect against +spoofing. + +This workflow is triggered with every push to a pull request. If this workflow +is already running and a new push is sent to the pull request, the workflow +running from the previous push will be cancelled and a new workflow run will be +started. + +The first step of this workflow is to check if it is valid (e.g. that no +workflow files have been modified). If there are workflow files that have been +modified, a comment is made that indicates that the workflow is not run. If +both a workflow file and lesson content is modified, an error will occurr. + +The second step (if valid) is to build the generated content from the pull +request. This builds the content and uploads three artifacts: + +1. The pull request number (pr) +2. A summary of changes after the rendering process (diff) +3. The rendered files (build) + +Because this workflow builds generated content, it follows the same general +process as the `sandpaper-main` workflow with the same caching mechanisms. + +The artifacts produced are used by the next workflow. + +### Comment on Pull Request (pr-comment.yaml) + +This workflow is triggered if the `pr-recieve.yaml` workflow is successful. +The steps in this workflow are: + +1. Test if the workflow is valid and comment the validity of the workflow to the + pull request. +2. If it is valid: create an orphan branch with two commits: the current state + of the repository and the proposed changes. +3. If it is valid: update the pull request comment with the summary of changes + +Importantly: if the pull request is invalid, the branch is not created so any +malicious code is not published. + +From here, the maintainer can request changes from the author and eventually +either merge or reject the PR. When this happens, if the PR was valid, the +preview branch needs to be deleted. + +### Send Close PR Signal (pr-close-signal.yaml) + +Triggered any time a pull request is closed. This emits an artifact that is the +pull request number for the next action + +### Remove Pull Request Branch (pr-post-remove-branch.yaml) + +Tiggered by `pr-close-signal.yaml`. This removes the temporary branch associated with +the pull request (if it was created). diff --git a/.github/workflows/pr-close-signal.yaml b/.github/workflows/pr-close-signal.yaml new file mode 100755 index 00000000..9b129d5d --- /dev/null +++ b/.github/workflows/pr-close-signal.yaml @@ -0,0 +1,23 @@ +name: "Bot: Send Close Pull Request Signal" + +on: + pull_request: + types: + [closed] + +jobs: + send-close-signal: + name: "Send closing signal" + runs-on: ubuntu-latest + if: ${{ github.event.action == 'closed' }} + steps: + - name: "Create PRtifact" + run: | + mkdir -p ./pr + printf ${{ github.event.number }} > ./pr/NUM + - name: Upload Diff + uses: actions/upload-artifact@v3 + with: + name: pr + path: ./pr + diff --git a/.github/workflows/pr-comment.yaml b/.github/workflows/pr-comment.yaml new file mode 100755 index 00000000..bb2eb03c --- /dev/null +++ b/.github/workflows/pr-comment.yaml @@ -0,0 +1,185 @@ +name: "Bot: Comment on the Pull Request" + +# read-write repo token +# access to secrets +on: + workflow_run: + workflows: ["Receive Pull Request"] + types: + - completed + +concurrency: + group: pr-${{ github.event.workflow_run.pull_requests[0].number }} + cancel-in-progress: true + + +jobs: + # Pull requests are valid if: + # - they match the sha of the workflow run head commit + # - they are open + # - no .github files were committed + test-pr: + name: "Test if pull request is valid" + runs-on: ubuntu-latest + if: > + github.event.workflow_run.event == 'pull_request' && + github.event.workflow_run.conclusion == 'success' + outputs: + is_valid: ${{ steps.check-pr.outputs.VALID }} + payload: ${{ steps.check-pr.outputs.payload }} + number: ${{ steps.get-pr.outputs.NUM }} + msg: ${{ steps.check-pr.outputs.MSG }} + steps: + - name: 'Download PR artifact' + id: dl + uses: carpentries/actions/download-workflow-artifact@main + with: + run: ${{ github.event.workflow_run.id }} + name: 'pr' + + - name: "Get PR Number" + if: ${{ steps.dl.outputs.success == 'true' }} + id: get-pr + run: | + unzip pr.zip + echo "NUM=$(<./NR)" >> $GITHUB_OUTPUT + + - name: "Fail if PR number was not present" + id: bad-pr + if: ${{ steps.dl.outputs.success != 'true' }} + run: | + echo '::error::A pull request number was not recorded. The pull request that triggered this workflow is likely malicious.' + exit 1 + - name: "Get Invalid Hashes File" + id: hash + run: | + echo "json<> $GITHUB_OUTPUT + - name: "Check PR" + id: check-pr + if: ${{ steps.dl.outputs.success == 'true' }} + uses: carpentries/actions/check-valid-pr@main + with: + pr: ${{ steps.get-pr.outputs.NUM }} + sha: ${{ github.event.workflow_run.head_sha }} + headroom: 3 # if it's within the last three commits, we can keep going, because it's likely rapid-fire + invalid: ${{ fromJSON(steps.hash.outputs.json)[github.repository] }} + fail_on_error: true + + # Create an orphan branch on this repository with two commits + # - the current HEAD of the md-outputs branch + # - the output from running the current HEAD of the pull request through + # the md generator + create-branch: + name: "Create Git Branch" + needs: test-pr + runs-on: ubuntu-latest + if: ${{ needs.test-pr.outputs.is_valid == 'true' }} + env: + NR: ${{ needs.test-pr.outputs.number }} + permissions: + contents: write + steps: + - name: 'Checkout md outputs' + uses: actions/checkout@v3 + with: + ref: md-outputs + path: built + fetch-depth: 1 + + - name: 'Download built markdown' + id: dl + uses: carpentries/actions/download-workflow-artifact@main + with: + run: ${{ github.event.workflow_run.id }} + name: 'built' + + - if: ${{ steps.dl.outputs.success == 'true' }} + run: unzip built.zip + + - name: "Create orphan and push" + if: ${{ steps.dl.outputs.success == 'true' }} + run: | + cd built/ + git config --local user.email "actions@github.com" + git config --local user.name "GitHub Actions" + CURR_HEAD=$(git rev-parse HEAD) + git checkout --orphan md-outputs-PR-${NR} + git add -A + git commit -m "source commit: ${CURR_HEAD}" + ls -A | grep -v '^.git$' | xargs -I _ rm -r '_' + cd .. + unzip -o -d built built.zip + cd built + git add -A + git commit --allow-empty -m "differences for PR #${NR}" + git push -u --force --set-upstream origin md-outputs-PR-${NR} + + # Comment on the Pull Request with a link to the branch and the diff + comment-pr: + name: "Comment on Pull Request" + needs: [test-pr, create-branch] + runs-on: ubuntu-latest + if: ${{ needs.test-pr.outputs.is_valid == 'true' }} + env: + NR: ${{ needs.test-pr.outputs.number }} + permissions: + pull-requests: write + steps: + - name: 'Download comment artifact' + id: dl + uses: carpentries/actions/download-workflow-artifact@main + with: + run: ${{ github.event.workflow_run.id }} + name: 'diff' + + - if: ${{ steps.dl.outputs.success == 'true' }} + run: unzip ${{ github.workspace }}/diff.zip + + - name: "Comment on PR" + id: comment-diff + if: ${{ steps.dl.outputs.success == 'true' }} + uses: carpentries/actions/comment-diff@main + with: + pr: ${{ env.NR }} + path: ${{ github.workspace }}/diff.md + + # Comment if the PR is open and matches the SHA, but the workflow files have + # changed + comment-changed-workflow: + name: "Comment if workflow files have changed" + needs: test-pr + runs-on: ubuntu-latest + if: ${{ always() && needs.test-pr.outputs.is_valid == 'false' }} + env: + NR: ${{ github.event.workflow_run.pull_requests[0].number }} + body: ${{ needs.test-pr.outputs.msg }} + permissions: + pull-requests: write + steps: + - name: 'Check for spoofing' + id: dl + uses: carpentries/actions/download-workflow-artifact@main + with: + run: ${{ github.event.workflow_run.id }} + name: 'built' + + - name: 'Alert if spoofed' + id: spoof + if: ${{ steps.dl.outputs.success == 'true' }} + run: | + echo 'body<> $GITHUB_ENV + echo '' >> $GITHUB_ENV + echo '## :x: DANGER :x:' >> $GITHUB_ENV + echo 'This pull request has modified workflows that created output. Close this now.' >> $GITHUB_ENV + echo '' >> $GITHUB_ENV + echo 'EOF' >> $GITHUB_ENV + + - name: "Comment on PR" + id: comment-diff + uses: carpentries/actions/comment-diff@main + with: + pr: ${{ env.NR }} + body: ${{ env.body }} + diff --git a/.github/workflows/pr-post-remove-branch.yaml b/.github/workflows/pr-post-remove-branch.yaml new file mode 100755 index 00000000..62c2e98d --- /dev/null +++ b/.github/workflows/pr-post-remove-branch.yaml @@ -0,0 +1,32 @@ +name: "Bot: Remove Temporary PR Branch" + +on: + workflow_run: + workflows: ["Bot: Send Close Pull Request Signal"] + types: + - completed + +jobs: + delete: + name: "Delete branch from Pull Request" + runs-on: ubuntu-latest + if: > + github.event.workflow_run.event == 'pull_request' && + github.event.workflow_run.conclusion == 'success' + permissions: + contents: write + steps: + - name: 'Download artifact' + uses: carpentries/actions/download-workflow-artifact@main + with: + run: ${{ github.event.workflow_run.id }} + name: pr + - name: "Get PR Number" + id: get-pr + run: | + unzip pr.zip + echo "NUM=$(<./NUM)" >> $GITHUB_OUTPUT + - name: 'Remove branch' + uses: carpentries/actions/remove-branch@main + with: + pr: ${{ steps.get-pr.outputs.NUM }} diff --git a/.github/workflows/pr-preflight.yaml b/.github/workflows/pr-preflight.yaml new file mode 100755 index 00000000..d0d7420d --- /dev/null +++ b/.github/workflows/pr-preflight.yaml @@ -0,0 +1,39 @@ +name: "Pull Request Preflight Check" + +on: + pull_request_target: + branches: + ["main"] + types: + ["opened", "synchronize", "reopened"] + +jobs: + test-pr: + name: "Test if pull request is valid" + if: ${{ github.event.action != 'closed' }} + runs-on: ubuntu-latest + outputs: + is_valid: ${{ steps.check-pr.outputs.VALID }} + permissions: + pull-requests: write + steps: + - name: "Get Invalid Hashes File" + id: hash + run: | + echo "json<> $GITHUB_OUTPUT + - name: "Check PR" + id: check-pr + uses: carpentries/actions/check-valid-pr@main + with: + pr: ${{ github.event.number }} + invalid: ${{ fromJSON(steps.hash.outputs.json)[github.repository] }} + fail_on_error: true + - name: "Comment result of validation" + id: comment-diff + if: ${{ always() }} + uses: carpentries/actions/comment-diff@main + with: + pr: ${{ github.event.number }} + body: ${{ steps.check-pr.outputs.MSG }} diff --git a/.github/workflows/pr-receive.yaml b/.github/workflows/pr-receive.yaml new file mode 100755 index 00000000..371ef542 --- /dev/null +++ b/.github/workflows/pr-receive.yaml @@ -0,0 +1,131 @@ +name: "Receive Pull Request" + +on: + pull_request: + types: + [opened, synchronize, reopened] + +concurrency: + group: ${{ github.ref }} + cancel-in-progress: true + +jobs: + test-pr: + name: "Record PR number" + if: ${{ github.event.action != 'closed' }} + runs-on: ubuntu-latest + outputs: + is_valid: ${{ steps.check-pr.outputs.VALID }} + steps: + - name: "Record PR number" + id: record + if: ${{ always() }} + run: | + echo ${{ github.event.number }} > ${{ github.workspace }}/NR # 2022-03-02: artifact name fixed to be NR + - name: "Upload PR number" + id: upload + if: ${{ always() }} + uses: actions/upload-artifact@v3 + with: + name: pr + path: ${{ github.workspace }}/NR + - name: "Get Invalid Hashes File" + id: hash + run: | + echo "json<> $GITHUB_OUTPUT + - name: "echo output" + run: | + echo "${{ steps.hash.outputs.json }}" + - name: "Check PR" + id: check-pr + uses: carpentries/actions/check-valid-pr@main + with: + pr: ${{ github.event.number }} + invalid: ${{ fromJSON(steps.hash.outputs.json)[github.repository] }} + + build-md-source: + name: "Build markdown source files if valid" + needs: test-pr + runs-on: ubuntu-latest + if: ${{ needs.test-pr.outputs.is_valid == 'true' }} + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + RENV_PATHS_ROOT: ~/.local/share/renv/ + CHIVE: ${{ github.workspace }}/site/chive + PR: ${{ github.workspace }}/site/pr + MD: ${{ github.workspace }}/site/built + steps: + - name: "Check Out Main Branch" + uses: actions/checkout@v3 + + - name: "Check Out Staging Branch" + uses: actions/checkout@v3 + with: + ref: md-outputs + path: ${{ env.MD }} + + - name: "Set up R" + uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + install-r: false + + - name: "Set up Pandoc" + uses: r-lib/actions/setup-pandoc@v2 + + - name: "Setup Lesson Engine" + uses: carpentries/actions/setup-sandpaper@main + with: + cache-version: ${{ secrets.CACHE_VERSION }} + + - name: "Setup Package Cache" + uses: carpentries/actions/setup-lesson-deps@main + with: + cache-version: ${{ secrets.CACHE_VERSION }} + + - name: "Validate and Build Markdown" + id: build-site + run: | + sandpaper::package_cache_trigger(TRUE) + sandpaper::validate_lesson(path = '${{ github.workspace }}') + sandpaper:::build_markdown(path = '${{ github.workspace }}', quiet = FALSE) + shell: Rscript {0} + + - name: "Generate Artifacts" + id: generate-artifacts + run: | + sandpaper:::ci_bundle_pr_artifacts( + repo = '${{ github.repository }}', + pr_number = '${{ github.event.number }}', + path_md = '${{ env.MD }}', + path_pr = '${{ env.PR }}', + path_archive = '${{ env.CHIVE }}', + branch = 'md-outputs' + ) + shell: Rscript {0} + + - name: "Upload PR" + uses: actions/upload-artifact@v3 + with: + name: pr + path: ${{ env.PR }} + + - name: "Upload Diff" + uses: actions/upload-artifact@v3 + with: + name: diff + path: ${{ env.CHIVE }} + retention-days: 1 + + - name: "Upload Build" + uses: actions/upload-artifact@v3 + with: + name: built + path: ${{ env.MD }} + retention-days: 1 + + - name: "Teardown" + run: sandpaper::reset_site() + shell: Rscript {0} diff --git a/.github/workflows/sandpaper-main.yaml b/.github/workflows/sandpaper-main.yaml new file mode 100755 index 00000000..e17707ac --- /dev/null +++ b/.github/workflows/sandpaper-main.yaml @@ -0,0 +1,61 @@ +name: "01 Build and Deploy Site" + +on: + push: + branches: + - main + - master + schedule: + - cron: '0 0 * * 2' + workflow_dispatch: + inputs: + name: + description: 'Who triggered this build?' + required: true + default: 'Maintainer (via GitHub)' + reset: + description: 'Reset cached markdown files' + required: false + default: false + type: boolean +jobs: + full-build: + name: "Build Full Site" + runs-on: ubuntu-latest + permissions: + checks: write + contents: write + pages: write + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + RENV_PATHS_ROOT: ~/.local/share/renv/ + steps: + + - name: "Checkout Lesson" + uses: actions/checkout@v3 + + - name: "Set up R" + uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + install-r: false + + - name: "Set up Pandoc" + uses: r-lib/actions/setup-pandoc@v2 + + - name: "Setup Lesson Engine" + uses: carpentries/actions/setup-sandpaper@main + with: + cache-version: ${{ secrets.CACHE_VERSION }} + + - name: "Setup Package Cache" + uses: carpentries/actions/setup-lesson-deps@main + with: + cache-version: ${{ secrets.CACHE_VERSION }} + + - name: "Deploy Site" + run: | + reset <- "${{ github.event.inputs.reset }}" == "true" + sandpaper::package_cache_trigger(TRUE) + sandpaper:::ci_deploy(reset = reset) + shell: Rscript {0} diff --git a/.github/workflows/sandpaper-version.txt b/.github/workflows/sandpaper-version.txt new file mode 100644 index 00000000..201a22c8 --- /dev/null +++ b/.github/workflows/sandpaper-version.txt @@ -0,0 +1 @@ +0.16.2 diff --git a/.github/workflows/update-cache.yaml b/.github/workflows/update-cache.yaml new file mode 100755 index 00000000..676d7424 --- /dev/null +++ b/.github/workflows/update-cache.yaml @@ -0,0 +1,125 @@ +name: "03 Maintain: Update Package Cache" + +on: + workflow_dispatch: + inputs: + name: + description: 'Who triggered this build (enter github username to tag yourself)?' + required: true + default: 'monthly run' + schedule: + # Run every tuesday + - cron: '0 0 * * 2' + +jobs: + preflight: + name: "Preflight Check" + runs-on: ubuntu-latest + outputs: + ok: ${{ steps.check.outputs.ok }} + steps: + - id: check + run: | + if [[ ${{ github.event_name }} == 'workflow_dispatch' ]]; then + echo "ok=true" >> $GITHUB_OUTPUT + echo "Running on request" + # using single brackets here to avoid 08 being interpreted as octal + # https://github.com/carpentries/sandpaper/issues/250 + elif [ `date +%d` -le 7 ]; then + # If the Tuesday lands in the first week of the month, run it + echo "ok=true" >> $GITHUB_OUTPUT + echo "Running on schedule" + else + echo "ok=false" >> $GITHUB_OUTPUT + echo "Not Running Today" + fi + + check_renv: + name: "Check if We Need {renv}" + runs-on: ubuntu-latest + needs: preflight + if: ${{ needs.preflight.outputs.ok == 'true'}} + outputs: + needed: ${{ steps.renv.outputs.exists }} + steps: + - name: "Checkout Lesson" + uses: actions/checkout@v3 + - id: renv + run: | + if [[ -d renv ]]; then + echo "exists=true" >> $GITHUB_OUTPUT + fi + + check_token: + name: "Check SANDPAPER_WORKFLOW token" + runs-on: ubuntu-latest + needs: check_renv + if: ${{ needs.check_renv.outputs.needed == 'true' }} + outputs: + workflow: ${{ steps.validate.outputs.wf }} + repo: ${{ steps.validate.outputs.repo }} + steps: + - name: "validate token" + id: validate + uses: carpentries/actions/check-valid-credentials@main + with: + token: ${{ secrets.SANDPAPER_WORKFLOW }} + + update_cache: + name: "Update Package Cache" + needs: check_token + if: ${{ needs.check_token.outputs.repo== 'true' }} + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + RENV_PATHS_ROOT: ~/.local/share/renv/ + steps: + + - name: "Checkout Lesson" + uses: actions/checkout@v3 + + - name: "Set up R" + uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + install-r: false + + - name: "Update {renv} deps and determine if a PR is needed" + id: update + uses: carpentries/actions/update-lockfile@main + with: + cache-version: ${{ secrets.CACHE_VERSION }} + + - name: Create Pull Request + id: cpr + if: ${{ steps.update.outputs.n > 0 }} + uses: carpentries/create-pull-request@main + with: + token: ${{ secrets.SANDPAPER_WORKFLOW }} + delete-branch: true + branch: "update/packages" + commit-message: "[actions] update ${{ steps.update.outputs.n }} packages" + title: "Update ${{ steps.update.outputs.n }} packages" + body: | + :robot: This is an automated build + + This will update ${{ steps.update.outputs.n }} packages in your lesson with the following versions: + + ``` + ${{ steps.update.outputs.report }} + ``` + + :stopwatch: In a few minutes, a comment will appear that will show you how the output has changed based on these updates. + + If you want to inspect these changes locally, you can use the following code to check out a new branch: + + ```bash + git fetch origin update/packages + git checkout update/packages + ``` + + - Auto-generated by [create-pull-request][1] on ${{ steps.update.outputs.date }} + + [1]: https://github.com/carpentries/create-pull-request/tree/main + labels: "type: package cache" + draft: false diff --git a/.github/workflows/update-workflows.yaml b/.github/workflows/update-workflows.yaml new file mode 100755 index 00000000..288bcd13 --- /dev/null +++ b/.github/workflows/update-workflows.yaml @@ -0,0 +1,66 @@ +name: "02 Maintain: Update Workflow Files" + +on: + workflow_dispatch: + inputs: + name: + description: 'Who triggered this build (enter github username to tag yourself)?' + required: true + default: 'weekly run' + clean: + description: 'Workflow files/file extensions to clean (no wildcards, enter "" for none)' + required: false + default: '.yaml' + schedule: + # Run every Tuesday + - cron: '0 0 * * 2' + +jobs: + check_token: + name: "Check SANDPAPER_WORKFLOW token" + runs-on: ubuntu-latest + outputs: + workflow: ${{ steps.validate.outputs.wf }} + repo: ${{ steps.validate.outputs.repo }} + steps: + - name: "validate token" + id: validate + uses: carpentries/actions/check-valid-credentials@main + with: + token: ${{ secrets.SANDPAPER_WORKFLOW }} + + update_workflow: + name: "Update Workflow" + runs-on: ubuntu-latest + needs: check_token + if: ${{ needs.check_token.outputs.workflow == 'true' }} + steps: + - name: "Checkout Repository" + uses: actions/checkout@v3 + + - name: Update Workflows + id: update + uses: carpentries/actions/update-workflows@main + with: + clean: ${{ github.event.inputs.clean }} + + - name: Create Pull Request + id: cpr + if: "${{ steps.update.outputs.new }}" + uses: carpentries/create-pull-request@main + with: + token: ${{ secrets.SANDPAPER_WORKFLOW }} + delete-branch: true + branch: "update/workflows" + commit-message: "[actions] update sandpaper workflow to version ${{ steps.update.outputs.new }}" + title: "Update Workflows to Version ${{ steps.update.outputs.new }}" + body: | + :robot: This is an automated build + + Update Workflows from sandpaper version ${{ steps.update.outputs.old }} -> ${{ steps.update.outputs.new }} + + - Auto-generated by [create-pull-request][1] on ${{ steps.update.outputs.date }} + + [1]: https://github.com/carpentries/create-pull-request/tree/main + labels: "type: template and tools" + draft: false diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..b8ab7062 --- /dev/null +++ b/.gitignore @@ -0,0 +1,55 @@ +# sandpaper files +episodes/*html +site/* +!site/README.md + +# History files +.Rhistory +.Rapp.history +# Session Data files +.RData +# User-specific files +.Ruserdata +# Example code in package build process +*-Ex.R +# Output files from R CMD build +/*.tar.gz +# Output files from R CMD check +/*.Rcheck/ +# RStudio files +.Rproj.user/ +# produced vignettes +vignettes/*.html +vignettes/*.pdf +# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3 +.httr-oauth +# knitr and R markdown default cache directories +*_cache/ +/cache/ +# Temporary files created by R markdown +*.utf8.md +*.knit.md +# R Environment Variables +.Renviron +# pkgdown site +docs/ +# translation temp files +po/*~ +# renv detritus +renv/sandbox/ +*.pyc +*~ +.DS_Store +.ipynb_checkpoints +.sass-cache +.jekyll-cache/ +.jekyll-metadata +__pycache__ +_site +.Rproj.user +.bundle/ +.vendor/ +vendor/ +.docker-vendor/ +Gemfile.lock +.*history diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md new file mode 100644 index 00000000..f19b8049 --- /dev/null +++ b/CODE_OF_CONDUCT.md @@ -0,0 +1,13 @@ +--- +title: "Contributor Code of Conduct" +--- + +As contributors and maintainers of this project, +we pledge to follow the [The Carpentries Code of Conduct][coc]. + +Instances of abusive, harassing, or otherwise unacceptable behavior +may be reported by following our [reporting guidelines][coc-reporting]. + + +[coc-reporting]: https://docs.carpentries.org/topic_folders/policies/incident-reporting.html +[coc]: https://docs.carpentries.org/topic_folders/policies/code-of-conduct.html diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 00000000..6c2b81c8 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,123 @@ +## Contributing + +[The Carpentries][cp-site] ([Software Carpentry][swc-site], [Data +Carpentry][dc-site], and [Library Carpentry][lc-site]) are open source +projects, and we welcome contributions of all kinds: new lessons, fixes to +existing material, bug reports, and reviews of proposed changes are all +welcome. + +### Contributor Agreement + +By contributing, you agree that we may redistribute your work under [our +license](LICENSE.md). In exchange, we will address your issues and/or assess +your change proposal as promptly as we can, and help you become a member of our +community. Everyone involved in [The Carpentries][cp-site] agrees to abide by +our [code of conduct](CODE_OF_CONDUCT.md). + +### How to Contribute + +The easiest way to get started is to file an issue to tell us about a spelling +mistake, some awkward wording, or a factual error. This is a good way to +introduce yourself and to meet some of our community members. + +1. If you do not have a [GitHub][github] account, you can [send us comments by + email][contact]. However, we will be able to respond more quickly if you use + one of the other methods described below. + +2. If you have a [GitHub][github] account, or are willing to [create + one][github-join], but do not know how to use Git, you can report problems + or suggest improvements by [creating an issue][repo-issues]. This allows us + to assign the item to someone and to respond to it in a threaded discussion. + +3. If you are comfortable with Git, and would like to add or change material, + you can submit a pull request (PR). Instructions for doing this are + [included below](#using-github). For inspiration about changes that need to + be made, check out the [list of open issues][issues] across the Carpentries. + +Note: if you want to build the website locally, please refer to [The Workbench +documentation][template-doc]. + +### Where to Contribute + +1. If you wish to change this lesson, add issues and pull requests here. +2. If you wish to change the template used for workshop websites, please refer + to [The Workbench documentation][template-doc]. + + +### What to Contribute + +There are many ways to contribute, from writing new exercises and improving +existing ones to updating or filling in the documentation and submitting [bug +reports][issues] about things that do not work, are not clear, or are missing. +If you are looking for ideas, please see [the list of issues for this +repository][repo-issues], or the issues for [Data Carpentry][dc-issues], +[Library Carpentry][lc-issues], and [Software Carpentry][swc-issues] projects. + +Comments on issues and reviews of pull requests are just as welcome: we are +smarter together than we are on our own. **Reviews from novices and newcomers +are particularly valuable**: it's easy for people who have been using these +lessons for a while to forget how impenetrable some of this material can be, so +fresh eyes are always welcome. + +### What *Not* to Contribute + +Our lessons already contain more material than we can cover in a typical +workshop, so we are usually *not* looking for more concepts or tools to add to +them. As a rule, if you want to introduce a new idea, you must (a) estimate how +long it will take to teach and (b) explain what you would take out to make room +for it. The first encourages contributors to be honest about requirements; the +second, to think hard about priorities. + +We are also not looking for exercises or other material that only run on one +platform. Our workshops typically contain a mixture of Windows, macOS, and +Linux users; in order to be usable, our lessons must run equally well on all +three. + +### Using GitHub + +If you choose to contribute via GitHub, you may want to look at [How to +Contribute to an Open Source Project on GitHub][how-contribute]. In brief, we +use [GitHub flow][github-flow] to manage changes: + +1. Create a new branch in your desktop copy of this repository for each + significant change. +2. Commit the change in that branch. +3. Push that branch to your fork of this repository on GitHub. +4. Submit a pull request from that branch to the [upstream repository][repo]. +5. If you receive feedback, make changes on your desktop and push to your + branch on GitHub: the pull request will update automatically. + +NB: The published copy of the lesson is usually in the `main` branch. + +Each lesson has a team of maintainers who review issues and pull requests or +encourage others to do so. The maintainers are community volunteers, and have +final say over what gets merged into the lesson. + +### Other Resources + +The Carpentries is a global organisation with volunteers and learners all over +the world. We share values of inclusivity and a passion for sharing knowledge, +teaching and learning. There are several ways to connect with The Carpentries +community listed at including via social +media, slack, newsletters, and email lists. You can also [reach us by +email][contact]. + +[repo]: https://example.com/FIXME +[repo-issues]: https://example.com/FIXME/issues +[contact]: mailto:team@carpentries.org +[cp-site]: https://carpentries.org/ +[dc-issues]: https://github.com/issues?q=user%3Adatacarpentry +[dc-lessons]: https://datacarpentry.org/lessons/ +[dc-site]: https://datacarpentry.org/ +[discuss-list]: https://carpentries.topicbox.com/groups/discuss +[github]: https://github.com +[github-flow]: https://guides.github.com/introduction/flow/ +[github-join]: https://github.com/join +[how-contribute]: https://egghead.io/courses/how-to-contribute-to-an-open-source-project-on-github +[issues]: https://carpentries.org/help-wanted-issues/ +[lc-issues]: https://github.com/issues?q=user%3ALibraryCarpentry +[swc-issues]: https://github.com/issues?q=user%3Aswcarpentry +[swc-lessons]: https://software-carpentry.org/lessons/ +[swc-site]: https://software-carpentry.org/ +[lc-site]: https://librarycarpentry.org/ +[template-doc]: https://carpentries.github.io/workbench/ diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 00000000..7632871f --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,79 @@ +--- +title: "Licenses" +--- + +## Instructional Material + +All Carpentries (Software Carpentry, Data Carpentry, and Library Carpentry) +instructional material is made available under the [Creative Commons +Attribution license][cc-by-human]. The following is a human-readable summary of +(and not a substitute for) the [full legal text of the CC BY 4.0 +license][cc-by-legal]. + +You are free: + +- to **Share**---copy and redistribute the material in any medium or format +- to **Adapt**---remix, transform, and build upon the material + +for any purpose, even commercially. + +The licensor cannot revoke these freedoms as long as you follow the license +terms. + +Under the following terms: + +- **Attribution**---You must give appropriate credit (mentioning that your work + is derived from work that is Copyright (c) The Carpentries and, where + practical, linking to ), provide a [link to the + license][cc-by-human], and indicate if changes were made. You may do so in + any reasonable manner, but not in any way that suggests the licensor endorses + you or your use. + +- **No additional restrictions**---You may not apply legal terms or + technological measures that legally restrict others from doing anything the + license permits. With the understanding that: + +Notices: + +* You do not have to comply with the license for elements of the material in + the public domain or where your use is permitted by an applicable exception + or limitation. +* No warranties are given. The license may not give you all of the permissions + necessary for your intended use. For example, other rights such as publicity, + privacy, or moral rights may limit how you use the material. + +## Software + +Except where otherwise noted, the example programs and other software provided +by The Carpentries are made available under the [OSI][osi]-approved [MIT +license][mit-license]. + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies +of the Software, and to permit persons to whom the Software is furnished to do +so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + +## Trademark + +"The Carpentries", "Software Carpentry", "Data Carpentry", and "Library +Carpentry" and their respective logos are registered trademarks of [Community +Initiatives][ci]. + +[cc-by-human]: https://creativecommons.org/licenses/by/4.0/ +[cc-by-legal]: https://creativecommons.org/licenses/by/4.0/legalcode +[mit-license]: https://opensource.org/licenses/mit-license.html +[ci]: https://communityin.org/ +[osi]: https://opensource.org diff --git a/README.md b/README.md index 59a90385..9956230d 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,6 @@ This is one sub-module within the [Neuroimaging cirriculum][neuro_cirriculum]. V **fMRI Analysis in Python** is a programme developed to facilitate reproducibility in functional neuroimaging analyses. Python is emerging as a standard language of data analysis, visualization, and workflow building. More recently, it has rapidly been adopted by the neuroimaging community as a means of developing powerful open-source tools in favour of historically used opaque software such as AFNI, FSL and SPM. In addition, the barrier to entry to Python is low - meaning that you as the user can easily develop your own packages and contribute to the open-source codebase of neuroimaging! - *** The **fMRI Analysis in Python** is a workshop series started up via a collaboration between researchers and staff at the Centre for Addiction and Mental Health (Toronto, ON), the University of Western Ontario (London, Ontario), and McGill University (Montreal, Quebec). @@ -23,18 +22,17 @@ This lesson covers fMRI imaging analysis from the basic steps of preprocessing a ### Episodes -| Time | Episode | Question(s) Answered | -| --- | --- | --- | -||Setup|Download files required for the lesson| -| 00:00 | 1. Course Overview and Introduction | What steps do I need to take before beginning to work with fMRI data? | -| 00:25 | 2. Exploring Preprocessed fMRI Data from fMRIPREP | How does fMRIPrep store preprocessed neuroimaging data? How do I access preprocessed neuroimaging data? | -| 00:50 | 3. Introduction to Image Manipulation using Nilearn | How can I perform arithmetic operations on MR images? | -| 01:35 | 4. Integrating Functional Data | How is fMRI data represented? How can I access fMRI data along spatial and temporal dimensions? How can I integrate fMRI and structural MRI together? | -| 02:20 | 6. Cleaning Confounders in your Data with Nilearn | How can I clean the data so that it more closely reflects BOLD instead of artifacts? | -| 02:50 | 7. Applying Parcellations to Resting State Data | How can I reduce amount of noise-related variance in my data? How can I frame my data as a set of meaningful features? | -| 03:30 | 8. Functional Connectivity Analysis | How can we estimate brain functional connectivity patterns from resting state data? | -| 04:15 | Finish | | - +| Time | Episode | Question(s) Answered | +| ----- | --------------------------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------- | +| | Setup | Download files required for the lesson | +| 00:00 | 1\. Course Overview and Introduction | What steps do I need to take before beginning to work with fMRI data? | +| 00:25 | 2\. Exploring Preprocessed fMRI Data from fMRIPREP | How does fMRIPrep store preprocessed neuroimaging data? How do I access preprocessed neuroimaging data? | +| 00:50 | 3\. Introduction to Image Manipulation using Nilearn | How can I perform arithmetic operations on MR images? | +| 01:35 | 4\. Integrating Functional Data | How is fMRI data represented? How can I access fMRI data along spatial and temporal dimensions? How can I integrate fMRI and structural MRI together? | +| 02:20 | 6\. Cleaning Confounders in your Data with Nilearn | How can I clean the data so that it more closely reflects BOLD instead of artifacts? | +| 02:50 | 7\. Applying Parcellations to Resting State Data | How can I reduce amount of noise-related variance in my data? How can I frame my data as a set of meaningful features? | +| 03:30 | 8\. Functional Connectivity Analysis | How can we estimate brain functional connectivity patterns from resting state data? | +| 04:15 | Finish | | ## Contributing @@ -47,18 +45,17 @@ how to write new episodes. Please see the current list of [issues](https://github.com/carpentries-incubator/SDC-BIDS-fMRI/issues) for ideas for contributing to this repository. For making your contribution, we use the GitHub flow, which is -nicely explained in the chapter [Contributing to a Project](http://git-scm.com/book/en/v2/GitHub-Contributing-to-a-Project) in Pro Git +nicely explained in the chapter [Contributing to a Project](https://git-scm.com/book/en/v2/GitHub-Contributing-to-a-Project) in Pro Git by Scott Chacon. -Look for the tag ![good_first_issue](https://img.shields.io/badge/-good%20first%20issue-gold.svg). This indicates that the mantainers will welcome a pull request fixing this issue. - +Look for the tag ![good\_first\_issue](https://img.shields.io/badge/-good%20first%20issue-gold.svg). This indicates that the mantainers will welcome a pull request fixing this issue. ## Maintainer(s) -* [Jerrold Jeyachandra][jerrold_jeyachandra] -* [Michael Joseph][michael_joseph] -* [Olivia Stanley][olivia_stanley] -* [Jason Kai][jason_kai] -* [Erin Dickie][erin_dickie] +- [Jerrold Jeyachandra][jerrold_jeyachandra] +- [Michael Joseph][michael_joseph] +- [Olivia Stanley][olivia_stanley] +- [Jason Kai][jason_kai] +- [Erin Dickie][erin_dickie] ## Authors @@ -68,10 +65,13 @@ A list of contributors to the lesson can be found in [AUTHORS](AUTHORS) To cite this lesson, please consult with [CITATION](CITATION) +[neuro_cirriculum]: https://carpentries.org/community-lessons/#neuroimaging [lesson-example]: https://carpentries.github.io/lesson-example [jerrold_jeyachandra]: https://github.com/jerdra -[olivia_stanley]: https://github.com/ostanley [michael_joseph]: https://github.com/josephmje +[olivia_stanley]: https://github.com/ostanley [jason_kai]: https://github.com/kaitj [erin_dickie]: https://github.com/edickie -[neuro_cirriculum]: https://carpentries.org/community-lessons/#neuroimaging + + + diff --git a/config.yaml b/config.yaml new file mode 100644 index 00000000..5c5fc38a --- /dev/null +++ b/config.yaml @@ -0,0 +1,88 @@ +#------------------------------------------------------------ +# Values for this lesson. +#------------------------------------------------------------ + +# Which carpentry is this (swc, dc, lc, or cp)? +# swc: Software Carpentry +# dc: Data Carpentry +# lc: Library Carpentry +# cp: Carpentries (to use for instructor training for instance) +# incubator: The Carpentries Incubator +carpentry: 'incubator' + +# Overall title for pages. +title: 'Functional Neuroimaging Analysis in Python' + +# Date the lesson was created (YYYY-MM-DD, this is empty by default) +created: + +# Comma-separated list of keywords for the lesson +keywords: 'software, data, lesson, The Carpentries' + +# Life cycle stage of the lesson +# possible values: pre-alpha, alpha, beta, stable +life_cycle: 'alpha' + +# License of the lesson materials (recommended CC-BY 4.0) +license: 'CC-BY 4.0' + +# Link to the source repository for this lesson +source: 'https://github.com/fishtree-attempt/SDC-BIDS-fMRI/' + +# Default branch of your lesson +branch: 'main' + +# Who to contact if there are any issues +contact: 'team@carpentries.org' + +# Navigation ------------------------------------------------ +# +# Use the following menu items to specify the order of +# individual pages in each dropdown section. Leave blank to +# include all pages in the folder. +# +# Example ------------- +# +# episodes: +# - introduction.md +# - first-steps.md +# +# learners: +# - setup.md +# +# instructors: +# - instructor-notes.md +# +# profiles: +# - one-learner.md +# - another-learner.md + +# Order of episodes in your lesson +episodes: +- 01-intro-and-preprocessing.md +- 02-exploring-fmriprep.md +- 03-basic_image_manipulation.md +- 04-integrating_functional_data.md +- 05-data-cleaning-with-nilearn.md +- 06-apply-a-parcellation.md +- 07-functional-connectivity-analysis.md + +# Information for Learners +learners: + +# Information for Instructors +instructors: + +# Learner Profiles +profiles: + +# Customisation --------------------------------------------- +# +# This space below is where custom yaml items (e.g. pinning +# sandpaper and varnish versions) should live + + +url: https://preview.carpentries.org/SDC-BIDS-fMRI +analytics: carpentries +lang: en +workbench-beta: yes diff --git a/episodes/01-intro-and-preprocessing.md b/episodes/01-intro-and-preprocessing.md index b6fc778e..569074b9 100644 --- a/episodes/01-intro-and-preprocessing.md +++ b/episodes/01-intro-and-preprocessing.md @@ -1,16 +1,18 @@ --- -title: "Course Overview and Introduction" +title: Course Overview and Introduction teaching: 25 exercises: 0 -questions: -- "What steps do I need to take before beginning to work with fMRI data?" -- "Learn about fMRIPrep derivatives" -keypoints: -- "fMRI data is dirty and needs to be cleaned, aligned, and transformed before being able to use" -- "There are standards in place which will allow you to effortlessly pull the data that you need for analysis" --- +:::::::::::::::::::::::::::::::::::::::: questions + +- What steps do I need to take before beginning to work with fMRI data? +- Learn about fMRIPrep derivatives + +:::::::::::::::::::::::::::::::::::::::::::::::::: + This document is based on the slides here [Google Slides](https://docs.google.com/presentation/d/1er6dQcERL-Yeb5-7A29tJnmqgHNaLpTLXM3e-SmpjDg/edit#slide=id.g484812a0c7_6_1) + ## Functional Neuroimaging in Python Welcome to the **Functional Neuroimaging Analysis in Python** workshop! In this workshpo we'll get you up to speed with the current tools and techniques used in the analysis of functional MRI (fMRI) data. The primary goals of this workshop are: @@ -58,62 +60,69 @@ First we'll want to deal with our structural data; this is called the **T1 image 1. Brain extraction - we want to analyze brains, not skulls 2. Normalization - since brains are different across people, we need a method to make them look more alike so we can perform group analysis. This is achieved using a non-linear warp which *squishes and pulls* the brain to look like a *template image* -![T1 Normalization](../fig/animated_t1_mni.gif){:class="img-responsive"} +![](fig/animated_t1_mni.gif){alt='T1 Normalization' class="img-responsive"} ### Visual Guide to Pre-processing fMRI data With fMRI data things are a bit more complicated since you have to deal with: 1. Motion across time -2. Distortion artifacts due to magnetic field inhomogeneities +2. Distortion artifacts due to magnetic field inhomogeneities The following steps are required (at minimum!): 1. Brain extraction - again we're only interested in the brain 2. Motion correction - we need to align the fMRI data *across time* 3. Susceptibility Distortion Correction (SDC) - we need to correct for magnetic field inhomogeneities -4. Alignment to the T1 image - aligning to the T1 image allows us to perform the *squishy/pully" normalization to make everyone's brain more alike +4. Alignment to the T1 image - aligning to the T1 image allows us to perform the \*squishy/pully" normalization to make everyone's brain more alike 5. Confound regression - not only does motion *misalign brains* it also corrupts the signal with *motion signal artifacts*, this also needs to be cleaned! - -![fMRI Preprocessing Steps](../fig/animated_fmri_preproc.gif){:class="img-responsive"} +![](fig/animated_fmri_preproc.gif){alt='fMRI Preprocessing Steps' class="img-responsive"} So **how does one begin to even accomplish this**? Traditionally, neuroimagers used a plethora of tools like, but not limited to: **FSL**, **AFNI**, **FREESURFER**, **ANTS**, **SPM**. Each with their own quirks and file format requirements. Unfortunately this is difficult to navigate, and each tool develops new techniques to better peform each of these pre-processing steps. Luckily, if your data is in a **BIDS Format**, there exists a tool, [**fMRIPrep**](https://fmriprep.org), which does this all for you while using the best methods across *most of these tools*!. An image below from their website depicts the processing steps they use: - -![fMRIPrep's Workflow](https://github.com/oesteban/fmriprep/raw/38a63e9504ab67812b63813c5fe9af882109408e/docs/_static/fmriprep-workflow-all.png){:class="img-responsive"} - +![](https://github.com/oesteban/fmriprep/raw/38a63e9504ab67812b63813c5fe9af882109408e/docs/_static/fmriprep-workflow-all.png){alt="fMRIPrep's Workflow" class="img-responsive"} Ultimately, fMRIPrep is an end-to-end pipeline - meaning that you feed it your raw organized data and it'll produce a bunch of outputs that you can use for analysis! What follows below are explanations of what those outputs are: #### fMRIPrep anatomical outputs -![fMRIPrep Anatomical Outputs](../fig/fmriprep_anat_out.png){:class="img-static"} +![](fig/fmriprep_anat_out.png){alt='fMRIPrep Anatomical Outputs' class="img-static"} ##### Native Space -1. **sub-xxxxx_T1w_brainmask.nii** - a binary mask which can be used to pull out just the brain -2. **sub-xxxxx_T1w_preproc.nii** - the fully cleaned T1 image which *has not been normalized* + +1. **sub-xxxxx\_T1w\_brainmask.nii** - a binary mask which can be used to pull out just the brain +2. **sub-xxxxx\_T1w\_preproc.nii** - the fully cleaned T1 image which *has not been normalized* ##### MNI (Normalized) Space -1. **sub-xxxxx_T1w_space-MNI152NLin2009cAsym_brainmask.nii.gz** - also a brain mask, but warped to fit a template brain (the template is MNI152NLin2009cAsym) -2. **sub-xxxxx_T1w_space-MNI152NLin2009cAsym_class-(CSF/GM/WM)_probtissue.nii.gz** - probability values for each of the tissue types. We won't get into too much detail with this one -3. **sub-xxxxx_T1w_space-MNI152NLin2009cAsym_preproc..nii** - the cleaned up T1 image that has been squished and warped into the MNI152NLin2009cAsym template space -#### fMRIPrep functional outputs +1. **sub-xxxxx\_T1w\_space-MNI152NLin2009cAsym\_brainmask.nii.gz** - also a brain mask, but warped to fit a template brain (the template is MNI152NLin2009cAsym) +2. **sub-xxxxx\_T1w\_space-MNI152NLin2009cAsym\_class-(CSF/GM/WM)\_probtissue.nii.gz** - probability values for each of the tissue types. We won't get into too much detail with this one +3. **sub-xxxxx\_T1w\_space-MNI152NLin2009cAsym\_preproc..nii** - the cleaned up T1 image that has been squished and warped into the MNI152NLin2009cAsym template space +#### fMRIPrep functional outputs -![fMRIPrep Functional Outputs](../fig/fmriprep_func_out.png){:class="img-static"} +![](fig/fmriprep_func_out.png){alt='fMRIPrep Functional Outputs' class="img-static"} -As above we have both **Native** and **Normalized** versions of the fMRI brain, we have a mask of each one as well as the preprocessed fMRI brain. +As above we have both **Native** and **Normalized** versions of the fMRI brain, we have a mask of each one as well as the preprocessed fMRI brain. **Note**: These have *not* been cleaned of motion artifacts. They have only been *aligned*, *distortion corrected*, and *skull-stripped*. -fMRIPrep also outputs a **sub-xxxxx_task-..._confounds.tsv** tab-delimited spreadsheet which contains a set of **nuisance regressors** that you can use to clean the *signal itself of motion artifacts*. We'll explore this in a later section. +fMRIPrep also outputs a **sub-xxxxx\_task-...\_confounds.tsv** tab-delimited spreadsheet which contains a set of **nuisance regressors** that you can use to clean the *signal itself of motion artifacts*. We'll explore this in a later section. ## End In the next section we'll start exploring the outputs generated by fMRIPrep to get a better handle of how to use them to manipulate images, clean motion signals, and perform analysis! -{% include links.md %} + + +:::::::::::::::::::::::::::::::::::::::: keypoints + +- fMRI data is dirty and needs to be cleaned, aligned, and transformed before being able to use +- There are standards in place which will allow you to effortlessly pull the data that you need for analysis + +:::::::::::::::::::::::::::::::::::::::::::::::::: + + diff --git a/episodes/02-exploring-fmriprep.md b/episodes/02-exploring-fmriprep.md index 739cf858..91176821 100644 --- a/episodes/02-exploring-fmriprep.md +++ b/episodes/02-exploring-fmriprep.md @@ -1,21 +1,26 @@ --- -title: "Exploring Preprocessed fMRI Data from fMRIPREP" +title: Exploring Preprocessed fMRI Data from fMRIPREP teaching: 20 exercises: 5 -questions: -- "How does fMRIPrep store preprocessed neuroimaging data" -- "How do I access preprocessed neuroimaging data" -objectives: -- "Learn about fMRIPrep derivatives" -- "Understand how preprocessed data is stored and how you can access key files for analysis" -keypoints: -- "fMRIPrep stores preprocessed data in a 'BIDS-like' fashion" -- "You can pull files using pyBIDS much like how you can navigate raw BIDS data" --- +::::::::::::::::::::::::::::::::::::::: objectives + +- Learn about fMRIPrep derivatives +- Understand how preprocessed data is stored and how you can access key files for analysis + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::: questions + +- How does fMRIPrep store preprocessed neuroimaging data +- How do I access preprocessed neuroimaging data + +:::::::::::::::::::::::::::::::::::::::::::::::::: + ## Exploring Preprocessed fMRI Data from fMRIPREP -BIDS applications such as fMRIPREP output data into a full data structure with strong similarity to BIDS organization principals. In fact, there is a specification for derivatives (outputs derived from) BIDS datasets; although this is a current work in progress, details can be found in: [BIDS Derivatives](https://bids-specification.readthedocs.io/en/latest/06-extensions.html). +BIDS applications such as fMRIPREP output data into a full data structure with strong similarity to BIDS organization principals. In fact, there is a specification for derivatives (outputs derived from) BIDS datasets; although this is a current work in progress, details can be found in: [BIDS Derivatives](https://bids-specification.readthedocs.io/en/latest/06-extensions.html). In this tutorial, we'll explore the outputs generated by fMRIPREP and get a handle of how the data is organized from this preprocessing pipeline @@ -25,12 +30,11 @@ Luckily the semi-standardized output for fMRIPrep is organized in such a way tha First let's take a quick look at the fMRIPrep data structure: -~~~ +```bash tree -L 1 '../data/ds000030/derivatives/fmriprep' -~~~ -{: .language-bash} +``` -~~~ +```output ../data/ds000030/derivatives/fmriprep/ ├── sub-10171 @@ -53,17 +57,15 @@ tree -L 1 '../data/ds000030/derivatives/fmriprep' ├── sub-50077 ├── sub-50081 └── sub-50083 -~~~ -{: .output} +``` First note that inside the fMRIPrep folder, we have a folder per-subject. Let's take a quick look at a single subject folder: -~~~ +```bash tree '../data/ds000030/derivatives/fmriprep/sub-10788' -~~~ -{: .language-bash} +``` -~~~ +```output ../data/ds000030/derivatives/fmriprep/sub-10788/ ├── anat │ ├── sub-10788_desc-aparcaseg_dseg.nii.gz @@ -106,18 +108,21 @@ tree '../data/ds000030/derivatives/fmriprep/sub-10788' │ ├── sub-10788_task-rest_space-T1w_desc-preproc_bold.json │ └── sub-10788_task-rest_space-T1w_desc-preproc_bold.nii.gz ... -~~~ -{: .output} +``` -As you can see above, each subject folder is organized into an anat and func sub-folder. +As you can see above, each subject folder is organized into an anat and func sub-folder. Specifically: - the anat folder contains the preprocessed anatomical data. If multiple T1 files are available (all T1s even across sessions), then these data are merged - you will always have one anat folder under the subject folder -- the func folder contains the preprocessed functional data. All tasks are dumped into the same folder and like the BIDS convention are indicated by the use of their filenames (task-[task_here]) +- the func folder contains the preprocessed functional data. All tasks are dumped into the same folder and like the BIDS convention are indicated by the use of their filenames (task-[task\_here]) + +::::::::::::::::::::::::::::::::::::::::: callout + +This data is single-session, so a session folder is missing here - but with multiple sessions you will see anat and ses-[insert\_session\_here] folders where each session folder contain a func folder. -> This data is single-session, so a session folder is missing here - but with multiple sessions you will see anat and ses-[insert_session_here] folders where each session folder contain a func folder. -{: .callout} + +:::::::::::::::::::::::::::::::::::::::::::::::::: Hopefully you're now convinced that the outputs of fMRIPREP roughly follows BIDS organization principles. The filenames themselves give you a full description of what each file is (check the [slides](https://docs.google.com/presentation/d/1er6dQcERL-Yeb5-7A29tJnmqgHNaLpTLXM3e-SmpjDg/edit?usp=sharing) to get an idea of what each file means! @@ -125,35 +130,31 @@ Now let's see how we can pull data in using pyBIDS! Let's import pyBIDS through the bids module first: -~~~ +```python import bids -~~~ -{: .language-python} +``` We can make a bids.BIDSLayout object as usual by just feeding in the fmriprep directory! However, there is one caveat... note that fMRIPrep doesn't exactly adhere to the *standard BIDS convention*. It uses fields such as desc- which are not part of the original BIDS specification. I.e: -~~~ +```bash sub-10788_desc-preproc_T1w.nii.gz -~~~ -{: .language-bash} +``` In fact, BIDS allows for *extensions* which enable you to add additional fields to the standard BIDS convention (such as desc-!). fMRIprep uses the derivatives extension of the BIDS standard. pyBIDS can handle standard extensions to the BIDS specification quite easily: -~~~ +```python layout = bids.BIDSLayout('../data/ds000030/derivatives/fmriprep/', config=['bids','derivatives']) -~~~ -{: .language-python} +``` Now that we have a layout object, we can pretend like we're working with a BIDS dataset! Let's try some common commands that you would've used with a BIDS dataset: First, we'll demonstrate that we can grab a list of pre-processed subjects much like in the way we would grab subjects from a raw BIDS dataset: -~~~ +```python layout.get_subjects() -~~~ -{: .language-python} +``` -~~~ +```output ['10171', '10292', '10365', @@ -174,31 +175,27 @@ layout.get_subjects() '50077', '50081', '50083'] - ~~~ - {: .output} +``` - We can also do the same for tasks +We can also do the same for tasks - ~~~ - layout.get_tasks() - ~~~ - {: .language-python} +```python +layout.get_tasks() +``` - ~~~ - ['rest'] - ~~~ - {: .output} +```output +['rest'] +``` - Now let's try fetching specific files. Similar to how you would fetch BIDS data using pyBIDS, the exact same syntax will work for fMRIPREP derivatives. Let's try pulling just the preprocessed anatomical data. +Now let's try fetching specific files. Similar to how you would fetch BIDS data using pyBIDS, the exact same syntax will work for fMRIPREP derivatives. Let's try pulling just the preprocessed anatomical data. Recall that the anatomical folder is organized as follows: -~~~ +```bash tree '../data/ds000030/derivatives/fmriprep/sub-10788/anat'^ -~~~ -{: .language-bash} +``` -~~~ +```output ../data/ds000030/derivatives/fmriprep/sub-10788/anat ├── sub-10788_desc-aparcaseg_dseg.nii.gz ├── sub-10788_desc-aseg_dseg.nii.gz @@ -222,84 +219,97 @@ tree '../data/ds000030/derivatives/fmriprep/sub-10788/anat'^ └── sub-10788_space-MNI152NLin2009cAsym_label-WM_probseg.nii.gz 0 directories, 20 files -~~~ -{: .output} +``` -The file that we're interested in is of form sub-[subject]_desc-preproc_T1w.nii.gz. Now we can construct a pyBIDS call to pull these types of files specifically: +The file that we're interested in is of form sub-[subject]\_desc-preproc\_T1w.nii.gz. Now we can construct a pyBIDS call to pull these types of files specifically: -~~~ +```python preproc_T1 = layout.get(datatype='anat', desc='preproc', extension=".nii.gz") preproc_T1 -~~~ -{: .language-python} +``` -~~~ +```output [, , , , -~~~ -{: .output} +``` + +::::::::::::::::::::::::::::::::::::::::: callout + +If we didn't configure pyBIDS with config=['bids','derivatives'] then the desc keyword would not work! + -> If we didn't configure pyBIDS with config=['bids','derivatives'] then the desc keyword would not work! -{: .callout} +:::::::::::::::::::::::::::::::::::::::::::::::::: -Note that we also pulled in MNI152NLin2009cAsym_preproc.nii.gz data as well. This is data that has been transformed into MNI152NLin2009cAsym template space. We can pull this data out by further specifying our layout.get using the space argument: +Note that we also pulled in MNI152NLin2009cAsym\_preproc.nii.gz data as well. This is data that has been transformed into MNI152NLin2009cAsym template space. We can pull this data out by further specifying our layout.get using the space argument: -~~~ +```python mni_preproc_T1 = layout.get(datatype='anat',desc='preproc',extension='.nii.gz',space='MNI152NLin2009cAsym') mni_preproc_T1 -~~~ -{: .language-python} +``` -~~~ +```output [, , ... -~~~ -{: .output} +``` What if we wanted to pull out the data in T1 "native space" (it really is a template space, since it is merged T1s)? Unfortunately for this isn't directly possible using layout.get. Instead we'll use a bit of python magic to pull the data that we want: -~~~ +```python native_preproc_T1 = [f for f in preproc_T1 if f not in mni_preproc_T1] native_preproc_T1 -~~~ -{: .language-python} - +``` Similarily fMRI data can be pulled by specifying datatype='func' and using the desc argument as appropriate: -> ## Exercise 1 -> 1. Get the list of **all** preprocessed functional data -> 2. Get the list of functional data in MNI152NLin2009cAsym space -> 3. Get the list of functional data in T1w space (native) -> Note that T1 space fMRI data can be pulled using space="T1w" (this is unlike the T1w data which required you to do some filtering) -> -> > ## Solution -> > *All the functional data* -> > ~~~ -> > func_data = layout.get(datatype='func', desc='preproc') -> > ~~~ -> > {: .language-python} -> > -> > *MNI152NLin2009cAsym Functional Data* -> > ~~~ -> > mni_func_data = layout.get(datatype='func', desc='preproc', space='MNI152NLin2009cAsym') -> > mni_func_data -> > ~~~ -> > {: .language-python} -> > -> > *Native/T1w space functional data* -> > ~~~ -> > t1w_func_data = layout.get(datatype='func', desc='preproc', space='T1w') -> > t1w_func_data -> > ~~~ -> > {: .language-python} -> > -> {: .solution} -{: .challenge} +::::::::::::::::::::::::::::::::::::::: challenge + +## Exercise 1 + +1. Get the list of **all** preprocessed functional data +2. Get the list of functional data in MNI152NLin2009cAsym space +3. Get the list of functional data in T1w space (native) + Note that T1 space fMRI data can be pulled using space="T1w" (this is unlike the T1w data which required you to do some filtering) + +::::::::::::::: solution + +## Solution + +*All the functional data* + +```python +func_data = layout.get(datatype='func', desc='preproc') +``` + +*MNI152NLin2009cAsym Functional Data* + +```python +mni_func_data = layout.get(datatype='func', desc='preproc', space='MNI152NLin2009cAsym') +mni_func_data +``` + +*Native/T1w space functional data* + +```python +t1w_func_data = layout.get(datatype='func', desc='preproc', space='T1w') +t1w_func_data +``` + +::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::::::::::::: Now that we have a handle on how fMRIPREP preprocessed data is organized and how we can pull this data. Let's start working with the actual data itself! -{% include links.md %} + + +:::::::::::::::::::::::::::::::::::::::: keypoints + +- fMRIPrep stores preprocessed data in a 'BIDS-like' fashion +- You can pull files using pyBIDS much like how you can navigate raw BIDS data + +:::::::::::::::::::::::::::::::::::::::::::::::::: + + diff --git a/episodes/03-basic_image_manipulation.md b/episodes/03-basic_image_manipulation.md index 6ad43b9d..315325a2 100644 --- a/episodes/03-basic_image_manipulation.md +++ b/episodes/03-basic_image_manipulation.md @@ -1,26 +1,31 @@ --- -title: "Introduction to Image Manipulation using Nilearn" +title: Introduction to Image Manipulation using Nilearn teaching: 30 exercises: 15 -questions: -- "How can be perform arithmetic operations on MR images" -objectives: -- "Use Nilearn to perform masking and mathematical operations" -- "Learn how to resample across modalities for image viewing and manipulation" -keypoints: -- "MR images are essentially 3D arrays where each voxel is represented by an (i,j,k) index" -- "Nilearn is Nibabel under the hood, knowing how Nibabel works is key to understanding Nilearn" --- -# Introduction +::::::::::::::::::::::::::::::::::::::: objectives + +- Use Nilearn to perform masking and mathematical operations +- Learn how to resample across modalities for image viewing and manipulation + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::: questions + +- How can be perform arithmetic operations on MR images + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +## Introduction Nilearn is a functional neuroimaging analysis and visualization library that wraps up a whole bunch of high-level operations (machine learning, statistical analysis, data cleaning, etc...) in easy-to-use commands. The neat thing about Nilearn is that it implements Nibabel under the hood. What this means is that everything you do in Nilearn can be represented by performing a set of operations on Nibabel objects. This has the important consequence of allowing you, yourself to perform high-level operations (like resampling) using Nilearn then dropping into Nibabel for more custom data processing then jumping back up to Nilearn for interactive image viewing. Pretty cool! -# Setting up +## Setting up The first thing we'll do is to important some Python modules that will allow us to use Nilearn: -~~~ +```python import os import matplotlib.pyplot as plt from nilearn import image as nimg @@ -29,15 +34,16 @@ from bids import BIDSLayout #for inline visualization in jupyter notebook %matplotlib inline -~~~ -{: .language-python} +``` Notice that we imported two things: + 1. `image as nimg` - allows us to load NIFTI images using nibabel under the hood -2. `plotting as nplot`- allows us to using Nilearn's plotting library for easy visualization +2. `plotting as nplot`\- allows us to using Nilearn's plotting library for easy visualization + +First let's grab some data from where we downloaded our **FMRIPREP** outputs. Note that we're using the argument return\_type='file' so that pyBIDS gives us file paths directly rather than the standard BIDSFile objects -First let's grab some data from where we downloaded our **FMRIPREP** outputs. Note that we're using the argument return_type='file' so that pyBIDS gives us file paths directly rather than the standard BIDSFile objects -~~~ +```python #Base directory for fmriprep output fmriprep_dir = '../data/ds000030/derivatives/fmriprep/' layout= BIDSLayout(fmriprep_dir, config=['bids','derivatives']) @@ -49,103 +55,99 @@ T1w_files = layout.get(subject='10788', datatype='anat', brainmask_files = layout.get(subject='10788', datatype='anat', suffix="mask", desc='brain', extension='.nii.gz', return_type='file') -~~~ -{: .language-python} +``` Here we used pyBIDS (as introduced in earlier sections) to pull a single participant. Specifically, we pulled all preprocessed (MNI152NLin2009cAsym, and T1w) anatomical files as well as their respective masks. Let's quickly view these files for review: -~~~ +```python #Display preprocessed files inside of anatomy folder T1w_files -~~~ -{: .language-python} +``` -~~~ +```output ['/home/jerry/projects/workshops/SDC-BIDS-fMRI/data/ds000030/derivatives/fmriprep/sub-10788/anat/sub-10788_desc-preproc_T1w.nii.gz', '/home/jerry/projects/workshops/SDC-BIDS-fMRI/data/ds000030/derivatives/fmriprep/sub-10788/anat/sub-10788_space-MNI152NLin2009cAsym_desc-preproc_T1w.nii.gz'] -~~~ -{: .output} +``` Now that we have our files set up, let's start performing some basic image operations. -## Basic Image Operations +### Basic Image Operations In this section we're going to deal with the following files: -1. sub-10171_desc-preproc_T1w.nii.gz - the T1 image in native space -2. sub-10171_desc-brain_mask.nii.gz - a mask with 1's representing the brain and 0's elsewhere. +1. sub-10171\_desc-preproc\_T1w.nii.gz - the T1 image in native space +2. sub-10171\_desc-brain\_mask.nii.gz - a mask with 1's representing the brain and 0's elsewhere. -~~~ +```python t1 = T1w_files[0] bm = brainmask_files[0] t1_img = nimg.load_img(t1) bm_img = nimg.load_img(bm) -~~~ -{: .language-python} +``` Using the plotting module (which we've aliased as nplot), we can view our MR image: -~~~ +```python nplot.plot_anat(t1_img) -~~~ -{: .language-python} +``` -![Nilearn antomical plotting](../fig/t1_img.png){:class="img-responsive"} +![](fig/t1_img.png){alt='Nilearn antomical plotting' class="img-responsive"} This gives just a still image of the brain. We can also view the brain more interactively using the `view_img` function. It will require some additional settings however: -~~~ +```python nplot.view_img(t1_img, bg_img=False, # Disable using a standard image as the background cmap='Greys_r', # Set color scale so white matter appears lighter than grey symmetric_cmap=False, # We don't have negative values threshold="auto", # Clears out the background ) -~~~ -{: .language-python} +``` Try clicking and dragging the image in each of the views that are generated! - Try viewing the mask as well! -~~~ +```python nplot.plot_anat(bm_img) -~~~ -{: .language-python} +``` -### Arithmetic Operations +#### Arithmetic Operations Let's start performing some image operations. The simplest operations we can perform is **element-wise**, what this means is that we want to perform some sort of mathematical operation on each **voxel** of the MR image. Since *voxels are represented in a 3D array, this is equivalent to performing an operation on each element (i,j,k) of a 3D array*. Let's try inverting the image, that is, flip the colour scale such that all blacks appear white and vice-versa. To do this, we'll use the method `nimg.math_img(formula, **imgs)` Where: + - `formula` is a mathematical expression such as `'a+1'` - `**imgs` is a set of key-value pairs linking variable names to images. For example `a=T1` In order to invert the image, we can simply flip the sign which will set the most positive elements (white) to the most negatve elements (black), and the least positives elements (black) to the least negative elements (white). This effectively flips the colour-scale: -~~~ +```python invert_img = nimg.math_img('-a', a=t1_img) nplot.plot_anat(invert_img) -~~~ -{: .language-python} +``` -![Nilearn image math example output](../fig/invert_img.png){:class="img-responsive"} +![](fig/invert_img.png){alt='Nilearn image math example output' class="img-responsive"} ->> Alternatively we don't need to first load in our t1_img using img.load_img. Instead we can feed in a path to img.math_img: +> > Alternatively we don't need to first load in our t1\_img using img.load\_img. Instead we can feed in a path to img.math\_img: -> ~~~ -> invert_img = nimg.math_img('-a', a=t1) -> nplot.plot_anat(invert_img) -> ~~~ -> -> This will yield the same result! -{: .callout} +::::::::::::::::::::::::::::::::::::::::: callout +``` +invert_img = nimg.math_img('-a', a=t1) +nplot.plot_anat(invert_img) +``` + +This will yield the same result! + + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +#### Applying a Mask -### Applying a Mask Let's extend this idea of applying operations to each element of an image to multiple images. Instead of specifying just one image like the following: `nimg.math_img('a+1',a=img_a)` @@ -158,87 +160,100 @@ The key requirement here is that when dealing with multiple images, that the *si We can take advantage of this property when masking our data using multiplication. Masking works by multipling a raw image (our `T1`), with some mask image (our `bm`). Whichever voxel (i,j,k) has a value of 0 in the mask multiplies with voxel (i,j,k) in the raw image resulting in a product of 0. Conversely, any voxel (i,j,k) in the mask with a value of 1 multiplies with voxel (i,j,k) in the raw image resulting in the same value. Let's try this out in practice and see what the result is: -~~~ +```python masked_t1 = nimg.math_img('a*b', a=t1, b=bm) nplot.plot_anat(masked_t1) -~~~ -{: .language-python} +``` -![Nilearn image masking output](../fig/masked_t1.png){:class="img-responsive"} +![](fig/masked_t1.png){alt='Nilearn image masking output' class="img-responsive"} As you can see areas where the mask image had a value of 1 were retained, everything else was set to 0 -> ## Exercise #1 -> Try applying the mask such that the brain is removed, but the rest of the head is intact! -> *Hint:* -> Remember that a mask is composed of 0's and 1's, where parts of the data labelled 1 are regions to keep, and parts of the data that are 0, are to throw away. -> You can do this in 2 steps: -> -> 1. Switch the 0's and 1's using an equation (simple addition/substraction) or condition (like x == 0). -> 2. Apply the mask -> > ## Solution -> > ~~~ -> > # Invert the mask -> > inverted_mask = nimg.math_img('1-x', x=bm) -> > nplot.plot_anat(inverted_mask) -> > ~~~ -> > {: .language-python} -> > ~~~ -> > # Apply the mask -> > inverted_mask_t1 = nimg.math_img('a*b', a=t1, b=inverted_mask) -> > nplot.plot_anat(inverted_mask_t1) -> > ~~~ -> > {: .language-python} -> > ![Episode 03 Exercise 1 inverted mask](../fig/inverted_mask_t1.png){:class="img-responsive"} -> {: .solution} -{: .challenge} - -## Slicing +::::::::::::::::::::::::::::::::::::::: challenge + +### Exercise #1 + +Try applying the mask such that the brain is removed, but the rest of the head is intact! +*Hint:* +Remember that a mask is composed of 0's and 1's, where parts of the data labelled 1 are regions to keep, and parts of the data that are 0, are to throw away. +You can do this in 2 steps: + +1. Switch the 0's and 1's using an equation (simple addition/substraction) or condition (like x == 0). +2. Apply the mask + +::::::::::::::: solution + +### Solution + +```python +# Invert the mask +inverted_mask = nimg.math_img('1-x', x=bm) +nplot.plot_anat(inverted_mask) +``` + +```python +# Apply the mask +inverted_mask_t1 = nimg.math_img('a*b', a=t1, b=inverted_mask) +nplot.plot_anat(inverted_mask_t1) +``` + +![](fig/inverted_mask_t1.png){alt='Episode 03 Exercise 1 inverted mask' class="img-responsive"} + + + +::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +### Slicing Recall that our data matrix is organized in the following manner: -### -![3D Array Representation](../fig/slicing.png){:class="img-responsive"} +#### + +![](fig/slicing.png){alt='3D Array Representation' class="img-responsive"} Slicing does exactly what it seems to imply. Given our 3D volume, we can pull out 2D subsets (called "slices"). Here's an example of slicing moving from left to right via an animated GIF: -### -![Animated Slicing of T1](../fig/animated_slicing.gif){:class="img-responsive"} +#### + +![](fig/animated_slicing.gif){alt='Animated Slicing of T1' class="img-responsive"} What you see here is a series of 2D images that start from the left, and move toward the right. Each frame of this GIF is a slice - a 2D subset of a 3D volume. Slicing can be useful for cases in which you'd want to loop through each MR slice and perform a computation; importantly in functional imaging data slicing is useful for pulling out timepoints as we'll see later! -> Sourced from: https://en.wikipedia.org/wiki/Neuroimaging#/media/File:Parasagittal_MRI_of_human_head_in_patient_with_benign_familial_macrocephaly_prior_to_brain_injury_(ANIMATED).gif -{: .callout} +::::::::::::::::::::::::::::::::::::::::: callout + +Sourced from: [https://en.wikipedia.org/wiki/Neuroimaging#/media/File:Parasagittal\_MRI\_of\_human\_head\_in\_patient\_with\_benign\_familial\_macrocephaly\_prior\_to\_brain\_injury\_(ANIMATED).gif](https://en.wikipedia.org/wiki/Neuroimaging#/media/File:Parasagittal_MRI_of_human_head_in_patient_with_benign_familial_macrocephaly_prior_to_brain_injury_\(ANIMATED\).gif) + + +:::::::::::::::::::::::::::::::::::::::::::::::::: Slicing is done easily on an image file using the attribute .slicer of a Nilearn image object. For example we can grab the $10^{\text{th}}$ slice along the x axis as follows: -~~~ +```python x_slice = t1_img.slicer[10:11,:,:] -~~~ -{: .language-python} +``` The statement $10:11$ is intentional and is required by .slicer. Alternatively we can slice along the x-axis using the data matrix itself: -~~~ +```python t1_data = t1_img.get_data() x_slice = t1_data[10,:,:] -~~~ -{: .language-python} +``` -This will yield the same result as above. Notice that when using the t1_data array we can just specify which slice to grab instead of using :. We can use slicing in order to modify visualizations. For example, when viewing the T1 image, we may want to specify at which slice we'd like to view the image. This can be done by specifying which coordinates to *cut* the image at: +This will yield the same result as above. Notice that when using the t1\_data array we can just specify which slice to grab instead of using :. We can use slicing in order to modify visualizations. For example, when viewing the T1 image, we may want to specify at which slice we'd like to view the image. This can be done by specifying which coordinates to *cut* the image at: -~~~ +```python nplot.plot_anat(t1_img,cut_coords=(50,30,20)) -~~~ -{: .language-python} +``` + +The cut\_coords option specifies 3 numbers: -The cut_coords option specifies 3 numbers: - The first number says cut the X coordinate at slice 50 and display (sagittal view in this case!) - The second number says cut the Y coordinate at slice 30 and display (coronal view) - The third number says cut the Z coordinate at slice 20 and display (axial view) -Remember nplot.plot_anat yields 3 images, therefore cut_coords allows you to display where to take cross-sections of the brain from different perspectives (axial, sagittal, coronal) - +Remember nplot.plot\_anat yields 3 images, therefore cut\_coords allows you to display where to take cross-sections of the brain from different perspectives (axial, sagittal, coronal) This covers the basics of image manipulation using T1 images. To review in this section we covered: @@ -248,4 +263,13 @@ This covers the basics of image manipulation using T1 images. To review in this In the next section we will cover how to integrate additional modalities (functional data) to what we've done so far using Nilearn. Then we can start using what we've learned in order to perform analysis and visualization! -{% include links.md %} + + +:::::::::::::::::::::::::::::::::::::::: keypoints + +- MR images are essentially 3D arrays where each voxel is represented by an (i,j,k) index +- Nilearn is Nibabel under the hood, knowing how Nibabel works is key to understanding Nilearn + +:::::::::::::::::::::::::::::::::::::::::::::::::: + + diff --git a/episodes/04-integrating_functional_data.md b/episodes/04-integrating_functional_data.md index 0a94c521..06ef83bd 100644 --- a/episodes/04-integrating_functional_data.md +++ b/episodes/04-integrating_functional_data.md @@ -1,37 +1,41 @@ --- -title: "Integrating Functional Data" +title: Integrating Functional Data teaching: 30 exercises: 15 -questions: -- "How is fMRI data represented" -- "How can we access fMRI data along spatial and temporal dimensions" -- "How can we integrate fMRI and structural MRI together" -objectives: -- "Extend the idea of slicing to 4 dimensions" -- "Apply resampling to T1 images to combine them with fMRI data" -keypoints: -- "fMRI data is represented by spatial (x,y,z) and temporal (t) dimensions, totalling 4 dimensions" -- "fMRI data is at a lower resolution than structural data. To be able to combine data requires resampling your data" --- -# Integrating Functional Data +::::::::::::::::::::::::::::::::::::::: objectives -So far most of our work has been examining anatomical images - the reason being is that it provides a nice visual way of exploring the effects of data manipulation and visualization is easy. In practice, you will most likely not analyze anatomical data using nilearn since there are other tools that are better suited for that kind of analysis (freesurfer, connectome-workbench, mindboggle, etc...). +- Extend the idea of slicing to 4 dimensions +- Apply resampling to T1 images to combine them with fMRI data + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::: questions + +- How is fMRI data represented +- How can we access fMRI data along spatial and temporal dimensions +- How can we integrate fMRI and structural MRI together + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +## Integrating Functional Data + +So far most of our work has been examining anatomical images - the reason being is that it provides a nice visual way of exploring the effects of data manipulation and visualization is easy. In practice, you will most likely not analyze anatomical data using nilearn since there are other tools that are better suited for that kind of analysis (freesurfer, connectome-workbench, mindboggle, etc...). In this notebook we'll finally start working with functional MR data - the modality of interest in this workshop. First we'll cover some basics about how the data is organized (similar to T1s but slightly more complex), and then how we can integrate our anatomical and functional data together using tools provided by nilearn Functional data consists of full 3D brain volumes that are *sampled* at multiple time points. Therefore you have a sequence of 3D brain volumes, stepping through sequences is stepping through time and therefore time is our 4th dimension! Here's a visualization to make this concept more clear: -![4D Array Representation](../fig/4D_array.png){:class="img-responsive"} +![](fig/4D_array.png){alt='4D Array Representation' class="img-responsive"} Each index along the 4th dimensions (called TR for "Repetition Time", or Sample) is a full 3D scan of the brain. Pulling out volumes from 4-dimensional images is similar to that of 3-dimensional images except you're now dealing with: - nimg.slicer[x,y,z,time] ! Let's try a couple of examples to familiarize ourselves with dealing with 4D images. But first, let's pull some functional data using PyBIDS! -~~~ +```python import os import matplotlib.pyplot as plt #to enable plotting within notebook from nilearn import image as nimg @@ -39,14 +43,13 @@ from nilearn import plotting as nplot from bids.layout import BIDSLayout import numpy as np %matplotlib inline -~~~ -{: .language-python} +``` These are the usual imports. Let's now pull some structural *and* functional data using pyBIDS. We'll be using functional files in MNI space rather than T1w space. Recall, that MNI space data is data that was been warped into standard space. These are the files you would typically use for a group-level functional imaging analysis! -~~~ +```python fmriprep_dir = '../data/ds000030/derivatives/fmriprep/' layout=BIDSLayout(fmriprep_dir, validate=False, config=['bids','derivatives']) @@ -76,107 +79,101 @@ func_mask_files = layout.get(subject='10788', space='MNI152NLin2009cAsym', extension="nii.gz", return_type='file') -~~~ -{: .language-python} - +``` -~~~ +```python func_mni = func_files[0] func_mni_img = nimg.load_img(func_mni) -~~~ -{: .language-python} +``` ## fMRI Dimensions First note that fMRI data contains both spatial dimensions (x,y,z) and a temporal dimension (t). This would mean that we require 4 dimensions in order to represent our data. Let's take a look at the shape of our data matrix to confirm this intuition: -~~~ +```python func_mni_img.shape -~~~ -{: .language-python} +``` -~~~ +```output (65, 77, 49, 152) -~~~ -{: .output} +``` -Notice that the Functional MR scan contains *4 dimensions*. This is in the form of $(x,y,z,t)$, where $t$ is time. -We can use slicer as usual where instead of using 3 dimensions we use 4. +Notice that the Functional MR scan contains *4 dimensions*. This is in the form of $(x,y,z,t)$, where $t$ is time. +We can use slicer as usual where instead of using 3 dimensions we use 4. For example: - func.slicer[x,y,z] + func.slicer[x,y,z] vs. func.slicer[x,y,z,t] +::::::::::::::::::::::::::::::::::::::: challenge + +## Exercise + +Try pulling out the 5th TR and visualizing it using nplot.view\_img. -> ## Exercise -> -> Try pulling out the 5th TR and visualizing it using nplot.view_img. -> -> Remember that nplot.view_img provides an interactive view of the brain, try scrolling around! -> -> > ## Solution -> > ~~~ -> > #Pull the 5th TR -> > func_vol5 = func_mni_img.slicer[:,:,:,4] -> > nplot.view_img(func_vol5) -> > ~~~ -> > {: .language-python} -> {: .solution} -{: .challenge} - -You may also use nplot.plot_epi. plot_epi is exactly the same as plot_anat except it displays using colors that make more sense for functional images... - -~~~ +Remember that nplot.view\_img provides an interactive view of the brain, try scrolling around! + +::::::::::::::: solution + +## Solution + +```python +#Pull the 5th TR +func_vol5 = func_mni_img.slicer[:,:,:,4] +nplot.view_img(func_vol5) +``` + +::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +You may also use nplot.plot\_epi. plot\_epi is exactly the same as plot\_anat except it displays using colors that make more sense for functional images... + +```python #Pull the 5th TR nplot.plot_epi(func_vol5) -~~~ -{: .language-python} +``` -![Visual of fMRI EPI Data](../fig/fmri_data.png){:class="img-responsive"} +![](fig/fmri_data.png){alt='Visual of fMRI EPI Data' class="img-responsive"} ## What fMRI actually represents -We've represented fMRI as a snapshot of MR signal over multiple timepoints. This is a useful way of understanding the organization of fMRI, however it isn't typically how we think about the data when we analyze fMRI data. fMRI is typically thought of as **time-series** data. We can think of each voxel (x,y,z coordinate) as having a time-series of length T. The length T represents the number of volumes/timepoints in the data. Let's pick an example voxel and examine its time-series using func_mni_img.slicer: +We've represented fMRI as a snapshot of MR signal over multiple timepoints. This is a useful way of understanding the organization of fMRI, however it isn't typically how we think about the data when we analyze fMRI data. fMRI is typically thought of as **time-series** data. We can think of each voxel (x,y,z coordinate) as having a time-series of length T. The length T represents the number of volumes/timepoints in the data. Let's pick an example voxel and examine its time-series using func\_mni\_img.slicer: -~~~ +```python #Pick one voxel at coordinate (60,45,88) single_vox = func_mni_img.slicer[59:60,45:46,30:31,:].get_data() single_vox.shape -~~~ -{: .language-python} +``` -~~~ +```output (1, 1, 1, 152) -~~~ -{: .output} +``` As you can see we have 1 element in (x,y,z) dimension representing a single voxel. In addition, we have 152 elements in the fourth dimension. In totality, this means we have a single voxel with 152 timepoints. Dealing with 4 dimensional arrays are difficult to work with - since we have a single element across the first 3 dimensions we can squish this down to a 1 dimensional array with 152 time-points. We no longer need the first 3 spatial dimensions since we're only looking at one voxel and don't need (x,y,z) anymore: -~~~ +```python single_vox = single_vox.flatten() single_vox.shape -~~~ -{: .language-python} +``` -~~~ +```output (152,) -~~~ -{: .output} +``` Here we've pulled out a voxel at a specific coordinate at every single time-point. This voxel has a single value for each timepoint and therefore is a time-series. We can visualize this time-series signal by using a standard python plotting library. We won't go into too much detail about python plotting, the intuition about what the data looks like is what is most important: First let's import the standard python plotting library matplotlib: -~~~ +```python import matplotlib.pyplot as plt -~~~ -{: .language-python} +``` -~~~ +```python # Make an array counting from 0 --> 152, this will be our x-axis x_axis = np.arange(0, single_vox.shape[0]) @@ -186,51 +183,52 @@ plt.plot( x_axis, single_vox, 'k') # Label our axes plt.xlabel('Timepoint') plt.ylabel('Signal Value') -~~~ -{: .language-python} +``` -![Example fMRI Timeseries](../fig/timeseries.png){:class="img-responsive"} +![](fig/timeseries.png){alt='Example fMRI Timeseries' class="img-responsive"} As you can see from the image above, fMRI data really is just a signal per voxel over time! ## Resampling + Recall from our introductory exploration of neuroimaging data: - T1 images are typically composed of voxels that are 1x1x1 in dimension - Functional images are typically composed of voxels that are 4x4x4 in dimension -If we'd like to overlay our functional on top of our T1 (for visualization purposes, or analyses), then we need to match the size of the voxels! +If we'd like to overlay our functional on top of our T1 (for visualization purposes, or analyses), then we need to match the size of the voxels! Think of this like trying to overlay a 10x10 JPEG and a 20x20 JPEG on top of each other. To get perfect overlay we need to resize (or more accurately *resample*) our JPEGs to match! ->Resampling is a method of interpolating in between data-points. When we stretch an image we need to figure out what goes in the spaces that are created via stretching - resampling does just that. In fact, resizing any type of image is actually just resampling to new dimensions. -{: .callout} +::::::::::::::::::::::::::::::::::::::::: callout + +Resampling is a method of interpolating in between data-points. When we stretch an image we need to figure out what goes in the spaces that are created via stretching - resampling does just that. In fact, resizing any type of image is actually just resampling to new dimensions. -Let's resampling some MRI data using nilearn. + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +Let's resampling some MRI data using nilearn. **Goal**: Match the dimensions of the structural image to that of the functional image -~~~ +```python # Files we'll be using (Notice that we're using _space-MNI... # which means they are normalized brains) T1_mni = T1w_files[0] T1_mni_img = nimg.load_img(T1_mni) -~~~ -{: .language-python} +``` Let's take a quick look at the sizes of both our functional and structural files: -~~~ +```python print(T1_mni_img.shape) print(func_mni_img.shape) -~~~ -{: .language-python} +``` -~~~ +```output (193, 229, 193) (60, 77, 49, 152) -~~~ -{: .output} +``` Taking a look at the spatial dimensions (first three dimensions), we can see that the number of voxels in the T1 image does not match that of the fMRI image. This is because the fMRI data (which has less voxels) is a *lower resolution image*. We either need to *upsample* our fMRI image to match that of the T1 image, or we need to *downsample* our T1 image to match that of the fMRI image. Typically, since the fMRI data is the one we'd like to ultimately use for analysis, we would leave it alone and downsample our T1 image. The reason being is that *resampling* requires interpolating values which may contaminate our data with artifacts. We don't mind having artifacts in our T1 data (for visualization purposes) since the fMRI data is the one actually being analyzed. @@ -238,192 +236,205 @@ Resampling in nilearn is as easy as telling it which image you want to sample an Structure of function: -nimg.resample_to_img(source_img,target_img,interpolation) -- source_img = the image you want to sample -- target_img = the image you wish to *resample to* +nimg.resample\_to\_img(source\_img,target\_img,interpolation) + +- source\_img = the image you want to sample +- target\_img = the image you wish to *resample to* - interpolation = the method of interpolation -> A note on **interpolation** -> nilearn supports 3 types of interpolation, the one you'll use depends on the type of data you're resampling! -> 1. **continuous** - Interpolate but maintain some edge features. Ideal for structural images where edges are well-defined. Uses $3^\text{rd}$-order spline interpolation. -> 2. **linear (default)** - Interpolate uses a combination of neighbouring voxels - will blur. Uses trilinear interpolation. -> 3. **nearest** - matches value of closest voxel (majority vote from neighbours). This is ideal for masks which are binary since it will preserve the 0's and 1's and will not produce in-between values (ex: 0.342). Also ideal for numeric labels where values are 0,1,2,3... (parcellations). Uses nearest-neighbours interpolation with majority vote. -{: .callout} +::::::::::::::::::::::::::::::::::::::::: callout + +A note on **interpolation** +nilearn supports 3 types of interpolation, the one you'll use depends on the type of data you're resampling! + +1. **continuous** - Interpolate but maintain some edge features. Ideal for structural images where edges are well-defined. Uses $3^\\text{rd}$-order spline interpolation. +2. **linear (default)** - Interpolate uses a combination of neighbouring voxels - will blur. Uses trilinear interpolation. +3. **nearest** - matches value of closest voxel (majority vote from neighbours). This is ideal for masks which are binary since it will preserve the 0's and 1's and will not produce in-between values (ex: 0.342). Also ideal for numeric labels where values are 0,1,2,3... (parcellations). Uses nearest-neighbours interpolation with majority vote. + -~~~ +:::::::::::::::::::::::::::::::::::::::::::::::::: + +```python #Try playing around with methods of interpolation #options: 'linear','continuous','nearest' resamp_t1 = nimg.resample_to_img(source_img=T1_mni_img,target_img=func_mni_img,interpolation='continuous') print(resamp_t1.shape) print(func_mni_img.shape) nplot.plot_anat(resamp_t1) -~~~ -{: .language-python} +``` -![Downsampled T1](../fig/downsample_t1.png){:class="img-responsive"} +![](fig/downsample_t1.png){alt='Downsampled T1' class="img-responsive"} Now that we've explored the idea of resampling let's do a cumulative exercise bringing together ideas from resampling and basic image operations. -> ## Exercise -> -> Using **Native** T1 and **T1w** resting state functional do the following: -> 1. Resample the T1 image to resting state size -> 2. Replace the brain in the T1 image with the first frame of the resting state brain -> -> ### Files we'll need -> -> -> #### Structural Files -> -> ~~~ -> #T1 image -> ex_t1 = nimg.load_img(T1w_files[0]) -> -> #mask file -> ex_t1_bm = nimg.load_img(brainmask_files[0]) -> ~~~ -> {: .language-python} -> -> #### Functional Files -> -> ~~~ -> #This is the pre-processed resting state data that hasn't been standardized -> ex_func = nimg.load_img(func_files[0]) -> -> #This is the associated mask for the resting state image. -> ex_func_bm = nimg.load_img(func_mask_files[0]) -> ~~~ -> {: .language-python} -> -> The first step is to remove the brain from the T1 image so that we're left with a hollow skull. This can be broken down into 2 steps: -> -> 1. Invert the mask so that all 1's become 0's and all 0's become 1's -> ~~~ -> # Invert the T1 mask -> invert_mask = nimg.math_img('??', a=??) -> nplot.plot_anat(??) -> ~~~ -> {: .language-python} -> -> ![Episode 04 Exercise Inverted Mask](../fig/exercise_inverted_mask.png){:class="img-responsive"} -> -> 2. Apply the mask onto the T1 image, this will effectively remove the brain -> -> ~~~ -> # Apply the mask onto the T1 image -> hollow_skull = nimg.math_img("??", a=??, b=??) -> nplot.plot_anat(??) -> ~~~ -> {: .language-python} -> -> ![Episode 04 Exercise Hollow Skull](../fig/exercise_hollow_skull.png){:class="img-responsive"} -> -> Our brain is now missing! -> -> Next we need to *resize* the hollow skull image to the dimensions of our resting state image. This can be done using resampling as we've done earlier in this episode. -> -> What kind of interpolation would we need to perform here? Recall that: -> - **Continuous**: Tries to maintain the edges of the image -> - **Linear**: Resizes the image but also blurs it a bit -> - **Nearest**: Sets values to the closest neighbouring values -> -> ~~~ -> #Resample the T1 to the size of the functional image! -> resamp_skull = nimg.resample_to_img(source_img=??, -> target_img=??, -> interpolation='??') -> nplot.plot_anat(resamp_skull) -> print(resamp_skull.shape) -> print(ex_func.shape) -> ~~~ -> {: .language-python} -> -> ![Episode 04 Exercise Resampled Hollow Skull](../fig/exercise_resampled_hollow_skull.png){:class="img-responsive"} -> -> We now have a skull missing the structural T1 brain that is resized to match the dimensions of the EPI image. -> -> The final steps are to: -> 1. Pull the first volume from the functional image -> 2. Place the functional image head into the hollow skull that we've created -> -> Since a functional image is 4-Dimensional, we'll need to pull the first volume to work with. This is because the structural image is 3-dimensional and operations will fail if we try to mix 3D and 4D data. -> -> ~~~ -> #Let's visualize the first volume of the functional image: -> first_vol = ex_func.slicer[??, ??, ??, ??] -> nplot.plot_epi(first_vol) -> ~~~ -> {: .language-python} -> -> -> ![Episode 04 Exercise fMRI](../fig/exercise_fmri.png){:class="img-responsive"} -> -> As shown in the figure above, the image has some "signal" outside of the brain. In order to place this within the now brainless head we made earlier, we need to mask out the functional MR data as well! -> -> ~~~ -> #Mask the first volume using ex_func_bm -> masked_func = nimg.math_img('??', a=??, b=??) -> nplot.plot_epi(masked_func) -> ~~~ -> {: .language-python} -> -> ![Episode 04 Exercise Masked fMRI](../fig/exercise_masked_fmri.png){:class="img-responsive"} -> -> The final step is to stick this data into the head of the T1 data. Since the hole in the T1 data is represented as $0$'s. We can add the two images together to place the functional data into the void: -> -> ~~~ -> #Now overlay the functional image on top of the anatomical -> combined_img = nimg.math_img(??) -> nplot.plot_anat(combined_img) -> ~~~ -> {: .language-python} -> -> ![Episode 04 Exercise Complete](../fig/exercise_complete.png){:class="img-responsive"} -> -> > ## Solution -> > -> > ~~~ -> > # Invert the mask -> > invert_mask = nimg.math_img('1-a', a=ex_t1_bm) -> > nplot.plot_anat(invert_mask) -> > -> > # Apply the mask onto the T1 image -> > hollow_skull = nimg.math_img("a*b", a=ex_t1, b=invert_mask) -> > nplot.plot_anat(hollow_skull) -> > -> > #Resample the T1 to the size of the functional image! -> > resamp_skull = nimg.resample_to_img(source_img=hollow_skull, -> > target_img=ex_func, -> > interpolation='continuous') -> > nplot.plot_anat(resamp_skull) -> > print(resamp_skull.shape) -> > print(ex_func.shape) -> > -> > #Let's visualize the first volume of the functional image: -> > first_vol = ex_func.slicer[:,:,:,0] -> > nplot.plot_epi(first_vol) -> > -> > masked_func = nimg.math_img('a*b', a=first_vol, b=ex_func_bm) -> > nplot.plot_epi(masked_func) -> > -> > #Now overlay the functional image on top of the anatomical -> > combined_img = nimg.math_img('a+b', -> > a=resamp_skull, -> > b=masked_func) -> > nplot.plot_anat(combined_img) -> > ~~~ -> > {: .language-python} -> {: .solution} -{: .challenge} +::::::::::::::::::::::::::::::::::::::: challenge + +## Exercise + +Using **Native** T1 and **T1w** resting state functional do the following: + +1. Resample the T1 image to resting state size +2. Replace the brain in the T1 image with the first frame of the resting state brain + +### Files we'll need + +#### Structural Files + +```python +#T1 image +ex_t1 = nimg.load_img(T1w_files[0]) + +#mask file +ex_t1_bm = nimg.load_img(brainmask_files[0]) +``` + +#### Functional Files + +```python +#This is the pre-processed resting state data that hasn't been standardized +ex_func = nimg.load_img(func_files[0]) + +#This is the associated mask for the resting state image. +ex_func_bm = nimg.load_img(func_mask_files[0]) +``` + +The first step is to remove the brain from the T1 image so that we're left with a hollow skull. This can be broken down into 2 steps: + +1. Invert the mask so that all 1's become 0's and all 0's become 1's + +```python +# Invert the T1 mask +invert_mask = nimg.math_img('??', a=??) +nplot.plot_anat(??) +``` + +![](fig/exercise_inverted_mask.png){alt='Episode 04 Exercise Inverted Mask' class="img-responsive"} + +2. Apply the mask onto the T1 image, this will effectively remove the brain + +```python +# Apply the mask onto the T1 image +hollow_skull = nimg.math_img("??", a=??, b=??) +nplot.plot_anat(??) +``` + +![](fig/exercise_hollow_skull.png){alt='Episode 04 Exercise Hollow Skull' class="img-responsive"} + +Our brain is now missing! + +Next we need to *resize* the hollow skull image to the dimensions of our resting state image. This can be done using resampling as we've done earlier in this episode. + +What kind of interpolation would we need to perform here? Recall that: + +- **Continuous**: Tries to maintain the edges of the image +- **Linear**: Resizes the image but also blurs it a bit +- **Nearest**: Sets values to the closest neighbouring values + +```python +#Resample the T1 to the size of the functional image! +resamp_skull = nimg.resample_to_img(source_img=??, + target_img=??, + interpolation='??') +nplot.plot_anat(resamp_skull) +print(resamp_skull.shape) +print(ex_func.shape) +``` + +![](fig/exercise_resampled_hollow_skull.png){alt='Episode 04 Exercise Resampled Hollow Skull' class="img-responsive"} + +We now have a skull missing the structural T1 brain that is resized to match the dimensions of the EPI image. + +The final steps are to: + +1. Pull the first volume from the functional image +2. Place the functional image head into the hollow skull that we've created + +Since a functional image is 4-Dimensional, we'll need to pull the first volume to work with. This is because the structural image is 3-dimensional and operations will fail if we try to mix 3D and 4D data. + +```python +#Let's visualize the first volume of the functional image: +first_vol = ex_func.slicer[??, ??, ??, ??] +nplot.plot_epi(first_vol) +``` + +![](fig/exercise_fmri.png){alt='Episode 04 Exercise fMRI' class="img-responsive"} + +As shown in the figure above, the image has some "signal" outside of the brain. In order to place this within the now brainless head we made earlier, we need to mask out the functional MR data as well! + +```python +#Mask the first volume using ex_func_bm +masked_func = nimg.math_img('??', a=??, b=??) +nplot.plot_epi(masked_func) +``` + +![](fig/exercise_masked_fmri.png){alt='Episode 04 Exercise Masked fMRI' class="img-responsive"} + +The final step is to stick this data into the head of the T1 data. Since the hole in the T1 data is represented as $0$'s. We can add the two images together to place the functional data into the void: + +```python +#Now overlay the functional image on top of the anatomical +combined_img = nimg.math_img(??) +nplot.plot_anat(combined_img) +``` + +![](fig/exercise_complete.png){alt='Episode 04 Exercise Complete' class="img-responsive"} + +::::::::::::::: solution + +## Solution + +```python +# Invert the mask +invert_mask = nimg.math_img('1-a', a=ex_t1_bm) +nplot.plot_anat(invert_mask) + +# Apply the mask onto the T1 image +hollow_skull = nimg.math_img("a*b", a=ex_t1, b=invert_mask) +nplot.plot_anat(hollow_skull) + +#Resample the T1 to the size of the functional image! +resamp_skull = nimg.resample_to_img(source_img=hollow_skull, + target_img=ex_func, + interpolation='continuous') +nplot.plot_anat(resamp_skull) +print(resamp_skull.shape) +print(ex_func.shape) + +#Let's visualize the first volume of the functional image: +first_vol = ex_func.slicer[:,:,:,0] +nplot.plot_epi(first_vol) + +masked_func = nimg.math_img('a*b', a=first_vol, b=ex_func_bm) +nplot.plot_epi(masked_func) + +#Now overlay the functional image on top of the anatomical +combined_img = nimg.math_img('a+b', + a=resamp_skull, + b=masked_func) +nplot.plot_anat(combined_img) +``` + +::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::::::::::::: This doesn't actually achieve anything useful in practice. However it has hopefully served to get you more comfortable with the idea of resampling and performing manipulations on MRI data! In this section we explored functional MR imaging. Specifically we covered: - + 1. How the data in a fMRI scan is organized - with the additional dimension of timepoints 2. How we can integrate functional MR images to our structural image using resampling 3. How we can just as easily manipulate functional images using nilearn -Now that we've covered all the basics, it's time to start working on data processing using the tools that we've picked up. +Now that we've covered all the basics, it's time to start working on data processing using the tools that we've picked up. + + + +:::::::::::::::::::::::::::::::::::::::: keypoints + +- fMRI data is represented by spatial (x,y,z) and temporal (t) dimensions, totalling 4 dimensions +- fMRI data is at a lower resolution than structural data. To be able to combine data requires resampling your data + +:::::::::::::::::::::::::::::::::::::::::::::::::: + -{% include links.md %} diff --git a/episodes/05-data-cleaning-with-nilearn.md b/episodes/05-data-cleaning-with-nilearn.md index 7f687079..e33dd5cb 100644 --- a/episodes/05-data-cleaning-with-nilearn.md +++ b/episodes/05-data-cleaning-with-nilearn.md @@ -1,17 +1,23 @@ --- -title: "Cleaning Confounders in your Data with Nilearn" +title: Cleaning Confounders in your Data with Nilearn teaching: 30 exercises: 0 -questions: -- "How can we clean the data so that it more closely reflects BOLD instead of artifacts" -objectives: -- "Understand the motivation behind confound/nuisance regression" -- "Learn how to implement cleaning using nilearn and fmriprep" -keypoints: -- "Nuisance regression is an attempt to make sure your results aren't driven by non-brain signals" -- "With resting state, we don't actually ever know the true signal - we can only attempt to estimate it" --- -# Introduction + +::::::::::::::::::::::::::::::::::::::: objectives + +- Understand the motivation behind confound/nuisance regression +- Learn how to implement cleaning using nilearn and fmriprep + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::: questions + +- How can we clean the data so that it more closely reflects BOLD instead of artifacts + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +## Introduction **Movement is the enemy of Neuroimagers** @@ -25,11 +31,11 @@ This process of removing motion-based artifacts from our data is called **confou In this section we'll implement confound regression for resting-state data using nilearn's high-level functionality. -# Setting up +## Setting up Let's load in some modules as we've done before -~~~ +```python import os from nilearn import image as nimg from nilearn import plotting as nplot @@ -37,35 +43,32 @@ import matplotlib.pyplot as plt import numpy as np import nibabel as nib %matplotlib inline -~~~ -{: .language-python} +``` -## Setting up our Motion Estimates +### Setting up our Motion Estimates The beauty of FMRIPREP is that it estimates a number of motion-related signals for you and outputs it into: -**sub-xxxx_task-xxxx_desc-confounds_timeseries.tsv** +**sub-xxxx\_task-xxxx\_desc-confounds\_timeseries.tsv** This is basically a spreadsheet that has columns related to each motion estimate type and rows for timepoints. We can view these using a language-python package called `pandas`. -~~~ +```python import pandas as pd -~~~ -{: .language-python} +``` Let's pick an fMRI file to clean and pull out the confound tsv that FMRIPREP computed for us: -~~~ +```python import bids # assuming pip install pybids was covered earlier sub = '10788' fmriprep_dir = '../data/ds000030/derivatives/fmriprep/' layout = bids.BIDSLayout(fmriprep_dir,validate=False, config=['bids','derivatives']) -~~~ -{: .language-python} +``` -~~~ +```python func_files = layout.get(subject=sub, datatype='func', task='rest', desc='preproc', @@ -86,209 +89,218 @@ confound_files = layout.get(subject=sub, desc='confounds', extension="tsv", return_type='file') -~~~ -{: .language-python} +``` -~~~ +```python func_file = func_files[0] mask_file = mask_files[0] confound_file = confound_files[0] -~~~ -{: .language-python} +``` Using `pandas` we can read in the confounds.tsv file as a spreadsheet and display some rows: -~~~ +```python #Delimiter is \t --> tsv is a tab-separated spreadsheet confound_df = pd.read_csv(confound_file, delimiter='\t') print(f"Read {len(confound_df.columns)} confounder items, each {len(confound_df)} TRs.") print(confound_df.head) -~~~ -{: .language-python} +``` -Each column in this DataFrame confound_df represents a specific confound variable that is either estimated directly from head motion during the functional scan or other noise characteristics that may capture noise (non grey-matter signal for example). Each row represents values from a TR/sample. So the number of rows in your confound_df should match the number of TRs you have in the functional MR data. +Each column in this DataFrame confound\_df represents a specific confound variable that is either estimated directly from head motion during the functional scan or other noise characteristics that may capture noise (non grey-matter signal for example). Each row represents values from a TR/sample. So the number of rows in your confound\_df should match the number of TRs you have in the functional MR data. -> ## Picking your Confounds -> The choice of which confounds to use in functional imaging analysis is a source of large debate. We recommend that you check out these sources for a start: -> 1. https://www.sciencedirect.com/science/article/pii/S1053811917302288#f0005 -> 2. https://www.sciencedirect.com/science/article/pii/S1053811917302288 -> For now we're going to replicate the pre-processing (mostly) from the seminal Yeo1000 17-networks paper: -> https://www.ncbi.nlm.nih.gov/pubmed/21653723 -{: .callout} +::::::::::::::::::::::::::::::::::::::::: callout ---- +### Picking your Confounds + +The choice of which confounds to use in functional imaging analysis is a source of large debate. We recommend that you check out these sources for a start: -### The Yeo 2011 Pre-processing schema +1. [https://www.sciencedirect.com/science/article/pii/S1053811917302288#f0005](https://www.sciencedirect.com/science/article/pii/S1053811917302288#f0005) +2. [https://www.sciencedirect.com/science/article/pii/S1053811917302288](https://www.sciencedirect.com/science/article/pii/S1053811917302288) + For now we're going to replicate the pre-processing (mostly) from the seminal Yeo1000 17-networks paper: + [https://www.ncbi.nlm.nih.gov/pubmed/21653723](https://www.ncbi.nlm.nih.gov/pubmed/21653723) + -#### Confound regressors -1. 6 motion parameters (trans_x, trans_y, trans_z, rot_x, rot_y, rot_z) -2. Global signal (global_signal) +:::::::::::::::::::::::::::::::::::::::::::::::::: + +*** + +#### The Yeo 2011 Pre-processing schema + +##### Confound regressors + +1. 6 motion parameters (trans\_x, trans\_y, trans\_z, rot\_x, rot\_y, rot\_z) +2. Global signal (global\_signal) 3. Cerebral spinal fluid signal (csf) -4. White matter signal (white_matter) +4. White matter signal (white\_matter) This is a total of 9 base confound regressor variables. Finally we add temporal derivatives of each of these signals as well (1 temporal derivative for each), the result is 18 confound regressors. *** -**Temporal Derivatives** are the changes in values across 2 consecutive samples. It represents change in signal over time. For example, when dealing with the confound variable "X", which represents motion along the "trans_x" direction, the temporal derivative represents *velocity in the X direction*. + +**Temporal Derivatives** are the changes in values across 2 consecutive samples. It represents change in signal over time. For example, when dealing with the confound variable "X", which represents motion along the "trans\_x" direction, the temporal derivative represents *velocity in the X direction*. *** -#### Low/High pass filtering -1. Low pass filtering cutoff: 0.08 +##### Low/High pass filtering + +1. Low pass filtering cutoff: 0.08 2. High pass filtering cutoff: 0.009 -Low pass filters out high frequency signals from our data. fMRI signals are slow evolving processes, any high frequency signals are likely due to noise +Low pass filters out high frequency signals from our data. fMRI signals are slow evolving processes, any high frequency signals are likely due to noise High pass filters out any very low frequency signals (below 0.009Hz), which may be due to intrinsic scanner instabilities -#### Drop dummy TRs -During the initial stages of a functional scan there is a strong signal decay artifact, thus the first 4ish or so -TRs are very high intensity signals that don't reflect the rest of the scan. Therefore we drop these timepoints. +##### Drop dummy TRs -#### Censoring + Interpolation (leaving out) -Censoring involves removal and interpolation of high-movement frames from the fMRI data. Interpolation is typically done using sophisticated algorithms much like [Power et al. 2014](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3849338/). +During the initial stages of a functional scan there is a strong signal decay artifact, thus the first 4ish or so +TRs are very high intensity signals that don't reflect the rest of the scan. Therefore we drop these timepoints. -**We won't be using censoring + interpolation since its fairly complicated and would take up too much time** +##### Censoring + Interpolation (leaving out) ---- +Censoring involves removal and interpolation of high-movement frames from the fMRI data. Interpolation is typically done using sophisticated algorithms much like [Power et al. 2014](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3849338/). +**We won't be using censoring + interpolation since its fairly complicated and would take up too much time** +*** -### Setting up Confound variables for regression +#### Setting up Confound variables for regression -#### Computing temporal derivatives for confound variables +##### Computing temporal derivatives for confound variables First we'll select our confound variables from our dataframe. You can do this by specifying a list of confounds, then using that list to pull out the associated columns -~~~ + +```python # Select confounds confound_vars = ['trans_x','trans_y','trans_z', 'rot_x','rot_y','rot_z', 'global_signal', 'csf', 'white_matter'] -~~~ -{: .language-python} +``` Next we need to get derivatives for each of these columns. Luckily fMRIPrep provides this for us. Derivative columns are denoted as the following: -- {NAME_OF_COLUMN}_derivative1 +- {NAME\_OF\_COLUMN}\_derivative1 -Since typing is alot of work, we'll use a for-loop instead to pick the derivatives for our confound_vars: -~~~ +Since typing is alot of work, we'll use a for-loop instead to pick the derivatives for our confound\_vars: + +```python # Get derivative column names derivative_columns = ['{}_derivative1'.format(c) for c in confound_vars] print(derivative_columns) -~~~ -{: .language-python} +``` Now we'll join these two lists together: -~~~ + +```python final_confounds = confound_vars + derivative_columns print(final_confounds) -~~~ -{: .language-python} +``` Finally we'll use this list to pick columns from our confounds table -~~~ +```python confound_df = confound_df[final_confounds] confound_df.head() -~~~ -{: .language-python} +``` + +::::::::::::::::::::::::::::::::::::::::: callout + +### What the NaN??? + +As you might have noticed, we have NaN's in our confound dataframe. This happens because there is no prior value to the first index to take a difference with, but this isn't a problem since we're going to be dropping 4 timepoints from our data and confounders anyway! + -> ## What the NaN??? -> As you might have noticed, we have NaN's in our confound dataframe. This happens because there is no prior value to the first index to take a difference with, but this isn't a problem since we're going to be dropping 4 timepoints from our data and confounders anyway! -{: .callout} +:::::::::::::::::::::::::::::::::::::::::::::::::: + +##### Dummy TR Drop -#### Dummy TR Drop Now we'll implement our **Dummy TR Drop**. Remember this means that we are removing the first 4 timepoints from our functional image (we'll also have to do this for our first 4 confound timepoints!): -~~~ +```python #First we'll load in our data and check the shape raw_func_img = nimg.load_img(func_file) raw_func_img.shape -~~~ -{: .language-python} +``` Recall that the fourth dimension represents frames/TRs(timepoints). We want to drop the first four timepoints entirely, to do so we use nibabel's slicer feature. We'll also drop the first 4 confound variable timepoints to match the functional scan -~~~ +```python func_img = raw_func_img.slicer[:,:,:,4:] func_img.shape -~~~ -{: .language-python} +``` -~~~ +```output (65, 77, 49, 148) -~~~ -{: .output} +``` -~~~ +```python #Drop confound dummy TRs drop_confound_df = confound_df.loc[4:] print(drop_confound_df.shape) #number of rows should match that of the functional image drop_confound_df.head() -~~~ -{: .language-python} - +``` -### Applying confound regression +#### Applying confound regression Now we'd like to clean our data of our selected confound variables. There are two ways to go about this: -1. If you have nilearn version 0.5.0 or higher use nilearn.image.clean_img(image,confounds,...) +1. If you have nilearn version 0.5.0 or higher use nilearn.image.clean\_img(image,confounds,...) 2. If you want full control over specific parts of the image you're cleaning use nilearn.signal.clean(signals,confounds,...) The first method is probably most practical and can be done in one line given what we've already set-up. However, in cases of very large datasets (HCP-style), the second method might be preferable for optimizing memory usage. First note that both methods take an argument confounds. This is a matrix: -![Confounds Matrix](../fig/matrix.png){:class="img-responsive"} +![](fig/matrix.png){alt='Confounds Matrix' class="img-responsive"} Therefore our goal is to take our confound matrix and work it into a matrix of the form above. The end goal is a matrix with 147 rows, and columns matching the number of confound variables (9x2=18) Luckily this is a one-liner! -~~~ +```python confounds_matrix = drop_confound_df.values #Confirm matrix size is correct confounds_matrix.shape -~~~ -{: .language-python} +``` Let's clean our image! -## Using nilearn.image.clean_img +### Using nilearn.image.clean\_img First we'll describe a couple of this function's important arguments. Any argument enclosed in [arg] is optional -nilearn.image.clean_img(image,confounds,[low_pass],[high_pass],[t_r],[mask_img],[detrend],[standardize]) +nilearn.image.clean\_img(image,confounds,[low\_pass],[high\_pass],[t\_r],[mask\_img],[detrend],[standardize]) **Required**: -- image: The functional image ( func_img ) -- confounds: The confound matrix ( confounds ) + +- image: The functional image ( func\_img ) +- confounds: The confound matrix ( confounds ) **Optional**: -- low_pass: A low pass filter cut-off -- high_pass A high pass filter cut-off -- t_r: This is required if using low/high pass, the repetition time of acquisition (imaging metadata) -- mask_img Apply a mask when performing confound regression, will speed up regression + +- low\_pass: A low pass filter cut-off +- high\_pass A high pass filter cut-off +- t\_r: This is required if using low/high pass, the repetition time of acquisition (imaging metadata) +- mask\_img Apply a mask when performing confound regression, will speed up regression - detrend: Remove drift from the data (useful for removing scanner instability artifacts) [default=True] - standardize: Set mean to 0, and variance to 1 --> sets up data for statistical analysis [default=True] -*** -**What we're using**: -The Repetition Time of our data is 2 seconds, in addition since we're replicating (mostly) Yeo 2011's analysis: -- high_pass = 0.009 -- low_pass = 0.08 +*** + +**What we're using**: + +The Repetition Time of our data is 2 seconds, in addition since we're replicating (mostly) Yeo 2011's analysis: + +- high\_pass = 0.009 +- low\_pass = 0.08 - detrend = True - standardize = True -In addition we'll use a mask of our MNI transformed functional image ( mask ) to speed up cleaning +In addition we'll use a mask of our MNI transformed functional image ( mask ) to speed up cleaning - -~~~ +```python #Set some constants high_pass= 0.009 low_pass = 0.08 @@ -302,23 +314,31 @@ clean_img = nimg.clean_img(func_img,confounds=confounds_matrix,detrend=True,stan nplot.plot_epi(clean_img.slicer[:,:,:,50]) before_figure = nplot.plot_carpet(clean_img, mask, t_r=t_r) after_figure = nplot.plot_carpet(func_img, mask, t_r=t_r) -~~~ -{: .language-python} +``` -### Done! +#### Done! -Hopefully by now you've learned what confound regression is, and how to perform it in nilearn using 2 different methods. We'd like to note that there are many more methods to perform confound regression (simultaneous signal extraction + confound regression for example) but all those methods fundamentally rely on what you've done here. +Hopefully by now you've learned what confound regression is, and how to perform it in nilearn using 2 different methods. We'd like to note that there are many more methods to perform confound regression (simultaneous signal extraction + confound regression for example) but all those methods fundamentally rely on what you've done here. -In addition, performing confound regression on *functional volumes*, is also not the only way to do data cleaning. More modern methods involve applying confound regression on *functional surfaces*, however, those methods are too advanced for an introductory course to functional data analysis and involve tools outside of python. +In addition, performing confound regression on *functional volumes*, is also not the only way to do data cleaning. More modern methods involve applying confound regression on *functional surfaces*, however, those methods are too advanced for an introductory course to functional data analysis and involve tools outside of python. If you're interested in surface-based analysis we recommend that you check out the following sources: -1. https://edickie.github.io/ciftify/#/ -2. https://www.humanconnectome.org/software/connectome-workbench +1. [https://edickie.github.io/ciftify/#/](https://edickie.github.io/ciftify/#/) +2. [https://www.humanconnectome.org/software/connectome-workbench](https://www.humanconnectome.org/software/connectome-workbench) 3. [The minimal preprocessing pipelines for the Human Connectome Project](https://www.ncbi.nlm.nih.gov/pubmed/23668970) *** -The section below is **optional** and is a more advanced dive into the underlying mechanics of how nilearn.clean_img works: +The section below is **optional** and is a more advanced dive into the underlying mechanics of how nilearn.clean\_img works: + + + +:::::::::::::::::::::::::::::::::::::::: keypoints + +- Nuisance regression is an attempt to make sure your results aren't driven by non-brain signals +- With resting state, we don't actually ever know the true signal - we can only attempt to estimate it + +:::::::::::::::::::::::::::::::::::::::::::::::::: + -{% include links.md %} diff --git a/episodes/06-apply-a-parcellation.md b/episodes/06-apply-a-parcellation.md index 9d82f0dd..1a0977df 100644 --- a/episodes/06-apply-a-parcellation.md +++ b/episodes/06-apply-a-parcellation.md @@ -1,22 +1,27 @@ --- -title: "Applying Parcellations to Resting State Data" +title: Applying Parcellations to Resting State Data teaching: 25 exercises: 15 -questions: -- "How can we reduce amount of noise-related variance in our data?" -- "How can we frame our data as a set of meaningful features?" -objectives: -- "Learn about the utility of parcellations as a data dimensionalty reduction tool" -- "Understand what the tradeoffs are when using parcellations to analyze your data" -keypoints: -- "Parcellations group voxels based on criteria such as similarities, orthogonality or some other criteria" -- "Nilearn stores several standard parcellations that can be applied to your data" -- "Parcellations are defined by assigning each voxel a parcel 'membership' value telling you which group the parcel belongs to" -- "Parcellations provide an interpretative framework for understanding resting state data. But beware, some of the techniques used to form parcellations may not represent actual brain functional units!" --- -# Introduction +::::::::::::::::::::::::::::::::::::::: objectives + +- Learn about the utility of parcellations as a data dimensionalty reduction tool +- Understand what the tradeoffs are when using parcellations to analyze your data + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::: questions + +- How can we reduce amount of noise-related variance in our data? +- How can we frame our data as a set of meaningful features? + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +## Introduction + ## What is a Brain Atlas or Parcellation? + A brain atlas/parcellation is a voxel-based labelling of your data into "structural or functional units". In a parcellation schema each voxel is assigned a numeric (integer) label corresponding to the structural/functional unit that the particular voxel is thought to belong to based on some criteria. You might wonder why someone would simply *average together a bunch of voxels* in a way that would reduce the richness of the data. This boils down to a few problems inherit to functional brain imaging: 1. Resting state data is noisy, averaging groups of "similar" voxels reduces the effect of random noise effects @@ -24,50 +29,46 @@ A brain atlas/parcellation is a voxel-based labelling of your data into "structu 3. Limit the number of statistical tests thereby reducing potential Type I errors without resorting to strong statistical correction techniques that might reduce statistical power. 4. A simpler way to visualize your data, instead of 40x40x40=6400 data points, you might have 17 or up to 200; this is still significantly less data to deal with! - ## Applying a Parcellation to your Data -Since the parcellation of a brain is defined (currently) by spatial locations, application of an parcellation to fMRI data only concerns the first 3 dimensions; the last dimension (time) is retained. Thus a parcellation assigns every voxel (x,y,z) to a particular parcel ID (an integer). +Since the parcellation of a brain is defined (currently) by spatial locations, application of an parcellation to fMRI data only concerns the first 3 dimensions; the last dimension (time) is retained. Thus a parcellation assigns every voxel (x,y,z) to a particular parcel ID (an integer). -Nilearn supports a large selection of different atlases that can be found [here](http://nilearn.github.io/modules/reference.html#module-nilearn.datasets). For information about how to select which parcellation to use for analysis of your data we refer you to Arslan et al. 2018. +Nilearn supports a large selection of different atlases that can be found [here](https://nilearn.github.io/modules/reference.html#module-nilearn.datasets). For information about how to select which parcellation to use for analysis of your data we refer you to Arslan et al. 2018. ### Retrieving the Atlas + For this tutorial we'll be using a set of parcellation from [Yeo et al. 2011](link). This atlas was generated from fMRI data from 1000 healthy control participants. First we'll load in our packages as usual: -~~~ +```python from nilearn import datasets from nilearn import image as nimg from nilearn import plotting as nplot %matplotlib inline -~~~ -{: .language-python} +``` To retrieve the Yeo atlas we'll use the `fetch_atlas_*` family of functions provided for by nilearn.datasets and download it into a local directory: -~~~ +```python parcel_dir = '../resources/rois/' atlas_yeo_2011 = datasets.fetch_atlas_yeo_2011(parcel_dir) -~~~ -{: .language-python} +``` The method `datasets.fetch_atlas_yeo_2011()` returns a `dict` object. Examining the keys of the dictionary yields the following: -~~~ +```python atlas_yeo_2011.keys() -~~~ -{: .language-python} +``` -~~~ +```python output -~~~ -{: .language-python} +``` Each of the values associated with a key in `atlas_yeo_2011` is a `.nii.gz` image which contains a 3D NIFTI volume with a label for a given (x,y,z) voxel. Since these images are 3D volumes (sort of like structural images), we can view them using nilearn's plotting utilities: -~~~ +```python #Define where to slice the image cut_coords(8, -4, 9) #Show a colorbar @@ -80,19 +81,18 @@ nplot.plot_roi(atlas_yeo_2011['thin_7'], cut_coords=cut_coords, colorbar=colorba nplot.plot_roi(atlas_yeo_2011['thin_17'], cut_coords=cut_coords, colorbar=colorbar, cmap=cmap, title='thin_17') nplot.plot_roi(atlas_yeo_2011['thick_7'], cut_coords=cut_coords, colorbar=colorbar, cmap=cmap, title='thick_7') nplot.plot_roi(atlas_yeo_2011['thick_17'], cut_coords=cut_coords, colorbar=colorbar, cmap=cmap, title='thick_17') -~~~ -{: .language-python} +``` -![Yeo Thin 7](../fig/thin_7.png){:class="img-responsive"} -![Yeo Thin 17](../fig/thin_17.png){:class="img-responsive"} -![Yeo Thick 7](../fig/thick_7.png){:class="img-responsive"} -![Yeo Thick 17](../fig/thick_17.png){:class="img-responsive"} +![](fig/thin_7.png){alt='Yeo Thin 7' class="img-responsive"} +![](fig/thin_17.png){alt='Yeo Thin 17' class="img-responsive"} +![](fig/thick_7.png){alt='Yeo Thick 7' class="img-responsive"} +![](fig/thick_17.png){alt='Yeo Thick 17' class="img-responsive"} You'll notice that the colour bar on the right shows the number of labels in each atlas and which colour corresponds to which network The 7 and 17 network parcellations correspond to the two most stable clustering solutions from the algorithm used by the authors. The thin/thick designation refer to how strict the voxel inclusion is (thick might include white matter/CSF, thin might exclude some regions of grey matter due to partial voluming effects). -For simplicity we'll use the thick_7 variation which includes the following networks: +For simplicity we'll use the thick\_7 variation which includes the following networks: 1. Visual 2. Somatosensory @@ -104,15 +104,15 @@ For simplicity we'll use the thick_7 variation which includes the following netw The parcel areas labelled with 0 are background voxels not associated with a particular network. -~~~ +```python atlas_yeo = atlas_yeo_2011['thick_7'] -~~~ -{: .language-python} +``` ### Spatial Separation of Network + A key feature of the Yeo2011 networks is that they are *spatially distributed*, meaning that the locations of two voxels in the same network need not be part of the same region. However, there could be some cases in which you might want to examine voxels belonging to a network within a particular region. To do this, we can separate parcels belonging to the same network based on spatial continuity. If there is a gap between two sets of voxels belonging to the same parcel group, we can assign new labels to separate them out. Nilearn has a feature to handle this: -~~~ +```python from nilearn.regions import connected_label_regions region_labels = connected_label_regions(atlas_yeo) nplot.plot_roi(region_labels, @@ -121,121 +121,134 @@ nplot.plot_roi(region_labels, colorbar=True, cmap='Paired', title='Relabeled Yeo Atlas') -~~~ -{: .language-python} +``` -![Separated Yeo Labels](../fig/yeo_sep.png){:class="img-responsive"} +![](fig/yeo_sep.png){alt='Separated Yeo Labels' class="img-responsive"} ### Resampling the Atlas Let's store the separated version of the atlas into a NIFTI file so that we can work with it later: -~~~ +```python region_labels.to_filename('../resources/rois/yeo_2011/Yeo_JNeurophysiol11_MNI152/relabeled_yeo_atlas.nii.gz') -~~~ -{: .language-python} - -> ## Resampling Exercise -> -> Our goal is to match the parcellation atlas dimensions to our functional file so that we can use it to extract the mean time series of each parcel region. Using `Nilearn`'s resampling capabilities match the dimensions of the atlas file to the functional file -> First let's pick our functional file. Atlases are typically defined in standard space so we will use the MNI152NLin2009cAsym version of the functional file: -> -> ~~~ -> func_file = '../data/ds000030/derivatives/fmriprep/sub-10788/func/sub-10788_task-rest_bold_space-MNI152NLin2009cAsym_preproc.nii.gz' -> func_img = nib.load(func_file) -> ~~~ -> {: .language-python} -> -> First examine the size of both files, if they match we are done: -> ~~~ -> print('Size of functional image:', func_img.shape) -> print('Size of atlas image:', ??) -> ~~~ -> {: .language-python} -> -> Looks like they don't match. To resolve this, we can use nimg.resample_to_img to resize the *atlas image* to match that of the *functional image*. Think about what kind of interpolation we'd like to use. Recall that the atlas contains integer values (i.e 0, 1, 2, 3,...), we *do not want any in-between values!* -> -> ~~~ -> resampled_yeo = nimg.resample_to_img(??, ??, interpolation = '??') -> ~~~ -> {: .language-python} -> -> -> > ## Solution -> > -> > ~~~ -> > # Print dimensions of functional image and atlas image -> > -> > print("Size of functional image:", func_img.shape) -> > print("Size of atlas image:", region_labels.shape) -> > -> > resampled_yeo = nimg.resample_to_img(region_labels, func_img, interpolation = 'nearest') -> > ~~~ -> > {: .language-python} -> {: .solution} -{: .challenge} +``` + +::::::::::::::::::::::::::::::::::::::: challenge + +## Resampling Exercise + +Our goal is to match the parcellation atlas dimensions to our functional file so that we can use it to extract the mean time series of each parcel region. Using `Nilearn`'s resampling capabilities match the dimensions of the atlas file to the functional file +First let's pick our functional file. Atlases are typically defined in standard space so we will use the MNI152NLin2009cAsym version of the functional file: + +```python +func_file = '../data/ds000030/derivatives/fmriprep/sub-10788/func/sub-10788_task-rest_bold_space-MNI152NLin2009cAsym_preproc.nii.gz' +func_img = nib.load(func_file) +``` + +First examine the size of both files, if they match we are done: + +```python +print('Size of functional image:', func_img.shape) +print('Size of atlas image:', ??) +``` + +Looks like they don't match. To resolve this, we can use nimg.resample\_to\_img to resize the *atlas image* to match that of the *functional image*. Think about what kind of interpolation we'd like to use. Recall that the atlas contains integer values (i.e 0, 1, 2, 3,...), we *do not want any in-between values!* + +```python +resampled_yeo = nimg.resample_to_img(??, ??, interpolation = '??') +``` + +::::::::::::::: solution + +## Solution + +```python +# Print dimensions of functional image and atlas image + +print("Size of functional image:", func_img.shape) +print("Size of atlas image:", region_labels.shape) + +resampled_yeo = nimg.resample_to_img(region_labels, func_img, interpolation = 'nearest') +``` + +::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::::::::::::: Let's see what the resampled atlas looks like overlayed on a slice of our NifTI file -~~~ +```python # Note that we're pulling a random timepoint from the fMRI data nplot.plot_roi(resampled_yeo, func_img.slicer[:, :, :, 54]) -~~~ -{: .language-python} - -![Episode 06 Exercise Resampled Yeo Labels](../fig/resampled_yeo.png){:class="img-responsive"} +``` +![](fig/resampled_yeo.png){alt='Episode 06 Exercise Resampled Yeo Labels' class="img-responsive"} ## Visualizing ROIs For the next section, we'll be performing an analysis using the Yeo parcellation on our functional data. Specifically, we'll be using two ROIs: 44 and 46. > ## Exercise + Visualize ROIs 44 and 46 in the Yeo atlas. We'll be looking at these 2 ROIs in more detail during our analysis -> -> ~~~ -> roi = 44 -> -> # Make a mask for ROI 44 -> roi_mask_44 = nimg.math_img('a == ??', a=resampled_yeo) -> -> # Visualize ROI -> nplot.plot_roi(??) -> ~~~ -> {: .language-python} -> -> ![Episode 06 Exercise Yeo ROI 44](../fig/roi_44.png){:class="img-responsive"} -> -> ~~~ -> roi = 46 -> -> # Make a mask for ROI 44 -> roi_mask_46 = nimg.math_img(??) -> -> # Visualize ROI -> ?? -> ~~~ -> {: .language-python} -> -> ![Episode 06 Exercise Yeo ROI 46](../fig/roi_46.png){:class="img-responsive"} -> -> > ## Solution -> > -> > ~~~ -> > # Make a mask for ROI 44 -> > roi_mask = nimg.math_img('a == 44', a=resampled_yeo) -> > -> > # Visualize ROI -> > nplot.plot_roi(masked_resamp_yeo) -> > -> > # Make a mask for ROI 44 -> > roi_mask = nimg.math_img('a == 46', a=resampled_yeo) -> > -> > # Visualize ROI -> > nplot.plot_roi(masked_resamp_yeo) -> > ~~~ -> > {: .language-python} -> {: .solution} -{: .challenge} - -{% include links.md %} + +::::::::::::::::::::::::::::::::::::::: challenge + +```python +roi = 44 + +# Make a mask for ROI 44 +roi_mask_44 = nimg.math_img('a == ??', a=resampled_yeo) + +# Visualize ROI +nplot.plot_roi(??) +``` + +![](fig/roi_44.png){alt='Episode 06 Exercise Yeo ROI 44' class="img-responsive"} + +```python +roi = 46 + +# Make a mask for ROI 44 +roi_mask_46 = nimg.math_img(??) + +# Visualize ROI +?? +``` + +![](fig/roi_46.png){alt='Episode 06 Exercise Yeo ROI 46' class="img-responsive"} + +::::::::::::::: solution + +## Solution + +```python +# Make a mask for ROI 44 +roi_mask = nimg.math_img('a == 44', a=resampled_yeo) + +# Visualize ROI +nplot.plot_roi(masked_resamp_yeo) + +# Make a mask for ROI 44 +roi_mask = nimg.math_img('a == 46', a=resampled_yeo) + +# Visualize ROI +nplot.plot_roi(masked_resamp_yeo) +``` + +::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::::::::::::: + + + +:::::::::::::::::::::::::::::::::::::::: keypoints + +- Parcellations group voxels based on criteria such as similarities, orthogonality or some other criteria +- Nilearn stores several standard parcellations that can be applied to your data +- Parcellations are defined by assigning each voxel a parcel 'membership' value telling you which group the parcel belongs to +- Parcellations provide an interpretative framework for understanding resting state data. But beware, some of the techniques used to form parcellations may not represent actual brain functional units! + +:::::::::::::::::::::::::::::::::::::::::::::::::: + + diff --git a/episodes/07-functional-connectivity-analysis.md b/episodes/07-functional-connectivity-analysis.md index a50b0ea3..ef514cbb 100644 --- a/episodes/07-functional-connectivity-analysis.md +++ b/episodes/07-functional-connectivity-analysis.md @@ -1,17 +1,22 @@ --- -title: "Functional Connectivity Analysis" +title: Functional Connectivity Analysis teaching: 45 exercises: 15 -questions: -- "How can we estimate brain functional connectivity patterns from resting state data?" -objectives: -- "Use parcellations to reduce fMRI noise and speed up computation of functional connectivity" -keypoints: -- "MR images are essentially 3D arrays where each voxel is represented by an (i,j,k) index" -- "Nilearn is Nibabel under the hood, knowing how Nibabel works is key to understanding Nilearn" --- -# Introduction +::::::::::::::::::::::::::::::::::::::: objectives + +- Use parcellations to reduce fMRI noise and speed up computation of functional connectivity + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::: questions + +- How can we estimate brain functional connectivity patterns from resting state data? + +:::::::::::::::::::::::::::::::::::::::::::::::::: + +## Introduction Now we have an idea of three important components to analyzing neuroimaging data: @@ -30,26 +35,24 @@ ROI-based correlational analysis forms the basis of many more sophisticated kind Nilearn has a built in function for extracting timeseries from functional files and doing all the extra signal processing at the same time. Let's walk through how this is done First we'll grab our imports as usual -~~~ + +```python from nilearn import image as nimg from nilearn import plotting as nplot import numpy as np import pandas as pd from bids import BIDSLayout -~~~ -{: .language-python} - +``` Let's grab the data that we want to perform our connectivity analysis on using PyBIDS: -~~~ +```python #Use PyBIDS to parse BIDS data structure layout = BIDSLayout(fmriprep_dir, config=['bids','derivatives']) -~~~ -{: .language-python} +``` -~~~ +```python #Get resting state data (preprocessed, mask, and confounds file) func_files = layout.get(subject=sub, datatype='func', task='rest', @@ -72,31 +75,30 @@ confound_files = layout.get(subject=sub, desc='confounds', extension='tsv', return_type='file') -~~~ -{: .language-python} - +``` Now that we have a list of subjects to peform our analysis on, let's load up our parcellation template file -~~~ + +```python #Load separated parcellation parcel_file = '../resources/rois/yeo_2011/Yeo_JNeurophysiol11_MNI152/relabeled_yeo_atlas.nii.gz' yeo_7 = nimg.load_img(parcel_file) -~~~ -{: .language-python} +``` -Now we'll import a package from nilearn, called input_data which allows us to pull data using the parcellation file, and at the same time applying data cleaning! +Now we'll import a package from nilearn, called input\_data which allows us to pull data using the parcellation file, and at the same time applying data cleaning! -We first create an object using the parcellation file yeo_7 and our cleaning settings which are the following: +We first create an object using the parcellation file yeo\_7 and our cleaning settings which are the following: Settings to use: -- Confounds: trans_x, trans_y, trans_z, rot_x, rot_y, rot_z, white_matter, csf, global_signal + +- Confounds: trans\_x, trans\_y, trans\_z, rot\_x, rot\_y, rot\_z, white\_matter, csf, global\_signal - Temporal Derivatives: Yes -- high_pass = 0.009 -- low_pass = 0.08 +- high\_pass = 0.009 +- low\_pass = 0.08 - detrend = True - standardize = True -~~~ +```python from nilearn import input_data @@ -108,26 +110,25 @@ masker = input_data.NiftiLabelsMasker(labels_img=yeo_7, low_pass = 0.08, high_pass = 0.009, t_r=2) -~~~ -{: .language-python} +``` The object masker is now able to be used on *any functional image of the same size*. The `input_data.NiftiLabelsMasker` object is a wrapper that applies parcellation, cleaning and averaging to an functional image. For example let's apply this to our first subject: -~~~ +```python # Pull the first subject's data func_file = func_files[0] mask_file = mask_files[0] confound_file = confound_files[0] -~~~ -{: .language-python} +``` Before we go ahead and start using the masker that we've created, we have to do some preparatory steps. The following should be done prior to use the masker object: + 1. Make your confounds matrix (as we've done in Episode 06) 2. Drop Dummy TRs that are to be excluded from our cleaning, parcellation, and averaging step To help us with the first part, let's define a function to help extract our confound regressors from the .tsv file for us. Note that we've handled pulling the appropriate `{confounds}_derivative1` columns for you! You just need to supply the base regressors! -~~~ +```python #Refer to part_06 for code + explanation def extract_confounds(confound_tsv,confounds,dt=True): ''' @@ -154,12 +155,11 @@ def extract_confounds(confound_tsv,confounds,dt=True): #Return confound matrix return confound_mat -~~~ -{: .language-python} +``` Finally we'll set up our image file for confound regression (as we did in Episode 6). To do this we'll drop 4 TRs from *both our functional image and our confounds file*. Note that our masker object will not do this for us! -~~~ +```python #Load functional image tr_drop = 4 func_img = nimg.load_img(func_file) @@ -175,49 +175,47 @@ confounds = extract_confounds(confound_file, 'white_matter','csf']) #Drop the first 4 rows of the confounds matrix confounds = confounds[tr_drop:,:] -~~~ -{: .language-python} +``` ### Using the masker + Finally with everything set up, we can now use the masker to perform our: + 1. Confounds cleaning 2. Parcellation 3. Averaging within a parcel -All in one step! + All in one step! -~~~ +```python #Apply cleaning, parcellation and extraction to functional data cleaned_and_averaged_time_series = masker.fit_transform(func_img,confounds) cleaned_and_averaged_time_series.shape -~~~ -{: .language-python} +``` -~~~ +```output (147,43) -~~~ -{: .output} +``` Just to be clear, this data is *automatically parcellated for you, and, in addition, is cleaned using the confounds you've specified already!* -The result of running masker.fit_transform is a matrix that has: +The result of running masker.fit\_transform is a matrix that has: + - Rows matching the number of timepoints (148) - Columns, each for one of the ROIs that are extracted (43) **But wait!** -We originally had **50 ROIs**, what happened to 3 of them? It turns out that masker drops ROIs that are empty (i.e contain no brain voxels inside of them), this means that 3 of our atlas' parcels did not correspond to any region with signal! To see which ROIs are kept after computing a parcellation you can look at the labels_ property of masker: +We originally had **50 ROIs**, what happened to 3 of them? It turns out that masker drops ROIs that are empty (i.e contain no brain voxels inside of them), this means that 3 of our atlas' parcels did not correspond to any region with signal! To see which ROIs are kept after computing a parcellation you can look at the labels\_ property of masker: -~~~ +```python print(masker.labels_) print("Number of labels", len(masker.labels_)) -~~~ -{: .language-python} +``` -~~~ +```output [1, 2, 4, 5, 6, 7, 8, 10, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 33, 34, 35, 36, 37, 38, 40, 41, 42, 43, 44, 45, 46, 47, 49] Number of labels 43 -~~~ -{: .output} +``` This means that our ROIs of interest (44 and 46) *cannot be accessed using the 44th and 46th columns directly*! @@ -225,7 +223,7 @@ There are many strategies to deal with this weirdness. What we're going to do is First we'll identify all ROIs from the original atlas. We're going to use the numpy package which will provide us with functions to work with our image arrays: -~~~ +```python import numpy as np # Get the label numbers from the atlas @@ -234,19 +232,18 @@ atlas_labels = np.unique(yeo_7.get_fdata().astype(int)) # Get number of labels that we have NUM_LABELS = len(atlas_labels) print(NUM_LABELS) -~~~ -{: .language-python} +``` -~~~ +```output 50 -~~~ -{: .output} +``` Now we're going to create an array that contains: + - A number of rows matching the number of timepoints - A number of columns matching the total number of regions -~~~ +```python # Remember fMRI images are of size (x,y,z,t) # where t is the number of timepoints num_timepoints = func_img.shape[3] @@ -261,69 +258,67 @@ regions_kept = np.array(masker.labels_) final_signal[:, regions_kept] = cleaned_and_averaged_time_series print(final_signal.shape) -~~~ -{: .language-python} +``` It's a bit of work, but now we have an array where: 1. The column number matches the ROI label number -2. Any column that is lost during the masker.fit_transform is filled with 0 values! - +2. Any column that is lost during the masker.fit\_transform is filled with 0 values! -To get the columns corresponding to the regions that we've kept, we can simply use the regions_kept variable to select columns corresponding to the regions that weren't removed: +To get the columns corresponding to the regions that we've kept, we can simply use the regions\_kept variable to select columns corresponding to the regions that weren't removed: -~~~ +```python valid_regions_signal = final_signal[:, regions_kept] print(valid_regions_signal.shape) -~~~ -{: .language-python} +``` -This is identical to the original output of masker.fit_transform +This is identical to the original output of masker.fit\_transform -~~~ +```python np.array_equal( valid_regions_signal, cleaned_and_averaged_time_series ) -~~~ -{: .language-python} +``` This might seem unnecessary for now, but as you'll see in a bit, it'll come in handy when we deal with multiple subjects!d - ### Calculating Connectivity -In fMRI imaging, connectivity typically refers to the *correlation of the timeseries of 2 ROIs*. Therefore we can calculate a *full connectivity matrix* by computing the correlation between *all pairs of ROIs* in our parcellation scheme! +In fMRI imaging, connectivity typically refers to the *correlation of the timeseries of 2 ROIs*. Therefore we can calculate a *full connectivity matrix* by computing the correlation between *all pairs of ROIs* in our parcellation scheme! We'll use another nilearn tool called ConnectivityMeasure from nilearn.connectome. This tool will perform the full set of pairwise correlations for us -~~~ +```python from nilearn.connectome import ConnectivityMeasure -~~~ -{: .language-python} +``` Like the masker, we need to make an object that will calculate connectivity for us. -~~~ +```python correlation_measure = ConnectivityMeasure(kind='correlation') -~~~ -{: .language-python} +``` + +::::::::::::::::::::::::::::::::::::::::: callout + +Try using SHIFT-TAB to see what options you can put into the kind argument of ConnectivityMeasure -> Try using SHIFT-TAB to see what options you can put into the kind argument of ConnectivityMeasure -{: .callout} -Then we use correlation_measure.fit_transform() in order to calculate the full correlation matrix for our parcellated data! +:::::::::::::::::::::::::::::::::::::::::::::::::: +Then we use correlation\_measure.fit\_transform() in order to calculate the full correlation matrix for our parcellated data! -~~~ +```python full_correlation_matrix = correlation_measure.fit_transform([cleaned_and_averaged_time_series]) full_correlation_matrix.shape -~~~ -{: .language-python} +``` -> Note that we're using a list [final_signal], this is because correlation_measure works on a *list of subjects*. We'll take advantage of this later! -{: .callout} +::::::::::::::::::::::::::::::::::::::::: callout +Note that we're using a list [final\_signal], this is because correlation\_measure works on a *list of subjects*. We'll take advantage of this later! + + +:::::::::::::::::::::::::::::::::::::::::::::::::: The result is a matrix which has: @@ -337,187 +332,192 @@ Suppose we wanted to know the correlation between ROI 30 and ROI 40 - Then Row 30, Column 40 gives us this correlation. - Row 40, Column 40 can also give us this correlation - This is because the correlation of $A --> B = B --> A$ -> **NOTE**: Remember we were supposed to lose 7 regions from the masker.fit_transform step. The correlations for these regions will be 0! -{: .callout} +::::::::::::::::::::::::::::::::::::::::: callout + +**NOTE**: Remember we were supposed to lose 7 regions from the masker.fit\_transform step. The correlations for these regions will be 0! + + +:::::::::::::::::::::::::::::::::::::::::::::::::: Let's try pulling the correlation for ROI 44 and 46! -~~~ +```python full_correlation_matrix[0, 43, 45] -~~~ -{: .language-python} +``` Note that it'll be the same if we swap the rows and columns! -~~~ +```python full_correlation_matrix[0, 45, 43] -~~~ -{: .language-python} - - -> ## Exercise -> Apply the data extract process shown above to all subjects in our subject list and collect the results. Your job is to fill in the blanks! -> -> ~~~ -> # First we're going to create some empty lists to store all our data in! -> ctrl_subjects = [] -> schz_subjects = [] -> -> -> # We're going to keep track of each of our subjects labels here -> # pulled from masker.labels_ -> labels_list = [] -> -> # Get the number of unique labels in our parcellation -> # We'll use this to figure out how many columns to make (as we did earlier) -> atlas_labels = np.unique(yeo_7.get_fdata().astype(int)) -> NUM_LABELS = len(atlas_labels) -> -> # Set the list of confound variables we'll be using -> confound_variables = ['trans_x','trans_y','trans_z', -> 'rot_x','rot_y','rot_z', -> 'global_signal', -> 'white_matter','csf'] -> -> # Number of TRs we should drop -> TR_DROP=4 -> -> # Lets get all the subjects we have -> subjects = layout.get_subjects() -> for sub in subjects: -> -> #Get the functional file for the subject (MNI space) -> func_file = layout.get(subject=??, -> datatype='??', task='rest', -> desc='??', -> space='??' -> extension="nii.gz", -> return_type='file')[0] -> -> #Get the confounds file for the subject (MNI space) -> confound_file=layout.get(subject=??, datatype='??', -> task='rest', -> desc='??', -> extension='tsv', -> return_type='file')[0] -> -> #Load the functional file in -> func_img = nimg.load_img(??) -> -> #Drop the first 4 TRs -> func_img = func_img.slicer[??,??,??,??] -> -> -> #Extract the confound variables using the function -> confounds = extract_confounds(confound_file, -> confound_variables) -> -> #Drop the first 4 rows from the confound matrix -> confounds = confounds[??] -> -> # Make our array of zeros to fill out -> # Number of rows should match number of timepoints -> # Number of columns should match the total number of regions -> fill_array = np.zeros((func_img.shape[??], ??)) -> -> #Apply the parcellation + cleaning to our data -> #What function of masker is used to clean and average data? -> time_series = masker.fit_transform(??,??) -> -> # Get the regions that were kept for this scan -> regions_kept = np.array(masker.labels_) -> -> # Fill the array, this is what we'll use -> # to make sure that all our array are of the same size -> fill_array[:, ??] = time_series -> -> #If the subject ID starts with a "1" then they are control -> if sub.startswith('1'): -> ctrl_subjects.append(fill_array) -> #If the subject ID starts with a "5" then they are case (case of schizophrenia) -> if sub.startswith('5'): -> schz_subjects.append(fill_array) -> ~~~ -> {: .language-python} -> > ## Solution -> > -> > ~~~ -> > # First we're going to create some empty lists to store all our data in! -> > pooled_subjects = [] -> > ctrl_subjects = [] -> > schz_subjects = [] -> > -> > #Which confound variables should we use? -> > confound_variables = ['trans_x','trans_y','trans_z', -> > 'rot_x','rot_y','rot_z', -> > 'global_signal', -> > 'white_matter','csf'] -> > for sub in subjects: -> > -> > #Get the functional file for the subject (MNI space) -> > func_file = layout.get(subject=sub, -> > datatype='func', task='rest', -> > desc='preproc', -> > extension="nii.gz", -> > return_type='file')[0] -> > -> > #Get the confounds file for the subject (MNI space) -> > confound_file=layout.get(subject=sub, datatype='func', -> > task='rest', -> > desc='confounds', -> > extension='tsv', -> > return_type='file')[0] -> > -> > #Load the functional file in -> > func_img = nimg.load_img(func_file) -> > -> > #Drop the first 4 TRs -> > func_img = func_img.slicer[:,:,:,tr_drop:] -> > -> > -> > #Extract the confound variables using the function -> > confounds = extract_confounds(confound_file, -> > confound_variables) -> > -> > #Drop the first 4 rows from the confound matrix -> > #Which rows and columns should we keep? -> > confounds = confounds[tr_drop:,:] -> > -> > #Apply the parcellation + cleaning to our data -> > #What function of masker is used to clean and average data? -> > time_series = masker.fit_transform(func_img,confounds) -> > -> > #This collects a list of all subjects -> > pooled_subjects.append(time_series) -> > -> > #If the subject ID starts with a "1" then they are control -> > if sub.startswith('1'): -> > ctrl_subjects.append(time_series) -> > #If the subject ID starts with a "5" then they are case (case of schizophrenia) -> > if sub.startswith('5'): -> > schz_subjects.append(time_series) -> > ~~~ -> > {: .language-python} -> {: .solution} -{: .challenge} +``` + +::::::::::::::::::::::::::::::::::::::: challenge + +## Exercise + +Apply the data extract process shown above to all subjects in our subject list and collect the results. Your job is to fill in the blanks! + +```python +# First we're going to create some empty lists to store all our data in! +ctrl_subjects = [] +schz_subjects = [] + + +# We're going to keep track of each of our subjects labels here +# pulled from masker.labels_ +labels_list = [] + +# Get the number of unique labels in our parcellation +# We'll use this to figure out how many columns to make (as we did earlier) +atlas_labels = np.unique(yeo_7.get_fdata().astype(int)) +NUM_LABELS = len(atlas_labels) + +# Set the list of confound variables we'll be using +confound_variables = ['trans_x','trans_y','trans_z', + 'rot_x','rot_y','rot_z', + 'global_signal', + 'white_matter','csf'] + +# Number of TRs we should drop +TR_DROP=4 + +# Lets get all the subjects we have +subjects = layout.get_subjects() +for sub in subjects: + + #Get the functional file for the subject (MNI space) + func_file = layout.get(subject=??, + datatype='??', task='rest', + desc='??', + space='??' + extension="nii.gz", + return_type='file')[0] + + #Get the confounds file for the subject (MNI space) + confound_file=layout.get(subject=??, datatype='??', + task='rest', + desc='??', + extension='tsv', + return_type='file')[0] + + #Load the functional file in + func_img = nimg.load_img(??) + + #Drop the first 4 TRs + func_img = func_img.slicer[??,??,??,??] + + + #Extract the confound variables using the function + confounds = extract_confounds(confound_file, + confound_variables) + + #Drop the first 4 rows from the confound matrix + confounds = confounds[??] + + # Make our array of zeros to fill out + # Number of rows should match number of timepoints + # Number of columns should match the total number of regions + fill_array = np.zeros((func_img.shape[??], ??)) + + #Apply the parcellation + cleaning to our data + #What function of masker is used to clean and average data? + time_series = masker.fit_transform(??,??) + + # Get the regions that were kept for this scan + regions_kept = np.array(masker.labels_) + + # Fill the array, this is what we'll use + # to make sure that all our array are of the same size + fill_array[:, ??] = time_series + + #If the subject ID starts with a "1" then they are control + if sub.startswith('1'): + ctrl_subjects.append(fill_array) + #If the subject ID starts with a "5" then they are case (case of schizophrenia) + if sub.startswith('5'): + schz_subjects.append(fill_array) +``` + +::::::::::::::: solution + +## Solution + +```python +# First we're going to create some empty lists to store all our data in! +pooled_subjects = [] +ctrl_subjects = [] +schz_subjects = [] + +#Which confound variables should we use? +confound_variables = ['trans_x','trans_y','trans_z', + 'rot_x','rot_y','rot_z', + 'global_signal', + 'white_matter','csf'] +for sub in subjects: + + #Get the functional file for the subject (MNI space) + func_file = layout.get(subject=sub, + datatype='func', task='rest', + desc='preproc', + extension="nii.gz", + return_type='file')[0] + + #Get the confounds file for the subject (MNI space) + confound_file=layout.get(subject=sub, datatype='func', + task='rest', + desc='confounds', + extension='tsv', + return_type='file')[0] + + #Load the functional file in + func_img = nimg.load_img(func_file) + + #Drop the first 4 TRs + func_img = func_img.slicer[:,:,:,tr_drop:] + + + #Extract the confound variables using the function + confounds = extract_confounds(confound_file, + confound_variables) + + #Drop the first 4 rows from the confound matrix + #Which rows and columns should we keep? + confounds = confounds[tr_drop:,:] + + #Apply the parcellation + cleaning to our data + #What function of masker is used to clean and average data? + time_series = masker.fit_transform(func_img,confounds) + + #This collects a list of all subjects + pooled_subjects.append(time_series) + + #If the subject ID starts with a "1" then they are control + if sub.startswith('1'): + ctrl_subjects.append(time_series) + #If the subject ID starts with a "5" then they are case (case of schizophrenia) + if sub.startswith('5'): + schz_subjects.append(time_series) +``` + +::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::::::::::::::: The result of all of this code is that: 1. Subjects who start with a "1" in their ID, are controls, and are placed into the `ctrl_subjects` list 2. Subjects who start with a "2" in their ID, have schizophrenia, and are placed into the `schz_subjects` list -What's actually being placed into the list? The cleaned, parcellated time series data for each subject (the output of masker.fit_transform)! +What's actually being placed into the list? The cleaned, parcellated time series data for each subject (the output of masker.fit\_transform)! -A helpful trick is that we can re-use the correlation_measure object we made earlier and apply it to a *list of subject data*! +A helpful trick is that we can re-use the correlation\_measure object we made earlier and apply it to a *list of subject data*! -~~~ +```python ctrl_correlation_matrices = correlation_measure.fit_transform(ctrl_subjects) schz_correlation_matrices = correlation_measure.fit_transform(schz_subjects) -~~~ -{: .language-python} +``` At this point, we have correlation matrices for each subject across two populations. The final step is to examine the differences between these groups in their correlation between ROI 43 and ROI 45. @@ -527,72 +527,64 @@ An important step in any analysis is visualizing the data that we have. We've cl To visualize data we'll be using a python package called seaborn which will allow us to create statistical visualizations with not much effort. -~~~ +```python import seaborn as sns import matplotlib.pyplot as plt -~~~ -{: .language-python} +``` We can view a single subject's correlation matrix by using seaborn's heatmap function: -~~~ +```python sns.heatmap(ctrl_correlation_matrices[0], cmap='RdBu_r') -~~~ -{: .language-python} +``` -![Connectivity Matrix Heatmap](../fig/heatmap.png){:class="img-responsive"} +![](fig/heatmap.png){alt='Connectivity Matrix Heatmap' class="img-responsive"} -Recall that cleaning and parcellating the data causes some ROIs to get dropped. We dealt with this by filling an array of zeros (fill_array) only for columns where the regions are kept (regions_kept). This means that we'll have some correlation values that are 0! +Recall that cleaning and parcellating the data causes some ROIs to get dropped. We dealt with this by filling an array of zeros (fill\_array) only for columns where the regions are kept (regions\_kept). This means that we'll have some correlation values that are 0! This is more apparent if we plot the data slightly differently. For demonstrative purposes we've: + - Taken the absolute value of our correlations so that the 0's are the darkest color - Used a different color scheme -~~~ +```python sns.heatmap(np.abs(ctrl_correlation_matrices[0]), cmap='viridis') -~~~ -{: .language-python} +``` -![Absolute Matrix Heatmap Zeros Highlighted](../fig/heatmap_show_zeroes.png){:class="img-responsive"} +![](fig/heatmap_show_zeroes.png){alt='Absolute Matrix Heatmap Zeros Highlighted' class="img-responsive"} The dark lines in the correlation matrix correspond to regions that were dropped and therefore have 0 correlation! We can now pull our ROI 44 and 46 by indexing our list of correlation matrices as if it were a 3D array (kind of like an MR volume). Take a look at the shape: -~~~ +```python print(ctrl_correlation_matrices.shape) -~~~ -{: .language-python} +``` This is of form: -ctrl_correlation_matrices[subject_index, row_index, column_index] +ctrl\_correlation\_matrices[subject\_index, row\_index, column\_index] Now we're going to pull out just the correlation values between ROI 43 and 45 *across all our subjects*. This can be done using standard array indexing: - -~~~ +```python ctrl_roi_vec = ctrl_correlation_matrices[:,43,45] schz_roi_vec = schz_correlation_matrices[:,43,45] -~~~ -{: .language-python} +``` Next we're going to arrange this data into a table. We'll create two tables (called dataframes in the python package we're using, `pandas`) - -~~~ +```python #Create control dataframe ctrl_df = pd.DataFrame(data={'dmn_corr':ctrl_roi_vec, 'group':'control'}) ctrl_df.head() -~~~ -{: .language-python} +``` -~~~ +```python # Create the schizophrenia dataframe scz_df = pd.DataFrame(data={'dmn_corr':schz_roi_vec, 'group' : 'schizophrenia'}) scz_df.head() -~~~ -{: .language-python} +``` The result is: @@ -602,20 +594,17 @@ The result is: For visualization we're going to stack the two tables together... - -~~~ +```python #Stack the two dataframes together df = pd.concat([ctrl_df,scz_df], ignore_index=True) # Show some random samples from dataframe df.sample(n=3) -~~~ -{: .language-python} +``` Finally we're going to visualize the results using the python package `seaborn`! - -~~~ +```python #Visualize results # Create a figure canvas of equal width and height @@ -633,12 +622,9 @@ ax.set_title('DMN Intra-network Connectivity') ax.set_ylabel(r'Intra-network connectivity $\mu_\rho$') plt.show() -~~~ -{: .language-python} - - -Although the results here aren't significant they seem to indicate that there might be three subclasses in our schizophrenia group - of course we'd need *a lot* more data to confirm this! The interpretation of these results should ideally be based on some *a priori* hypothesis! +``` +Although the results here aren't significant they seem to indicate that there might be three subclasses in our schizophrenia group - of course we'd need *a lot* more data to confirm this! The interpretation of these results should ideally be based on some *a priori* hypothesis! ## Congratulations! @@ -650,4 +636,13 @@ Hopefully now you understand that: 4. Parcellating is a method of simplifying and "averaging" data. The type of parcellation reflect assumptions you make about the structure of your data 5. Functional Connectivity is really just time-series correlations between two signals! -{% include links.md %} + + +:::::::::::::::::::::::::::::::::::::::: keypoints + +- MR images are essentially 3D arrays where each voxel is represented by an (i,j,k) index +- Nilearn is Nibabel under the hood, knowing how Nibabel works is key to understanding Nilearn + +:::::::::::::::::::::::::::::::::::::::::::::::::: + + diff --git a/data/download_list b/episodes/data/download_list similarity index 100% rename from data/download_list rename to episodes/data/download_list diff --git a/data/setup_download b/episodes/data/setup_download similarity index 100% rename from data/setup_download rename to episodes/data/setup_download diff --git a/data/test_download b/episodes/data/test_download similarity index 100% rename from data/test_download rename to episodes/data/test_download diff --git a/fig/4D_array.png b/episodes/fig/4D_array.png similarity index 100% rename from fig/4D_array.png rename to episodes/fig/4D_array.png diff --git a/fig/animated_fmri_preproc.gif b/episodes/fig/animated_fmri_preproc.gif similarity index 100% rename from fig/animated_fmri_preproc.gif rename to episodes/fig/animated_fmri_preproc.gif diff --git a/fig/animated_slicing.gif b/episodes/fig/animated_slicing.gif similarity index 100% rename from fig/animated_slicing.gif rename to episodes/fig/animated_slicing.gif diff --git a/fig/animated_t1_mni.gif b/episodes/fig/animated_t1_mni.gif similarity index 100% rename from fig/animated_t1_mni.gif rename to episodes/fig/animated_t1_mni.gif diff --git a/fig/combined_img.png b/episodes/fig/combined_img.png similarity index 100% rename from fig/combined_img.png rename to episodes/fig/combined_img.png diff --git a/fig/console_filled.png b/episodes/fig/console_filled.png similarity index 100% rename from fig/console_filled.png rename to episodes/fig/console_filled.png diff --git a/fig/ctrl_r.png b/episodes/fig/ctrl_r.png similarity index 100% rename from fig/ctrl_r.png rename to episodes/fig/ctrl_r.png diff --git a/fig/dicom_to_nifti.png b/episodes/fig/dicom_to_nifti.png similarity index 100% rename from fig/dicom_to_nifti.png rename to episodes/fig/dicom_to_nifti.png diff --git a/fig/downsample_t1.png b/episodes/fig/downsample_t1.png similarity index 100% rename from fig/downsample_t1.png rename to episodes/fig/downsample_t1.png diff --git a/fig/exercise_complete.png b/episodes/fig/exercise_complete.png similarity index 100% rename from fig/exercise_complete.png rename to episodes/fig/exercise_complete.png diff --git a/fig/exercise_fmri.png b/episodes/fig/exercise_fmri.png similarity index 100% rename from fig/exercise_fmri.png rename to episodes/fig/exercise_fmri.png diff --git a/fig/exercise_hollow_skull.png b/episodes/fig/exercise_hollow_skull.png similarity index 100% rename from fig/exercise_hollow_skull.png rename to episodes/fig/exercise_hollow_skull.png diff --git a/fig/exercise_inverted_mask.png b/episodes/fig/exercise_inverted_mask.png similarity index 100% rename from fig/exercise_inverted_mask.png rename to episodes/fig/exercise_inverted_mask.png diff --git a/fig/exercise_masked_fmri.png b/episodes/fig/exercise_masked_fmri.png similarity index 100% rename from fig/exercise_masked_fmri.png rename to episodes/fig/exercise_masked_fmri.png diff --git a/fig/exercise_resampled_hollow_skull.png b/episodes/fig/exercise_resampled_hollow_skull.png similarity index 100% rename from fig/exercise_resampled_hollow_skull.png rename to episodes/fig/exercise_resampled_hollow_skull.png diff --git a/fig/exercise_t1.png b/episodes/fig/exercise_t1.png similarity index 100% rename from fig/exercise_t1.png rename to episodes/fig/exercise_t1.png diff --git a/fig/fmri_data.png b/episodes/fig/fmri_data.png similarity index 100% rename from fig/fmri_data.png rename to episodes/fig/fmri_data.png diff --git a/fig/fmriprep_anat_out.png b/episodes/fig/fmriprep_anat_out.png similarity index 100% rename from fig/fmriprep_anat_out.png rename to episodes/fig/fmriprep_anat_out.png diff --git a/fig/fmriprep_func_out.png b/episodes/fig/fmriprep_func_out.png similarity index 100% rename from fig/fmriprep_func_out.png rename to episodes/fig/fmriprep_func_out.png diff --git a/fig/fmriprep_surf_out.png b/episodes/fig/fmriprep_surf_out.png similarity index 100% rename from fig/fmriprep_surf_out.png rename to episodes/fig/fmriprep_surf_out.png diff --git a/fig/group_compare.png b/episodes/fig/group_compare.png similarity index 100% rename from fig/group_compare.png rename to episodes/fig/group_compare.png diff --git a/fig/heatmap.png b/episodes/fig/heatmap.png similarity index 100% rename from fig/heatmap.png rename to episodes/fig/heatmap.png diff --git a/fig/heatmap_show_zeroes.png b/episodes/fig/heatmap_show_zeroes.png similarity index 100% rename from fig/heatmap_show_zeroes.png rename to episodes/fig/heatmap_show_zeroes.png diff --git a/fig/invert_img.png b/episodes/fig/invert_img.png similarity index 100% rename from fig/invert_img.png rename to episodes/fig/invert_img.png diff --git a/fig/inverted_mask_t1.png b/episodes/fig/inverted_mask_t1.png similarity index 100% rename from fig/inverted_mask_t1.png rename to episodes/fig/inverted_mask_t1.png diff --git a/fig/jupyterhub.png b/episodes/fig/jupyterhub.png similarity index 100% rename from fig/jupyterhub.png rename to episodes/fig/jupyterhub.png diff --git a/fig/masked_func.png b/episodes/fig/masked_func.png similarity index 100% rename from fig/masked_func.png rename to episodes/fig/masked_func.png diff --git a/fig/masked_t1.png b/episodes/fig/masked_t1.png similarity index 100% rename from fig/masked_t1.png rename to episodes/fig/masked_t1.png diff --git a/fig/matrix.png b/episodes/fig/matrix.png similarity index 100% rename from fig/matrix.png rename to episodes/fig/matrix.png diff --git a/fig/mr_scan_types.png b/episodes/fig/mr_scan_types.png similarity index 100% rename from fig/mr_scan_types.png rename to episodes/fig/mr_scan_types.png diff --git a/fig/numpy_arrays.png b/episodes/fig/numpy_arrays.png similarity index 100% rename from fig/numpy_arrays.png rename to episodes/fig/numpy_arrays.png diff --git a/fig/removed_t1.png b/episodes/fig/removed_t1.png similarity index 100% rename from fig/removed_t1.png rename to episodes/fig/removed_t1.png diff --git a/fig/resamp_t1.png b/episodes/fig/resamp_t1.png similarity index 100% rename from fig/resamp_t1.png rename to episodes/fig/resamp_t1.png diff --git a/fig/resampled_yeo.png b/episodes/fig/resampled_yeo.png similarity index 100% rename from fig/resampled_yeo.png rename to episodes/fig/resampled_yeo.png diff --git a/fig/roi_44.png b/episodes/fig/roi_44.png similarity index 100% rename from fig/roi_44.png rename to episodes/fig/roi_44.png diff --git a/fig/roi_46.png b/episodes/fig/roi_46.png similarity index 100% rename from fig/roi_46.png rename to episodes/fig/roi_46.png diff --git a/fig/schz_r.png b/episodes/fig/schz_r.png similarity index 100% rename from fig/schz_r.png rename to episodes/fig/schz_r.png diff --git a/fig/slicing.png b/episodes/fig/slicing.png similarity index 100% rename from fig/slicing.png rename to episodes/fig/slicing.png diff --git a/fig/t1_img.png b/episodes/fig/t1_img.png similarity index 100% rename from fig/t1_img.png rename to episodes/fig/t1_img.png diff --git a/fig/terminal.png b/episodes/fig/terminal.png similarity index 100% rename from fig/terminal.png rename to episodes/fig/terminal.png diff --git a/fig/thick_17.png b/episodes/fig/thick_17.png similarity index 100% rename from fig/thick_17.png rename to episodes/fig/thick_17.png diff --git a/fig/thick_7.png b/episodes/fig/thick_7.png similarity index 100% rename from fig/thick_7.png rename to episodes/fig/thick_7.png diff --git a/fig/thin_17.png b/episodes/fig/thin_17.png similarity index 100% rename from fig/thin_17.png rename to episodes/fig/thin_17.png diff --git a/fig/thin_7.png b/episodes/fig/thin_7.png similarity index 100% rename from fig/thin_7.png rename to episodes/fig/thin_7.png diff --git a/fig/timeseries.png b/episodes/fig/timeseries.png similarity index 100% rename from fig/timeseries.png rename to episodes/fig/timeseries.png diff --git a/fig/yeo_sep.png b/episodes/fig/yeo_sep.png similarity index 100% rename from fig/yeo_sep.png rename to episodes/fig/yeo_sep.png diff --git a/files/scwg_neuroimaging_workshop.pptx b/episodes/files/scwg_neuroimaging_workshop.pptx similarity index 100% rename from files/scwg_neuroimaging_workshop.pptx rename to episodes/files/scwg_neuroimaging_workshop.pptx diff --git a/index.md b/index.md index 06a7c1b8..1cccaae4 100644 --- a/index.md +++ b/index.md @@ -1,14 +1,19 @@ --- -layout: lesson -root: . # Is the only page that doesn't follow the pattern /:path/index.html -permalink: index.html # Is the only page that doesn't follow the pattern /:path/index.html +permalink: index.html +site: sandpaper::sandpaper_site --- + FIXME: Add information about fMRI analysis in Neuroimaging +:::::::::::::::::::::::::::::::::::::::::: prereq + +## Prerequisites + +Attendees must have some base familiarity with Python in order to comfortably progress through the lesson + + +:::::::::::::::::::::::::::::::::::::::::::::::::: + + -> ## Prerequisites -> -> Attendees must have some base familiarity with Python in order to comfortably progress through the lesson -{: .prereq} -{% include links.md %} diff --git a/_extras/guide.md b/instructors/instructor-notes.md similarity index 93% rename from _extras/guide.md rename to instructors/instructor-notes.md index 6d07f16c..b6ca67a6 100644 --- a/_extras/guide.md +++ b/instructors/instructor-notes.md @@ -1,5 +1,5 @@ --- -title: "Instructor Notes" +title: Instructor Notes --- ## Scope of Lesson @@ -13,6 +13,7 @@ Although the primary learning path presented is focused in on Resting State conn Intro to fMRI analyses is meant to be taught as part of a broader Neuroimaging carpentries lesson. As such, it makes some basic assumptions (although we provide very quick overviews) of concepts taught in the pre-requisite workshop: [Intro to MRI](https://carpentries-incubator.github.io/SDC-BIDS-IntroMRI/). In addition, we expect learners to have some basic familiarity with Python (Software Carpentries provides such a [workshop](https://swcarpentry.github.io/python-novice-inflammation/)). In summary the core pre-requisites that we expect learners to have some familiarity with are the following: + - Basic Python programming knowledge - [DataLad](datalad.org) - [BIDS](bids.neuroimaging.io) and the pyBIDS for querying BIDS datasets @@ -22,7 +23,7 @@ In summary the core pre-requisites that we expect learners to have some familiar ### Teaching Style -This workshop as originally designed to be a live-coding session where learners are walked through the process of cleaning and analyzing resting state fMRI data. To this end, all lessons in the carpentries have a mirror Jupyter Notebook codebook. Jupyter Notebooks can be found in [the code directory of the git repo](https://github.com/carpentries-incubator/SDC-BIDS-fMRI/tree/gh-pages/code). Each episode in `code/` mirrors the corresponding episode in the workshop's carpentries site. +This workshop as originally designed to be a live-coding session where learners are walked through the process of cleaning and analyzing resting state fMRI data. To this end, all lessons in the carpentries have a mirror Jupyter Notebook codebook. Jupyter Notebooks can be found in [the code directory of the git repo](https://github.com/carpentries-incubator/SDC-BIDS-fMRI/tree/gh-pages/code). Each episode in `code/` mirrors the corresponding episode in the workshop's carpentries site. In addition, we provide two versions of each notebook i.e: @@ -38,34 +39,36 @@ The notebook without the `_solutions` suffix ("workshop notebook") is mostly bla Companion [slides](https://docs.google.com/presentation/d/1er6dQcERL-Yeb5-7A29tJnmqgHNaLpTLXM3e-SmpjDg/edit#slide=id.g484812a0c7_6_1) are provided for instructors to provide visual aid during certain components of the workshop. We recommend that slides are used for the following episodes: ##### Slides 1-5 - Introduction + **May be used in lieu of Episode 1** -Introductory slides that present the goals of the workshop as well as give an overview of the material to be presented to the learners. +Introductory slides that present the goals of the workshop as well as give an overview of the material to be presented to the learners. #### Slides 10-18 - Preprocessing + **To be presented prior to Episode 2** These slides present an overview of fMRI pre-processing steps as well as motivation for its use prior to statistical analysis of fMRI. It also introduces the fMRIPrep pipeline as the recommended pre-processing pipeline - this is typically of interest to many learners from our experience. In addition, the data used in the workshop are outputs from the fMRIPrep pipeline. - #### Slides 19-20 - Confounds + **To be presented prior to Episode 5** This slide presents a brief motivation behind the need for confound/nuisance regression in fMRI analysis. - ### Tips for hosting a successful workshop -- From experience, the most difficult component of the workshop is the set-up process. It is highly recommended that learners use the [workshop's Binder image](https://mybinder.org/v2/gh/carpentries-incubator/SDC-BIDS-fMRI/gh-pages). It may be useful to set up a separate time-slot to focus on local set-up after the main workshop. +- From experience, the most difficult component of the workshop is the set-up process. It is highly recommended that learners use the [workshop's Binder image](https://mybinder.org/v2/gh/carpentries-incubator/SDC-BIDS-fMRI/gh-pages). It may be useful to set up a separate time-slot to focus on local set-up after the main workshop. - It is often helpful (wherever possible) to have multiple screens on display during the workshop: - - one with the workshop live-coding notebook - - one with a solutions notebook - - This is helpful for learners falling behind -- This workshop is _sequentially dependent_, if learners fail to run a preceding code block they will not be able to run subsequent blocks causing frustration. Try to ensure that any issues from learners get resolved as soon as it is brought up + - one with the workshop live-coding notebook + - one with a solutions notebook + - This is helpful for learners falling behind +- This workshop is *sequentially dependent*, if learners fail to run a preceding code block they will not be able to run subsequent blocks causing frustration. Try to ensure that any issues from learners get resolved as soon as it is brought up - We recommend additional helpers be available during teaching of the workshop. It is often easier to keep the pace/timing of the workshop when a helper is available to address minor issues learners may run into ### Frequent Issues -- **Python Kernel Dies in Binder**: This is a result of the limited memory available (2GB) available in each Binder instance. This can be prevented by **shutting down a notebook** once it is done with. *Simply closing a notebook does not shut it down and will continue to take up memory!*. - - Another solution is to request [additional resources](https://github.com/jupyterhub/mybinder.org-deploy/issues/new?labels=impact&template=request_resources.md) from the Binder maintainers +- **Python Kernel Dies in Binder**: This is a result of the limited memory available (2GB) available in each Binder instance. This can be prevented by **shutting down a notebook** once it is done with. *Simply closing a notebook does not shut it down and will continue to take up memory!*. + - Another solution is to request [additional resources](https://github.com/jupyterhub/mybinder.org-deploy/issues/new?labels=impact&template=request_resources.md) from the Binder maintainers - The request form above can also be done to extend the timeout duration so that the connection isn't lost during a break -{% include links.md %} + + diff --git a/_extras/discuss.md b/learners/discuss.md similarity index 58% rename from _extras/discuss.md rename to learners/discuss.md index bfc33c50..515e3baf 100644 --- a/_extras/discuss.md +++ b/learners/discuss.md @@ -1,6 +1,9 @@ --- title: Discussion --- + FIXME -{% include links.md %} + + + diff --git a/learners/reference.md b/learners/reference.md new file mode 100644 index 00000000..cf6f94d5 --- /dev/null +++ b/learners/reference.md @@ -0,0 +1,11 @@ +--- +title: 'Glossary' +--- + +## Glossary + +FIXME + + + + diff --git a/setup.md b/learners/setup.md similarity index 92% rename from setup.md rename to learners/setup.md index b79d36ed..f22f0e8b 100644 --- a/setup.md +++ b/learners/setup.md @@ -1,159 +1,172 @@ ---- -title: Setup ---- -# Setting up the tutorial environment - -## Binder - -Using Binder is the easiest and fastest way to get started with the workshop. Binder is a virtual environment containing the full computing environment required in order to go through the workshop. Note that Binder hosts the environment on a cloud server and therefore internet access is required to launch it. - -### Setting up Binder - -#### **Step 1**: - -Click the link here to spin up the workshop environment: [Binder Workshop](https://mybinder.org/v2/gh/carpentries-incubator/SDC-BIDS-fMRI/gh-pages?urlpath=lab/tree/code) -You will see an interface that looks like the following: - -![JupyterHub](./fig/jupyterhub.png){:class="img-static"} - -The left-hand pane shows a list of workshop notebooks that contain the content of the workshop itself. Before jumping into the workshop notebooks we need to perform a setup step to pull the neuroimaging data that will be used in the workshop... - -#### **Step 2**: - -Once the environment is launched click the "Terminal" button under "Other". - -![JupyterHub Launch Terminal](./fig/terminal.png){:class="img-static"} - -Clicking on this will launch a console. From here, copy and paste the following lines - -```{bash} -# Move into the code directory -cd code - -# Run the setup_workshop script -./setup_workshop -``` - -Hit enter once pasted and you should see the following - -![Console Filled](./fig/console_filled.png){:class="img-static"} - -This will begin downloading the data required for the workshop onto your Binder instance so that it is usable for the workshop. Once started *do not close the tab by pressing the "x" button.* Instead, you may now open and begin working through the workshop notebooks. - - -## Getting workshop material locally - -### Method 1: Downloading directly from the repository - -On the GitHub repo (this page), click the green button that says "Clone or download", then click **Download ZIP**. Once downloaded, extract the ZIP file. - -### Method 2: Using Git - -Using this method requires a (very) useful piece of software called git. The process of installing git depends heavily on whether you're using MacOS, Windows or Linux. Follow the instructions in the link below to set up git on your PC: - -[Installing Git](https://git-scm.com/book/en/v2/Getting-Started-Installing-Git) - -Once you've installed git, open up your terminal and do the following: - -``` -git clone https://github.com/carpentries-incubator/SDC-BIDS-fMRI.git -``` - -This will download the repository directly into your current directory. - -### Setting up Python environment -We use python version 3.6.0, but any newer version should also work (Python 2 versions haven't been tested). There are many methods to setting up a python environment but we'd recommend using some sort of virtual environment as to not break your system python install. Two methods (of many) are listed below: - -### Method 1: Setting up conda environment (easiest) [Windows, Linux, MacOS] -For easy set-up we recommend [Anaconda](https://www.anaconda.com/download/) to manage python packages for scientific computing. Once installed, setting up the python environment can be done quite easily: - -#### Windows -1. Install Anaconda Python version 3.7 -2. Open **Anaconda Navigator** -3. Click on **Environments** on the left pane -4. Click **Create** then type in sdc-bids-fmri -5. In the sdc-bids-fmri entry click the play button then click **Open Terminal** -6. In terminal type: -``` -conda install -y numpy pandas scipy scikit-learn matplotlib jupyter ipykernel nb_conda -conda install -y -c conda-forge awscli -pip install nilearn nibabel -``` -7. Close the terminal, click on the play button again and open **Jupyter Notebook** -8. Navigate to sdc-bids-fmri folder you downloaded earlier. -9. Done! - -#### [Linux](Linux) and MacOS - -After installing Anaconda, open terminal and type: - -``` -cd sdc-bids-fmri -conda create -p ./sdc-fmri -source activate $(pwd)/sdc-fmri -conda install numpy pandas scipy scikit-learn matplotlib jupyter ipykernel nb_conda -conda install -c conda-forge awscli -pip install nilearn nibabel - -``` -##### Method 2: Using pyenv [Linux, MacOS] -An alternative method uses [pyenv](https://github.com/pyenv/pyenv) with [pyenv virtualenv](https://github.com/pyenv/pyenv-virtualenv). This is a favourite because it seamlessly integrates multiple python versions and environments into your system while maintaining use of pip (instead of conda). -``` -cd sdc-bids-fmri -pyenv virtualenv 3.6.0 sdc-fmri -echo sdc-fmri > .python-version -pip install --requirement requirements.txt -``` - -## Acquiring the data -This tutorial uses data derived from the **UCLA Consortium for Neuropsychiatric Phenomics LA5c Study [1]**. - -To download (**warning: large download size!**) the subset of the data used for the tutorial: - -``` -cd sdc-bids-fmri - -# download T1w scans -cat download_list | \ - xargs -I '{}' aws s3 sync --no-sign-request \ - s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/{}/anat \ - ./data/ds000030/{}/anat - -# download resting state fMRI scans -cat download_list | \ - xargs -I '{}' aws s3 sync --no-sign-request \ - s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/{}/func \ - ./data/ds000030/{}/func \ - --exclude '*' \ - --include '*task-rest_bold*' - -# download fmriprep preprocessed anat data -cat download_list | \ - xargs -I '{}' aws s3 sync --no-sign-request \ - s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/derivatives/fmriprep/{}/anat \ - ./data/ds000030/derivatives/fmriprep/{}/anat - -# download fmriprep preprocessed func data -cat download_list | \ - xargs -I '{}' aws s3 sync --no-sign-request \ - s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/derivatives/fmriprep/{}/func \ - ./data/ds000030/derivatives/fmriprep/{}/func \ - --exclude '*' \ - --include '*task-rest_bold*' -``` -Finally open up the jupyter notebook to explore the tutorials: -``` -cd sdc-bids-fmri - -#Include below line if using anaconda environment -source activate $(pwd)/sdc-fmri - -jupyter notebook -``` - -**Reference** - -[1] Gorgolewski KJ, Durnez J and Poldrack RA. Preprocessed Consortium for Neuropsychiatric Phenomics dataset [version 2; referees: 2 approved]. F1000Research 2017, 6:1262 -(https://doi.org/10.12688/f1000research.11964.2) - -{% include links.md %} +--- +title: Setup +--- + +# Setting up the tutorial environment + +## Binder + +Using Binder is the easiest and fastest way to get started with the workshop. Binder is a virtual environment containing the full computing environment required in order to go through the workshop. Note that Binder hosts the environment on a cloud server and therefore internet access is required to launch it. + +### Setting up Binder + +#### **Step 1**: + +Click the link here to spin up the workshop environment: [Binder Workshop](https://mybinder.org/v2/gh/carpentries-incubator/SDC-BIDS-fMRI/gh-pages?urlpath=lab/tree/code) +You will see an interface that looks like the following: + +![](./fig/jupyterhub.png){alt='JupyterHub' class="img-static"} + +The left-hand pane shows a list of workshop notebooks that contain the content of the workshop itself. Before jumping into the workshop notebooks we need to perform a setup step to pull the neuroimaging data that will be used in the workshop... + +#### **Step 2**: + +Once the environment is launched click the "Terminal" button under "Other". + +![](./fig/terminal.png){alt='JupyterHub Launch Terminal' class="img-static"} + +Clicking on this will launch a console. From here, copy and paste the following lines + +```{bash} +# Move into the code directory +cd code + +# Run the setup_workshop script +./setup_workshop +``` + +Hit enter once pasted and you should see the following + +![](./fig/console_filled.png){alt='Console Filled' class="img-static"} + +This will begin downloading the data required for the workshop onto your Binder instance so that it is usable for the workshop. Once started *do not close the tab by pressing the "x" button.* Instead, you may now open and begin working through the workshop notebooks. + +## Getting workshop material locally + +### Method 1: Downloading directly from the repository + +On the GitHub repo (this page), click the green button that says "Clone or download", then click **Download ZIP**. Once downloaded, extract the ZIP file. + +### Method 2: Using Git + +Using this method requires a (very) useful piece of software called git. The process of installing git depends heavily on whether you're using MacOS, Windows or Linux. Follow the instructions in the link below to set up git on your PC: + +[Installing Git](https://git-scm.com/book/en/v2/Getting-Started-Installing-Git) + +Once you've installed git, open up your terminal and do the following: + +``` +git clone https://github.com/carpentries-incubator/SDC-BIDS-fMRI.git +``` + +This will download the repository directly into your current directory. + +### Setting up Python environment + +We use python version 3.6.0, but any newer version should also work (Python 2 versions haven't been tested). There are many methods to setting up a python environment but we'd recommend using some sort of virtual environment as to not break your system python install. Two methods (of many) are listed below: + +### Method 1: Setting up conda environment (easiest) [Windows, Linux, MacOS] + +For easy set-up we recommend [Anaconda](https://www.anaconda.com/download/) to manage python packages for scientific computing. Once installed, setting up the python environment can be done quite easily: + +#### Windows + +1. Install Anaconda Python version 3.7 +2. Open **Anaconda Navigator** +3. Click on **Environments** on the left pane +4. Click **Create** then type in sdc-bids-fmri +5. In the sdc-bids-fmri entry click the play button then click **Open Terminal** +6. In terminal type: + +``` +conda install -y numpy pandas scipy scikit-learn matplotlib jupyter ipykernel nb_conda +conda install -y -c conda-forge awscli +pip install nilearn nibabel +``` + +7. Close the terminal, click on the play button again and open **Jupyter Notebook** +8. Navigate to sdc-bids-fmri folder you downloaded earlier. +9. Done! + +#### [Linux](Linux) and MacOS + +After installing Anaconda, open terminal and type: + +``` +cd sdc-bids-fmri +conda create -p ./sdc-fmri +source activate $(pwd)/sdc-fmri +conda install numpy pandas scipy scikit-learn matplotlib jupyter ipykernel nb_conda +conda install -c conda-forge awscli +pip install nilearn nibabel + +``` + +##### Method 2: Using pyenv [Linux, MacOS] + +An alternative method uses [pyenv](https://github.com/pyenv/pyenv) with [pyenv virtualenv](https://github.com/pyenv/pyenv-virtualenv). This is a favourite because it seamlessly integrates multiple python versions and environments into your system while maintaining use of pip (instead of conda). + +``` +cd sdc-bids-fmri +pyenv virtualenv 3.6.0 sdc-fmri +echo sdc-fmri > .python-version +pip install --requirement requirements.txt +``` + +## Acquiring the data + +This tutorial uses data derived from the **UCLA Consortium for Neuropsychiatric Phenomics LA5c Study [1]**. + +To download (**warning: large download size!**) the subset of the data used for the tutorial: + +``` +cd sdc-bids-fmri + +# download T1w scans +cat download_list | \ + xargs -I '{}' aws s3 sync --no-sign-request \ + s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/{}/anat \ + ./data/ds000030/{}/anat + +# download resting state fMRI scans +cat download_list | \ + xargs -I '{}' aws s3 sync --no-sign-request \ + s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/{}/func \ + ./data/ds000030/{}/func \ + --exclude '*' \ + --include '*task-rest_bold*' + +# download fmriprep preprocessed anat data +cat download_list | \ + xargs -I '{}' aws s3 sync --no-sign-request \ + s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/derivatives/fmriprep/{}/anat \ + ./data/ds000030/derivatives/fmriprep/{}/anat + +# download fmriprep preprocessed func data +cat download_list | \ + xargs -I '{}' aws s3 sync --no-sign-request \ + s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/derivatives/fmriprep/{}/func \ + ./data/ds000030/derivatives/fmriprep/{}/func \ + --exclude '*' \ + --include '*task-rest_bold*' +``` + +Finally open up the jupyter notebook to explore the tutorials: + +``` +cd sdc-bids-fmri + +#Include below line if using anaconda environment +source activate $(pwd)/sdc-fmri + +jupyter notebook +``` + +**Reference** + +[1] Gorgolewski KJ, Durnez J and Poldrack RA. Preprocessed Consortium for Neuropsychiatric Phenomics dataset [version 2; referees: 2 approved]. F1000Research 2017, 6:1262 +([https://doi.org/10.12688/f1000research.11964.2](https://doi.org/10.12688/f1000research.11964.2)) + + + + diff --git a/profiles/learner-profiles.md b/profiles/learner-profiles.md new file mode 100644 index 00000000..434e335a --- /dev/null +++ b/profiles/learner-profiles.md @@ -0,0 +1,5 @@ +--- +title: FIXME +--- + +This is a placeholder file. Please add content here. diff --git a/reference.md b/reference.md deleted file mode 100644 index 8c826167..00000000 --- a/reference.md +++ /dev/null @@ -1,9 +0,0 @@ ---- -layout: reference ---- - -## Glossary - -FIXME - -{% include links.md %} diff --git a/site/README.md b/site/README.md new file mode 100644 index 00000000..42997e3d --- /dev/null +++ b/site/README.md @@ -0,0 +1,2 @@ +This directory contains rendered lesson materials. Please do not edit files +here.