As Bogumil says, it's not ideal to ask many questions in one post on StackOverflow - your question should be well defined and targeted, ideally with a minimum working example to make it easiest for people to help you.
Let's therefore answer what from the title of the post I take as your main question: how can I fit a linear regression model with GLM which includes many response columns. That question is almost a duplicate of this question, so a very similar answer applies: broadcast the term function over the names in your DataFrame which you want to include on the right hand side like this:
julia> using DataFrames, GLM
julia> df = hcat(DataFrame(RESP = rand(100)), DataFrame(rand(100, 10), :auto));
julia> mymodel = lm(term(:RESP) ~ sum(term.(names(df[!, Not(:RESP)]))), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}}}}, Matrix{Float64}}
RESP ~ 1 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
(...)
You also ask about plotting, but that's probably best dealt with in a separate question. You can access the estimated coefficients of your linear model using the coef function, e.g.
julia> coef(mymodel)
11-element Vector{Float64}:
0.504236533528822
0.11712812154185266
-0.0206810430546413
-0.15693089456050294
-0.011916514466331067
-0.1030171434361648
0.10378957999147352
-0.09447743618381275
-0.08860977078650123
0.0816071818033377
0.09939548661830626
and the full output with coeftable.
Finally note that you won't necessarily see from that which columns "impact most" your model as you say in the comment, unless you have standardized your regressors.
