<p>tl:dr use limma (not voom, that's for raw counts) and refer to previous questions:</p>
<ul>
<li><a href="https://support.bioconductor.org/p/98820/#98875" rel="nofollow noreferrer">https://support.bioconductor.org/p/98820/#98875</a></li>
<li><a href="https://support.bioconductor.org/p/125144/#125156" rel="nofollow noreferrer">https://support.bioconductor.org/p/125144/#125156</a></li>
</ul>
<p>The problem you have is not so much the threshold but the method to infer any metric to be used for this threshold. The arguably most common methods in R for DEG analysis are DESeq2, edgeR and limma, based on number of citations. The former two require raw counts which you apparently do not have, but if you do then use this with any of the mentioned tools.</p>
<p>limma is what is left and could be used here with TPM. Note that TPM has the problem that counts are divided by a length factor. A brief background: Any statistical test aims to decide whether the given variance between two groups (basically the logFC) is greater than expected by chance, which is the variance within group. In RNA-seq the variance is not equal across the range of counts (low, middle, high, for simplicity). On log2 count scale variance is greater for genes with low(er) counts than for higher counts. Usually, using raw counts the method (DESeq2, edgeR, limma, and others) learn from the data what is the expected variance at the given average expression level, and incorporate that into the model to make adequate inference even at small sample size. That is by the way where the common assumption &quot;most genes do not change&quot; comes from. Assuming that most genes are not DE, then the mean-variance trend represents the expected variance at the given average count range, and this can then be used to help decide whether a gene is DE or not. With TPM you basically throw this out of the window because the actual counts (as they come from the sequencer) are manipulated by dividing by length and the entire mean-variance model gets distorted. Another factor is that TPM, which corrects by total library size can fail to properly normalize counts, I wrote a summary here:</p>
<ul>
<li><a href="https://www.biostars.org/p/9465851/#9465854" rel="nofollow noreferrer">https://www.biostars.org/p/9465851/#9465854</a></li>
</ul>
<p><strong>Long story short: TPM is not good for DEG.</strong></p>
<p>If it's all you have then use what is called &quot;limma-trend&quot;, see the limma manual at Bioconductor, which will directly model the log2-TPMs and give some logFC and FDR. You can then use standard cutoffs, such as FDR &lt; 0.05 &amp; logFC &gt;/&lt; 1.5. These cutoffs are not perfect, but common, and given that TPMs are bad anyway you can and should stick to standard cutoffs as much as possible, rather than infering your own ones. This at least gives some standardization to this suboptimal analysis.</p>

CC BY-SA 4.0
2023-06-11T22:23:32.780Z
2023-06-11T22:23:32.780Z
ATpoint
<p>Two solutions - the underlying issue is called Type 1 error (below).</p>
<ol>
<li>If you are performing studies on humans then the Benjamini-Hochberg method is without question the method of choice, because this includes the <em>false detection</em> rate for EdgeR and DeSeq2. Thats been independently calculated.</li>
</ol>
<p>Thus, the paper that caused a stir was <a href="https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02648-4" rel="nofollow noreferrer">Li <em>et al</em> (2022)</a> who wrote,</p>
<blockquote>
<p>Particularly, the actual FDRs of DESeq2 and edgeR sometimes exceed 20% when the target FDR is 5%.</p>
</blockquote>
<p><em><strong>Benjamini-Hochberg (BH) method solves that for humans (EdgeR and DeSeq2)</strong></em>, because the problem with the BH method is you need a false discovery rate, Li <em>et al</em> 2022, calculated that its about 20% ... so for EdgeR and DeSeq2 the problem is solved. The BH tests will simply use the estimate for the FDR by Li <em>et al</em> (2022) instead of the FDR by EdgeR and DeSeq2. Cool? It is :-) it wasn't what Li <em>et al</em> (2022) intended but it is the solution to a difficult conundrum (concerning a statistical issue or parametrics vs. non-parametrics). <strong>So yeah it's cool</strong> (needs count data).</p>
<p>To perform BH you simply rank the observed probabilities (e.g. from Limma), and there are significant if each is below (or equal to) the BH rank. The BH rank is the total number of tests as denominator, the rank position as numerator and this is multiplied by the false detection rate as a frequency.</p>
<p><em>HOWEVER</em>, you are using Limma and I don't how the model works. I know the underlying stats behind EdgeR, DeSeq2 and NOISeq.</p>
<hr />
<p>Okay so the false discovery rate is not known ... in the way Li <em>et al</em> (2022) calculated it. What now? Limma is likely to calculate an FDR, you could use that via BH. My system is so far from humans there's no way Li <em>et al</em> (2022) FDR is relevant.</p>
<ol start="2">
<li>Bonferoni's approach - I don't know how mainstream this is, but its what I use (I've a background in classical stats), because I've no idea what the true FDR is. This approach can be exemplified via Volcano plots is to visualise the Type 1 error issues.</li>
</ol>
<p>Here's a volcano plot which permits online viewing via links <a href="https://www.htgmolecular.com/blog/2022-08-25/understanding-volcano-plots" rel="nofollow noreferrer">obtained here</a>. I can't change the size BTW (Bioconductor and loads more BTW don't permit this),</p>
<p><a href="https://i.sstatic.net/mFeN2.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/mFeN2.jpg" alt="enter image description here" /></a></p>
<p>This is using a flat probability - which in the picture is very low - it's the horizontal threshold (<em>y</em>-axis). So the P value is one parameter approach to significance ....</p>
<p>The other parameter is the log<sub>2</sub> threshold ...</p>
<p>Here's a two-parameter approach obtained from <a href="https://www.genekitr.fun/plot-others-1.html" rel="nofollow noreferrer">GeneKitR's site</a>, who do permit images to be link viewed (thank you genekitr).</p>
<p><a href="https://i.sstatic.net/gnq9l.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gnq9l.jpg" alt="enter image description here" /></a></p>
<p>Here the authors have used log<sub>2</sub> as the threshold (<em>x</em>-axis). The threshold is probably log<sub>2</sub> &gt; 1 and log<sub>2</sub> &lt; -1, which is reasonable. They have set the P-value threshold to be really low, looks to be about 5%. The plot reveals the effective P-value is 1E-2.5, which is low, but nothing really starts happening until 1E-3, which ... okay. Although this is a two-parameter approach the guys really just relied on log<sub>2</sub>.</p>
<p><strong>Type 1 error and a 1-parameter method threshold</strong></p>
<p>The massive problem is the shear number of tests performed and this will result in significant values PURELY BY CHANCE ...artefacts of critical probability. BH is one approach, Bonferoni's is an alternative way in.</p>
<p>Put simply for your data to calculate the Bonferoni threshold if a 5% critical probability is sought then its ...</p>
<blockquote>
<p>0.05 / 19000 = 2.6e-6</p>
</blockquote>
<p>Which is okay. So any single result &lt;2.6E-6 is significant. This okay, but you will lose results so its vulnerable to Type 2 error. Why do you lose stuff ...</p>
<p><strong>Type 2 error and 1-parameter method threshold</strong></p>
<p>The biological question is e.g. is 0.1log<sub>2</sub> biological sensible? What does it actually mean? Honestly I don't have a clue ... and the problem with a 1-parameter approach, it's very conservative but infers weird stuff. An expression difference of 10% isn't really serious, but thats what it's measuring. ... because its measuring this it means the probability via Bonferoni's is set high - <em>for meaningless stuff</em> and you can be missing important shifts in log<sub>2</sub>, e.g. &gt;1 because they are falling under the high threshold.</p>
<p><strong>Type 1 error and a 2-parameter method thresholds</strong></p>
<p>This is where it all starts getting interesting. As soon as a +/- log<sub>2</sub> threshold is set the number of total tests will drop A LOT. That means the P-value threshold raises. For example, if selecting &gt; (or &lt;) +/- 1log<sub>2</sub> as used in the above example, supposing the number of observations above/below the +/- 1log<sub>2</sub> was 5000. The P-value now raises via Bonferoni, e.g.</p>
<blockquote>
<p>0.05/5000 = 1e-05</p>
</blockquote>
<p>That looks better there'll be more results. Supposing we set a threshold of  +/- 1.5log<sub>2</sub> and obtained 500 results ... then the new Bonferoni derived threshold is,</p>
<blockquote>
<p>0.05/500 = 1e-04 (0.0001)</p>
</blockquote>
<p>That looks fabulous giving loads of significant results ... Wait isn't this all getting a bit crazy, because it depends on the biological assumption of what a reasonable log<sub>2</sub> is. Yep, indeed, agreed.</p>
<p><strong>Type 1 error and a 2-parameter method thresholds: weaknesses</strong></p>
<p>The issue is that a reasonable log<sub>2</sub> is sensible, because there's a point where log<sub>2</sub> expression differences becoming utterly meaningless. However there's a point up to -/+ 2.5 log<sub>2</sub> (which is very high) where you could choose anything - that shifts the critical probability and that shifts how many results are positive.</p>
<p>So how do you ...</p>
<p><strong>Set a log<sub>2</sub> within a 2-parameter threshold?</strong></p>
<p>I've just solved that. I was really happy. The practical answer is it's complicated. My contractual basis is I can't reveal that because public disclosure isn't permitted without top-level approval.</p>
<p>The approach is to use an independent method to calculation DEG  (if you check my posts you'll see I do AI). You can then perform parameter tuning of DEG against the independent method and statistically assess the association. In my case I was surprised by the answer - intuitively I wouldn't have gone there, but my system ain't humans.</p>
<p><strong>Conclusion</strong></p>
<p>In your case,</p>
<ul>
<li>simply check the literature to see what other log<sub>2</sub>  thresholds are in use, use that and adjust the P-value threshold accordingly as described above and cite that work.</li>
<li>look at the biological reality of the results at different log<sub>2</sub> thresholds. For example, an operon should behave like an operon. If the log<sub>2</sub> threshold is smashing-up that operon - that doesn't make sense and adjust log<sub>2</sub> threshold by lowering.</li>
<li>alternatively use BH, this however assumes Limma's FDR calculation is correct and the problem is there no independent study examining that.</li>
</ul>
<p>Either of the above 3 approaches will be fine your context. This is not cutting-edge statistical methodology it simply needs to be reasonable.</p>
<p><strong>Ultra-conservative parameters and Type 2 error</strong></p>
<p>To avoid all the complexities some investigators simply use &gt; (or &lt;) +/- 2.5log<sub>2</sub> and P-value = 1e-6. These values in ANY system are extreme when used together, but there is no chance of Type 1 error. The danger is Type 2 error , where significant values are either omitted or considered non-significant. This could happen here from two fronts:</p>
<ul>
<li>If this results e.g. 100 observations, and DEGs with P-values &lt; 5e-4 but &gt; 1e-6 are omitted</li>
<li>DEGs &lt; (or &gt;) +/- 2.5 log<sub>2</sub> and above the lower threshold are omitted.</li>
</ul>
<p>However, it will depend on the question.</p>
<p><strong>Note</strong> Much of the theory in <em>point 2</em>, would also apply to the BH-method. However, I've worked this theory exclusively for Bonferoni's correction and I've not given much thought to BH. This is because it's not relevant to my purpose.</p>
<hr />
<p><em>See comments below your post about your TPM approach</em>. This post is cut-off criteria.</p>
<hr />
<p>Just to finish this post, here's a volcano with thanks to <a href="https://www.britannica.com/science/volcano" rel="nofollow noreferrer">Britannica</a>  for permitting the link-view. I think the idea of a volcano plot refers to the &quot;plume&quot; of ash, dust and occasionally lava (which some of us have real-world experience).</p>
<p><a href="https://i.sstatic.net/I74jU.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/I74jU.jpg" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2023-06-11T18:08:05.207Z
2023-06-12T16:06:01.410Z
2023-06-12T16:06:01.410Z
M__
<p>So for context, I have a set of TPM values for multiple genes for different samples, and I need to calculate the differential expression for RNA-seq values. I've been using this website as a guide: <a href="https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2019_March_UCSF_mRNAseq_Workshop/master/differential_expression/DE_Analysis.Rmd" rel="nofollow noreferrer">https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2019_March_UCSF_mRNAseq_Workshop/master/differential_expression/DE_Analysis.Rmd</a>. The overall aim of this calculation is so that I can use PANDA to map TM motifs/combine gene expression data to create networks representing interactions between transcription factors and genes, and GSEA to analyze genes ranked by fold change or differential expression p-value.</p>
<p>Note that I am using log2(TPM+1) values. How do I calculate what the filtering cutoff should be? Note that I already excluded all genes that had any values of 0 for any of the samples, and I will use the Benjamini-Hochberg method to adjust p-values for multiple testing afterwards.</p>

CC BY-SA 4.0
2023-06-08T18:36:53.370Z
The number of observations is little less than 19,000
Feynman
CC BY-SA 4.0
2023-06-12T14:50:05.717Z
TPM + Limma looks fine.  Cut-offs are cut-offs, end of story.  The DEG issue is the variance is much greater than the mean: DESeq2, edgeR and voom attempted to "hardwire this relationship" hence criticisms from Li et al (2022) and others. DESeq2 and edgeR need tests to assess the fit between data and model (residual). That's common in other stats.
M__
CC BY-SA 4.0
2023-06-12T15:22:06.623Z
TPM *per se* has been criticised in some published comparative studies  recently. **I wouldn't discard your approach** because: 1) *it's a recognised approach*; 2) transformations are complicated and the alternative (counts) are back to edgeR, DESeq2, NOIseq ... not fun; 3) the weaknesses are not formalised like 20+% FDR (theoretically always anticipated for DESeq2 etc.). You may be requested to do an alternative method post peer-review, or you could go into comparative methodology (which is good analytics). *Your approach is ok stick with it*, complimenting it is your call
M__
<differential-expression>
How do I calculate DEG cutoffs?
# How do I calculate DEG cutoffs?

(Asked by: username_0 on 2023-06-08)

So for context, I have a set of TPM values for multiple genes for different samples, and I need to calculate the differential expression for RNA-seq values. I've been using this website as a guide: <https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2019_March_UCSF_mRNAseq_Workshop/master/differential_expression/DE_Analysis.Rmd>. The overall aim of this calculation is so that I can use PANDA to map TM motifs/combine gene expression data to create networks representing interactions between transcription factors and genes, and GSEA to analyze genes ranked by fold change or differential expression p-value.

Note that I am using log2(TPM+1) values. How do I calculate what the filtering cutoff should be? Note that I already excluded all genes that had any values of 0 for any of the samples, and I will use the Benjamini-Hochberg method to adjust p-values for multiple testing afterwards.

**Tags:** differential-expression
> **Comments:**
> **username_0**: The number of observations is little less than 19,000  
> **username_1**: TPM + Limma looks fine.  Cut-offs are cut-offs, end of story.  The DEG issue is the variance is much greater than the mean: DESeq2, edgeR and voom attempted to "hardwire this relationship" hence criticisms from Li et al (2022) and others. DESeq2 and edgeR need tests to assess the fit between data and model (residual). That's common in other stats.  
> **username_1**: TPM *per se* has been criticised in some published comparative studies  recently. **I wouldn't discard your approach** because: 1) *it's a recognised approach*; 2) transformations are complicated and the alternative (counts) are back to edgeR, DESeq2, NOIseq ... not fun; 3) the weaknesses are not formalised like 20+% FDR (theoretically always anticipated for DESeq2 etc.). You may be requested to do an alternative method post peer-review, or you could go into comparative methodology (which is good analytics). *Your approach is ok stick with it*, complimenting it is your call  

