3
$\begingroup$

I am new to glmnet and other ML techniques, so apologise in advance if this sounds a trivial question. However, your guidance is very much appreciated.

So, I have nearly 700 genes with their expression profile across 4 conditions. I want to identify signature genes for each of the 4 conditions. I am using glmnet but I am not getting anywhere near the results.

so my data looks like this

           Subtype     HES4   FNDC10    PRKCZ PRKCZ.AS1
Sample1     NRT 5.452955 6.241387 8.546059  6.260245
Sample2     RT 5.125350 7.038719 8.900024  6.267325
Sample3     NRT 5.689560 5.903811 8.029434  4.213996
Sample4     NRP 5.834468 5.722589 8.061339  5.072274
Sample5     RP 5.401623 7.272111 9.108342  6.599629

purely in terms of numbers, the subtypes have following number of samples.

NRT 52
RT 42
NRP 30
RP 7

Q1. May be the groups are imbalanced and caution should be taken?

I have applied glmnet as following

modellinglasso=read.table("LRMatrix.txt", header=T, row.names = 1, sep="\t" )
modellinglasso$Subtype=as.factor(modellinglasso$Subtype)
modellinglasso$Subtype
#build a single pass (single fold) lasso-penalised model
lassoModel <- glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  family="multinomial")
plot(lassoModel, xvar="lambda")
Warning message:
In lognet(xd, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :
  one multinomial or binomial class has fewer than 8  observations; dangerous ground

enter image description here

#perform 10-fold cross validation
cv.lassoModel <- cv.glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  nfolds=10,
  family="multinomial",
  parallel=TRUE)
# plot variable deviances vs. shrinkage parameter, λ (lambda)
plot(cv.lassoModel)

enter image description here

idealLambda <- cv.lassoModel$lambda.min
idealLambda1se <- cv.lassoModel$lambda.1se
print(idealLambda); print(idealLambda1se)

[1] 0.1186125
[1] 0.1186125

# derive coefficients for each gene
co <- coef(cv.lassoModel, s=idealLambda, exact=TRUE)
head(co)

co.se <- coef(cv.lassoModel, s=idealLambda1se, exact=TRUE)
head(co.se)

# identify predictors for each sub-type
rownames(co.se$PR)[which(co.se$NRT != 0)]

[1] "(Intercept)"

output for

co.se$RT!=0
6 x 1 sparse Matrix of class "lgCMatrix"
            1
(Intercept) |
HES4        .
FNDC10      .
PRKCZ       .
PRKCZ.AS1   .
MMEL1       .

It doesn't matter which subtype I put, I get the same output "(Intercept)" instead of gene names.

Any help to get the gene signatures for the subtypes is sincerely appreciated.

Thank you

$\endgroup$
0

1 Answer 1

2
$\begingroup$

Final answer. Keep in mind I work in Python so I'm trying to translate here.

  1. The output of coef(cv.lassoModel) ... those are your genes of interest thats your feature selection. head I assume gives you the top values ... I think thats right (I need to check its not coef(lassoModel)) ... Basically thats where your results are.

This might be over simplistic but co.se without the head is your answer. Not every gene will carry a weight BTW. If the entirety of that is (700) . ... try point 2 below and then point 4. I don't think that will happen (see Whats happened)

  1. I would redo this with alpha=0 (ridge regression) and alpha=0.5 elasticnet, but you don't need to
  2. Generate an accuracy score
    y_predicted <- predict(best_model, s = best_lambda, newx = x)
    sst <- sum((y - mean(y))^2)
    sse <- sum((y_predicted - y)^2)
    rsq <- 1 - sse/sst
    rsq

Hopefully this is large e.g. around 0.8 (in Python this automatic BTW). You could ignore that, but its what I'd do

  1. The intercept is simply the weight when your gene expression score is 0.
  2. Lasso regression handles sparse data - so no worries about ".". I don't whether that Python's version of "nan" (not a number), or if its zero. ... @haci knows what this means and his answer is authoritative for the second part of your question.
  3. In Python transformations are common, but that doesn't seem an option here, so find providing your accuracy score is reasonable, personally I'm happy.

