We propose a new data-driven method to solve multiscale elliptic PDEs with high-dimensional random coefficients, where both the local problem and global problem are considered. Our method consists of offline and online stages. In the offline stage, we construct a small number of data-driven basis functions either on subdomains or on the whole computational domain, which can be used to approximate the multiscale random solutions. In the online stage, with the help of the data-driven basis, we can quickly compute the multiscale solutions for each new coefficient. We provide some analysis for the proposed method, which provides some guidance on how to determine the number of learning samples and the truncation threshold in building the data-driven basis. Finally, we present numerical examples to demonstrate the accuracy and efficiency of the proposed method.