Answers.Body: <blockquote>
    <p>I'd like to find genes that were not expressed in a group of samples and were expressed in another group.</p>
  </blockquote>
  
  <p>This is, fundamentally, a differential expression analysis, with a twist. To solve this, you’d first use a differential expression library of your choice (e.g. DESeq2) and perform a one-tailed test of differential expression.</p>
  
  <p>Briefly, you’d perform the normal setup and then use</p>
  
  <pre><code>results(dds, altHypothesis = 'greater')
  </code></pre>
  
  <p>To perform a one-tailed test. This will give you only those genes that are significantly upregulated in one group. Check chapter 3.9 of the vignette for details.</p>
  
  <p>Of course this won’t tell you that the genes are <em>unexpressed</em> in the other group. Unfortunately I don’t know of a good value to threshold the results; I would start by plotting a histogram of the (variance stabilised) expression values in your first group, and then visually choose an expression threshold that cleanly separates genes that are clearly expressed from zeros:</p>
  
  <pre><code>vst_counts = assay(vst(dds))
  dens = density(vst_counts[, replicate])
  plot(dens, log = 'y')
  </code></pre>
  
  <p>(This merges the replicates in the group, which should be fine.)</p>
  
  <p>Counts follow a multimodal distribution, with one mode for unexpressed and one or more for expressed genes. The expression threshold can be set somewhere between the clearly unexpressed and expressed peaks:</p>
  
  <p><a href="https://i.sstatic.net/MF0vy.png" rel="noreferrer"><img src="https://i.sstatic.net/MF0vy.png" alt="enter image description here"></a></p>
  
  <p>Here I used <code>identify(dens)</code> to identify the threshold interactively but you could also use an analytical method:</p>
  
  <pre><code>threshold = identify(dens)
  quantile = sum(dens$x &lt; dens$x[threshold]) / length(dens$x)
  
  # Using just one replicate here; more robust would be to use a mean value.
  nonzero_counts = counts(dds, normalized = TRUE)[, replicates[1]]
  nonzero_counts = nonzero_counts[nonzero_counts &gt; 0]
  
  (expression_threshold = quantile(nonzero_counts, probs = quantile))
  </code></pre>
  
  
  
  <pre><code>26.5625%
  4.112033
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T10:35:30.620Z
Answers.Text: I like this approach but I have a few concerns. In the DE step, how does the differential expression test perform with zero count values? If one group has zero values and the other group has near zero, or clearly expressed, but with high variability etc.
  In the cutoff step, (assuming that the blue line is the threshold) it seems you have a very low cutoff value (~1), so you will basically accept all DE genes (or accept all genes above count 0)?
Answers.UserDisplayName: Peter
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T10:37:51.827Z
Answers.Text: Also the graph is not bimodal, but has 4 local maxima. The second peak corresponds to roughly 1-2 counts, and when you have 20 million reads, they are probably not not expressed. Is it not better to have the cutoff after the 3rd peak, around rlog=3?
Answers.UserDisplayName: Peter
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T10:48:57.443Z
Answers.Text: @Peter Your observation about the different modes is a good one that I glossed over. In fact, the distinction should be between a mode for unexpressed genes, and “all the rest”. In order to determine which mode corresponds to truly unexpressed genes, it might be necessary to plot the density with different bandwidths and see how the modes change. My plot isn’t perfect in this regard since I used simulated data based on averaged counts. Real counts would behave better.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T12:24:31.943Z
Answers.Text: @Peter Your comments made many good points so I redid the example with actual read count data, and used VST instead of rlog, which appears to work somewhat better here.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T11:13:08.053Z
Answers.LastActivityDate: 2017-06-14T12:23:37.687Z
Answers.LastEditDate: 2017-06-14T12:23:37.687Z
Answers.UserDisplayName: Konrad Rudolph
Answers.Body: <p>A common method is to use <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3870982/" rel="noreferrer">zFPKMs</a>, which you can find implemented in R <a href="https://github.com/severinEvo/gene_expression/blob/master/zFPKM.R" rel="noreferrer">here</a>.</p>
  
  <p>Having said that, there's an inherent problem in declaring a difference between two things on either side of a given threshold. Given that, I would encourage you to use more than "expressed in one and not in another", likely adding at least a "with a minimal difference of X between them" metric. You may also find the tau metric useful, which is implemented in R <a href="https://github.com/severinEvo/gene_expression/blob/master/tau.R" rel="noreferrer">here</a>. This is meant to measure tissue-specificity, which is more akin to what you're probably interested in doing.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T11:09:16.817Z
