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
#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