<p>Final answer. Keep in mind I work in Python so I'm trying to translate here.</p>
<ol>
<li>The output of <code>coef(cv.lassoModel)</code> ... those are your genes of interest thats your <strong>feature selection</strong>. <code>head</code> I assume gives you the top values ... I think thats right (I need to check its not <code>coef(lassoModel)</code>) ... Basically thats where your results are.</li>
</ol>
<p>This might be over simplistic but <code>co.se</code> without the <code>head</code> is your answer. Not every gene will carry a weight BTW. If the entirety of that is (700) <code>.</code> ... try point 2 below and then point 4. I don't think that will happen (see <strong>Whats happened</strong>)</p>
<ol start="2">
<li>I would redo this with <code>alpha=0</code> (ridge regression) and <code>alpha=0.5</code> elasticnet, but you don't need to</li>
<li>Generate an accuracy score</li>
</ol>
<pre><code>    y_predicted &lt;- predict(best_model, s = best_lambda, newx = x)
    sst &lt;- sum((y - mean(y))^2)
    sse &lt;- sum((y_predicted - y)^2)
    rsq &lt;- 1 - sse/sst
    rsq
</code></pre>
<p>Hopefully this is large e.g. around 0.8 (in Python this automatic BTW). You could ignore that, but its what I'd do</p>
<ol start="4">
<li>The intercept is simply the weight when your gene expression score is 0.</li>
<li>Lasso regression handles sparse data - so no worries about <code>&quot;.&quot;</code>. I don't whether that Python's version of &quot;nan&quot; (not a number), or if its zero.  ... @haci knows what this means and his answer is authoritative for the second part of your question.</li>
<li>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.</li>
</ol>
<p>In Python what you want is feature selection
<code>print(pd.DataFrame({'Features': X_train.columns, 'Coefficients': lasso.coef_}))</code></p>
<p><strong>Whats happened?</strong> The coefficients for <code>head(co.se)</code> have been shrunk to zero, thats only the top 10? results, you've 700 weights. The full list <em>should</em> (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.</p>
<hr />
<p>Background</p>
<p><strong>Lasso regression</strong>, 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.</p>
<p><strong>Imbalance</strong>, 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.</p>
<p><strong>Imbalance</strong> is when it's 1 in 10 and over. You've got 1 in 6, 1 in 7, not huge. However, - <em>you might be right</em> - and if you're worried simply repeat the analysis excluding <code>RP</code> (value 7). Then check whether the parameterisation plot (lambda against weight) is more stable or more &quot;homogeneous&quot; (multinomial plot), i.e. the genes behaviour in a more uniform way. If that is does work - it gives you nicer results AND <code>RP</code> 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.</p>
<p><strong>Transformation</strong> 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 <em>max-min</em> transformation - thats where I might start.</p>
<p><strong>Overall good results</strong> 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 <em>red</em> (top), <em>blue</em> (bottom) and cyan (negative y-axis) are real results showing that within the 700 genes <strong>3 genes</strong> 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. <strong>From 700 genes down to 5 genes</strong>. Whats wrong with that?</p>
<p><strong>Whats the answer to the data set?</strong></p>
<p>Its ...</p>
<blockquote>
<p>NRT
RT
NRP
RP</p>
</blockquote>
<p>... is due to 'red', 'blue', 'cyan', 'black' and 'deep pink' genes. Thats it.</p>
<p><strong>Frank assessment</strong> To be perfectly honest, looks good. No-one is going to say &quot;Wooaaa lambda ain't stable, destroy! ;-)&quot;, they'll say &quot;Wow thats cool and clear&quot;.</p>
<hr />
<p><strong>How I'd do it..</strong> 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.</p>
<p><strong>Note</strong> I use Python not R for ML, both work well however. <code>glmnet</code> within ML is new to me BTW.</p>

