Makes samples comparable by correcting for sequencing depth and, for some methods, RNA composition. The genes and samples stay the same; only the values change. A message reports the method and the per-sample factors.
Arguments
- x
The gene by sample table, as returned by
rsem_read(): gene names as row names (or in the first column) and one numeric column per sample.- expr
"count"or"tpm": what the values are. Required.- method
For counts:
"cpm","tmm","upper_quartile","deseq2"or"quantile". For TPM:"quantile". Required.
Value
A data frame of normalized values in the same layout as x.
The per-sample factors are in the attribute "factors": library size
for "cpm", effective library size (library size times the
normalization factor) for "tmm" and "upper_quartile", and size
factors for "deseq2".
Methods
For counts (expr = "count"):
"cpm": counts per million,count / library size * 1e6. Corrects for depth only."tmm": CPM with TMM normalization factors (trimmed mean of M values; Robinson and Oshlack 2010), asedgeR::calcNormFactors()thenedgeR::cpm(). Also corrects for differences in RNA composition, for example when a few genes take up a large part of one sample."upper_quartile": CPM with upper-quartile factors (Bullard et al. 2010), asedgeR::calcNormFactors(method = "upperquartile")."deseq2": counts divided by the DESeq2 median-of-ratios size factors (Anders and Huber 2010), the values ofDESeq2::counts(dds, normalized = TRUE). Needs genes that are nonzero in every sample."quantile": quantile normalization, aslimma::normalizeQuantiles().
For TPM (expr = "tpm"), TPM already corrects for depth and gene length,
so only "quantile" is offered. It makes every sample have the same
distribution of values, which is a strong assumption: use it only when the
samples are expected to differ in a few genes.
The results are the same as those of the packages named above (checked in
the tests), computed here without them. Low-count genes affect TMM and the
size factors, so filter first (rsem_filter()). Methods for differential
expression (edgeR, DESeq2) need the raw counts, not these values. For a log
scale, pass the result to rsem_transform().
References
Robinson MD, Oshlack A (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology 11:R25.
Anders S, Huber W (2010). Differential expression analysis for sequence count data. Genome Biology 11:R106.
Bullard JH, Purdom E, Hansen KD, Dudoit S (2010). Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics 11:94.
Bolstad BM, Irizarry RA, Astrand M, Speed TP (2003). A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19:185-193.
Examples
if (requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
ex <- system.file("extdata", "example", package = "rsemverse")
counts <- rsem_read(file.path(ex, "data"), file.path(ex, "metadata.csv"),
species = "human", expr = "count")
tmm <- rsem_normalize(counts, expr = "count", method = "tmm")
attr(tmm, "factors")
rsem_transform(tmm, "log2") # log2(TMM-CPM + 1)
}
#> rsem_read: 3 samples, 44 genes -> 43 rows (species = human, expr = count, org.Hs.eg.db 3.23.1)
#> with symbol: 42 genes (1 symbols merged from several IDs); kept as original IDs: 2 genes
#> rsem_normalize (tmm, count): 43 genes x 3 samples; effective library sizes: S1 7478, S2 7861, S3 8821
#> rsem_transform: 43 genes x 3 samples; log2(value + 1)
#> S1 S2 S3
#> TP53 14.56699 14.87684 14.50084
#> GAPDH 15.14795 14.51079 14.58975
#> ACTB 12.74981 12.97996 12.96455
#> MYC 14.15594 13.69856 14.06808
#> EGFR 10.70782 10.24022 10.34949
#> BRCA1 14.26774 14.45468 13.67668
#> BRCA2 13.61011 13.67069 13.20867
#> KRAS 12.93365 12.25818 13.02466
#> PTEN 14.55902 14.74937 14.16029
#> CDKN2A 13.41196 13.37494 13.72579
#> RB1 14.65935 14.09447 14.38331
#> CD3E 14.27749 14.81451 14.09636
#> CD4 15.13728 15.37267 14.40608
#> CD8A 15.58080 15.21272 15.48841
#> IL2 13.65567 13.65654 14.42851
#> IFNG 13.79791 13.32215 14.32074
#> TNF 14.04618 14.17605 14.26368
#> VEGFA 13.20281 13.42585 13.25959
#> STAT3 14.59064 14.65293 14.81639
#> JAK2 15.10480 15.29719 14.66722
#> MTOR 13.49785 13.19092 13.54604
#> AKT1 12.08580 12.67782 12.76754
#> PIK3CA 14.77396 14.95981 15.26157
#> BRAF 13.95193 13.56864 14.18684
#> NRAS 12.90878 12.67782 13.04416
#> ERBB2 13.85110 13.91600 13.77329
#> ALB 13.82475 13.80498 13.68911
#> INS 13.20281 12.78572 13.17370
#> HBB 14.16647 13.95119 13.63871
#> CD19 15.41407 15.42367 15.49196
#> MS4A1 12.47275 13.04655 13.00490
#> PDCD1 14.56699 15.03826 14.80501
#> CD274 14.85431 15.42784 14.69215
#> CTLA4 13.32063 13.64226 13.11960
#> LAG3 15.43600 15.41109 15.62413
#> TIGIT 15.62728 15.24143 15.28633
#> FOXP3 12.08580 11.60627 11.84758
#> GZMB 13.56307 11.92220 13.37185
#> PRF1 16.40408 16.70964 16.29753
#> IL6 13.59460 13.83040 14.02947
#> AAAS 15.31106 15.21995 15.58307
#> ERCC-00002 15.15854 15.56288 15.41169
#> 10904 16.08132 15.70708 16.04037