Geometric variation causes functional and aesthetic problems in the assemblies. The challenge is predicting the moments of the distribution of geometric deviations of assemblies to evaluate compliance with the set requirements. Variation analysis of non-rigid sheet metal assemblies with respect to the geometric quality is traditionally performed using Method of Influence Coefficient (MIC), based on Finite Element Method (FEM). The joining operation is one of the most crucial steps in the assembly process, imposing forces on the parts and causes bending and deformation during the assembly, consequently contributing considerably to the final geometric outcome of the assembly. To model the behavior of the assembly realistically and achieve accurate simulation results, considering the sequence of joining is essential. In a digital twin of the assembly process, where the scanned geometry of parts are available and the assembly parameters are optimized for a batch of fabricated components, joining sequences need to be provided for the optimal geometric outcome of the batch of assemblies. The sequence optimization of the joining processes is a time-consuming combinatorial problem to solve. Variation analysis of non-rigid assemblies with stochastic part inputs, including optimal joining sequences, requires an extensive amount of computational effort. More efficient approaches for evaluating assembly geometric variation are desired. In this paper, a computationally efficient approach is proposed for geometric variation analysis and optimization of non-rigid assemblies with stochastic part inputs with respect to the joining sequences. A clustering approach is proposed categorizing the incoming parts based on the part variation. Sequence optimization is performed, and geometric variation is analyzed for each cluster. The results show that the proposed method drastically reduces the computation time needed for sequence optimization compared to individualized optimization for each assembly.