CC BY-SA 4.0
2023-07-22T09:02:50.540Z
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.
Angelo
CC BY-SA 4.0
2023-07-22T09:43:27.533Z
Another puzzling aspect is when I do co.se$RT!=0, I get a gene list with "." next to it :(
Angelo
CC BY-SA 4.0
2023-07-22T11:14:58.627Z
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)`.
haci
CC BY-SA 4.0
2023-07-22T11:54:27.233Z
All of them :( I tried to normalise the values by getting zscores. HES4 is a gene name and its reading in different samples.
Angelo
CC BY-SA 4.0
2023-07-22T13:59:23.700Z
@Angelo I've addressed the comments in the first part of the answer.
M__
CC BY-SA 4.0
2023-07-22T04:39:07.903Z
2023-07-22T16:23:48.990Z
2023-07-22T16:23:48.990Z
M__
<p>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.</p>
<p>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.</p>
<p>so my data looks like this</p>
<pre class="lang-r prettyprint-override"><code>           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
</code></pre>
<p>purely in terms of numbers, the subtypes have following number of samples.</p>
<pre><code>NRT 52
RT 42
NRP 30
RP 7
</code></pre>
<p>Q1. May be the groups are imbalanced and caution should be taken?</p>
<p>I have applied glmnet as following</p>
<pre class="lang-r prettyprint-override"><code>modellinglasso=read.table(&quot;LRMatrix.txt&quot;, header=T, row.names = 1, sep=&quot;\t&quot; )
modellinglasso$Subtype=as.factor(modellinglasso$Subtype)
modellinglasso$Subtype
#build a single pass (single fold) lasso-penalised model
lassoModel &lt;- glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  family=&quot;multinomial&quot;)
plot(lassoModel, xvar=&quot;lambda&quot;)
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
</code></pre>
<p><a href="https://i.sstatic.net/0lUP0.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/0lUP0.png" alt="enter image description here" /></a></p>
<pre class="lang-r prettyprint-override"><code>#perform 10-fold cross validation
cv.lassoModel &lt;- cv.glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  nfolds=10,
  family=&quot;multinomial&quot;,
  parallel=TRUE)
# plot variable deviances vs. shrinkage parameter, λ (lambda)
plot(cv.lassoModel)
</code></pre>
<p><a href="https://i.sstatic.net/E02wY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/E02wY.png" alt="enter image description here" /></a></p>
<pre class="lang-r prettyprint-override"><code>idealLambda &lt;- cv.lassoModel<span class="math-container">$lambda.min
idealLambda1se &lt;- cv.lassoModel$</span>lambda.1se
print(idealLambda); print(idealLambda1se)

[1] 0.1186125
[1] 0.1186125

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

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

# identify predictors for each sub-type
rownames(co.se<span class="math-container">$PR)[which(co.se$</span>NRT != 0)]

[1] &quot;(Intercept)&quot;
</code></pre>
<p>output for</p>
<pre><code>co.se$RT!=0
6 x 1 sparse Matrix of class &quot;lgCMatrix&quot;
            1
(Intercept) |
HES4        .
FNDC10      .
PRKCZ       .
PRKCZ.AS1   .
MMEL1       .
</code></pre>
<p>It doesn't matter which subtype I put, I get the same output &quot;(Intercept)&quot; instead of gene names.</p>
<p>Any help to get the gene signatures for the subtypes is sincerely appreciated.</p>
<p>Thank  you</p>

<r><rna-seq><machine-learning>
Applying glmnet to identify predictors for subtypes
# Applying glmnet to identify predictors for subtypes

(Asked by: username_0 on 2023-07-21)

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


#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)


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

**Tags:** r, rna-seq, machine-learning

---

## 1 answers

### Answer by username_1 on 2023-07-22. Score: 2 | Accepted answer

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

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)

I would redo this with alpha=0 (ridge regression) and alpha=0.5 elasticnet, but you don't need to
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

The intercept is simply the weight when your gene expression score is 0.
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.
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.

> **Comments:**
> **username_0**: 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.  
> **username_0**: Another puzzling aspect is when I do co.se$RT!=0, I get a gene list with "." next to it :(  
> **username_2**: 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)`.  
> **username_0**: All of them :( I tried to normalise the values by getting zscores. HES4 is a gene name and its reading in different samples.  
> **username_1**: @Angelo I've addressed the comments in the first part of the answer.  

---

