With the advance of new sequencing technologies producing massive short reads data, metagenomics is rapidly growing, especially in the fields of environmental biology and medical science. The metagenomic data are not only high-dimensional with a large number of features and limited number of samples, but also complex with skewed distribution. Efficient computational and statistical tools are needed to deal with these unique characteristics of metagenomic data. In metagenomic studies, one main aim is to assess whether and how the multiple microbial communities differ under various environmental conditions. In this research we propose a two-stage statistical procedure for selecting informative features and identifying differentially abundant features between two or more microbial communities. In the functional analysis of metagenomes the features may refer to the pathways, subsystems, functional roles, and so on. Compared with other available methods the proposed approach demonstrates better performance via comprehensive simulation studies. Applications to real data will also be presented.