text: user3225: I want to see the expression of a gene in a group of patient amongst the entire cohort using my RNA-Seq data. While I can do a differential expression analysis with limma or DESeq2, I want to see how much each sample from my cohort expresses the gene. 
  The plan is to plot a waterfall plot (as defined in this paper). However, when I used TPM and rlog(CPM) they gave me very different looking graphs (Not just the shape, which is expected, but the ranking of expression). My concerns are, for TPM, the batch effect with two or three lots in one cohort, and, for rlog(CPM), the results don't match with prediction but is adjusted for batch effect. Which one would be a right or better choice?
  Upper one is with plot with. Lower one is plot with voom transformed counts (log(CPM)). I didn't include all cases since there are too many. This look more or less what it is like in the one including every case. Target is the condition group in which I am interested. Differential expression analysis definitely showed that it is upregulated.
                       logFC       CI.L       CI.R    AveExpr         t      P.Value    adj.P.Val         B
  GeneA  1.6576257  0.9143490  2.4009024  5.8713219  4.391796 1.646627e-05 0.0002581676  2.484413
  
  Working steps
  TPM is basically from kallisto
  log(CPM) calculation
  #Import data from tximport
  txi = tximport(files = read.hd5s,
                       type = "kallisto",
                       tx2gene = tx2gene,
                       countsFromAbundance = "lengthScaledTPM",
                       ignoreTxVersion = TRUE)
  raw_data = DGEList(txi$counts)
  # Filter low count genes
  keep <- filterByExpr(raw_data)
  raw_data <- raw_data[keep,,keep.lib.sizes=FALSE]
  # Normalise counts
  raw_data = calcNormFactors(raw_data)
  # log(CPM) by voom
  voom_data = voom(raw_data, design = NULL, normalize="none", plot = TRUE)
  # Batch effect correction
  voom_expr_psva = psva(as.matrix(voom_data$E),Batch)
  
  Then I plot graph with TPM and log(CPM) respectively
  
  
  
  user5077: To add to what @gringer, when you do use TPM, the normalization done is for both library size and gene length. When you use rlog, the normalization is done via median normalization (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/). See below for a quick example about the rlog, and size factor normalization.
  dds=makeExampleDESeqDataSet(betaSD=0.5,sizeFactors =runif(12))
  rld = rlog(dds)
  ## you can see the data is variance stablized, i.e lower counts are converted
  ### see the size factor calculated 
  rld$sizeFactor
  ###
  boxplot(assay(rld))
  ## look at unnormalized distributions
  boxplot(log2(counts(dds)+1))
  ## now we normalize and see
  dds = estimateSizeFactors(dds)
  # the size factors should be similar
  table(sizeFactors(dds) == rld$sizeFactor)
  ### here you can see the medians are similar but its a bit different from rlog
  boxplot(log2(counts(dds,normalize=TRUE)+1))
  ###
  
  To answer your question specifically, the sizeFactor is unaware of any batch effect. It assumes median expression of all genes to be similar and normalize using that. If you think there is a batch effect, it has to be present across many genes, as opposed to one gene. I would check it with a PCA plot (https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/plotPCA) and then test for differential expression by including the batch effect in the model.
  
  user73: For comparing the counts of different samples from DESeq2, Michael Love recommends using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation:
  ## [assuming "dds <- DESeq(dds)" has been run]
  dds.counts <- assay(vst(dds, blind=FALSE));
  ## [This makes the scaled counts resemble more closely the actual read counts]
  dds.quantile99 <- quantile(dds.counts[dds.counts > min(dds.counts)], 0.99);
  dds.counts <- ((dds.counts - min(dds.counts)) / (dds.quantile99-min(dds.counts))) * (dds.quantile99);
  
  The matrix dds.counts should now represent a model-normalised approximation of the log of counts, with column names as the sample names, and row names as the gene names. The default behaviour of the VST changes zero counts to an arbitrary positive number, which is why I've modified the output so that count values of zero will be changed to 0 in log space (i.e. a count of 1). To plot results for a single gene, subset based on the row name:
  gene.VST <- sort(dds.counts["geneOfInterest",]);
  plot(gene.VST);
  
  
