Answers.Body: <p>edgeR functions do not do any filtering, i.e., they do not remove any genes from the results or from the plots. The only filtering was that which you did explicitly yourself by:</p>
  <pre><code>keep &lt;- filterByExpr(R10LB_0vs0.5)
  R10LB_0vs0.5_filter &lt;- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]
  </code></pre>
  <p>For your data, <code>filterByExpr</code> will keep any gene that is expressed at a minimal level in at least 7 samples.
  You could relax that requirement, to allow genes to be expressed in only 5 or 6 samples say, by reducing the <code>large.n</code> argument to <code>filterByExpr</code>.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-12-24T03:52:00.557Z
Answers.LastActivityDate: 2021-12-24T04:04:15.170Z
Answers.LastEditDate: 2021-12-24T04:04:15.170Z
Answers.UserDisplayName: Gordon Smyth
Answers.Body: <p><em>[disclaimer: I understand DESeq2, not EdgeR]</em></p>
  <p>Most likely, low expression values (as in raw counts, prior to any normalisation such as TPM). Expression that has high dispersion across all samples is less likely to be biologically significant.</p>
  <p>It's easier to work out what's going on with an MA plot, with average gene expression on the X axis, and log2 fold change on the Y axis.</p>
  <p>I'm used to DESeq2 results, so I'd guess the construction would be something like this (assuming that the 'baseMean' field exists in EdgeR results):</p>
  <pre><code>R10LB_0vs0.5lrt %&gt;%
      ggplot() +
      aes(y = logFC, col = padj, x=log(baseMean)) +
      geom_point()
  </code></pre>
  <p>It's possible that EdgeR doesn't have anything like baseMean in its default results. As a proxy for this, you could use the count totals from the input matrix (e.g. <code>R10LB_0vs0.5 %&gt;% rowSums() %&gt;% log2()</code>).</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-03-12T07:22:35.613Z
