Skip to contents

Remove lowly expressed or uninformative genes

Source: R/rsem_filter.R
rsem_filter.Rd

Keeps the genes that pass one of four common filters that need no sample groups. The table keeps its values and row order; only rows are removed. A message reports the rule and how many genes were kept.

Usage

rsem_filter(
  x,
  expr,
  method = NULL,
  min_value = NULL,
  min_samples = 0.2,
  n_top = 2000
)

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

"edger", "detected", "nonzero" or "variable". The default is "edger" for counts and "detected" for TPM.

min_value

For "edger", the count needed in a sample of median size (default 10). For "detected", the value a gene must reach (default 1 for TPM, 10 for counts).

min_samples

For "detected": in how many samples the gene must reach min_value. A number below 1 is read as a proportion of the samples (0.5 = half). Default 0.2.

n_top

For "variable": how many genes to keep. Default 2000.

Value

A data frame with the kept genes, in the same layout as x (genes as row names, samples as columns).

Methods

  • "edger" (counts only; the default for counts): the rule of edgeR::filterByExpr() with all samples treated as one group. A gene is kept when its count per million (CPM) reaches min_value / median library size * 1e6 in at least k samples and its total count is at least 15. k is the number of samples, or 10 + (n - 10) * 0.7 when there are more than 10. min_value is the count a gene needs in a sample of median size (default 10, the edgeR default).

  • "detected" (the default for TPM): a gene is kept when its value is at least min_value in at least min_samples samples. The defaults are TPM >= 1 or count >= 10, in at least 20% of the samples (at least one sample). This is the form used for bulk RNA-seq (for example, GTEx keeps genes above a TPM and a count threshold in at least 20% of samples) and for single-cell data ("expressed in at least n cells").

  • "nonzero": removes genes that are zero in every sample.

  • "variable": keeps the n_top genes with the largest variance of log2(value + 1) across samples (default 2000), as used before PCA, clustering or heatmaps. Genes that are zero everywhere are never kept. Rows stay in their original order, not in the order of variance.

Filters can be chained, for example "edger" and then "variable".

References

Chen Y, Lun ATL, Smyth GK (2016). From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research 5:1438.

GTEx Consortium (2020). The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science 369:1318-1330.

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")
  kept <- rsem_filter(counts, expr = "count")                 # edgeR rule
  tpm <- rsem_read(file.path(ex, "data"), file.path(ex, "metadata.csv"),
                   species = "human", expr = "tpm")
  rsem_filter(tpm, expr = "tpm", min_value = 1, min_samples = 2)
  rsem_filter(tpm, expr = "tpm", method = "variable", n_top = 10)
}
#> 
#> 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_filter (edger): kept 42 of 43 genes (CPM >= 1260 in >= 3 of 3 samples and total count >= 15); removed 1
#> 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_filter (detected): kept 43 of 43 genes (TPM >= 1 in >= 2 of 3 samples); removed 0
#> rsem_filter (variable): kept 10 of 43 genes (the 10 genes with the largest variance of log2(tpm + 1)); removed 33
#>              S1        S2        S3
#> GAPDH 18999.156 12154.876 12726.933
#> BRCA1 27363.142 30993.360 17916.614
#> KRAS   4751.162  2959.748  4991.503
#> CD3E   9980.064 14408.709  8682.311
#> CD4   16529.815 19362.216  9821.504
#> IL2    6139.603  6112.643 10347.432
#> IFNG  18732.716 13402.846 26547.789
#> AKT1   5231.922  7847.568  8278.526
#> CD274 32646.205 48340.620 28777.024
#> GZMB  15091.937  4814.445 13037.912