Feature selection can disclose biomarkers of mental disorders that have unclear biological mechanisms. Although neighborhood rough set (NRS) has been applied to discover important sparse features, it has hardly ever been utilized in neuroimaging-based biomarker identification, probably due to the inadequate feature evaluation metric and incomplete information provided under a single-granularity. Here, we propose a new NRS-based feature selection method and successfully identify brain functional connectivity biomarkers of schizophrenia (SZ) using functional magnetic resonance imaging (fMRI) data. Specifically, we develop a new weighted metric based on NRS combined with information entropy to evaluate the capacity of features in distinguishing different groups. Inspired by multi-granularity information maximization theory, we further take advantage of the complementary information from different neighborhood sizes via a multi-granularity fusion to obtain the most discriminative and stable features. For validation, we compare our method with six popular feature selection methods using three public omics datasets as well as resting-state fMRI data of 393 SZ patients and 429 healthy controls. Results show that our method obtained higher classification accuracies on both omics data (100.0%, 88.6%, and 72.2% for three omics datasets, respectively) and fMRI data (93.9% for main dataset, and 76.3% and 83.8% for two independent datasets, respectively). Moreover, our findings reveal biologically meaningful substrates of SZ, notably involving the connectivity between the thalamus and superior temporal gyrus as well as between the postcentral gyrus and calcarine gyrus. Taken together, we propose a new NRS-based feature selection method that shows the potential of exploring effective and sparse neuroimaging-based biomarkers of mental disorders.