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
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 iscountsFromAbundance="no"
. $\endgroup$raw_data = calcNormFactors(raw_data)
was intended to renormalise to library size after filtering due tokeep.lib.sizes = FALSE
. Still trying to work this around my head. $\endgroup$