Keep or remove groups of genes — rsem

Selects genes by built-in gene groups, so that no gene list has to be prepared. Gene types and chromosomes come from the same annotation package as rsem_read() (org.Hs.eg.db or org.Mm.eg.db); the other groups are fixed symbol patterns. Values are not changed. A message reports how many genes each group matched.

Usage

rsem_subset(x, species, keep = NULL, remove = NULL)

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.
species: "human" or "mouse". Required.
keep, remove: Gene groups: one or more of "protein_coding", "mito", "ribo", "hb", "sex" and "ercc". At least one of the two is required.

Value

A data frame with the selected genes, in the same layout and order as x.

Gene groups

"protein_coding": genes whose type is protein-coding in the annotation package (NCBI Gene). Many analyses use only these genes, for example to compare with proteomics, to reduce noise from lowly expressed non-coding RNAs and pseudogenes, or because a signature was built on them.
"mito": the 37 genes of the mitochondrial genome (chromosome MT; symbols MT- in human and mt- in mouse). The rRNA and tRNA genes have no Ensembl link in the OrgDb, so rsem_read() keeps them under their Ensembl IDs; a built-in list of the 37 Ensembl IDs (Ensembl release 116) finds them as well, together with the symbols the OrgDb gives to these IDs (the ID of MT-TS1 maps to MIR12136, a microRNA inside it). A high share of these reads is a sign of damaged cells or degraded RNA.
"ribo": cytosolic ribosomal protein genes (RPL*, RPS*, RPLP0-2, RPSA; for mouse Rpl*, Rps*), protein-coding only. Mitochondrial ribosomal genes (MRPL*, MRPS*), kinases such as RPS6KA1 and pseudogenes such as RPL10P1 are not included. They are highly expressed and often removed before clustering.
"hb": hemoglobin genes (HBA1, HBA2, HBB, HBD, HBE1, HBG1, HBG2, HBM, HBQ1, HBZ; for mouse Hba-*, Hbb-*, Hbq1*). protein-coding only. They signal blood contamination. HBEGF, hemoglobin pseudogenes and the locus control regions are not included.
"sex": genes on chromosome Y and XIST. They separate samples by sex and are often removed when sex is not of interest. Pseudoautosomal genes (on both X and Y, expressed in both sexes) are not included.
"ercc": ERCC spike-in controls (IDs starting with ERCC-).

keep is applied first (a gene is kept when it is in any of the listed groups), then remove. Genes are matched by symbol; rows that kept their Ensembl ID in rsem_read() are matched by that ID. Genes that the annotation package does not know are not in "protein_coding", "ribo", "hb" or "sex"; for example, chromosome Y genes that rsem_read() kept under their Ensembl ID are not found by "sex".

Examples

if (requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
  ex <- system.file("extdata", "example", package = "rsemverse")
  tpm <- rsem_read(file.path(ex, "data"), file.path(ex, "metadata.csv"),
                   species = "human", expr = "tpm")
  pc <- rsem_subset(tpm, "human", keep = "protein_coding")
  rsem_subset(tpm, "human", keep = "protein_coding", remove = c("mito", "ribo"))
}
#> rsem_read: 3 samples, 44 genes -> 43 rows (species = human, expr = tpm, org.Hs.eg.db 3.23.1)
#>   with symbol: 42 genes (1 symbols merged from several IDs); kept as original IDs: 2 genes
#> rsem_subset: kept 41 of 43 genes; keep protein_coding (genes per group: protein_coding 41)
#> rsem_subset: kept 41 of 43 genes; keep protein_coding; remove mito, ribo (genes per group: protein_coding 41, mito 0, ribo 0)
#>               S1         S2          S3
#> TP53   13599.421 16773.5477  12812.5954
#> GAPDH  18999.156 12154.8762  12726.9327
#> ACTB    3421.496  3993.3057   3916.4985
#> MYC    10197.296  7389.5612   9463.7881
#> EGFR    1129.441   812.5069    868.8631
#> BRCA1  27363.142 30993.3597  17916.6144
#> BRCA2  33191.910 34442.4861  24786.0383
#> KRAS    4751.162  2959.7477   4991.5028
#> PTEN   30402.758 34517.8455  22746.1721
#> CDKN2A 10066.129  9762.0579  12341.6385
#> RB1    30467.761 20493.3175  24818.2378
#> CD3E    9980.064 14408.7088   8682.3114
#> CD4    16529.815 19362.2158   9821.5043
#> CD8A   27778.200 21415.3146  25699.2824
#> IL2     6139.603  6112.6427  10347.4321
#> IFNG   18732.716 13402.8457  26547.7890
#> TNF    50822.554 55332.3292  58286.1612
#> VEGFA   4783.798  5555.8414   4907.9469
#> STAT3  17983.739 18683.5181  20742.9935
#> JAK2   19563.871 22243.2315  14248.0973
#> MTOR    9272.595  7458.0025   9456.6576
#> AKT1    5231.922  7847.5683   8278.5258
#> PIK3CA 85303.690 96547.6080 117974.2004
#> BRAF   20916.558 15956.0444  24279.3906
#> NRAS   13462.988 11413.8149  14585.7876
#> ERBB2   8271.326  8608.7315   7730.0064
#> ALB    35978.639 35311.6097  32302.8570
#> INS    12764.924  9512.0053  12339.1921
#> HBB    30291.348 25961.8272  20723.7546
#> CD19   41323.654 41391.7296  43020.5629
#> MS4A1  17662.341 26159.5923  25194.0129
#> PDCD1  20855.896 28769.2073  24261.2748
#> CD274  32646.205 48340.6198  28777.0238
#> CTLA4   5654.803  7031.8906   4852.0550
#> LAG3   38936.855 38079.2084  43754.7351
#> TIGIT  31889.919 24284.2576  24833.9997
#> FOXP3   2357.916  1682.5235   1971.6436
#> GZMB   15091.937  4814.4453  13037.9119
#> PRF1   55255.115 67948.8978  50620.5202
#> IL6    11415.408 13375.3517  15220.9500
#> AAAS   31996.429 26097.8198  32555.7921