The brain activity is organized in clusters/modules that have different roles in our behavior (Tononi et al., 1999). Alterations in the clustering pattern can be associated with neurologic disorders (Grossberg, 2000; Sato et al., 2016). Thus, it is important to systematically discriminate the clustering structures among different populations. This leads to the problem of how to statistically test the equality of clustering structures of two or more populations and how to identify the features that contribute to the differential clustering structure. These statistical problems were recently solved for a large class of clustering algorithms by using the Analysis of Cluster Variability—ANOCVA (Fujita et al., 2014a).

Here, we provide an implementation of ANOCVA in R for a better dissemination of this technique in the scientific community. ANOCVA was designed to test whether the clustering structures of several populations are equal. Briefly, ANOCVA uses the silhouette statistic (Rousseeuw, 1987) as a measure of variability of the clustering structure of each population and then compares the variability among populations using an idea similar to the classical analysis of variance (ANOVA). To calculate the statistical significance value, we use a bootstrap procedure that was previously shown to control the type I error.

We illustrate the step-by-step application of ANOCVA by analyzing a large functional magnetic resonance imaging (fMRI) data acquired under a resting-state protocol (ABIDE—The Autism Brain Imaging Data Exchange Consortium) composed of 529 controls and 285 patients diagnosed with autism. Subjects with Autism Spectrum Disorders (ASD) have significant differences in the resting state functional connectivity when compared to healthy subjects (for review, see Kana et al., 2011), suggesting that ASD is as a neural systems disorder with disruptions in several distributed neurocognitive networks of brain regions (Ecker et al., 2015). However, most studies describe integration (Washington et al., 2014; Sporns and Betzel, 2016) and segregation (Assaf et al., 2013) as separate processes. Instead, in this study we consider both processes simultaneously using the idea of clusters, where structures within are integrated and structures between are segregated.

To formalize ANOCVA, we will first describe the silhouette statistic to define “clustering variability” and then we introduce the ANOCVA. Finally, we describe its implementation and application to ABIDE dataset.

The problem that we want to solve is the following. Given k populations T₁, T₂, …, T_k where each population T_j (j = 1, …, k) is composed of n_j subjects, and each subject has N items that are clustered, we would like to verify whether the clustering structures of the brain networks of the k populations are equal and, if not, which ROIs are differently clustered. In our case, we have k = 2 populations with T₁ and T₂ as controls and ASD, respectively. The number of subjects in each population is n₁ = 529 and n₂ = 285, for T₁ and T₂, respectively. The number of ROIs (items) to be clustered is N = 316. Since head movement during magnetic resonance scanning may affect statistical analysis, ANOCVA was applied to both “scrubbed” and “not scrubbed” data (Power et al., 2012) with the number of bootstrap samples set to 1000.

The first step in ANOCVA analysis is the construction of the average dissimilarity matrix A= and its clustering. The estimated number of clusters by the silhouette criterion was five as depicted in Figure 2. Notice that the highest silhouette statistic was obtained when the number of clusters is five. The sub-networks obtained by applying the spectral clustering on the dissimilarity matrix A= can be visualized in Figure 3 where each color represents one sub-network (cluster).

Then, ANOCVA calculates the silhouette statistic for each ROI by using the labels obtained by clustering the dissimilarity matrix A= and performs the test. We verified that in fact the entire clustering structure of subjects diagnosed with ASD differs from controls (p < 0.001). Next, we tested each ROI to identify which ones significantly contribute to the differential clustering between controls and subjects diagnosed with ASD. ROIs that presented a difference in p > 5% between “scrubbed” and “not scrubbed” datasets were excluded for subsequent analysis. Remaining p-values were corrected for multiple comparisons by the Bonferroni method. Figure 4 illustrates the statistically significant ROIs at a p-value threshold of 0.05 after Bonferroni correction. The highlighted regions include portions of the cerebellum and middle frontal gyrus, pre- and post-central gyri, inferior temporal gyrus, and lateral occipital cortex.