In Python what you want is feature selection print(pd.DataFrame({'Features': X_train.columns, 'Coefficients': lasso.coef_}))

Whats happened? The coefficients for head(co.se) have been shrunk to zero, thats only the top 10? results, you've 700 weights. The full list should (famous last words) contain a minority of really big coefficients ... those are really important and those are your answer. So by printing only the top few results ... not the whole lot the results are still there you've just not seen them. Thats it, basically you were almost there and you could happily settle on that as the final answer. Basically, alot of genes ain't doing much a few genes are doing all the work - thats the answer. Below I stress repeatedly I don't think anything is wrong here and I personally like the result.


Background

Lasso regression, it's a cousin of ridge regression and key parameter is lambda (who'd guessed). What your plots are doing is tuning the parameter and finding the associated weight individually or collectively. Fairly straight forward. So like whats the value of lambda for your data set ('answered' below)? It's worth doing some reading on how lambda is calculated to understand the plot, if you're not sure how it works.

Imbalance, well to be honest I'm not convinced there's anything wrong with your analysis. Okay there's a lack of stationarity for some of your weights ... but some are fine and gist is pretty clear. Yeah okay there's a problem with stabilisation of lambda - thats not necessarily imbalance.

Imbalance is when it's 1 in 10 and over. You've got 1 in 6, 1 in 7, not huge. However, - you might be right - and if you're worried simply repeat the analysis excluding RP (value 7). Then check whether the parameterisation plot (lambda against weight) is more stable or more "homogeneous" (multinomial plot), i.e. the genes behaviour in a more uniform way. If that is does work - it gives you nicer results AND RP is essential to the analysis thats a separate question about augmentation. That is complicated. For example, the multinomial plot with a 3 target categorical model might give a clear value for lambda.

Transformation Here, I would suggest if you want to stabilise the multinomial plot (if it's possible) I would look at a transformation. The data you presented doesn't look like it needs a transformation, but we can't see the whole data. In Python there's a max-min transformation - thats where I might start.

Overall good results Personally, I'd say these were good results, except the weights ain't homogeneous for some genes (most are fine BTW), but the patterns clear. The red (top), blue (bottom) and cyan (negative y-axis) are real results showing that within the 700 genes 3 genes are having a massive impact. Sure you have positive and negative and thats how it works one gene is positively contributing to the classification, others are negatively contributing to it. Then you've the black and 'deep pink' lines which have stabilised around 5E4. From 700 genes down to 5 genes. Whats wrong with that?

Whats the answer to the data set?

Its ...

NRT RT NRP RP

... is due to 'red', 'blue', 'cyan', 'black' and 'deep pink' genes. Thats it.

Frank assessment To be perfectly honest, looks good. No-one is going to say "Wooaaa lambda ain't stable, destroy! ;-)", they'll say "Wow thats cool and clear".


How I'd do it.. I would look at basic ML robustness statistics and use this to assess the model rather than the parameterisation of lambda. Its horses for courses as they say. Upvotes.

Note I use Python not R for ML, both work well however. glmnet within ML is new to me BTW.

$\endgroup$
5
  • $\begingroup$ Thank you for your reply and detailed explanation. What I do not understand is why I am getting output as Intercept instead of gene names. $\endgroup$
    – Angelo
    Jul 22, 2023 at 9:02
  • 1
    $\begingroup$ Another puzzling aspect is when I do co.se$RT!=0, I get a gene list with "." next to it :( $\endgroup$
    – Angelo
    Jul 22, 2023 at 9:43
  • 2
    $\begingroup$ Technically you are getting a sparse Matrix with co.se$RT!=0 and the "dots" are representing zeros or in this specific case FALSE values. You can see yourself by as.matrix(co.se$RT!=0). $\endgroup$
    – haci
    Jul 22, 2023 at 11:14
  • $\begingroup$ All of them :( I tried to normalise the values by getting zscores. HES4 is a gene name and its reading in different samples. $\endgroup$
    – Angelo
    Jul 22, 2023 at 11:54
  • 1
    $\begingroup$ @Angelo I've addressed the comments in the first part of the answer. $\endgroup$
    – M__
    Jul 22, 2023 at 13:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.