Pathway scores#

gsva(), ssgsea(), combined_z_scores() and plage() turn the expression values of every sample into one score per gene set, as the four methods of the R package GSVA compute them, from a samples x genes matrix or a cell type of a cube. The scores, a PathwayScores, are samples x sets, so that fit_linear_model() tests them with a design built from their samples, and plot_pathway_scores() draws them as a heatmap.

Per-sample scores of gene sets: GSVA, ssGSEA, combined z-scores and PLAGE, as the R package GSVA computes them.

Each function turns the expression values of every sample into one score per gene set, so that pathways can be compared across samples, drawn as a heatmap (plot_pathway_scores()) and tested with the designs and contrasts of the linear models, as GSVA’s guide tests its scores with limma. They compute the four methods of GSVA 2.0.0 (Haenzelmann, Castelo and Guinney 2013) with its options and defaults: gsva(), the method of that paper; ssgsea(), single sample GSEA (Barbie et al. 2009); combined_z_scores() (Lee et al. 2008); and plage() (Tomfohr, Lu and Kepler 2005). validation/pathway_scores.py compares them with GSVA run as a black box. GSVA is licensed under the GPL: the methods are implemented from the papers and GSVA’s documentation, and what these leave open was found by running GSVA on chosen inputs, as marked “found by probing” below.

The data:

Values are samples x genes: a 2-D array with one gene id per column, or a group (cell type) of a Cube, of which the samples with at least min_cells cells and some counts are scored, over the genes that edgeR’s filterByExpr keeps (filter_by_expression()), as log2 counts per million after TMM normalization (values="log_cpm", edgeR’s cpm(log = TRUE)) or as counts (values="counts"), the values GSVA’s guide passes. NaN marks a missing value; infinite values are refused. Each method’s docstring says how it treats missing values.

gsva(), combined_z_scores() and plage() leave out the genes whose values are all equal (constant), before the sets are matched, as GSVA does; ssgsea() keeps them, as GSVA does (found by probing).

Gene sets:

The sets are a GeneSets or a mapping of set name to gene ids, matched to the genes’ ids exactly. A set is scored when min_size to max_size of its genes are among the genes scored (GSVA’s minSize and maxSize, counted after matching); the others are left out, as GSVA leaves them out. The result keeps the order of the sets.

GSVA (gsva()):

Each gene’s values are turned into the probability of a value at most as high under the gene’s distribution across samples, estimated by a kernel (Haenzelmann, Castelo and Guinney 2013): for gene i in sample j,

F_ij = (1/n) sum_k K(x_ij, x_ik),

a sum over the n samples in their order, with K the Gaussian kernel Phi((x_ij - x_ik) / h_i), h_i a quarter of the gene’s standard deviation, for continuous values; the Poisson kernel P(Y <= x_ij) for Y ~ Poisson(x_ik + 1/2) for counts; or, without a kernel, the indicator of x_ik <= x_ij, the empirical distribution function. GSVA evaluates Phi from a table of its values at multiples of 0.001: at |t| rounded down to such a multiple, and as one minus that value for negative t (found by probing; it moves F by less than 4e-4). kernel="auto" chooses the Poisson kernel for values stored as integers (an integer dtype, as GSVA chooses it for R’s integer matrices), the Gaussian kernel otherwise, and no kernel from min_samples_without_kernel samples on (found by probing).
The genes of each sample are ranked by the statistic z_ij = -log((1 - F_ij) / F_ij), from 1 for the lowest to p for the highest of the p genes; genes with equal statistics take their ranks in reverse order of the genes (found by probing). Each gene weighs |p/2 - r_ij|^tau, so that both ends of the ranking weigh most (exponent, GSVA’s tau).
A walk down the genes from rank p rises at each gene of a set by its share of the set’s weight and falls at each other gene by 1/(p - k), for k genes in the set, as in gene set enrichment analysis. With E+ and E- the sizes of the walk’s largest deviations above and below zero, the score is E+ - E- (score_type="difference", the default, GSVA’s maxDiff = TRUE), which is large for a set whose genes are high and low for one whose genes are low; E+ + E- ("kuiper", absRanking = TRUE), Kuiper’s statistic, large for a set whose genes are at either end; or the larger of the two deviations, with its sign ("max_deviation", maxDiff = FALSE; -E- when they are equal, found by probing).