Answers.Text: Do you mean that the TPM values are too low? Because I believe the ones that have been filtered out have a similar value to other "0" tpm genes that were significant. Also if I use the average of two conditions, will the genes I pull out be less significant because the MA plot will not account for variance between replicates? Thanks for the answer
Answers.UserDisplayName: ry0-H
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-03-12T03:40:57.220Z
Answers.LastActivityDate: 2021-03-13T09:52:29.413Z
Answers.LastEditDate: 2021-03-13T09:52:29.413Z
Answers.UserDisplayName: gringer
Body: <p>I am using RNA seq data to analyze genes via a volcano plot comparing differential gene expression of bacteria with and without antibiotic in R. After having created my plot, I am unsure why some of my values which read &quot;0&quot; TPM in one condition and a TPM value that's not &quot;0&quot; in the other condition were not determined to be differentially expressed. Some of the genes on my volcano plot have this difference in TPM and show up as significant on my plot while others with this &quot;0&quot; value difference are not considered to be significant according to my plot.</p>
  <p><strong>Sample of my data:</strong>
  (UO1_D712_##### represents the locus number, the top column represents the different replicates where &quot;1&quot; is untreated and &quot;2&quot; is treated with antibiotic, the letters next to the numbers represent a different technical replicate):</p>
  <pre><code>structure(list(`1` = c(&quot;U01_D712_00001&quot;, 
  
  &quot;U01_D712_00002&quot;, &quot;U01_D712_00003&quot;, 
  &quot;U01_D712_00004&quot;, &quot;U01_D712_00005&quot;, &quot;U01_D712_00006&quot;, &quot;U01_D712_00007&quot;, 
  &quot;U01_D712_00008&quot;, &quot;U01_D712_00009&quot;, &quot;U01_D712_00010&quot;, &quot;U01_D712_00011&quot;, 
  &quot;U01_D712_00012&quot;, &quot;U01_D712_00013&quot;, &quot;U01_D712_00014&quot;, &quot;U01_D712_00015&quot;, 
  &quot;U01_D712_00016&quot;, &quot;U01_D712_00017&quot;, &quot;U01_D712_00018&quot;, &quot;U01_D712_00019&quot;, 
  &quot;U01_D712_00020&quot;, &quot;U01_D712_00021&quot;, &quot;U01_D712_00022&quot;, &quot;U01_D712_00023&quot;, 
  &quot;U01_D712_00024&quot;, &quot;U01_D712_00025&quot;, &quot;U01_D712_00026&quot;, &quot;U01_D712_00027&quot;, 
  &quot;U01_D712_00028&quot;, &quot;U01_D712_00029&quot;, &quot;U01_D712_00030&quot;), `1a` = c(&quot;6.590456502&quot;, 
  &quot;6.741694688&quot;, &quot;7.110342585&quot;, &quot;6.299482433&quot;, &quot;2.77173982&quot;, &quot;2.470330508&quot;, 
  &quot;3.827125186&quot;, &quot;6.267229842&quot;, &quot;5.524708663&quot;, &quot;2.657913228&quot;, &quot;2.87209859&quot;, 
  &quot;7.479820548&quot;, &quot;4.980185572&quot;, &quot;4.210955199&quot;, &quot;0&quot;, &quot;4.492492822&quot;, 
  &quot;3.611091371&quot;, &quot;3.714270433&quot;, &quot;7.455036914&quot;, &quot;7.045203025&quot;, &quot;6.061860857&quot;, 
  &quot;2.925268313&quot;, &quot;6.544077039&quot;, &quot;5.747318013&quot;, &quot;9.97083444&quot;, &quot;9.22523089&quot;, 
  &quot;4.20205383&quot;, &quot;6.097040679&quot;, &quot;2.621192351&quot;, &quot;0&quot;), `1b` = c(&quot;6.544454427&quot;, 
  &quot;6.601489488&quot;, &quot;7.134224619&quot;, &quot;5.814043553&quot;, &quot;3.280958379&quot;, &quot;2.649180803&quot;, 
  &quot;3.860542083&quot;, &quot;6.256648363&quot;, &quot;5.380427766&quot;, &quot;2.581027705&quot;, &quot;3.016132165&quot;, 
  &quot;7.405329447&quot;, &quot;4.701503289&quot;, &quot;4.073814818&quot;, &quot;0&quot;, &quot;4.304196924&quot;, 
  &quot;3.515977329&quot;, &quot;3.843535649&quot;, &quot;7.342972625&quot;, &quot;6.966606769&quot;, &quot;6.122878624&quot;, 
  &quot;3.007522306&quot;, &quot;6.495797641&quot;, &quot;5.965431621&quot;, &quot;9.828050269&quot;, &quot;9.219563915&quot;, 
  &quot;4.065778989&quot;, &quot;6.105331066&quot;, &quot;3.061209408&quot;, &quot;0&quot;), `1c` = c(&quot;6.608196006&quot;, 
  &quot;6.743010138&quot;, &quot;7.102600793&quot;, &quot;6.146518601&quot;, &quot;3.555184202&quot;, &quot;2.364971542&quot;, 
  &quot;4.034053983&quot;, &quot;6.158523627&quot;, &quot;5.656051812&quot;, &quot;2.658660735&quot;, &quot;3.054717455&quot;, 
  &quot;7.392164473&quot;, &quot;4.950953264&quot;, &quot;4.277770477&quot;, &quot;0&quot;, &quot;4.507936666&quot;, 
  &quot;3.794842979&quot;, &quot;3.610794578&quot;, &quot;7.471646548&quot;, &quot;7.104792624&quot;, &quot;6.484767016&quot;, 
  &quot;3.071205184&quot;, &quot;6.584425715&quot;, &quot;6.20466015&quot;, &quot;9.986342122&quot;, &quot;9.282943758&quot;, 
  &quot;4.179958213&quot;, &quot;6.219551653&quot;, &quot;2.984738345&quot;, &quot;0&quot;), `1d` = c(&quot;6.547155382&quot;, 
  &quot;6.558328892&quot;, &quot;6.992501615&quot;, &quot;6.449558793&quot;, &quot;3.059801464&quot;, &quot;2.418800257&quot;, 
  &quot;3.96498952&quot;, &quot;6.013208538&quot;, &quot;5.279919645&quot;, &quot;2.893295085&quot;, &quot;1.750510471&quot;, 
  &quot;7.408735671&quot;, &quot;4.9425624&quot;, &quot;3.804549986&quot;, &quot;0&quot;, &quot;4.174835979&quot;, 
  &quot;3.806006888&quot;, &quot;3.570390524&quot;, &quot;7.641137006&quot;, &quot;6.976672494&quot;, &quot;6.363030106&quot;, 
  &quot;3.083061726&quot;, &quot;6.300910093&quot;, &quot;6.007490342&quot;, &quot;9.926316442&quot;, &quot;9.09671588&quot;, 
  &quot;4.320556917&quot;, &quot;6.153860107&quot;, &quot;2.877230446&quot;, &quot;0&quot;), `1e` = c(&quot;6.626724417&quot;, 
  &quot;6.577345176&quot;, &quot;7.156821278&quot;, &quot;6.296873411&quot;, &quot;3.618089702&quot;, &quot;2.394444986&quot;, 
  &quot;4.129376392&quot;, &quot;6.011715246&quot;, &quot;5.31197869&quot;, &quot;2.00754706&quot;, &quot;2.695493528&quot;, 
  &quot;7.538910448&quot;, &quot;4.606060035&quot;, &quot;3.909472643&quot;, &quot;0&quot;, &quot;4.346616047&quot;, 
  &quot;3.468681284&quot;, &quot;3.338231445&quot;, &quot;7.559599613&quot;, &quot;7.1527452&quot;, &quot;6.232923513&quot;, 
  &quot;3.108624209&quot;, &quot;6.535435309&quot;, &quot;6.12922864&quot;, &quot;10.13108497&quot;, &quot;9.310331313&quot;, 
  &quot;3.959568571&quot;, &quot;6.182335537&quot;, &quot;2.902736258&quot;, &quot;0&quot;), `1f` = c(&quot;6.419219179&quot;, 
  &quot;6.650459302&quot;, &quot;7.319125725&quot;, &quot;5.570869357&quot;, &quot;2.962845933&quot;, &quot;2.55903176&quot;, 
  &quot;4.087597573&quot;, &quot;5.995610538&quot;, &quot;5.386268651&quot;, &quot;2.750800859&quot;, &quot;2.416572678&quot;, 
  &quot;7.579148955&quot;, &quot;3.952633067&quot;, &quot;3.615227674&quot;, &quot;0&quot;, &quot;4.562838935&quot;, 
  &quot;3.76942104&quot;, &quot;3.920096905&quot;, &quot;7.935320749&quot;, &quot;7.501470652&quot;, &quot;6.13700099&quot;, 
  &quot;3.123910608&quot;, &quot;6.035971952&quot;, &quot;5.706235015&quot;, &quot;9.254751395&quot;, &quot;8.379630979&quot;, 
  &quot;4.51391973&quot;, &quot;5.6890651&quot;, &quot;2.43285316&quot;, &quot;0&quot;), `1g` = c(&quot;6.553221221&quot;, 
  &quot;6.633949852&quot;, &quot;7.182305386&quot;, &quot;5.769365973&quot;, &quot;2.721354972&quot;, &quot;1.668390466&quot;, 
  &quot;4.148367057&quot;, &quot;6.240883382&quot;, &quot;5.458877133&quot;, &quot;1.733842637&quot;, &quot;2.723803355&quot;, 
  &quot;7.522249899&quot;, &quot;4.149567197&quot;, &quot;3.780763096&quot;, &quot;0&quot;, &quot;4.496306813&quot;, 
  &quot;3.645643535&quot;, &quot;3.851001768&quot;, &quot;7.678552875&quot;, &quot;7.283411279&quot;, &quot;6.591585956&quot;, 
  &quot;2.879345378&quot;, &quot;6.389427003&quot;, &quot;5.911222165&quot;, &quot;9.851084493&quot;, &quot;9.084575304&quot;, 
  &quot;4.272587776&quot;, &quot;5.974762147&quot;, &quot;2.98852705&quot;, &quot;0&quot;), `2a` = c(&quot;9.769887737&quot;, 
  &quot;4.550226652&quot;, &quot;1.869021464&quot;, &quot;4.944848987&quot;, &quot;7.9678549&quot;, &quot;2.682865013&quot;, 
  &quot;8.495575559&quot;, &quot;2.234521659&quot;, &quot;3.667196316&quot;, &quot;9.180445037&quot;, &quot;3.210107621&quot;, 
  &quot;7.21523691&quot;, &quot;5.714579923&quot;, &quot;5.423986751&quot;, &quot;9.118981459&quot;, &quot;8.635701597&quot;, 
  &quot;2.742889473&quot;, &quot;3.712618983&quot;, &quot;8.006057144&quot;, &quot;4.999541279&quot;, &quot;10.54351774&quot;, 
  &quot;5.880978085&quot;, &quot;7.145433526&quot;, &quot;6.982416661&quot;, &quot;9.339651188&quot;, &quot;5.360835327&quot;, 
  &quot;4.699680905&quot;, &quot;3.423826225&quot;, &quot;6.408271885&quot;, &quot;5.038170992&quot;), 
      `2b` = c(&quot;10.26397519&quot;, &quot;1.945664005&quot;, &quot;2.086763158&quot;, &quot;2.800904763&quot;, 
      &quot;8.583418657&quot;, &quot;1.536094563&quot;, &quot;9.32057547&quot;, &quot;3.10685839&quot;, 
      &quot;4.224502319&quot;, &quot;10.1252842&quot;, &quot;1.811175407&quot;, &quot;5.714439316&quot;, 
      &quot;6.039142559&quot;, &quot;3.833361174&quot;, &quot;9.757360286&quot;, &quot;9.565906731&quot;, 
      &quot;1.523640473&quot;, &quot;2.315033488&quot;, &quot;6.312524363&quot;, &quot;4.889986456&quot;, 
      &quot;10.23020108&quot;, &quot;4.848727685&quot;, &quot;7.533256999&quot;, &quot;7.138160378&quot;, 
      &quot;10.30380331&quot;, &quot;3.955469283&quot;, &quot;3.167940742&quot;, &quot;3.599655687&quot;, 
      &quot;4.828945262&quot;, &quot;3.701043054&quot;), `2c` = c(&quot;9.478481216&quot;, &quot;2.789289131&quot;, 
      &quot;3.393949527&quot;, &quot;3.754810933&quot;, &quot;8.710777154&quot;, &quot;1.806170784&quot;, 
      &quot;9.150005253&quot;, &quot;2.612275457&quot;, &quot;4.961073313&quot;, &quot;9.802701699&quot;, 
      &quot;2.933183115&quot;, &quot;6.532384958&quot;, &quot;6.919449225&quot;, &quot;4.432699799&quot;, 
      &quot;9.715063475&quot;, &quot;9.265691356&quot;, &quot;1.412064593&quot;, &quot;3.330131873&quot;, 
      &quot;6.665979896&quot;, &quot;5.158526421&quot;, &quot;9.417365584&quot;, &quot;4.899531204&quot;, 
      &quot;8.173459354&quot;, &quot;7.271400938&quot;, &quot;9.813068613&quot;, &quot;4.384622077&quot;, 
      &quot;3.700645365&quot;, &quot;4.457874829&quot;, &quot;5.649440022&quot;, &quot;3.531010379&quot;
      ), `2d` = c(&quot;8.199795497&quot;, &quot;2.565711524&quot;, &quot;2.202287889&quot;, 
      &quot;3.856354444&quot;, &quot;7.380224849&quot;, &quot;2.192476466&quot;, &quot;8.14446837&quot;, 
      &quot;1.144258862&quot;, &quot;3.31447122&quot;, &quot;8.713146629&quot;, &quot;2.697890381&quot;, 
      &quot;6.304428859&quot;, &quot;5.745291803&quot;, &quot;4.898396114&quot;, &quot;9.173362747&quot;, 
      &quot;8.339933849&quot;, &quot;3.159678152&quot;, &quot;4.094587234&quot;, &quot;7.608649692&quot;, 
      &quot;5.280206424&quot;, &quot;10.34630403&quot;, &quot;5.098585806&quot;, &quot;7.262400625&quot;, 
      &quot;6.150190905&quot;, &quot;9.316698845&quot;, &quot;5.073027993&quot;, &quot;4.695003229&quot;, 
      &quot;2.485847024&quot;, &quot;5.545300465&quot;, &quot;4.350571411&quot;), `2e` = c(&quot;8.74935033&quot;, 
      &quot;3.739489484&quot;, &quot;0.217205045&quot;, &quot;4.413657999&quot;, &quot;8.745525588&quot;, 
      &quot;3.657060086&quot;, &quot;9.279834921&quot;, &quot;2.898951179&quot;, &quot;4.282874018&quot;, 
      &quot;9.610485827&quot;, &quot;3.561102455&quot;, &quot;7.228334332&quot;, &quot;6.388491443&quot;, 
      &quot;4.389908652&quot;, &quot;8.781086564&quot;, &quot;9.178866581&quot;, &quot;2.374603596&quot;, 
      &quot;3.961037408&quot;, &quot;7.864369809&quot;, &quot;3.654728044&quot;, &quot;10.15284858&quot;, 
      &quot;5.894439123&quot;, &quot;7.68020282&quot;, &quot;7.12243523&quot;, &quot;9.998637438&quot;, 
      &quot;5.092174395&quot;, &quot;4.111530392&quot;, &quot;3.776835632&quot;, &quot;5.624523213&quot;, 
      &quot;4.095011377&quot;), `2f` = c(&quot;7.648310926&quot;, &quot;4.345557215&quot;, &quot;1.986576876&quot;, 
      &quot;4.99426288&quot;, &quot;7.087937177&quot;, &quot;2.810917253&quot;, &quot;7.77706637&quot;, 
      &quot;2.62822773&quot;, &quot;3.581811188&quot;, &quot;8.470225989&quot;, &quot;3.335757437&quot;, 
      &quot;7.416094847&quot;, &quot;5.208841128&quot;, &quot;5.536128034&quot;, &quot;8.255571138&quot;, 
      &quot;7.993319997&quot;, &quot;1.9209089&quot;, &quot;3.573861828&quot;, &quot;7.318814519&quot;, 
      &quot;5.233804806&quot;, &quot;11.05855833&quot;, &quot;6.247720809&quot;, &quot;6.673407583&quot;, 
      &quot;6.029960625&quot;, &quot;8.806867591&quot;, &quot;5.459208493&quot;, &quot;4.001428729&quot;, 
      &quot;3.609936979&quot;, &quot;5.876522973&quot;, &quot;4.652839671&quot;), `2g` = c(&quot;7.555235468&quot;, 
      &quot;3.892899549&quot;, &quot;1.726443458&quot;, &quot;5.304546796&quot;, &quot;7.039588042&quot;, 
      &quot;3.027235295&quot;, &quot;7.703852207&quot;, &quot;1.753183519&quot;, &quot;2.909288568&quot;, 
      &quot;8.385169315&quot;, &quot;3.902707541&quot;, &quot;7.523315081&quot;, &quot;4.978364017&quot;, 
      &quot;5.49103181&quot;, &quot;8.096218606&quot;, &quot;7.944822989&quot;, &quot;2.352609608&quot;, 
      &quot;4.155433517&quot;, &quot;7.227355741&quot;, &quot;5.532668321&quot;, &quot;11.24946953&quot;, 
      &quot;6.159185473&quot;, &quot;6.443203375&quot;, &quot;5.931761874&quot;, &quot;8.7421732&quot;, 
      &quot;5.502000205&quot;, &quot;4.652883503&quot;, &quot;3.458323017&quot;, &quot;6.566487449&quot;, 
      &quot;4.89790353&quot;)), row.names = c(NA, 30L), class = &quot;data.frame&quot;)
  </code></pre>
  <p><strong>Design matrix:</strong></p>
  <pre><code>structure(list(SampleID = c(&quot;1a &quot;, &quot;1b&quot;, &quot;1c&quot;, 
  
  &quot;1d&quot;, &quot;1e&quot;, &quot;1f&quot;, 
  &quot;1g&quot;, &quot;2a&quot;, &quot;2b&quot;, &quot;2c&quot;, &quot;2d&quot;, &quot;2e&quot;, &quot;2f&quot;, &quot;2g&quot;), Group = c(1L, 
  1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Replicate = c(&quot;a&quot;, 
  &quot;b&quot;, &quot;c&quot;, &quot;d&quot;, &quot;e&quot;, &quot;f&quot;, &quot;g&quot;, &quot;a&quot;, &quot;b&quot;, &quot;c&quot;, &quot;d&quot;, &quot;e&quot;, &quot;f&quot;, &quot;g&quot;
  )), row.names = c(NA, 14L), class = &quot;data.frame&quot;)
  </code></pre>
  <p>Here is how I wrote my code to sample the RNA seq data (note ~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv refers to the the sample of data above and ~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv refers to the above design matrix):</p>
  <pre><code>R10LB_0vs0.5 &lt;- read.csv2(&quot;~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv&quot;, sep=&quot;,&quot;, check.names = F)
  R10LB_0vs0.5 &lt;- janitor::remove_empty(R10LB_0vs0.5, which = &quot;cols&quot;)
  R10LB_0vs0.5_design &lt;- read.csv2(&quot;~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv&quot;, sep=&quot;,&quot;)
  
  rownames(R10LB_0vs0.5) &lt;- R10LB_0vs0.5$`1`
  R10LB_0vs0.5 &lt;- R10LB_0vs0.5[,-1]
  
  #colnames(R10LB_0ugvs0.5ug) &lt;- R10LB_0ugvs0.5ug[1,]
  #R10LB_0ugvs0.5ug &lt;- R10LB_0ugvs0.5ug[-1,]
  
  R10LB_0vs0.5 &lt;- data.matrix(R10LB_0vs0.5)
  
  R10LB_0vs0.5 &lt;- DGEList(counts = R10LB_0vs0.5, group = R10LB_0vs0.5_design$Group)
  str(R10LB_0vs0.5)
  
  Group &lt;- as.vector(as.character(R10LB_0vs0.5_design<span class="math-container">$Group))
  Replicate &lt;- as.vector((R10LB_0vs0.5_design$</span>Replicate))
  R10LB_0vs0.5_designmatrix &lt;- model.matrix(~0+Group+Replicate)
  
  keep &lt;- filterByExpr(R10LB_0vs0.5)
  R10LB_0vs0.5_filter &lt;- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]
  
  R10LB_0vs0.5_Disp &lt;- estimateDisp(R10LB_0vs0.5_filter, R10LB_0vs0.5_designmatrix)
  
  CONTRASTS &lt;- makeContrasts( Group1vs2 = Group1-Group2,
                              levels = R10LB_0vs0.5_designmatrix)
  
  R10LB_0vs0.5fit &lt;- glmFit(R10LB_0vs0.5_Disp, contrast = CONTRASTS[,1])
  R10LB_0vs0.5lrt &lt;- glmLRT(R10LB_0vs0.5fit)
  R10LB_0vs0.5TT &lt;- topTags(R10LB_0vs0.5lrt, n=nrow(R10LB_0vs0.5_Disp))
  write.csv(R10LB_0vs0.5TT, file = &quot;R10LB_0vs0.5_comparison.csv&quot;)
  saveRDS(R10LB_0vs0.5TT, file = &quot;R10LB_0vs0.5_comparison.RDS&quot;)
  </code></pre>
  <p><strong>Volcano plot:</strong></p>
  <pre><code>ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = PValue)) + geom_point()
  
  p4 &lt;- ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = -log10(PValue))) + geom_point() + theme_minimal()
  
  #Add vertical lines for logFC thresholds, and one horizontal line for the p-value threshold
  p5 &lt;- p4 + geom_vline(xintercept = c(-0.6, 0.6), col = &quot;red&quot;) + geom_hline(yintercept = -log10(0.05), col = &quot;red&quot;)
  
  
  # add a column of NAs
  R10LB_0vs0.5_comparison<span class="math-container">$diffexpressed &lt;- "NO"
  # if logFC &gt; 0.6 and PValue &lt; 0.05, set as "UP" 
  R10LB_0vs0.5_comparison$</span>diffexpressed[R10LB_0vs0.5_comparison<span class="math-container">$logFC &gt; 0.6 &amp; R10LB_0vs0.5_comparison$</span>PValue &lt; 0.05] &lt;- &quot;UP&quot;
  # if logFC &lt; -0.6 and PValue &lt; 0.05, set as &quot;DOWN&quot;
  R10LB_0vs0.5_comparison<span class="math-container">$diffexpressed[R10LB_0vs0.5_comparison$</span>logFC &lt; -0.6 &amp; R10LB_0vs0.5_comparison$PValue &lt; 0.05] &lt;- &quot;DOWN&quot;
  
  p4 &lt;- ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed)) + geom_point() + theme_minimal()
  
  p5 &lt;- p4 + geom_vline(xintercept=c(-0.6, 0.6), col=&quot;red&quot;) + geom_hline(yintercept=-log10(0.05), col=&quot;red&quot;)
  
  p6 &lt;- p5 + scale_color_manual(values=c(&quot;blue&quot;, &quot;black&quot;, &quot;red&quot;))
  p6 &lt;- p5 + scale_colour_manual(values = mycolors)
  
  R10LB_0vs0.5_comparison<span class="math-container">$delabel &lt;- NA
  R10LB_0vs0.5_comparison$</span>delabel[R10LB_0vs0.5_comparison<span class="math-container">$diffexpressed != "NO"] &lt;- R10LB_0vs0.5_comparison$</span>gene_symbol[R10LB_0vs0.5_comparison$diffexpressed != &quot;NO&quot;]
  
  ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) + geom_point() + theme_minimal() +geom_text()
  
  library(ggrepel)
  ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) +geom_point() + theme_minimal() +geom_text_repel() + scale_color_manual(values=c(&quot;blue&quot;, &quot;black&quot;, &quot;red&quot;)) + geom_vline(xintercept=c(-0.6, 0.6), col=&quot;red&quot;) +geom_hline(yintercept=-log10(0.05), col=&quot;red&quot;)
  </code></pre>
  <p>This is the volcano plot I currently have with all of my data:
  <img src="https://i.sstatic.net/tidIM.png" alt="volcano plot" /></p>
  <p>Can you explain how some of these functions filtered my data (i.e. filterByExpr, estimateDisp, makeContrasts, glmFit, glmLRT)? Why are some of my data points that change from 0 TPM in one condition to some value in another condition showing up on my plot while others are not? Or should I accept that my data analysis is valid or &quot;correct&quot;?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-03-12T20:57:41.677Z