---

## 2 answers

### Answer by username_2 on 2023-06-11. Score: 1

tl:dr use limma (not voom, that's for raw counts) and refer to previous questions:

* <https://support.bioconductor.org/p/98820/#98875>
* <https://support.bioconductor.org/p/125144/#125156>

The problem you have is not so much the threshold but the method to infer any metric to be used for this threshold. The arguably most common methods in R for DEG analysis are DESeq2, edgeR and limma, based on number of citations. The former two require raw counts which you apparently do not have, but if you do then use this with any of the mentioned tools.

limma is what is left and could be used here with TPM. Note that TPM has the problem that counts are divided by a length factor. A brief background: Any statistical test aims to decide whether the given variance between two groups (basically the logFC) is greater than expected by chance, which is the variance within group. In RNA-seq the variance is not equal across the range of counts (low, middle, high, for simplicity). On log2 count scale variance is greater for genes with low(er) counts than for higher counts. Usually, using raw counts the method (DESeq2, edgeR, limma, and others) learn from the data what is the expected variance at the given average expression level, and incorporate that into the model to make adequate inference even at small sample size. That is by the way where the common assumption "most genes do not change" comes from. Assuming that most genes are not DE, then the mean-variance trend represents the expected variance at the given average count range, and this can then be used to help decide whether a gene is DE or not. With TPM you basically throw this out of the window because the actual counts (as they come from the sequencer) are manipulated by dividing by length and the entire mean-variance model gets distorted. Another factor is that TPM, which corrects by total library size can fail to properly normalize counts, I wrote a summary here:

* <https://www.biostars.org/p/9465851/#9465854>

**Long story short: TPM is not good for DEG.**

If it's all you have then use what is called "limma-trend", see the limma manual at Bioconductor, which will directly model the log2-TPMs and give some logFC and FDR. You can then use standard cutoffs, such as FDR < 0.05 & logFC >/< 1.5. These cutoffs are not perfect, but common, and given that TPMs are bad anyway you can and should stick to standard cutoffs as much as possible, rather than infering your own ones. This at least gives some standardization to this suboptimal analysis.


---

### Answer by username_1 on 2023-06-11. Score: 0

Two solutions - the underlying issue is called Type 1 error (below).

1. If you are performing studies on humans then the Benjamini-Hochberg method is without question the method of choice, because this includes the *false detection* rate for EdgeR and DeSeq2. Thats been independently calculated.

Thus, the paper that caused a stir was [Li *et al* (2022)](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02648-4) who wrote,

> Particularly, the actual FDRs of DESeq2 and edgeR sometimes exceed 20% when the target FDR is 5%.

***Benjamini-Hochberg (BH) method solves that for humans (EdgeR and DeSeq2)***, because the problem with the BH method is you need a false discovery rate, Li *et al* 2022, calculated that its about 20% ... so for EdgeR and DeSeq2 the problem is solved. The BH tests will simply use the estimate for the FDR by Li *et al* (2022) instead of the FDR by EdgeR and DeSeq2. Cool? It is :-) it wasn't what Li *et al* (2022) intended but it is the solution to a difficult conundrum (concerning a statistical issue or parametrics vs. non-parametrics). **So yeah it's cool** (needs count data).