Answers.LastActivityDate: 2017-06-13T11:09:16.817Z
Answers.UserDisplayName: Devon Ryan
Answers.Body: <p>Depending on how much effort you wish to put into this, here is one suggestion I have see used before (uses more than just R). Steps 1-4 come from <a href="http://genome.cshlp.org/content/24/12/1918.long#ref-55" rel="nofollow noreferrer">this paper</a> (see supplimentary material section "Calculation of per gene Local FDR").</p>
  
  <ol>
  <li><p><em>Take your gene models and collapse the introns down to maybe 100bp</em><br>
  Thus if we have a gene with two exons and a 1kb intron eg. exons (1000, 1100) and (2000,2100), we reduce the intron size so that the exons are (1000,1100) and (1200,1300). This is because we need to find sufficient gene free space to fit the null gene in.    </p></li>
  <li><p><em>Shift each gene into the nearest genomic space at least 5 kp from a genomic annotation and free of ESTs.</em></p></li>
  <li><p><em>Calculate the FPKM distribution of this shifted null set</em>  </p></li>
  <li><p><em>From these two distributions it should be possible to calculate a <a href="https://cran.r-project.org/web/packages/locfdr/vignettes/locfdr-example.pdf" rel="nofollow noreferrer">local FDR</a> for expression.</em><br>
  You will need to use a "mixing proportion" to do this. The above reference used <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2660870/" rel="nofollow noreferrer">Qvality</a> to do this.</p></li>
  <li><p><em>Carry out differential expression analysis as described in <a href="https://bioinformatics.stackexchange.com/a/690/235">Konrad's answer</a></em></p></li>
  <li><p><em>Subset the differentially expressed genes to only consider those that are not expressed in the control condition</em></p></li>
  <li><p><em>You might want to also recalculate the FDRs for differential expression as you are only considering a subset of the tests you might have.</em> </p></li>
  </ol>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T09:53:18.767Z
Answers.Text: This sounds like a promising approach for calculating a threshold. However,  I don’t entirely understand your step 1, could you expand on that?
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-15T09:54:14.353Z
Answers.Text: Edited to include more detail on step 1.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-15T10:58:23.797Z
Answers.Text: Yes the calculation of local FDR is something that has come up lately and been supported both by Gordon Smith and Mike Love. In any case if you do not want to do that one can always start from a list of genes that are above a certain threshold of expression levels to restrict the universe of expressed genes and then perform your DE with DESeq2.
Answers.UserDisplayName: ivivek_ngs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T18:46:27.060Z
Answers.LastActivityDate: 2017-06-15T09:53:50.680Z
Answers.LastEditDate: 2017-06-15T09:53:50.680Z
Answers.UserDisplayName: Ian Sudbery
Answers.Body: <p>Unfortunatly that is not doable - there is now way of distinguishing between features not expressed and features not expressed for technical reasons. Furthermore a lot of genes might seem lowly expressed to to technical noise.</p>
  
  <p>This also means that you will have to choose a cutoff and simply use that. For choosing a cutoff looking at the replicate distribution of log10(counts) and log10(FPKM) are probably the best options.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T13:14:32.560Z
