In this big data era, interpretable machine learning models are strongly demanded for the comprehensive analytics of large-scale multiclass data. Characterizing all features from such data is a key but challenging step to understand the complexity. However, existing feature selection methods do not meet this need. In this paper, to address this problem, we propose a Bayesian multiclass nonnegative matrix factorization (MC-NMF) model with structured sparsity that is able to discover ubiquitous and class-specific features. Variational update rules were derived for efficient decomposition. In order to relieve the need of model selection and stably describe feature patterns, we further propose MC-NMF with stability selection, an ensemble method that collectively detects feature patterns from many runs of MC-NMF using different hyperparameter values and training subsets. We assessed our models on both simulated count data and multitumor ribonucleic acid-seq data. The experiments revealed that our models were able to recover predefined feature patterns from the simulated data and identify biologically meaningful patterns from the pan-cancer data.