One-way ANOVA for multiple features? I have 60 continuous features (i.e., lab results) for 35,000 patients, and the patients are divided into three groups. I want to evaluate each feature to find out what the differences are between the groups.
Here is my plan:

*

*For each feature, iteratively run a one-way ANOVA independently.

*After 60 ANOVA are run, correct for multiple comparisons using Bonferroni correction.

*Using only the significant features (as per the ANOVA with Bonferroni), run a post-hoc analysis (Tukey's HSD) for each one.

*Use the p-value results from Tukey's HSD to confirm any significant pairwise differences between the groups. Assume there is no need for additional Bonferroni correction at this stage.

Is this an acceptable approach? Or do I need to use a multi-way ANOVA for this? Or MANOVA?

I would suggest a different approach, that yields two features that separate the three groups best. The approach is called "Linear Discriminant Analysis" (LDA) and is a means to extract new features as linear combinations from the original features.
The number of features it extracts are the number of groups minus one, which is two in your case. This has the nice side effect that you can easily visualize the discrimination of the three groups by means of these new features in a two dimensional plot.
Here is an example R code with the iris dataset (3 classe, 4 original features):
> library(MASS)
> res <- lda(iris[,-5], iris$Species)
> scale.norm <- scale(res$scaling, center=FALSE, scale=sqrt(colSums(res$scaling^2)))
> scale.norm
                    LD1          LD2
Sepal.Length  0.2087418  0.006531964
Sepal.Width   0.3862037  0.586610553
Petal.Length -0.5540117 -0.252561540
Petal.Width  -0.7073504  0.769453092
> x.lda <- as.matrix(iris[,-5]) %*% as.matrix(scale.norm)
> plot(x.lda, col=iris$Species)