A set with every gene, or whose genes all weigh zero (one gene at rank p/2), has a NaN score, as in GSVA. With missing values, each gene’s distribution is estimated from its observed samples and each sample ranks its observed genes, as GSVA’s use = "na.rm" does; the walk then falls by 1/(p_j - k_j) at each of the sample’s observed genes outside the set, with p_j and k_j counting the genes observed in sample j.

ssGSEA (ssgsea()):

Each sample ranks its genes by their values, and the walk down the ranking rises at each gene of a set by its rank to the power alpha (exponent), as a share of the set’s total, and falls at each other gene by 1/(p - k) (Barbie et al. 2009). The score is the sum of the walk over all p genes, not its extreme:

ES = sum_{g in G} w_g a_g / sum_{g in G} w_g - sum_{g not in G} a_g / (p - k),

with a_g the gene’s rank from the bottom of the walk. Genes with equal values weigh their average rank truncated to an integer, and the walk passes them in the order of the genes (found by probing). normalize=True divides every score by the range of all the scores, over sets and samples, as Barbie et al. did and GSVA does by default. Missing values are left out of each sample’s ranking, as GSVA’s use = "na.rm" does.

Combined z-scores (combined_z_scores()):

Each gene is standardized across samples, z = (x - mean) / sd, and a set’s score in a sample is sum_{g in G} z_g / sqrt(k) (Lee et al. 2008). A missing value leaves its gene out of the mean and standard deviation, and makes the score of every set that holds the gene NaN in that sample, as in GSVA.

PLAGE (plage()):

The standardized values of a set’s k genes form a k x n matrix, and the set’s scores are the first right singular vector of that matrix: the samples’ coordinates on the direction along which the set’s genes vary most together (Tomfohr, Lu and Kepler 2005). Missing values are refused, as GSVA refuses them.

Deliberate differences from GSVA:

PLAGE’s sign. A singular vector has no sign of its own; GSVA returns the one its singular value decomposition gives (found by probing), which depends on the linear algebra library. Here the scores rise with the mean standardized expression of the set’s genes: the sum of the set’s loadings (the first left singular vector) is positive, or, when it is zero, the loading of largest size, the first of equal ones, is positive. A set of one gene then has the gene’s standardized values, scaled to length one, where GSVA gives 1 in every sample.
GSVA with missing values: GSVA’s use = "na.rm" lets the walk of a sample fall by 1/(p - k_j) at its observed genes outside the set, p counting the genes missing in the sample too, so that the walk ends above zero; here it falls by 1/(p_j - k_j) and ends at zero. GSVA’s default, use = "everything", gives meaningless scores with the GSVA method (found by probing), and a sample without any value stops GSVA with an error; here such a sample has NaN scores.
ssGSEA with a set of every gene: its score is NaN and the normalization uses the other scores, where GSVA stops with an error.
kernel="auto" with an integer dtype and min_samples_without_kernel samples or more uses no kernel, as GSVA documents; GSVA 2.0.0 stops with an error. The Poisson kernel refuses values that are not counts, where GSVA stops with an error on negative values and scores fractional ones (found by probing).

Agreement with GSVA:

The scores are GSVA’s to rounding when the genes of each sample are ranked as GSVA ranks them. Ranks are decided by comparing the kernel estimates F of different genes, and two genes whose estimates are equal in exact arithmetic, such as genes with the same counts in different samples, can be ordered by the last bit of R’s normal and Poisson probabilities, which CytePool computes accurately but not always to the same bit: such genes, mostly genes with low counts, can take each other’s ranks.

References

