2
$\begingroup$

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

Waterfall plot using TPM Waterfall plot using voom

Improve this question
$\endgroup$
4
  • $\begingroup$ in addition to TPM, I think kallisto also produces count-like data as output, which would be better for downstream analysis (presumably the default when tximport is used in 'kallisto' mode). The count-level data gives a better idea of uncertainty based on shot noise (i.e. the likelihood of missing a gene because its expression levels are low). $\endgroup$
    – gringer ♦
    Sep 17, 2019 at 21:40
  • $\begingroup$ I think there might be some misunderstanding. I actually used the count data as shown in this line: raw_data = DGEList(txi$counts), which is to generate a DGEList object from the count-like data you mentioned. $\endgroup$
    – Kent
    Sep 17, 2019 at 22:01
  • $\begingroup$ You're using the countsFromAbundance="lengthScaledTPM" argument on tximport, which looks suspiciously like a form of normalisation to me ("Whether to generate estimated counts using abundance estimates"). The import default for that, as used in the kallisto import example is countsFromAbundance="no". $\endgroup$
    – gringer ♦
    Sep 17, 2019 at 23:04
  • $\begingroup$ Thanks for pointing that out. I didn't notice since it was something shown to me by a bioinformatician when I first started learning R. Yes it looks like the reads are scaled to the library size and average transcript length as show in this tutorial. I guess raw_data = calcNormFactors(raw_data) was intended to renormalise to library size after filtering due to keep.lib.sizes = FALSE. Still trying to work this around my head. $\endgroup$
    – Kent
    Sep 18, 2019 at 13:16

2 Answers 2

Reset to default
3
$\begingroup$

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);
Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks! Sorry for not providing an example. I thought it was more about the practice of the community than a technical question so I omitted it. I use voom transformation followed by psva() in sva package but I assume a model-normalisation can also be done. The reason why I have this question is that using TPM shows a result closer to our expectation, while limma actually uses the voom transformation counts to detect differentially expressed genes. I am just not sure if I am misleading the reader if I use TPM just to show a better looking graph. $\endgroup$
    – Kent
    Sep 17, 2019 at 15:13
  • $\begingroup$ Specific examples are helpful in bioinformatics because there are a huge number of solutions available that are context-dependent. It'd be great if you could add this information to your question to improve the quality of the answers you get. $\endgroup$
    – gringer ♦
    Sep 17, 2019 at 20:18
  • $\begingroup$ Hi! I have added the graphs that I get and the results. I am not sure if that's an example that would help answering the question, but that's basically the problem I am facing: will using TPM misleading/exaggerating the upregulation and trick people in interpreting the data wrong? And is it better to use voom since they are variance stabled and batch effect corrected so that reader can better interpret the data e.g. observe any trend in clinical data rather than randomness that might lead to the difference in gene A expression level within the group? $\endgroup$
    – Kent
    Sep 17, 2019 at 21:22
  • $\begingroup$ It looks like you're doing a double normalization (i.e. using kallisto TPM values in addition to calcNormFactors(raw_data)). I'm not familiar with voom, so don't know if voom is attempting additional correction despite the use of the normalize="none" option. $\endgroup$
    – gringer ♦
    Sep 17, 2019 at 21:45
2
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.