Answers.Text: I think you’ve posted your answer on the wrong question: did you mean to post it here instead? https://bioinformatics.stackexchange.com/q/693/29
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-15T10:03:25.983Z
Answers.Text: @KonradRudolph Actually, the linked question did not exist when this answer was posted.
Answers.UserDisplayName: Peter
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T11:30:59.540Z
Answers.LastActivityDate: 2017-06-13T11:30:59.540Z
Answers.UserDisplayName: Kristoffer Vitting-Seerup
Body: <p>I have a gene expression count matrix produced from bulk RNA-seq data. I'd like to find genes that were <strong>not</strong> expressed in a group of samples and were expressed in another group.</p>
  
  <p>The problem of course is that not all <em>effectively</em> non-expressed genes will have 0 counts due to sequencing errors, or because they were expressed in a small subset of cells.</p>
  
  <p>I'm interested in solutions using <code>R</code>.</p>
  
Tags: <rna-seq><r>
Title: What methods are available to find a cutoff value for non-expressed genes in RNA-seq?
ThreadText: # What methods are available to find a cutoff value for non-expressed genes in RNA-seq?
  
  (Asked by: username_0 on 2017-06-13)
  
  I have a gene expression count matrix produced from bulk RNA-seq data. I'd like to find genes that were **not** expressed in a group of samples and were expressed in another group.
  
  The problem of course is that not all *effectively* non-expressed genes will have 0 counts due to sequencing errors, or because they were expressed in a small subset of cells.
  
  I'm interested in solutions using `R`.
  
  **Tags:** rna-seq, r
  
  
  ---
  
  ## 4 answers
  
  ### Answer by username_1 on 2017-06-13. Score: 13
  
  > I'd like to find genes that were not expressed in a group of samples and were expressed in another group.
  
  This is, fundamentally, a differential expression analysis, with a twist. To solve this, you’d first use a differential expression library of your choice (e.g. DESeq2) and perform a one-tailed test of differential expression.
  
  Briefly, you’d perform the normal setup and then use
  
  ```
  results(dds, altHypothesis = 'greater')
  
  ```
  
  To perform a one-tailed test. This will give you only those genes that are significantly upregulated in one group. Check chapter 3.9 of the vignette for details.
  
  Of course this won’t tell you that the genes are *unexpressed* in the other group. Unfortunately I don’t know of a good value to threshold the results; I would start by plotting a histogram of the (variance stabilised) expression values in your first group, and then visually choose an expression threshold that cleanly separates genes that are clearly expressed from zeros:
  
  ```
  vst_counts = assay(vst(dds))
  dens = density(vst_counts[, replicate])
  plot(dens, log = 'y')
  
  ```
  
  (This merges the replicates in the group, which should be fine.)
  
  Counts follow a multimodal distribution, with one mode for unexpressed and one or more for expressed genes. The expression threshold can be set somewhere between the clearly unexpressed and expressed peaks:
  
  [![enter image description here](https://i.sstatic.net/MF0vy.png)](https://i.sstatic.net/MF0vy.png)
  
  Here I used `identify(dens)` to identify the threshold interactively but you could also use an analytical method:
  
  ```
  threshold = identify(dens)
  quantile = sum(dens$x < dens$x[threshold]) / length(dens$x)
  
  # Using just one replicate here; more robust would be to use a mean value.
  nonzero_counts = counts(dds, normalized = TRUE)[, replicates[1]]
  nonzero_counts = nonzero_counts[nonzero_counts > 0]
  
  (expression_threshold = quantile(nonzero_counts, probs = quantile))
  
  ```
  
  ```
  26.5625%
  4.112033
  
  ```
  
  > **Comments:**
  > **username_0**: I like this approach but I have a few concerns. In the DE step, how does the differential expression test perform with zero count values? If one group has zero values and the other group has near zero, or clearly expressed, but with high variability etc. In the cutoff step, (assuming that the blue line is the threshold) it seems you have a very low cutoff value (~1), so you will basically accept all DE genes (or accept all genes above count 0)?  
  > **username_0**: Also the graph is not bimodal, but has 4 local maxima. The second peak corresponds to roughly 1-2 counts, and when you have 20 million reads, they are probably not not expressed. Is it not better to have the cutoff after the 3rd peak, around rlog=3?  
  > **username_1**: @username_0 Your observation about the different modes is a good one that I glossed over. In fact, the distinction should be between a mode for unexpressed genes, and “all the rest”. In order to determine which mode corresponds to truly unexpressed genes, it might be necessary to plot the density with different bandwidths and see how the modes change. My plot isn’t perfect in this regard since I used simulated data based on averaged counts. Real counts would behave better.  
  > **username_1**: @username_0 Your comments made many good points so I redid the example with actual read count data, and used VST instead of rlog, which appears to work somewhat better here.  
  
  ---
  
  ### Answer by username_2 on 2017-06-13. Score: 10
  
  A common method is to use [zFPKMs](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3870982/), which you can find implemented in R [here](https://github.com/severinEvo/gene_expression/blob/master/zFPKM.R).
  
  Having said that, there's an inherent problem in declaring a difference between two things on either side of a given threshold. Given that, I would encourage you to use more than "expressed in one and not in another", likely adding at least a "with a minimal difference of X between them" metric. You may also find the tau metric useful, which is implemented in R [here](https://github.com/severinEvo/gene_expression/blob/master/tau.R). This is meant to measure tissue-specificity, which is more akin to what you're probably interested in doing.
  
  
  ---
  
  ### Answer by username_3 on 2017-06-13. Score: 5
  
  Depending on how much effort you wish to put into this, here is one suggestion I have see used before (uses more than just R). Steps 1-4 come from [this paper](http://genome.cshlp.org/content/24/12/1918.long#ref-55) (see supplimentary material section "Calculation of per gene Local FDR").
  
  1. *Take your gene models and collapse the introns down to maybe 100bp*  
     Thus if we have a gene with two exons and a 1kb intron eg. exons (1000, 1100) and (2000,2100), we reduce the intron size so that the exons are (1000,1100) and (1200,1300). This is because we need to find sufficient gene free space to fit the null gene in.
  2. *Shift each gene into the nearest genomic space at least 5 kp from a genomic annotation and free of ESTs.*
  3. *Calculate the FPKM distribution of this shifted null set*
  4. *From these two distributions it should be possible to calculate a [local FDR](https://cran.r-project.org/web/packages/locfdr/vignettes/locfdr-example.pdf) for expression.*  
     You will need to use a "mixing proportion" to do this. The above reference used [Qvality](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2660870/) to do this.
  5. *Carry out differential expression analysis as described in [username_1's answer](https://bioinformatics.stackexchange.com/a/690/235)*
  6. *Subset the differentially expressed genes to only consider those that are not expressed in the control condition*
  7. *You might want to also recalculate the FDRs for differential expression as you are only considering a subset of the tests you might have.*
  
  > **Comments:**
  > **username_1**: This sounds like a promising approach for calculating a threshold. However,  I don’t entirely understand your step 1, could you expand on that?  
  > **username_3**: Edited to include more detail on step 1.  
  > **username_4**: Yes the calculation of local FDR is something that has come up lately and been supported both by Gordon Smith and Mike Love. In any case if you do not want to do that one can always start from a list of genes that are above a certain threshold of expression levels to restrict the universe of expressed genes and then perform your DE with DESeq2.  
  
  ---
  
  ### Answer by username_5 on 2017-06-13. Score: 0
  
  Unfortunatly that is not doable - there is now way of distinguishing between features not expressed and features not expressed for technical reasons. Furthermore a lot of genes might seem lowly expressed to to technical noise.
  
  This also means that you will have to choose a cutoff and simply use that. For choosing a cutoff looking at the replicate distribution of log10(counts) and log10(FPKM) are probably the best options.
  
  > **Comments:**
  > **username_1**: I think you’ve posted your answer on the wrong question: did you mean to post it here instead? https://bioinformatics.stackexchange.com/q/693/29  
  > **username_0**: @username_1 Actually, the linked question did not exist when this answer was posted.  
  
  ---
  
  