In the current study, we combined spectral clustering analysis with ANOCVA implemented in R to investigate which brain regions are clustered in a different way between controls and ASD groups. Our results suggest that several regions distributed across different neurocognitive systems significantly contributed to the different clustering network structure observed in ASD. First we demonstrated that the spectral clustering method yielded partitions that were well-characterized as functional modules of the brain that have been consistently identified in previous studies using different approaches (Damoiseaux et al., 2006; Power et al., 2011), including the fronto-temporal, sensorimotor, visual, and cerebellar systems. This is consistent with the hypothesis that the spectral clustering algorithm groups anatomically contiguous and also spatially distributed areas with common brain functionalities in the same cluster. Then, using ANOCVA we showed that the superior division of the lateral parietal cortex, precentral, and postcentral gyri, anterior dorsal middle frontal gyrus, and a medial portion of the cerebellum and of the brainstem have a distinct cluster organization between ASD and controls. All these brain regions have been previously identified as presenting ASD-related differences in studies using functional MRI. For example, the recruitment of portions of the precentral and postcentral gyri as well as the cerebellum across sensorimotor tasks are atypical in ASD, and may underlie deficits in fine motor sequencing and visual motor learning observed in autistic individuals (Müller et al., 2001; Mostofsky et al., 2009). Interestingly, these regions have also been implicated in cognitive process crucial for interpersonal interactions such as theory-of-mind (Martineau et al., 2010; Wang et al., 2014). This suggests that these areas are involved in the social communication deficits that are a core clinical feature of ASD. Moreover, the lateral parietal cortex is an important node of the default-mode network, and abnormalities in the connectivity between nodes of this network have been widely investigated in ASD (Kennedy and Courchesne, 2008; Assaf et al., 2010; Weng et al., 2010) giving its associations with social cognition (Buckner et al., 2008). The identification of these regions by our study therefore confirms that they are key brain structures in ASD that may have a role in the development of sub-networks organization in this population.

Head motion is one of the most challenging obstacles in functional connectivity studies involving clinical populations, which usually present high levels of movement. Our attempt to handle this problem was to apply the scrubbing method proposed by Power et al. (2012), which discards scans acquired under excessive head motion. However, although this approach may reduce the influence of movement artifacts, they may still be present in the scrubbed data. Thus, we opted for a more conservative approach, which consisted in excluding the regions where the p-values were more sensitive to scrubbing. We assumed that the analyses of these regions were more vulnerable to artifacts and thus they were removed. This approach is also helpful to reduce the number of multiple comparisons, by excluding the less reliable tests. Another important limitation to be mentioned is that the ABIDE data is multicentric with heterogeneous acquisition parameters across sites. We minimized the site effect by removing it in the pre-processing stage of the data. Finally, all analyses are based on the CC400 atlas (Craddock et al., 2012), obtained by using a functional parcellation. Since other atlases are different on ROIs size, number of ROIs and spatial location, the parcellation choice is expected to influence our findings. However, this variability does not invalidate the results obtained with CC400 because the procedures adopted here are conservative (regarding type I error control). Finally, an important future question for the presented results is whether the contribution of these specific brain regions to a differential network clustering in ASD is static or may exhibit dynamic changes during rest (Hutchison et al., 2013).

MV, JS, DT, and AF designed the work. MV pre-processed and analyzed the data. JS and JB interpreted the results. All authors drafted the work, read, and approved the final version of the manuscript.

MV was supported by CAPES and CNPq Fellowships. JS was supported by State of São Paulo Research Foundation—FAPESP (#2013/10498-6). DT was partially supported by Pew Latin American Fellowship and Ciência sem Fronteiras Fellowship (CNPq #246778/2012-1). AF was partially supported by FAPESP (#2014/09576-5, #2013/01715-3, #2015/01587-0, #2016/13422-9, and #2013/03447-6), CNPq (#304020/2013-3 and #473063/2013-1), and NAP eScience—PRP—USP.