original_text: 3225: <p>I want to see the expression of a gene in a group of patient amongst the entire cohort using my RNA-Seq data. While I can do a differential expression analysis with limma or DESeq2, I want to see how much each sample from my cohort expresses the gene. </p>
  
  <p>The plan is to plot a waterfall plot (as defined in <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4093310/" rel="nofollow noreferrer">this paper</a>). However, when I used TPM and rlog(CPM) they gave me very different looking graphs (Not just the shape, which is expected, but the ranking of expression). My concerns are, for TPM, the batch effect with two or three lots in one cohort, and, for rlog(CPM), the results don't match with prediction but is adjusted for batch effect. Which one would be a right or better choice?</p>
  
  <p>Upper one is with plot with. Lower one is plot with voom transformed counts (log(CPM)). I didn't include all cases since there are too many. This look more or less what it is like in the one including every case. Target is the condition group in which I am interested. Differential expression analysis definitely showed that it is upregulated.</p>
  
  <pre><code>                     logFC       CI.L       CI.R    AveExpr         t      P.Value    adj.P.Val         B
  GeneA  1.6576257  0.9143490  2.4009024  5.8713219  4.391796 1.646627e-05 0.0002581676  2.484413
  </code></pre>
  
  <p>Working steps
  TPM is basically from kallisto</p>
  
  <p>log(CPM) calculation</p>
  
  <pre><code>#Import data from tximport
  txi = tximport(files = read.hd5s,
                       type = "kallisto",
                       tx2gene = tx2gene,
                       countsFromAbundance = "lengthScaledTPM",
                       ignoreTxVersion = TRUE)
  raw_data = DGEList(txi$counts)
  # Filter low count genes
  keep &lt;- filterByExpr(raw_data)
  raw_data &lt;- raw_data[keep,,keep.lib.sizes=FALSE]
  # Normalise counts
  raw_data = calcNormFactors(raw_data)
  # log(CPM) by voom
  voom_data = voom(raw_data, design = NULL, normalize="none", plot = TRUE)
  # Batch effect correction
  voom_expr_psva = psva(as.matrix(voom_data$E),Batch)
  </code></pre>
  
  <p>Then I plot graph with TPM and log(CPM) respectively</p>
  
  <p><a href="https://i.stack.imgur.com/f1t1V.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/f1t1V.png" alt="Waterfall plot using TPM"></a>
  <a href="https://i.stack.imgur.com/jxRIp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/jxRIp.png" alt="Waterfall plot using voom"></a></p>
   5077: <p>To add to what @gringer, when you do use TPM, the normalization done is for both library size and gene length. When you use rlog, the normalization is done via median normalization (<a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/" rel="nofollow noreferrer">https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/</a>). See below for a quick example about the rlog, and size factor normalization.</p>
  
  <pre><code>dds=makeExampleDESeqDataSet(betaSD=0.5,sizeFactors =runif(12))
  rld = rlog(dds)
  ## you can see the data is variance stablized, i.e lower counts are converted
  ### see the size factor calculated 
  rld$sizeFactor
  ###
  boxplot(assay(rld))
  ## look at unnormalized distributions
  boxplot(log2(counts(dds)+1))
  ## now we normalize and see
  dds = estimateSizeFactors(dds)
  # the size factors should be similar
  table(sizeFactors(dds) == rld$sizeFactor)
  ### here you can see the medians are similar but its a bit different from rlog
  boxplot(log2(counts(dds,normalize=TRUE)+1))
  ###
  </code></pre>
  
  <p>To answer your question specifically, the sizeFactor is unaware of any batch effect. It assumes median expression of all genes to be similar and normalize using that. If you think there is a batch effect, it has to be present across many genes, as opposed to one gene. I would check it with a PCA plot (<a href="https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/plotPCA" rel="nofollow noreferrer">https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/plotPCA</a>) and then test for differential expression by including the batch effect in the model.</p>
   73: <p>For comparing the counts of different samples from DESeq2, Michael Love <a href="https://support.bioconductor.org/p/120764/#120843" rel="nofollow noreferrer">recommends</a> using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation:</p>
  
  <pre><code>## [assuming "dds &lt;- DESeq(dds)" has been run]
  dds.counts &lt;- assay(vst(dds, blind=FALSE));
  ## [This makes the scaled counts resemble more closely the actual read counts]
  dds.quantile99 &lt;- quantile(dds.counts[dds.counts &gt; min(dds.counts)], 0.99);
  dds.counts &lt;- ((dds.counts - min(dds.counts)) / (dds.quantile99-min(dds.counts))) * (dds.quantile99);
  </code></pre>
  
  <p>The matrix <em>dds.counts</em> should now represent a model-normalised approximation of the log of counts, with column names as the sample names, and row names as the gene names. The default behaviour of the VST changes zero counts to an arbitrary positive number, which is why I've modified the output so that count values of zero will be changed to 0 in log space (i.e. a count of 1). To plot results for a single gene, subset based on the row name:</p>
  
  <pre><code>gene.VST &lt;- sort(dds.counts["geneOfInterest",]);
  plot(gene.VST);
  </code></pre>
  