Tutorial 1: Statistics Module
In this tutorial you will set up your first linear model with BrainStat. To this end, we will load some sample data from the ABIDE dataset. Note that,contrary to the results shown in our manuscript, we are only using a few sites to reduce computation time in this tutorial. As such the results shown here differ from those reported in our manuscript. To get identical results, simply set `sites` to `nan`.
But before we get started, please note that this tutorial requires that you have BrainStat and all its dependencies installed. For more details see our installation guide.
sites = ["PITT", "OLIN", "OHSU"];
[cortical_thickness, demographics] = fetch_abide_data('sites', sites);
[surf_left, surf_right] = fetch_template_surface('civet41k');
mask = fetch_mask('civet41k');
Lets have a look at the cortical thickness data. To do this, we will use the surface plotter included with BrainSpace. Note that we use a recent version of BrainSpace in this tutorial (released on 16 July, 2021). If you've installed an older version please update before continuing. Lets plot mean thickness. Note that we've defined a function, pretty_plot, at the end of this file to make the plots a little prettier when they're embedded inside live scripts. If you run this tutorial outside of the live scripts you'll want to comment out the pretty_plot lines.
obj = plot_hemispheres(...
mean(cortical_thickness)', ...
{surf_left, surf_right}, ...
'labeltext', {{'Cortical', 'Thickness'}} ...
);
pretty_plot(obj);
Next, lets see whether cortical thickness is related to age in our sample data. To this end we can create a linear model with BrainStat. First we declare the behavioral variables as FixedEffects. Once, that's done we can create the model by simply adding the terms together. Lets set up a model with age and patient status as fixed effects.
% Change the encoding of patient/control from 1/2 to strings.
conversion = ["Patient"; "Control"];
demographics.DX_GROUP = conversion(demographics.DX_GROUP);
term_age = FixedEffect(demographics.AGE_AT_SCAN, 'Age');
term_patient = FixedEffect(demographics.DX_GROUP);
model = term_age + term_patient;
disp(model)
As you can see, model contains two properties: the names of each variable and a matrix containing the effects in a observation-by-variable matrix.
Beside simple fixed effects, we may also be interested in interaction effects. We can add these to the model by multiplying terms. Lets create a model with an interaction between age and patient status.
model_interaction = 1 + term_age + term_patient + ...
term_age * term_patient;
Now, imagine we have some cortical marker (e.g. cortical thickness) for each subject, and we want to evaluate whether this marker changes with age whilst correcting for effects of patient status. To do this, we can use the model we defined before, and a contrast in observations (here: age). Then we simply initialize an SLM model and fit it to the cortical thickness data. Note that the fitting may take a few minutes.
contrast_age = demographics.AGE_AT_SCAN;
BrainStatModel = SLM( ...
model, ...
contrast_age, ...
'surf', 'civet41k', ...
'correction', 'rft', ...
'mask', mask);
BrainStatModel.fit(cortical_thickness);
We can access the outcome of this model through its properties. Lets plot the t-values and the vertexwise p-values derived from a random field theory correction.
to_plot1 = [BrainStatModel.t', BrainStatModel.P.pval.P'];
to_plot1(~mask,:) = inf;
plot1 = plot_hemispheres(...
to_plot1, ...
{surf_left, surf_right}, ...
'labeltext', {'t-values', 'p-values'} ...
);
plot1.colormaps({[parula; .7 .7 .7], [flipud(hot); .7 .7 .7]});
plot1.colorlimits([-6, 4; 0, 0.05]);
pretty_plot(plot1);
BrainStat also allows for assessing significant clusters and peaks. The data on clusters are stored in tables inside BrainStatModel.P.clus and information on the peaks is stored in BrainStatModel.P.peak. If a two-tailed test is run (BrainStat defaults to two-tailed), a table is returned for each tail. The first table uses the contrast as provided, the second table uses the inverse contrast. If a one-tailed test is performed, then only a single table is returned. Lets print the first 15 rows of the inverted contrast cluster table.
disp(BrainStatModel.P.clus{2}(1:15, :))
Here, we see that cluster 1 contains 8738 vertices and is significant at a p-value of 7.45e-08. Clusters are sorted by p-value; later clusters will generally be smaller and have higher p-values. Lets now have a look at the peaks within these clusters.
disp(BrainStatModel.P.peak{2}(1:15, :))
Within cluster 1, we are able to detect several peaks. The most significant peak (t=6.3155) occurs at vertex 72883, which is inside the visual network as defined by the Yeo-7 networks. Note that the Yeo network membership is only provided if the surface is specified as a template name as we did here. For custom surfaces, or pre-loaded surfaces (as we will use below) this column is omitted.
By default BrainStat uses a two-tailed test. If you want to get a one-tailed test, simply specify it in the SLM model initialization with 'two_tailed', false. Note that the one-tailed test will test for positive t-values. If you want to test for negative t-values, simply invert the contrast. We may hypothesize based on prior research that cortical thickness decreases with age, so we could specify this as follows. Note the minus in front of contrast_age to test for decreasing thickness with age.
BrainStatModel_onetailed = SLM(...
model, ...
-contrast_age, ...
'correction', 'rft', ...
'surf', {surf_left, surf_right}, ...
'mask', mask, ...
'two_tailed', false ...
);
BrainStatModel_onetailed.fit(cortical_thickness);
to_plot2 = [BrainStatModel_onetailed.t', ...
BrainStatModel_onetailed.P.pval.P'];
to_plot2(~mask,:) = inf;
plot2 = plot_hemispheres(...
to_plot2, ...
{surf_left, surf_right}, ...
'labeltext', {'t-values', 'p-values'} ...
);
plot2.colormaps({[parula; .7 .7 .7], [flipud(hot); .7 .7 .7]});
plot2.colorlimits([-4, 6; 0, 0.05]);
pretty_plot(plot2)
Similarly, we could also use the patient/control status as a contrast, instead of age. This would work as follows:
contrast_patient = model.Patient - model.Control;
BrainStatModel_patient = SLM( ...
model, ...
contrast_patient, ...
'surf', 'civet41k', ...
'correction', 'rft', ...
'mask', mask);
BrainStatModel_patient.fit(cortical_thickness);
Now, imagine that instead of using a fixed effects model, you would prefer a mixed effects model wherein the scanning site is a random variable. This is simple to set up. All you need to do is initialize the site term with the MixedEffect class, all other code remains identical.
random_site = MixedEffect(demographics.SITE_ID);
model_random = term_age + term_patient + random_site;
BrainStatModel_random = SLM(model_random, contrast_age, ...
'surf', {surf_left, surf_right}, 'mask', mask, 'correction', 'rft');
BrainStatModel_random.fit(cortical_thickness)
to_plot3 = [BrainStatModel_random.t', BrainStatModel_random.P.pval.P'];
to_plot3(~mask,:) = inf;
plot3 = plot_hemispheres(...
to_plot3, ...
{surf_left, surf_right}, ...
'labeltext', {'t-values', 'p-values'} ...
);
plot3.colormaps({[parula; .7 .7 .7], [flipud(hot); .7 .7 .7]});
plot3.colorlimits([-4, 6; 0, 0.05]);
pretty_plot(plot3);
It appears that, after correction for site as a random variable, the significant results are far more localized. That concludes the basic usage of the BrainStat for statistical models. In the next tutorial we'll show you how to use the context decoding module.
function pretty_plot(obj)
% Makes the BrainSpace plots prettier inside live scripts.
obj.handles.figure.Units = 'pixels';
set(obj.handles.axes, 'Units', 'pixels')
set(obj.handles.colorbar, 'Units', 'pixels');
obj.handles.figure.Position = [1, 1, 700, size(obj.data, 2) * 95];
set(obj.handles.text, 'FontSize', 12);
for jj = 1:size(obj.data, 2) % variates
obj.handles.colorbar(jj).Position = [...
500, ...
30 + (size(obj.data, 2) * 95) - 10 - jj*90, ...
5, ...
50];
for ii = 1:4 % views
obj.handles.axes(jj, ii).Position = [...
75 + (ii-1)*105, ...
(size(obj.data, 2) * 95) - jj*95, ...
100, ...
100];
end
end
end
Tutorial 1: Statistics Module
In this tutorial you will set up your first linear model with BrainStat. To this end, we will load some sample data from the MICS dataset. But before we get started, please note that this tutorial requires that you have BrainStat and all its dependencies installed. For more details see our installation guide.
Lets have a look at the cortical thickness data and the demographics data. Note that we've defined a function, pretty_plot, at the end of this file to make the plots a little prettier when they're embedded inside live scripts. If you run this tutorial outside of a live script you'll want to comment out the pretty_plot lines.
[cortical_thickness, demographics] = fetch_mics_data();
[surf_left, surf_right] = fetch_template_surface('fsaverage5');
surface = fetch_template_surface('fsaverage5', 'join', true);
obj = plot_hemispheres(...
mean(cortical_thickness)', ...
{surf_left, surf_right}, ...
'labeltext', {{'Cortical', 'Thickness'}} ...
);
pretty_plot(obj);
disp(demographics(1:10,:))
for ii = 1:2
fprintf(['Visit %d: N=%d, %d females, mean subject age: %2.2f, ' ...
'standard deviation of age: %2.2f.\n'], ...
ii, ...
sum(demographics.VISIT==ii), ...
sum(demographics.SEX(demographics.VISIT==ii) == "F"), ...
mean(demographics.AGE_AT_SCAN(demographics.VISIT==ii)), ...
std(demographics.AGE_AT_SCAN(demographics.VISIT==ii)));
end
Demographics contains four variables: a subject ID, a visit number (some subjects visited multiple times), their age at the time of scanning and their sex. We've also printed the means and standard deviations for visits 1, 2, and both. Next, we will assess whether a subject's age is related to their cortical thickness. To this end we can create a linear model with BrainStat. For our first model, we will only consider the effect of age, i.e. we will disregard the effect of sex and that some subjects visit twice.
A basic model
term_age = FixedEffect(demographics.AGE_AT_SCAN, 'Age');
model = term_age;
disp(model)
As you can see, model contains two properties: the names of each variable and a matrix containing the effects in a observation-by-variable matrix. As you can see, an intercept is included in the model by default. We can also access the vectors related to each effect by their name i.e. model.intercept and model.Age will return the vectors of the intercept and age, respectively.
Next, we must define a contrast in observations (here: age). With the model, contrast, and cortical thickness data, we can then assess whether there is an effect of age on cortical thickness.
contrast_age = demographics.AGE_AT_SCAN;
mask = fetch_mask('fsaverage5');
slm = SLM( ...
model, ...
contrast_age, ...
'surf', 'fsaverage5', ...
'correction', {'rft', 'fdr'}, ...
'cluster_threshold', 0.01, ...
'mask', mask);
slm.fit(cortical_thickness);
We can access the outcome of this model through its properties. Lets plot the t-values, as well as the peak and cluster p-values derived from a random field theory correction. We've defined a function for slm results at the end of this file as we'll be doing this a lot.
plot_slm(slm, {surf_left, surf_right}, mask, true, true);
to_plot = [slm.t', slm.P.pval.C', slm.P.pval.P', slm.Q'];
to_plot(~mask, :) = inf;
obj = plot_hemispheres(to_plot, {surf_left, surf_right}, ...
'LabelText', {'t-values', 'Cluster p-values', ...
'Peak p-values', 'Vertex p-values'});
obj.colorlimits([-6, 4; 0, 0.05; 0, 0.05; 0, 0.05])
obj.colormaps({[parula; 0.7, 0.7, 0.7], ...
[flipud(autumn); 0.7, 0.7, 0.7], ...
[flipud(autumn); 0.7, 0.7, 0.7], ...
[flipud(autumn); 0.7, 0.7, 0.7]});
Only clusters are significant, and not peaks. This suggests that the age effect covers large regions, rather than local foci. Furthermore, at the vertex level we only find few small groups of significant vertices. Lets have a closer look at the clusters and their peaks. The data on clusters are stored in tables inside slm.P.clus and information on the peaks is stored in slm.P.peak. If a two-tailed test is run (BrainStat defaults to two-tailed), a table is returned for each tail. The first table uses the contrast as provided, the second table uses the inverse contrast. If a one-tailed test is performed, then only a single table is returned. Lets print the first 15 rows of the inverted contrast cluster table.
disp(slm.P.clus{2}(1:15, :))
Here, we see that cluster 1 contains 141 vertices. Furthermore, we see that clusters 1-9 are significant P<0.05. Lets now have a look at the peaks within these clusters.
disp(slm.P.peak{2}(1:15, :))
We are able to detect several peaks. The peak with the highest t-statistic (t=5.6954) occurs at vertex 18720, which is inside the ventral attention network as defined by the Yeo-7 networks. However, this peak occurs inside a non-significant cluster (11). Note that the Yeo network membership is only provided if the surface is specified as a template name as we did here. For custom surfaces, or pre-loaded surfaces (as we will use below) this column is omitted.
Interaction models
Similarly to age, we can also test for the effect of sex on cortical thickness.
term_sex = FixedEffect(demographics.SEX);
model_sex = term_sex;
contrast_sex = (demographics.SEX == "M") - (demographics.SEX == "F");
slm_sex = SLM( ...
model_sex, ...
contrast_sex, ...
'surf', surface, ...
'correction', {'rft', 'fdr'}, ...
'cluster_threshold', 0.01, ...
'mask', mask);
slm_sex.fit(cortical_thickness);
plot_slm(slm_sex, {surf_left, surf_right}, mask);
Here, we find no significant effects of sex on cortical thickness. However, as we've already established, age has an effect on cortical thickness. So we may want to correct for this effect before evaluating whether sex has an effect on cortical thickenss. Lets make a new model that includes the effect of age.
model_sexage = term_age + term_sex;
slm_sexage = SLM( ...
model_sexage, ...
contrast_sex, ...
'surf', surface, ...
'correction', {'rft', 'fdr'}, ...
'cluster_threshold', 0.01, ...
'mask', mask);
slm_sexage.fit(cortical_thickness);
plot_slm(slm_sexage, {surf_left, surf_right}, mask);
After accounting for the effect of age, we still don't find significant clusters of effect of sex on cortical thickness. However, it could be that age affects men and women differently. To account for this, we could include an interaction effect into the model. Lets run the model again with an interaction effect.
model_sexage_int = term_age + term_sex + term_age * term_sex;
slm_sexage = SLM( ...
model_sexage_int, ...
contrast_sex, ...
'surf', surface, ...
'correction', {'rft', 'fdr'}, ...
'cluster_threshold', 0.01, ...
'mask', mask);
slm_sexage.fit(cortical_thickness);
plot_slm(slm_sexage, {surf_left, surf_right}, mask);
After including the interaction effect, we now do find significant effects of sex on cortical thickness in several clusters.
Instead of looking at the effects of age or sex, we could also look at whether the cortex of men and women changes differently with age by comparing their interaction effects. To do this, we have to modify the contrast.
contrast_sex_int = (demographics.AGE_AT_SCAN .* ...
(demographics.SEX == "M")) - ...
(demographics.AGE_AT_SCAN .* ...
(demographics.SEX == "F"));
slm_sexage_int = SLM( ...
model_sexage_int, ...
contrast_sex_int, ...
'surf', surface, ...
'correction', 'rft', ...
'cluster_threshold', 0.01, ...
'mask', mask);
slm_sexage_int.fit(cortical_thickness);
obj = plot_slm(slm_sexage_int, {surf_left, surf_right}, mask);
obj.colorlimits([-3, 3; 0, 0.05]);
Indeed, it appears that the interaction effect between sex and age is quite different across men and women, with stronger effects occuring in women.
One-tailed models
Imagine that, based on prior research, we hypothesize that men have higher cortical thickness than women. In that case, we could run this same model with a one-tailed test, rather than a two-tailed test. By default BrainStat uses a two-tailed test. If you want to get a one-tailed test, simply specify it in the SLM model initialization with 'two_tailed', false. Note that the one-tailed test will test for the significance of positive t-values. If you want to test for the significance of negative t-values, simply change the sign of the contrast. We may hypothesize based on prior research that cortical thickness decreases with age, so we could specify this as follows. Note the minus in front of contrast_age to test for decreasing thickness with age.
slm_onetail = SLM( ...
model_sexage_int, ...
contrast_sex, ...
'surf', surface, ...
'correction', {'rft', 'fdr'}, ...
'cluster_threshold', 0.01, ...
'two_tailed', false, ...
'mask', mask);
slm_onetail.fit(cortical_thickness);
plot_slm(slm_onetail, {surf_left, surf_right}, mask);
Notice the additional clusters that we find when using a one-tailed test.
Mixed Effects Models
So far, we've considered multiple visits of the same subject as two separate, independent measurements. Clearly, however, such measurements are not independent of each other. To account for this, we could add subject ID as a random effect. Lets do this and test the effect of age on cortical thickness again.
term_subject = MixedEffect(demographics.SUB_ID);
model_mixed = term_age + term_sex + term_age * term_sex + term_subject;
slm_mixed = SLM( ...
model_mixed, ...
-contrast_age, ...
'surf', surface, ...
'correction', 'rft', ...
'cluster_threshold', 0.01, ...
'two_tailed', false, ...
'mask', mask);
slm_mixed.fit(cortical_thickness);
obj = plot_slm(slm_mixed, {surf_left, surf_right}, mask);
obj.colorlimits([-6, 4; 0, 0.05]);
After inclusion of subject as a random variable into the one-tailed model shown earlier, we find fewer and smaller clusters, indicating that by not accounting for the repeated measures structure of the data we were overestimating the significance of effects.
That concludes the basic usage of the BrainStat for statistical models. In the next tutorial we'll show you how to use the context decoding module.
Below are the support functions used in this live script.
function pretty_plot(obj)
% Makes the BrainSpace plots prettier inside live scripts.
obj.handles.figure.Units = 'pixels';
set(obj.handles.axes, 'Units', 'pixels')
set(obj.handles.colorbar, 'Units', 'pixels');
obj.handles.figure.Position = [1, 1, 700, size(obj.data, 2) * 95];
set(obj.handles.text, 'FontSize', 12);
for jj = 1:size(obj.data, 2) % variates
obj.handles.colorbar(jj).Position = [...
500, ...
30 + (size(obj.data, 2) * 95) - 10 - jj*90, ...
5, ...
50];
for ii = 1:4 % views
obj.handles.axes(jj, ii).Position = [...
75 + (ii-1)*105, ...
(size(obj.data, 2) * 95) - jj*95, ...
100, ...
100];
end
end
end
function obj = plot_slm(slm, surfaces, mask, plot_peak, plot_fdr)
% Plot the t-values and p-values of an SLM.
arguments
slm
surfaces
mask
plot_peak (1,1) logical = false
plot_fdr (1,1) logical = false
end
to_plot = [slm.t(:), slm.P.pval.C(:)];
labels = {'t-values', {'Cluster', 'p-values (RFT)'}};
colormaps = {[parula; .7 .7 .7];
[flipud(autumn); .7 .7 .7]};
colorlimits = [-6, 4; 0, 0.05];
if plot_peak
to_plot = [to_plot, slm.P.pval.P(:)];
labels = [labels, {{'Peak', 'p-values (RFT)'}}];
colormaps = [colormaps; {[flipud(autumn); .7 .7 .7]}];
colorlimits = [colorlimits; 0, 0.05];
end
if plot_fdr
to_plot = [to_plot, slm.Q(:)];
labels = [labels, {{'Vertexwise', 'p-values (FDR)'}}];
colormaps = [colormaps; {[flipud(autumn); .7 .7 .7]}];
colorlimits = [colorlimits; 0, 0.05];
end
to_plot(~mask, :) = inf;
obj = plot_hemispheres(...
to_plot, ...
{surfaces{1}, surfaces{2}}, ...
'labeltext', labels ...
);
obj.colormaps(colormaps);
obj.colorlimits(colorlimits);
pretty_plot(obj);
end
Tutorial 2: Context Module
In this tutorial you will learn about the context decoding tools included with BrainStat. The context decoding module consists of three parts: genetic decoding, meta-analytic decoding and histological comparisons. First, we'll consider how to run the genetic decoding analysis. But before we get started, please note that this tutorial, alike the previous one, requires that you have BrainStat and all its dependencies installed. For more details see our installation guide.
Genetics
For genetic decoding we use the Allen Human Brain Atlas through the abagen toolbox. Note that, as abagen is only available in Python, we simply supply the genetic expression of pre-computed regions of interest with default parameters in the MATLAB toolbox. This should suffice for most use-cases, but for more specific use-cases we advise you to use the Python implementation of BrainStat. For a full list of supported parcellations see help fetch_genetic_expression.
[expression, gene_names] = fetch_genetic_expression('schaefer', 100);
img = imagesc(expression);
xlabel('Genes');
ylabel('Regions');
colorbar;
The variable expression is a matrix containing the genetic expression of genes within each region of the atlas. These values will fall in the range [0, 1] where higher values represent higher expression. Some regions may return NaN values for all genes (e.g. see region 60 in the plot). This occurs when there are no samples within this region across all donors. The variable gene_names contains the names of each gene in the same order as the columns of expression.
Meta-Analytic Decoding
To perform meta-analytic decoding, BrainStat uses precomputed Neurosynth maps. Here we test which terms are most associated with a map of cortical thickness. A simple example analysis can be run as follows. First, we will load some cortical thickness data and two cortical surfaces. The ABIDE dataset provides this data on the CIVET surface, so we will also load those surfaces. The surface decoder interpolates the data from the surface to the voxels in the volume that are in between the two input surfaces. Please be aware that this analysis may take several minutes to download and run.
mask = fetch_mask('civet41k');
thickness = fetch_abide_data('sites', 'PITT');
[correlation, feature] = surface_decoder(mean(thickness), ...
'template', 'civet41k');
disp(correlation(1:3));
disp(feature(1:3));
The variable correlation now contains the correlation values between the input map mean(thickness) and each of the neurosynth features. The corresponding feature names are stored in the variable feature.
Histology
For histological decoding we use microstructural profile covariance gradients, as first shown by (Paquola et al, 2019, Plos Biology), computed from the BigBrain dataset. Firstly, lets download the MPC data and compute its gradients.
schaefer_400 = fetch_parcellation('fsaverage5' ,'schaefer', 400);
Reading from version 2
+.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
+.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
+.S10 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: 400; text-align: left; }
+.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
+.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
Tutorial 2: Context Module
In this tutorial you will learn about the context decoding tools included with BrainStat. The context decoding module consists of three parts: genetic decoding, meta-analytic decoding and histological comparisons. First, we'll consider how to run the genetic decoding analysis. But before we get started, please note that this tutorial, alike the previous one, requires that you have BrainStat and all its dependencies installed. For more details see our installation guide. Also, we'll run a simple linear model again like we did in the previous tutorial, as we'll be using the results of this model.
[cortical_thickness, demographics] = fetch_mics_data();
[surf_left, surf_right] = fetch_template_surface('fsaverage5');
mask = fetch_mask('fsaverage5');
term_age = FixedEffect(demographics.AGE_AT_SCAN, 'Age');
term_sex = FixedEffect(demographics.SEX);
term_subject = MixedEffect(demographics.SUB_ID);
model = term_age + term_sex + term_age * term_sex + term_subject;
contrast_age = -demographics.AGE_AT_SCAN;
slm = SLM( ...
model, ...
contrast_age, ...
'surf', 'fsaverage5', ...
'correction', {'rft', 'fdr'}, ...
'mask', mask, ...
'two_tailed', false, ...
'cluster_threshold', 0.01);
slm.fit(cortical_thickness);
Genetics
For genetic decoding we use the Allen Human Brain Atlas through the abagen toolbox. Note that, as abagen is only available in Python, we simply supply the genetic expression of pre-computed regions of interest with default parameters in the MATLAB toolbox. This should suffice for most use-cases, but for more specific use-cases we advise you to use the Python implementation of BrainStat. For a full list of supported parcellations see help fetch_genetic_expression.
[expression, gene_names] = fetch_genetic_expression('schaefer', 100);
img = imagesc(expression);
xlabel('Genes');
ylabel('Regions');
colorbar;
The variable expression is a matrix containing the genetic expression of genes within each region of the atlas. These values will fall in the range [0, 1] where higher values represent higher expression. Some regions may return NaN values for all genes (e.g. see region 60 in the plot). This occurs when there are no samples within this region across all donors. The variable gene_names contains the names of each gene in the same order as the columns of expression. We can also correlate these genes to the t-statistic map derived earlier. Lets do this for the highest correlating gene, WFDC1, and create a scatter plot as well as a surface map.
schaefer_100 = fetch_parcellation('fsaverage5', 'schaefer', 100);
t_schaefer_100 = full2parcel(slm.t, schaefer_100);
WFDC1_expression = expression(:, gene_names == "WFDC1");
no_nan = ~isnan(WFDC1_expression);
t_WFDC1_corr = corr(t_schaefer_100(no_nan), WFDC1_expression(no_nan));
bestfit_coefficients = polyfit(t_schaefer_100(no_nan), WFDC1_expression(no_nan), 1);
figure();
scatter(t_schaefer_100, WFDC1_expression, 200, 'k', '.');
hold on
plot(t_schaefer_100, polyval(bestfit_coefficients, t_schaefer_100), ...
'k', 'LineWidth', 3);
xlabel('t-statistic')
ylabel('WFDC1 Expression')
set(gca, 'FontSize', 16, 'YLim', [0, 1], 'YTick', [0, 1], 'Xtick', [-1.5, 2])
text(-1.4, 0.9, sprintf('r = %0.2f', t_WFDC1_corr), 'FontSize', 16);
WFDC1_expression(isnan(WFDC1_expression)) = -inf;
obj = plot_hemispheres(WFDC1_expression, {surf_left, surf_right}, ...
'parcellation', schaefer_100);
obj.colormaps([.7 .7 .7; parula])
obj.colorlimits([0, 0.85])
pretty_plot(obj)
Meta-Analytic Decoding
To perform meta-analytic decoding, BrainStat uses precomputed Neurosynth maps. Here we test which terms are most associated with our t-statistic map. The surface decoder interpolates the data from the surface to the voxels in the volume that are in between pial and white matter surface. Next, using all voxels that are non-zero both in this volume and in the Neurosynth database, a Pearson correlation is computed for every feature. Please be aware that this analysis may take several minutes to download and run.
[correlation, feature] = meta_analytic_decoder(slm.t, ...
'template', 'fsaverage5');
disp(correlation(1:3));
disp(feature(1:3));
The variable correlation now contains the correlation values between the t-statistic map and each of the neurosynth features. The corresponding feature names are stored in the variable feature. An easy way of visualizing this data is by using a word cloud. These are simple to build with MATLAB:
figure
wc = wordcloud(feature(correlation>0), correlation(correlation>0));
Here, larger words represent higher correlations. It's common to see the largest words be related to anatomy. Looking past anatomical terms, we see several words related to motivation and risk taking behavior e.g. reward, motivation, additction, gambling, punishment.
Histology
For histological decoding we use microstructural profile covariance gradients, as first shown by (Paquola et al, 2019, Plos Biology), computed from the BigBrain dataset. Firstly, lets download the MPC data and compute its gradients.
schaefer_400 = fetch_parcellation('fsaverage5' ,'schaefer', 400);
profiles = read_histology_profile('template', 'fsaverage5');
mpc = compute_mpc(profiles, schaefer_400);
gm = compute_histology_gradients(mpc);
The variable profiles now contains histological profiles sampled at 50 different depths across the cortex, mpc contains the covariance of these profiles, and gm contains their gradients. Depending on your use-case, each of these variables could be of interest, but for purposes of this tutorial we'll plot the gradients to the surface with BrainSpace. For details on what the GradientMaps class, gm, contains please consult the BrainSpace documentation.
[surface_left, surface_right] = fetch_template_surface('fsaverage5');
vertexwise_gradients = parcel2full(gm.gradients{1}(:,1:2), schaefer_400);
vertexwise_gradients(isnan(vertexwise_gradients)) = inf;
obj = plot_hemispheres(vertexwise_gradients, ...
{surface_left, surface_right}, ...
'labeltext', {'MPC Gradient 1', 'MPC Gradient 2'});
obj.colormaps([parula; .7 .7 .7])
pretty_plot(obj);
Note that we no longer use the y-axis regression used in (Paquola et al., 2019, Plos Biology), as such the first gradient becomes an anterior-posterior- gradient.
Resting-State Contextualization
Lastly, BrainStat provides contextualization using resting-state fMRI markers: specifically, with the Yeo functional networks (Yeo et al., 2011, Journal of Neurophysiology), a clustering of resting-state connectivity, and the functional gradients (Margulies et al., 2016, PNAS), a lower dimensional manifold of resting-state connectivity.
Lets first have a look at contextualization of cortical thickness using the Yeo networks. We'll use some of the sample cortical thickness data included with BrainSpace, and see what its mean is within each Yeo network.
% Load cortical thickness
[thickness_left, thickness_right] = load_marker('thickness');
thickness = [thickness_left; thickness_right];
% Load Yeo networks
yeo_7_networks = fetch_parcellation('fslr32k', 'yeo', 7);
[network_names, colormap] = fetch_yeo_networks_metadata(7);
% Get mean intensity within yeo networks
mean_thickness = full2parcel(thickness', yeo_7_networks)';
% Plot
figure;
s = spider_plot_class(mean_thickness, ...
'AxesLabels', cellstr(network_names)', ...;
'AxesLimits', [2; 3.2] .* ones(2,7), ...
'AxesFontSize', 14, ...
'LabelFontSize', 14, ...
'LegendLabels', {'Thickness'});
Here we can see that, on average, the somatomotor/visual cortices have low cortical thickness whereas the default/limbic cortices have high thickness.
Next, lets have a look at how cortical thickness relates to the first functional gradient which describes a sensory-transmodal axis in the brain. First lets plot the first few gradients.
functional_gradients = fetch_gradients('fslr32k', 'margulies2016');
gradients_plot = functional_gradients;
gradients_plot(isnan(gradients_plot)) = -inf;
[surface_left, surface_right] = fetch_template_surface('fslr32k');
obj = plot_hemispheres(...
gradients_plot(:,1:3), ...
{surface_left, surface_right}, ...
'LabelText', "Gradient " + (1:3));
obj.colormaps([.7 .7 .7; parula])
obj.colorlimits([-10, 11; -8, 5; -5, 3])
pretty_plot(obj)
There are many ways to compare these gradients to cortical markers such as cortical thickness. In general, we recommend using corrections for spatial autocorrelation which are implemented in BrainSpace. We'll show a correction with spin test in this tutorial; for other methods and further details please consult the BrainSpace tutorials.
In a spin test we compare the empirical correlation between the gradient and the cortical marker to a distribution of correlations derived from data rotated across the cortical surface. The p-value then depends on the percentile of the empirical correlation within the permuted distribution.
[sphere_left, sphere_right] = load_conte69('spheres');
% Note: we generally recommend running >=1000 permutations. For purposes
% of this tutorial we will only run 100.
thickness_permuted = spin_permutations( ...
{thickness_left, thickness_right}, ...
{sphere_left, sphere_right}, 100, ...
'random_state', 2021);
thickness_permuted_full = squeeze( ...
[thickness_permuted{1}; thickness_permuted{2}]);
r_empirical = corr(thickness, functional_gradients(:, 1), ...
'rows', 'pairwise');
r_permuted = corr(thickness_permuted_full, functional_gradients(:, 1), ...
'rows', 'pairwise');
% Significance depends on whether we do a one-tailed or two-tailed test.
% If one-tailed it depends on in which direction the test is.
p_value_right_tailed = mean(r_empirical > r_permuted);
p_value_left_tailed = mean(r_empirical < r_permuted);
p_value_two_tailed = min(p_value_right_tailed, p_value_left_tailed) * 2;
fprintf(['The two-tailed p-value for the correlation between ' ...
'cortical thickness and \n gradient 1 is %0.2g.\n'], ...
p_value_two_tailed);
That concludes the tutorials of BrainStat. If anything is unclear, or if you think you've found a bug, please post it to the Issues page of our Github.
Happy BrainStating!
function pretty_plot(obj)
% Makes the BrainSpace plots prettier inside live scripts.
obj.handles.figure.Units = 'pixels';
set(obj.handles.axes, 'Units', 'pixels')
set(obj.handles.colorbar, 'Units', 'pixels');
obj.handles.figure.Position = [1, 1, 700, size(obj.data, 2) * 95];
set(obj.handles.text, 'FontSize', 12);
for jj = 1:size(obj.data, 2) % variates
obj.handles.colorbar(jj).Position = [...
500, ...
30 + (size(obj.data, 2) * 95) - 10 - jj*90, ...
5, ...
50];
for ii = 1:4 % views
obj.handles.axes(jj, ii).Position = [...
75 + (ii-1)*105, ...
(size(obj.data, 2) * 95) - jj*95, ...
100, ...
100];
end
end
end
+colortable with 201 entries read (originally Schaefer2018_400Parcels_7Networks)
profiles = read_histology_profile('template', 'fsaverage5');
mpc = compute_mpc(profiles, schaefer_400);
gm = compute_histology_gradients(mpc);
The variable profiles now contains histological profiles sampled at 50 different depths across the cortex, mpc contains the covariance of these profiles, and gm contains their gradients. Depending on your use-case, each of these variables could be of interest, but for purposes of this tutorial we'll plot the gradients to the surface with BrainSpace. For details on what the GradientMaps class, gm, contains please consult the BrainSpace documentation.
vertexwise_gradients = parcel2full(gm.gradients{1}(:,1:2), schaefer_400);
vertexwise_gradients(isnan(vertexwise_gradients)) = inf;
obj = plot_hemispheres(vertexwise_gradients, ...
{surf_left, surf_right}, ...
'labeltext', {'MPC Gradient 1', 'MPC Gradient 2'});
obj.colormaps([parula; .7 .7 .7])
pretty_plot(obj);
Note that we no longer use the y-axis regression used in (Paquola et al., 2019, Plos Biology), as such the first gradient becomes an anterior-posterior- gradient. Next, lets plot a scatter plot between our t-statistic map and the first gradient of MPC.
t_schaefer_400 = full2parcel(slm.t, schaefer_400);
t_mpc_corr = corr(t_schaefer_400, gm.gradients{1}(:,1));
bestfit_coefficients = polyfit(t_schaefer_400, gm.gradients{1}(:,1), 1);
figure();
scatter(t_schaefer_400, gm.gradients{1}(:,1), 200, 'k', '.');
hold on
plot(t_schaefer_400, polyval(bestfit_coefficients, t_schaefer_400), ...
'k', 'LineWidth', 3);
xlabel('t-statistic')
ylabel('MPC Gradient 1')
set(gca, 'FontSize', 16, 'YLim', [-1, 1], 'YTick', [-1, 1], ...
'Xtick', [-2, 3], 'XLim', [-2, 3])
text(-1.9, 0.8, sprintf('r = %0.2f', t_mpc_corr), 'FontSize', 16)
Resting-State Contextualization
Lastly, BrainStat provides contextualization using resting-state fMRI markers: specifically, with the Yeo functional networks (Yeo et al., 2011, Journal of Neurophysiology), a clustering of resting-state connectivity, and the functional gradients (Margulies et al., 2016, PNAS), a lower dimensional manifold of resting-state connectivity.
Lets first have a look at contextualization of our t-statistic map using their mean within the Yeo networks.
% Get mean/sem t-stat within yeo networks
compute_sem = @(x) std(x, 'omitnan') / sqrt(sum(~isnan(x)));
mean_t_stat = yeo_networks_association(slm.t, ...
'template', 'fsaverage5');
sem_t_stat = yeo_networks_association(slm.t, ...
'template', 'fsaverage5', 'reduction_operation', compute_sem);
[network_names, colormap] = fetch_yeo_networks_metadata(7);
% Create a t-statistic bar plot for the Yeo networks.
figure;
axes('XTick', 1:7, 'XTickLabel', network_names, 'FontSize', 16);
hold on
bar(mean_t_stat, 'FaceColor', 'flat', 'CData', colormap);
errorbar(1:numel(mean_t_stat), mean_t_stat, sem_t_stat, sem_t_stat, ...
'Color', 'k', 'LineStyle', 'None');
Here we can see that the effects are largest in the ventral attention and limbic networks, and that the visual network is the only region with a negative t-statistic.
Next, lets have a look at how cortical thickness relates to the first functional gradient which describes a sensory-transmodal axis in the brain. First lets plot the first few gradients and a scatter plot of the t-statistic map versus the first gradient.
functional_gradients = fetch_gradients('fsaverage5', 'margulies2016');
gradients_plot = functional_gradients;
gradients_plot(isnan(gradients_plot)) = -inf;
% Plot the gradients on the surface.
obj = plot_hemispheres(...
gradients_plot(:,1:3), ...
{surf_left, surf_right}, ...
'LabelText', "Gradient " + (1:3));
obj.colormaps([.7 .7 .7; parula])
obj.colorlimits([-10, 11; -8, 5; -5, 3])
pretty_plot(obj)
% Plot a scatter plot
figure();
coefficients = polyfit(slm.t(mask), functional_gradients(mask, 1), 1);
no_nan = ~isnan(slm.t)' & ~isnan(functional_gradients(:,1));
r = corr(slm.t(no_nan)', functional_gradients(no_nan, 1));
scatter(slm.t(no_nan), functional_gradients(no_nan, 1), 10, 'k', '.');
hold on
plot(slm.t(no_nan), polyval(coefficients, slm.t(no_nan)), ...
'k', 'LineWidth', 3);
xlabel('t-statistic')
ylabel('Functional Gradient 1')
set(gca, 'FontSize', 16, 'YLim', [-10, 15], 'YTick', [-10, 15], 'Xtick', [-6, 4])
text(-5.5, 13, sprintf('r = %0.2f', r), 'FontSize', 16)
There are many ways to compare these gradients to cortical markers such as cortical thickness. In general, we recommend using corrections for spatial autocorrelation which are implemented in BrainSpace. We'll show a correction with spin test in this tutorial; for other methods and further details please consult the BrainSpace tutorials.
In a spin test we compare the empirical correlation between the gradient and the cortical marker to a distribution of correlations derived from data rotated across the cortical surface. The p-value then depends on the percentile of the empirical correlation within the permuted distribution.
[sphere_left, sphere_right] = fetch_template_surface('fsaverage5','layer','sphere');
t_left = slm.t(1:end/2)';
t_right = slm.t(end/2+1:end)';
% Note: for the tutorial we run only 100 permutations.
% We generally recommend at least 1000.
t_permuted = spin_permutations( ...
{t_left, t_right}, ...
{sphere_left, sphere_right}, 100, ...
'random_state', 2021);
t_permuted_full = squeeze([t_permuted{1}; t_permuted{2}]);
r_empirical = corr(slm.t', functional_gradients(:, 1), ...
'rows', 'pairwise');
r_permuted = corr(t_permuted_full, functional_gradients(:, 1), ...
'rows', 'pairwise');
% Significance depends on whether we do a one-tailed or two-tailed test.
% If one-tailed it depends on in which direction the test is.
p_value_right_tailed = mean(r_empirical > r_permuted);
p_value_left_tailed = mean(r_empirical < r_permuted);
p_value_two_tailed = min(p_value_right_tailed, p_value_left_tailed) * 2;
fprintf(['The two-tailed p-value between our t-statistic map ' ...
'and \ngradient 1 is %0.2g.\n'], ...
p_value_two_tailed);
That concludes the tutorials of BrainStat. If anything is unclear, or if you think you've found a bug, please post it to the Issues page of our Github.
Happy BrainStating!
function pretty_plot(obj)
% Makes the BrainSpace plots prettier inside live scripts.
obj.handles.figure.Units = 'pixels';
set(obj.handles.axes, 'Units', 'pixels')
set(obj.handles.colorbar, 'Units', 'pixels');
obj.handles.figure.Position = [1, 1, 700, size(obj.data, 2) * 95];
set(obj.handles.text, 'FontSize', 12);
for jj = 1:size(obj.data, 2) % variates
obj.handles.colorbar(jj).Position = [...
500, ...
30 + (size(obj.data, 2) * 95) - 10 - jj*90, ...
5, ...
50];
for ii = 1:4 % views
obj.handles.axes(jj, ii).Position = [...
75 + (ii-1)*105, ...
(size(obj.data, 2) * 95) - jj*95, ...
100, ...
100];
end
end
end