Comments.Text: I heard that low expression genes are often filtered out in a first past due to a high chance of false positives, so they could be dropped from filterByExpr. High variance within a condition might undermine the statistical significance of the mean difference across conditions. Is there a way to look at the variance of my  data, and is there a threshold for variance between my replicates that causes them to be filtered out by filterByExpr or other functions?
Comments.UserDisplayName: ry0-H
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-01-27T22:23:53.053Z
Comments.Text: It is a long time that I don't use R. However, I had a sinilar problem once. Using not cartesian coordinates, ggplot cut several points out of the plot. If you plan to continue using R, search about this problem. Otherwise, if you can accede to it go with MATLAB that has the best tool for volcano plot, by my opinion. It never fails.
Comments.UserDisplayName: Spartan 117
Tags: <r><rna-seq><bioconductor><edger><limma>
Title: Determining what RNAseq data is filtered on volcano plot
ThreadText: # Determining what RNAseq data is filtered on volcano plot
  
  (Asked by: username_0 on 2021-03-12)
  
  I am using RNA seq data to analyze genes via a volcano plot comparing differential gene expression of bacteria with and without antibiotic in R. After having created my plot, I am unsure why some of my values which read "0" TPM in one condition and a TPM value that's not "0" in the other condition were not determined to be differentially expressed. Some of the genes on my volcano plot have this difference in TPM and show up as significant on my plot while others with this "0" value difference are not considered to be significant according to my plot.
  
  **Sample of my data:**
  (UO1\_D712\_##### represents the locus number, the top column represents the different replicates where "1" is untreated and "2" is treated with antibiotic, the letters next to the numbers represent a different technical replicate):
  
  ```
  structure(list(`1` = c("U01_D712_00001", 
  
  "U01_D712_00002", "U01_D712_00003", 
  "U01_D712_00004", "U01_D712_00005", "U01_D712_00006", "U01_D712_00007", 
  "U01_D712_00008", "U01_D712_00009", "U01_D712_00010", "U01_D712_00011", 
  "U01_D712_00012", "U01_D712_00013", "U01_D712_00014", "U01_D712_00015", 
  "U01_D712_00016", "U01_D712_00017", "U01_D712_00018", "U01_D712_00019", 
  "U01_D712_00020", "U01_D712_00021", "U01_D712_00022", "U01_D712_00023", 
  "U01_D712_00024", "U01_D712_00025", "U01_D712_00026", "U01_D712_00027", 
  "U01_D712_00028", "U01_D712_00029", "U01_D712_00030"), `1a` = c("6.590456502", 
  "6.741694688", "7.110342585", "6.299482433", "2.77173982", "2.470330508", 
  "3.827125186", "6.267229842", "5.524708663", "2.657913228", "2.87209859", 
  "7.479820548", "4.980185572", "4.210955199", "0", "4.492492822", 
  "3.611091371", "3.714270433", "7.455036914", "7.045203025", "6.061860857", 
  "2.925268313", "6.544077039", "5.747318013", "9.97083444", "9.22523089", 
  "4.20205383", "6.097040679", "2.621192351", "0"), `1b` = c("6.544454427", 
  "6.601489488", "7.134224619", "5.814043553", "3.280958379", "2.649180803", 
  "3.860542083", "6.256648363", "5.380427766", "2.581027705", "3.016132165", 
  "7.405329447", "4.701503289", "4.073814818", "0", "4.304196924", 
  "3.515977329", "3.843535649", "7.342972625", "6.966606769", "6.122878624", 
  "3.007522306", "6.495797641", "5.965431621", "9.828050269", "9.219563915", 
  "4.065778989", "6.105331066", "3.061209408", "0"), `1c` = c("6.608196006", 
  "6.743010138", "7.102600793", "6.146518601", "3.555184202", "2.364971542", 
  "4.034053983", "6.158523627", "5.656051812", "2.658660735", "3.054717455", 
  "7.392164473", "4.950953264", "4.277770477", "0", "4.507936666", 
  "3.794842979", "3.610794578", "7.471646548", "7.104792624", "6.484767016", 
  "3.071205184", "6.584425715", "6.20466015", "9.986342122", "9.282943758", 
  "4.179958213", "6.219551653", "2.984738345", "0"), `1d` = c("6.547155382", 
  "6.558328892", "6.992501615", "6.449558793", "3.059801464", "2.418800257", 
  "3.96498952", "6.013208538", "5.279919645", "2.893295085", "1.750510471", 
  "7.408735671", "4.9425624", "3.804549986", "0", "4.174835979", 
  "3.806006888", "3.570390524", "7.641137006", "6.976672494", "6.363030106", 
  "3.083061726", "6.300910093", "6.007490342", "9.926316442", "9.09671588", 
  "4.320556917", "6.153860107", "2.877230446", "0"), `1e` = c("6.626724417", 
  "6.577345176", "7.156821278", "6.296873411", "3.618089702", "2.394444986", 
  "4.129376392", "6.011715246", "5.31197869", "2.00754706", "2.695493528", 
  "7.538910448", "4.606060035", "3.909472643", "0", "4.346616047", 
  "3.468681284", "3.338231445", "7.559599613", "7.1527452", "6.232923513", 
  "3.108624209", "6.535435309", "6.12922864", "10.13108497", "9.310331313", 
  "3.959568571", "6.182335537", "2.902736258", "0"), `1f` = c("6.419219179", 
  "6.650459302", "7.319125725", "5.570869357", "2.962845933", "2.55903176", 
  "4.087597573", "5.995610538", "5.386268651", "2.750800859", "2.416572678", 
  "7.579148955", "3.952633067", "3.615227674", "0", "4.562838935", 
  "3.76942104", "3.920096905", "7.935320749", "7.501470652", "6.13700099", 
  "3.123910608", "6.035971952", "5.706235015", "9.254751395", "8.379630979", 
  "4.51391973", "5.6890651", "2.43285316", "0"), `1g` = c("6.553221221", 
  "6.633949852", "7.182305386", "5.769365973", "2.721354972", "1.668390466", 
  "4.148367057", "6.240883382", "5.458877133", "1.733842637", "2.723803355", 
  "7.522249899", "4.149567197", "3.780763096", "0", "4.496306813", 
  "3.645643535", "3.851001768", "7.678552875", "7.283411279", "6.591585956", 
  "2.879345378", "6.389427003", "5.911222165", "9.851084493", "9.084575304", 
  "4.272587776", "5.974762147", "2.98852705", "0"), `2a` = c("9.769887737", 
  "4.550226652", "1.869021464", "4.944848987", "7.9678549", "2.682865013", 
  "8.495575559", "2.234521659", "3.667196316", "9.180445037", "3.210107621", 
  "7.21523691", "5.714579923", "5.423986751", "9.118981459", "8.635701597", 
  "2.742889473", "3.712618983", "8.006057144", "4.999541279", "10.54351774", 
  "5.880978085", "7.145433526", "6.982416661", "9.339651188", "5.360835327", 
  "4.699680905", "3.423826225", "6.408271885", "5.038170992"), 
      `2b` = c("10.26397519", "1.945664005", "2.086763158", "2.800904763", 
      "8.583418657", "1.536094563", "9.32057547", "3.10685839", 
      "4.224502319", "10.1252842", "1.811175407", "5.714439316", 
      "6.039142559", "3.833361174", "9.757360286", "9.565906731", 
      "1.523640473", "2.315033488", "6.312524363", "4.889986456", 
      "10.23020108", "4.848727685", "7.533256999", "7.138160378", 
      "10.30380331", "3.955469283", "3.167940742", "3.599655687", 
      "4.828945262", "3.701043054"), `2c` = c("9.478481216", "2.789289131", 
      "3.393949527", "3.754810933", "8.710777154", "1.806170784", 
      "9.150005253", "2.612275457", "4.961073313", "9.802701699", 
      "2.933183115", "6.532384958", "6.919449225", "4.432699799", 
      "9.715063475", "9.265691356", "1.412064593", "3.330131873", 
      "6.665979896", "5.158526421", "9.417365584", "4.899531204", 
      "8.173459354", "7.271400938", "9.813068613", "4.384622077", 
      "3.700645365", "4.457874829", "5.649440022", "3.531010379"
      ), `2d` = c("8.199795497", "2.565711524", "2.202287889", 
      "3.856354444", "7.380224849", "2.192476466", "8.14446837", 
      "1.144258862", "3.31447122", "8.713146629", "2.697890381", 
      "6.304428859", "5.745291803", "4.898396114", "9.173362747", 
      "8.339933849", "3.159678152", "4.094587234", "7.608649692", 
      "5.280206424", "10.34630403", "5.098585806", "7.262400625", 
      "6.150190905", "9.316698845", "5.073027993", "4.695003229", 
      "2.485847024", "5.545300465", "4.350571411"), `2e` = c("8.74935033", 
      "3.739489484", "0.217205045", "4.413657999", "8.745525588", 
      "3.657060086", "9.279834921", "2.898951179", "4.282874018", 
      "9.610485827", "3.561102455", "7.228334332", "6.388491443", 
      "4.389908652", "8.781086564", "9.178866581", "2.374603596", 
      "3.961037408", "7.864369809", "3.654728044", "10.15284858", 
      "5.894439123", "7.68020282", "7.12243523", "9.998637438", 
      "5.092174395", "4.111530392", "3.776835632", "5.624523213", 
      "4.095011377"), `2f` = c("7.648310926", "4.345557215", "1.986576876", 
      "4.99426288", "7.087937177", "2.810917253", "7.77706637", 
      "2.62822773", "3.581811188", "8.470225989", "3.335757437", 
      "7.416094847", "5.208841128", "5.536128034", "8.255571138", 
      "7.993319997", "1.9209089", "3.573861828", "7.318814519", 
      "5.233804806", "11.05855833", "6.247720809", "6.673407583", 
      "6.029960625", "8.806867591", "5.459208493", "4.001428729", 
      "3.609936979", "5.876522973", "4.652839671"), `2g` = c("7.555235468", 
      "3.892899549", "1.726443458", "5.304546796", "7.039588042", 
      "3.027235295", "7.703852207", "1.753183519", "2.909288568", 
      "8.385169315", "3.902707541", "7.523315081", "4.978364017", 
      "5.49103181", "8.096218606", "7.944822989", "2.352609608", 
      "4.155433517", "7.227355741", "5.532668321", "11.24946953", 
      "6.159185473", "6.443203375", "5.931761874", "8.7421732", 
      "5.502000205", "4.652883503", "3.458323017", "6.566487449", 
      "4.89790353")), row.names = c(NA, 30L), class = "data.frame")
  
  ```
  
  **Design matrix:**
  
  ```
  structure(list(SampleID = c("1a ", "1b", "1c", 
  
  "1d", "1e", "1f", 
  "1g", "2a", "2b", "2c", "2d", "2e", "2f", "2g"), Group = c(1L, 
  1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Replicate = c("a", 
  "b", "c", "d", "e", "f", "g", "a", "b", "c", "d", "e", "f", "g"
  )), row.names = c(NA, 14L), class = "data.frame")
  
  ```
  
  Here is how I wrote my code to sample the RNA seq data (note ~/Documents/VOLCANO/R10LB\_0ugvs0.5ug.csv refers to the the sample of data above and ~/Documents/VOLCANO/R10LB\_0vs0.5\_designmatrix.csv refers to the above design matrix):
  
  ```
  R10LB_0vs0.5 <- read.csv2("~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv", sep=",", check.names = F)
  R10LB_0vs0.5 <- janitor::remove_empty(R10LB_0vs0.5, which = "cols")
  R10LB_0vs0.5_design <- read.csv2("~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv", sep=",")
  
  rownames(R10LB_0vs0.5) <- R10LB_0vs0.5$`1`
  R10LB_0vs0.5 <- R10LB_0vs0.5[,-1]
  
  #colnames(R10LB_0ugvs0.5ug) <- R10LB_0ugvs0.5ug[1,]
  #R10LB_0ugvs0.5ug <- R10LB_0ugvs0.5ug[-1,]
  
  R10LB_0vs0.5 <- data.matrix(R10LB_0vs0.5)
  
  R10LB_0vs0.5 <- DGEList(counts = R10LB_0vs0.5, group = R10LB_0vs0.5_design$Group)
  str(R10LB_0vs0.5)
  
  Group <- as.vector(as.character(R10LB_0vs0.5_design$Group))
  Replicate <- as.vector((R10LB_0vs0.5_design$Replicate))
  R10LB_0vs0.5_designmatrix <- model.matrix(~0+Group+Replicate)
  
  keep <- filterByExpr(R10LB_0vs0.5)
  R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]
  
  R10LB_0vs0.5_Disp <- estimateDisp(R10LB_0vs0.5_filter, R10LB_0vs0.5_designmatrix)
  
  CONTRASTS <- makeContrasts( Group1vs2 = Group1-Group2,
                              levels = R10LB_0vs0.5_designmatrix)
  
  R10LB_0vs0.5fit <- glmFit(R10LB_0vs0.5_Disp, contrast = CONTRASTS[,1])
  R10LB_0vs0.5lrt <- glmLRT(R10LB_0vs0.5fit)
  R10LB_0vs0.5TT <- topTags(R10LB_0vs0.5lrt, n=nrow(R10LB_0vs0.5_Disp))
  write.csv(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.csv")
  saveRDS(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.RDS")
  
  ```
  
  **Volcano plot:**
  
  ```
  ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = PValue)) + geom_point()
  
  p4 <- ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = -log10(PValue))) + geom_point() + theme_minimal()
  
  #Add vertical lines for logFC thresholds, and one horizontal line for the p-value threshold
  p5 <- p4 + geom_vline(xintercept = c(-0.6, 0.6), col = "red") + geom_hline(yintercept = -log10(0.05), col = "red")
  
  # add a column of NAs
  R10LB_0vs0.5_comparison$diffexpressed <- "NO"
  # if logFC > 0.6 and PValue < 0.05, set as "UP" 
  R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC > 0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "UP"
  # if logFC < -0.6 and PValue < 0.05, set as "DOWN"
  R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC < -0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "DOWN"
  
  p4 <- ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed)) + geom_point() + theme_minimal()
  
  p5 <- p4 + geom_vline(xintercept=c(-0.6, 0.6), col="red") + geom_hline(yintercept=-log10(0.05), col="red")
  
  p6 <- p5 + scale_color_manual(values=c("blue", "black", "red"))
  p6 <- p5 + scale_colour_manual(values = mycolors)
  
  R10LB_0vs0.5_comparison$delabel <- NA
  R10LB_0vs0.5_comparison$delabel[R10LB_0vs0.5_comparison$diffexpressed != "NO"] <- R10LB_0vs0.5_comparison$gene_symbol[R10LB_0vs0.5_comparison$diffexpressed != "NO"]
  
  ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) + geom_point() + theme_minimal() +geom_text()
  
  library(ggrepel)
  ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) +geom_point() + theme_minimal() +geom_text_repel() + scale_color_manual(values=c("blue", "black", "red")) + geom_vline(xintercept=c(-0.6, 0.6), col="red") +geom_hline(yintercept=-log10(0.05), col="red")
  
  ```
  
  This is the volcano plot I currently have with all of my data:
  ![volcano plot](https://i.sstatic.net/tidIM.png)
  
  Can you explain how some of these functions filtered my data (i.e. filterByExpr, estimateDisp, makeContrasts, glmFit, glmLRT)? Why are some of my data points that change from 0 TPM in one condition to some value in another condition showing up on my plot while others are not? Or should I accept that my data analysis is valid or "correct"?
  
  **Tags:** r, rna-seq, bioconductor, edger, limma
  
  > **Comments:**
  > **username_0**: I heard that low expression genes are often filtered out in a first past due to a high chance of false positives, so they could be dropped from filterByExpr. High variance within a condition might undermine the statistical significance of the mean difference across conditions. Is there a way to look at the variance of my  data, and is there a threshold for variance between my replicates that causes them to be filtered out by filterByExpr or other functions?  
  > **username_1**: It is a long time that I don't use R. However, I had a sinilar problem once. Using not cartesian coordinates, ggplot cut several points out of the plot. If you plan to continue using R, search about this problem. Otherwise, if you can accede to it go with MATLAB that has the best tool for volcano plot, by my opinion. It never fails.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2021-12-24. Score: 1
  
  edgeR functions do not do any filtering, i.e., they do not remove any genes from the results or from the plots. The only filtering was that which you did explicitly yourself by:
  
  ```
  keep <- filterByExpr(R10LB_0vs0.5)
  R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]
  
  ```
  
  For your data, `filterByExpr` will keep any gene that is expressed at a minimal level in at least 7 samples.
  You could relax that requirement, to allow genes to be expressed in only 5 or 6 samples say, by reducing the `large.n` argument to `filterByExpr`.
  
  
  ---
  
  ### Answer by username_3 on 2021-03-12. Score: 0
  
  *[disclaimer: I understand DESeq2, not EdgeR]*
  
  Most likely, low expression values (as in raw counts, prior to any normalisation such as TPM). Expression that has high dispersion across all samples is less likely to be biologically significant.
  
  It's easier to work out what's going on with an MA plot, with average gene expression on the X axis, and log2 fold change on the Y axis.
  
  I'm used to DESeq2 results, so I'd guess the construction would be something like this (assuming that the 'baseMean' field exists in EdgeR results):
  
  ```
  R10LB_0vs0.5lrt %>%
      ggplot() +
      aes(y = logFC, col = padj, x=log(baseMean)) +
      geom_point()
  
  ```
  
  It's possible that EdgeR doesn't have anything like baseMean in its default results. As a proxy for this, you could use the count totals from the input matrix (e.g. `R10LB_0vs0.5 %>% rowSums() %>% log2()`).
  
  > **Comments:**
  > **username_0**: Do you mean that the TPM values are too low? Because I believe the ones that have been filtered out have a similar value to other "0" tpm genes that were significant. Also if I use the average of two conditions, will the genes I pull out be less significant because the MA plot will not account for variance between replicates? Thanks for the answer  
  
  ---
  
  