To perform BH you simply rank the observed probabilities (e.g. from Limma), and there are significant if each is below (or equal to) the BH rank. The BH rank is the total number of tests as denominator, the rank position as numerator and this is multiplied by the false detection rate as a frequency.

*HOWEVER*, you are using Limma and I don't how the model works. I know the underlying stats behind EdgeR, DeSeq2 and NOISeq.

---

Okay so the false discovery rate is not known ... in the way Li *et al* (2022) calculated it. What now? Limma is likely to calculate an FDR, you could use that via BH. My system is so far from humans there's no way Li *et al* (2022) FDR is relevant.

2. Bonferoni's approach - I don't know how mainstream this is, but its what I use (I've a background in classical stats), because I've no idea what the true FDR is. This approach can be exemplified via Volcano plots is to visualise the Type 1 error issues.

Here's a volcano plot which permits online viewing via links [obtained here](https://www.htgmolecular.com/blog/2022-08-25/understanding-volcano-plots). I can't change the size BTW (Bioconductor and loads more BTW don't permit this),

[![enter image description here](https://i.sstatic.net/mFeN2.jpg)](https://i.sstatic.net/mFeN2.jpg)

This is using a flat probability - which in the picture is very low - it's the horizontal threshold (*y*-axis). So the P value is one parameter approach to significance ....

The other parameter is the log2 threshold ...

Here's a two-parameter approach obtained from [GeneKitR's site](https://www.genekitr.fun/plot-others-1.html), who do permit images to be link viewed (thank you genekitr).

[![enter image description here](https://i.sstatic.net/gnq9l.jpg)](https://i.sstatic.net/gnq9l.jpg)

Here the authors have used log2 as the threshold (*x*-axis). The threshold is probably log2 > 1 and log2 < -1, which is reasonable. They have set the P-value threshold to be really low, looks to be about 5%. The plot reveals the effective P-value is 1E-2.5, which is low, but nothing really starts happening until 1E-3, which ... okay. Although this is a two-parameter approach the guys really just relied on log2.

**Type 1 error and a 1-parameter method threshold**

The massive problem is the shear number of tests performed and this will result in significant values PURELY BY CHANCE ...artefacts of critical probability. BH is one approach, Bonferoni's is an alternative way in.

Put simply for your data to calculate the Bonferoni threshold if a 5% critical probability is sought then its ...

> 0.05 / 19000 = 2.6e-6

Which is okay. So any single result <2.6E-6 is significant. This okay, but you will lose results so its vulnerable to Type 2 error. Why do you lose stuff ...

**Type 2 error and 1-parameter method threshold**

The biological question is e.g. is 0.1log2 biological sensible? What does it actually mean? Honestly I don't have a clue ... and the problem with a 1-parameter approach, it's very conservative but infers weird stuff. An expression difference of 10% isn't really serious, but thats what it's measuring. ... because its measuring this it means the probability via Bonferoni's is set high - *for meaningless stuff* and you can be missing important shifts in log2, e.g. >1 because they are falling under the high threshold.

**Type 1 error and a 2-parameter method thresholds**

This is where it all starts getting interesting. As soon as a +/- log2 threshold is set the number of total tests will drop A LOT. That means the P-value threshold raises. For example, if selecting > (or <) +/- 1log2 as used in the above example, supposing the number of observations above/below the +/- 1log2 was 5000. The P-value now raises via Bonferoni, e.g.

> 0.05/5000 = 1e-05

That looks better there'll be more results. Supposing we set a threshold of +/- 1.5log2 and obtained 500 results ... then the new Bonferoni derived threshold is,

> 0.05/500 = 1e-04 (0.0001)

That looks fabulous giving loads of significant results ... Wait isn't this all getting a bit crazy, because it depends on the biological assumption of what a reasonable log2 is. Yep, indeed, agreed.

**Type 1 error and a 2-parameter method thresholds: weaknesses**

The issue is that a reasonable log2 is sensible, because there's a point where log2 expression differences becoming utterly meaningless. However there's a point up to -/+ 2.5 log2 (which is very high) where you could choose anything - that shifts the critical probability and that shifts how many results are positive.

So how do you ...

**Set a log2 within a 2-parameter threshold?**

I've just solved that. I was really happy. The practical answer is it's complicated. My contractual basis is I can't reveal that because public disclosure isn't permitted without top-level approval.

The approach is to use an independent method to calculation DEG (if you check my posts you'll see I do AI). You can then perform parameter tuning of DEG against the independent method and statistically assess the association. In my case I was surprised by the answer - intuitively I wouldn't have gone there, but my system ain't humans.

**Conclusion**

In your case,

* simply check the literature to see what other log2 thresholds are in use, use that and adjust the P-value threshold accordingly as described above and cite that work.
* look at the biological reality of the results at different log2 thresholds. For example, an operon should behave like an operon. If the log2 threshold is smashing-up that operon - that doesn't make sense and adjust log2 threshold by lowering.
* alternatively use BH, this however assumes Limma's FDR calculation is correct and the problem is there no independent study examining that.

Either of the above 3 approaches will be fine your context. This is not cutting-edge statistical methodology it simply needs to be reasonable.

**Ultra-conservative parameters and Type 2 error**

To avoid all the complexities some investigators simply use > (or <) +/- 2.5log2 and P-value = 1e-6. These values in ANY system are extreme when used together, but there is no chance of Type 1 error. The danger is Type 2 error , where significant values are either omitted or considered non-significant. This could happen here from two fronts:

* If this results e.g. 100 observations, and DEGs with P-values < 5e-4 but > 1e-6 are omitted
* DEGs < (or >) +/- 2.5 log2 and above the lower threshold are omitted.

However, it will depend on the question.

**Note** Much of the theory in *point 2*, would also apply to the BH-method. However, I've worked this theory exclusively for Bonferoni's correction and I've not given much thought to BH. This is because it's not relevant to my purpose.

---

*See comments below your post about your TPM approach*. This post is cut-off criteria.

---

Just to finish this post, here's a volcano with thanks to [Britannica](https://www.britannica.com/science/volcano) for permitting the link-view. I think the idea of a volcano plot refers to the "plume" of ash, dust and occasionally lava (which some of us have real-world experience).

[![enter image description here](https://i.sstatic.net/I74jU.jpg)](https://i.sstatic.net/I74jU.jpg)


---

