This work is concerned with feature screening in the situation when both the responses and predictors are high-dimensional with adjusting for the confounders. We develop a screening procedure and show that it has a sure screening property under some mild technical conditions. Three prominent merits of the proposed method include none requirement of model specification, high-dimensionality of the responses and non-euclidean data. The Monter Carlo results and application to mice pleiotropic genetic association study demonstrate that the proposed procedures work better than two existing methods.