Skip to contents

Read RSEM gene-level results as a gene symbol by sample table

Source: R/rsem_read.R
rsem_read.Rd

Reads the *.genes.results files listed in a metadata CSV and returns one expression value per gene symbol and sample, either the RSEM expected count or the TPM.

Usage

rsem_read(data = "data", metadata = "metadata.csv", species, expr)

Arguments

data

Directory that contains the RSEM files.

metadata

Path to the metadata CSV.

species

"human" or "mouse". Required. The function stops if the gene IDs are of the other species.

expr

"count" (expected_count) or "tpm" (TPM). Required.

Value

A data frame with one row per gene symbol (as row names) and one column per sample, in the order of the metadata. A message reports how many genes received a symbol.

Input

metadata is a CSV with the columns sample_id and ref; other columns are ignored. ref is the name of the sample's RSEM file inside data (that is, the file read is file.path(data, ref)).

Every file needs the columns gene_id and expected_count or TPM (whichever expr selects); other columns are ignored. All files must contain the same genes in the same order, without duplicates.

Gene symbols

Symbols come from org.Hs.eg.db (human) or org.Mm.eg.db (mouse). The version suffix of an Ensembl gene ID (ENSG00000141510.18) is removed before the lookup, because the OrgDb keys have no version. The rules are:

  • An Ensembl ID with several symbols gets the symbol of its lowest Entrez ID.

  • Rows that end up with the same symbol are summed (counts and TPM are both additive).

  • A gene without a symbol (not in the OrgDb, or not an Ensembl ID, such as a spike-in or an ID ending in _PAR_Y) keeps its gene_id from the file, so no expression is lost. If such an ID equals a gene symbol, the function stops instead of merging the two.

Rows appear in the order in which they first occur in the files. Values are not rounded and TPM is not renormalized.

Examples

if (requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
  ex <- system.file("extdata", "example", package = "rsemverse")
  tpm <- rsem_read(
    data = file.path(ex, "data"),
    metadata = file.path(ex, "metadata.csv"),
    species = "human",
    expr = "tpm"
  )
  head(tpm)
}
#> 
#> 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
#>              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