Haenzelmann S, Castelo R, Guinney J (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, Froehling S, Chan EM, Sos ML, Michel K, Mermel C, Silver SJ, Weir BA, Reiling JH, Sheng Q, Gupta PB, Wadlow RC, Le H, Hoersch S, Wittner BS, Ramaswamy S, Livingston DM, Sabatini DM, Meyerson M, Thomas RK, Lander ES, Mesirov JP, Root DE, Gilliland DG, Jacks T, Hahn WC (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462, 108-112. doi:10.1038/nature08460

Lee E, Chuang HY, Kim JW, Ideker T, Lee D (2008). Inferring pathway activity toward precise disease classification. PLoS Computational Biology 4, e1000217. doi:10.1371/journal.pcbi.1000217

Tomfohr J, Lu J, Kepler TB (2005). Pathway level analysis of gene expression using singular value decomposition. BMC Bioinformatics 6, 225. doi:10.1186/1471-2105-6-225

Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences 102, 15545-15550. doi:10.1073/pnas.0506580102

Robinson MD, Oshlack A (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology 11, R25. doi:10.1186/gb-2010-11-3-r25

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. doi:10.12688/f1000research.8987.2

class cytepool.pathway_scores.PathwayScores(scores, samples, gene_sets, method, settings)#

Scores of gene sets in every sample.

Rows of scores are samples, in the order of samples, and columns are the sets scored, in the order of gene_sets, so that fit_linear_model() takes scores as its values, with a design built from samples. The arrays are read-only.

Variables:

scores – (n_samples, n_sets) float64 scores; NaN where a set has none in a sample.
samples – One row per sample: sample_id and, for a cube’s group, the columns of the cube’s sample table.
gene_sets – The sets scored, each with the genes that were used, in the order given, with the description and source of the collection.
method – "gsva", "ssgsea", "combined_z_scores" or "plage".
settings – The options of the method, the kernel that "auto" chose included, and for a cube how its values were taken.

Explained in: Gene-set tests

property n_genes#: The number of genes of each set that were used, in the order of gene_sets.

table()#

Return the scores as a table with one row per set and sample.

Returns:: A table with gene_set, n_genes, sample_id and score, the rows of each set together, in the order of gene_sets, and its samples in the order of samples.

Explained in: Gene-set tests

cytepool.pathway_scores.combined_z_scores(data, gene_sets, *, gene_ids=None, sample_ids=None, group=None, values='log_cpm', min_cells=MIN_CELLS, expression_filter=True, min_size=1, max_size=None)#

Score gene sets in every sample by their combined z-scores.

The method of Lee et al. (2008) as GSVA 2.0.0’s gsva(zscoreParam(...)) computes it: each gene is standardized across samples, and a set’s score in a sample is the sum of its genes’ z-scores divided by the square root of their number, which is a standard normal score when the genes are independent standard normal values.

Parameters:

data – (n_samples, n_genes) expression values on a log scale; NaN marks a missing value. Or a Cube, whose group is scored.
gene_sets – A GeneSets or a mapping of set name to gene ids.
gene_ids – The genes’ ids, as for gsva().
sample_ids – For an array, one id per row, as for gsva().
group – The group of a cube, as for gsva().
values – For a cube, "log_cpm" or "counts", as for gsva().
min_cells – For a cube, as for gsva().
expression_filter – For a cube, as for gsva().
min_size – The fewest of a set’s genes among the genes for the set to be scored, at least 1.
max_size – The most of them, at least min_size; None for no limit.

Returns:

The scores, samples x sets; NaN where a gene of the set is missing in the sample.

Raises:

TypeError – If an argument has the wrong type, values are not real numbers, gene ids are not strings, or an argument is or holds a numpy masked array.
ValueError – If the values are not 2-D or hold an infinite value, the ids do not match the values, an id names two genes, an option is out of range, no gene varies across samples, no gene of any set is among the genes, or no set has min_size to max_size genes.
KeyError – If the cube has no such group or features column.

Used in: Pathway and regulator activity

Explained in: Choosing a method, Gene-set tests

cytepool.pathway_scores.gsva(data, gene_sets, *, gene_ids=None, sample_ids=None, group=None, values='log_cpm', min_cells=MIN_CELLS, expression_filter=True, kernel='auto', min_samples_without_kernel=200, exponent=1.0, score_type='difference', min_size=1, max_size=None, n_threads=None)#

Score gene sets in every sample by gene set variation analysis (GSVA).

The method of Haenzelmann, Castelo and Guinney (2013) as GSVA 2.0.0’s gsva(gsvaParam(...)) computes it: each gene’s values are turned into ranks of the probabilities of its distribution across samples, estimated by a kernel, and each sample scores a set by a walk down the ranking of its genes, as gene set enrichment analysis scores a ranking (Subramanian et al. 2005). Scores lie in [-1, 1], or in [0, 2] with score_type="kuiper". The module docstring gives the method, its missing values and the deliberate differences from GSVA.

Parameters:

data – (n_samples, n_genes) expression values, such as log-CPM or, for the Poisson kernel, counts; NaN marks a missing value. Or a Cube, whose group is scored.
gene_sets – A GeneSets, such as one read by read_gmt(), or a mapping of set name to gene ids.
gene_ids – The genes’ ids, which the sets name: one per column of an array (None, NaN or an empty string for a gene without one); for a cube, the name of a column of its features table, "feature_id" by default or "symbol", or one id per feature.
sample_ids – For an array, one id per row; by default "1", "2", …
group – The group (cell type) of a cube to score; it may be left out when the cube has one.
values – For a cube, the values scored: "log_cpm", log2 counts per million after TMM normalization, or "counts".
min_cells – For a cube, the fewest cells a sample needs in the group to be scored, at least 0, as for limma_voom().
expression_filter – For a cube, whether to score only the genes that filter_by_expression() keeps.
kernel – How each gene’s distribution is estimated (GSVA’s kcdf): "gaussian" for continuous values, "poisson" for counts (whole numbers of at least 0), "none" for the empirical distribution function, or "auto".
min_samples_without_kernel – The number of samples from which "auto" uses no kernel, at least 1 (GSVA’s kcdfNoneMinSampleSize).
exponent – The power tau of the weights |p/2 - r|^tau of the ranks, at least 0 (GSVA’s tau).
score_type – "difference", "kuiper" or "max_deviation", as the module docstring describes (GSVA’s maxDiff and absRanking).
min_size – The fewest of a set’s genes among the genes scored for the set to be scored, at least 1 (GSVA’s minSize).
max_size – The most of them, at least min_size; None for no limit (GSVA’s maxSize).
n_threads – The most threads to use; None uses numba’s current setting. The result does not depend on it.

Returns:

The scores, samples x sets, of the sets with min_size to max_size genes, in the order of gene_sets.

Raises:

TypeError – If an argument has the wrong type, values are not real numbers, gene ids are not strings, or an argument is or holds a numpy masked array.
ValueError – If the values are not 2-D or hold an infinite value, the ids do not match the values, an id names two genes, an option is out of range, the Poisson kernel is given values that are not counts, no gene varies across samples, no gene of any set is among the genes, or no set has min_size to max_size genes.
KeyError – If the cube has no such group or features column.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(1)
>>> y = rng.normal(size=(6, 40))
>>> y[3:, :5] += 2.0  # the first five genes are higher in the last samples
>>> genes = [f"g{j}" for j in range(40)]
>>> result = gsva(y, {"up": genes[:5], "other": genes[20:30]}, gene_ids=genes)
>>> result.scores[:, 0].round(2).tolist()
[-0.54, -0.82, -0.69, 0.78, 0.34, 0.77]

Used in: Pathway and regulator activity

Explained in: Choosing a method, Gene-set tests

cytepool.pathway_scores.plage(data, gene_sets, *, gene_ids=None, sample_ids=None, group=None, values='log_cpm', min_cells=MIN_CELLS, expression_filter=True, min_size=1, max_size=None)#

Score gene sets in every sample by pathway level analysis of gene expression (PLAGE).

The method of Tomfohr, Lu and Kepler (2005) as GSVA 2.0.0’s gsva(plageParam(...)) computes it: the standardized values of a set’s genes form a genes x samples matrix, and the set’s scores are its first right singular vector, the samples’ coordinates on the axis along which the genes vary most together. The scores of a set have mean zero and length one across samples, and rise with the mean standardized expression of its genes, a sign chosen here where GSVA’s is arbitrary (see the module docstring).

Parameters:

data – (n_samples, n_genes) expression values on a log scale, without missing values. Or a Cube, whose group is scored.
gene_sets – A GeneSets or a mapping of set name to gene ids.
gene_ids – The genes’ ids, as for gsva().
sample_ids – For an array, one id per row, as for gsva().
group – The group of a cube, as for gsva().
values – For a cube, "log_cpm" or "counts", as for gsva().
min_cells – For a cube, as for gsva().
expression_filter – For a cube, as for gsva().
min_size – The fewest of a set’s genes among the genes for the set to be scored, at least 1.
max_size – The most of them, at least min_size; None for no limit.

Returns:

The scores, samples x sets.

Raises:

TypeError – If an argument has the wrong type, values are not real numbers, gene ids are not strings, or an argument is or holds a numpy masked array.
ValueError – If the values are not 2-D, hold an infinite or a missing value, the ids do not match the values, an id names two genes, an option is out of range, no gene varies across samples, no gene of any set is among the genes, or no set has min_size to max_size genes.
KeyError – If the cube has no such group or features column.

Used in: Pathway and regulator activity

Explained in: Choosing a method, Gene-set tests

cytepool.pathway_scores.ssgsea(data, gene_sets, *, gene_ids=None, sample_ids=None, group=None, values='log_cpm', min_cells=MIN_CELLS, expression_filter=True, exponent=0.25, normalize=True, min_size=1, max_size=None, n_threads=None)#

Score gene sets in every sample by single-sample GSEA (ssGSEA).

The method of Barbie et al. (2009) as GSVA 2.0.0’s gsva(ssgseaParam(...)) computes it: each sample ranks its genes by their values, and a set’s score is the sum of the walk of gene set enrichment analysis down that ranking, with the genes of the set weighted by their ranks to the power exponent. Only the order of the values within each sample matters. The module docstring gives the method, its ties and its missing values.

Parameters:

data – (n_samples, n_genes) expression values; NaN marks a missing value, which is left out of its sample’s ranking. Or a Cube, whose group is scored.
gene_sets – A GeneSets or a mapping of set name to gene ids.
gene_ids – The genes’ ids, as for gsva().
sample_ids – For an array, one id per row, as for gsva().
group – The group of a cube, as for gsva().
values – For a cube, "log_cpm" or "counts", as for gsva(); both rank the genes of a sample alike.
min_cells – For a cube, as for gsva().
expression_filter – For a cube, as for gsva().
exponent – The power alpha of the weights of the ranks, at least 0 (GSVA’s alpha).
normalize – Whether to divide every score by the range of all the scores (GSVA’s normalize).
min_size – The fewest of a set’s genes among the genes for the set to be scored, at least 1.
max_size – The most of them, at least min_size; None for no limit.
n_threads – The most threads to use; None uses numba’s current setting. The result does not depend on it.

Returns:

The scores, samples x sets, of the sets with min_size to max_size genes, in the order of gene_sets.

Raises:

TypeError – If an argument has the wrong type, values are not real numbers, gene ids are not strings, or an argument is or holds a numpy masked array.
ValueError – If the values are not 2-D or hold an infinite value, the ids do not match the values, an id names two genes, an option is out of range, no gene of any set is among the genes, no set has min_size to max_size genes, or normalize would divide by a range of zero.
KeyError – If the cube has no such group or features column.

Used in: Pathway and regulator activity

Explained in: Choosing a method, Gene-set tests

Pathway scores

Contents

Pathway scores#