Contents

• Probability distribution of counts
  • Negative binomial distribution
    • Properties
    • Formulation as Gamma-Poisson
  • Count data
• Differential expression analysis
  • Dataset
    • Model (aka design) matrix
  • Normalization
    • Naive normalization
    • Median-of-ratios
    • Trimmed mean of M values (TMM)
    • Gene length-dependent metrics
  • Model
    • DESeq2
      • Sharing dispersion information across genes
• Knowledge base-driven pathway analysis
  • Over-representation analysis (ORA)
  • Functional class scoring (FCS)
    • GSEA
  • Pathway-topology (PT)-based approaches
    • SPIA
• References

Probability distribution of counts

Negative binomial distribution

\(K \sim \text{NB}(r, p)\): number of failures until \(r\) successes have occurred; \(p\) is the probability of success

Notes and References

1. The negative binomial may alternatively be formulated as \(K + r = Y \sim \text{NB}(r, p)\) = the number of the trial on which the \(r\)th success occurs. See Wackerly’s Mathematical Statistics with Applications, 7th Edition, or Rice’s Mathematical Statistics and Data Analysis, 3rd Edition.
2. Wikipedia uses “success” and “failure” oppositely of the formulation presented here. Swap \(p\) and \(1 - p\) for the equations to match.
3. For a reference that matches the formulation used here, see https://www.johndcook.com/negative_binomial.pdf.

Properties

where \(\alpha = \frac{1}{r}\) is called the dispersion parameter.

Then, we can re-express \(r\) and \(p\) in terms of \(\mu\) and \(\sigma\), or \(\mu\) and \(\alpha\):

Now, we can equivalently parameterize the negative binomial distribution as follows:

• Mean and variance: \(K \sim \text{NB}(\mu, \sigma^2)\)

  \[P(K = k) = {k + \frac{\mu^2}{\sigma^2 - \mu} - 1 \choose k} \left(\frac{\mu}{\sigma^2} \right)^{\frac{\mu^2}{\sigma^2 - \mu}} \left(\frac{\sigma^2 - \mu}{\sigma^2} \right)^k\]

• Mean and dispersion: \(K \sim \text{NB}(\mu, \alpha)\)

  \[P(K = k) = {k + \frac{1}{\alpha} - 1 \choose k} \left(\frac{1}{1 + \alpha \mu} \right)^{\frac{1}{\alpha}} \left(\frac{\alpha \mu}{1 + \alpha \mu} \right)^k\]

Formulation as Gamma-Poisson

Idea: \(K\) is a Poisson distribution where the mean of the Poisson distribution is proportional to a Gamma-distributed random variable.

\(K \sim \text{NB}(r, p)\) is equivalent to \(K \mid \Lambda = \lambda \sim \text{Poi}(s\lambda)\) where \(\Lambda \sim \text{Gamma}(a = r, \theta = \frac{1-p}{sp})\)

• \(p = \frac{1}{s\theta + 1}\)
• The \(a\) parameter in the Gamma distribution is not the same as the \(\alpha\) dispersion parameter.

Count data

Let \(N\) be the total number of RNA transcripts in a sample and \(p_i\) be the real proportion of those transcripts belonging to gene \(i\).

Consider a sequencing run that samples \(n\) of those transcripts (i.e., \(n\) is the total number of reads). Note that \(n\) is usually very large (\(n > 10^5\)) and \(p_i\) is usually very small (\(p_i < 0.01\)). (For example, a TPM value of \(10000\) corresponds to a \(p_i\) of \(0.01\).)

Let \(k_i\) be the observed read counts of gene \(i\).

Notes

1. For the variance under the Binomial model, the approximation \(n p_i (1 - p_i) \approx n p_i\) holds because \(p_i\) is small.
2. Since \(n\) is large and \(p_i\) is small, the Poisson distribution accurately approximates the Binomial distribution, and we see that the means and variance under both models are the same.
3. The symbol \(p_i\) used here is not the same as the symbol \(p\) used in the previous section describing the negative binomial distribution, hence the use of \(\phi_i\) in the table.

Negative binomial interpretation

• Gamma-Poisson interpretation: In the Binomial / Poisson models, we assume that \(p_i\) is fixed (constant) for gene \(i\) across all samples in the same condition. In a Gamma-Poisson model, we assume that \(p_i\) is the value of a random variable \(P_i\) whose distribution across samples in the given condition follows a Gamma distribution. [DESeq paper]
  • We can think of a fixed \(p_i\) as corresponding to technical replicates: samples are measurements of the same underlying population of \(N\) transcripts, of which a proportion \(p_i\) belong to gene \(i\). Then, there should be no overdispersion, and the counts should follow a Poisson distribution.
  • We can think of a varying \(p_i\) as corresponding to biological replicates: samples are measurements of different mice, etc., where the proportion of transcripts belonging to gene \(i\) varies over the population of mice, etc., according to a Gamma distribution. [Robinson & Oshlack] [Bioramble blog]

• Direct interpretation: Each read is a trial, and the read (trial) is “successful” if it does not come from gene \(i\). The negative binomial models the number of “failed” reads (reads from gene \(i\)) until \(r\) “successful” reads (reads not from gene \(i\)) have been observed, given that \(\phi_i = 1 - p_i\) is the probability of observing a “successful” read. Finally,

  \[\frac{r_i(1 - \phi_i)}{\phi_i} = \frac{r_i}{\phi_i} p_i = \frac{\text{# of reads not from gene $i$}}{\text{proportion of reads not from gene $i$}} p_i = (\text{# total reads}) p_i = np_i\]
• Empirical evidence: Empirically, the variability of read counts is larger than the Binomial and Poisson distributions allows and is better approximated by a Negative Binomial distribution. [DESeq paper] [Bioramble blog]

Differential expression analysis

Dataset

Dataset size

• \(n\): number of genes
• \(m\): number of samples
• \(c\): number of conditions (including intercept)

Count matrix: \(K \in \mathbb{N}^{n \times m}\)

• Rows: genes
• Columns: samples
• \(K_{ij}\): number of sequencing reads mapped to gene \(i\) in sample \(j\)
  • When referring to the count matrix as a random variable (e.g., see the DESeq2 model), a capital letter \(K\) is used. When referring to the observed count matrix (e.g., a real dataset), a lower-case \(k\) is used.
  • Let \(N_j = \sum_{i=1}^n k_{ij}\) be the total number of mapped reads in sample \(j\).

Model (aka design) matrix

\(X \in \{0,1\}^{m \times c}\)

• Rows: samples
• Columns: conditions
• Example: Each sample is a patient
  • Intercept: The intercept parameter is used to model gene expression for a healthy, non-drugged patient.
  • Drug: The patient was administered the drug (1) or not (0).
  • Diseased: The patient is diseased (1) or healthy (0).

If the normalized count \(q_{ij}\) of gene \(i\) is modeled as

then the \(\beta_{il}\) parameters give the \(\log_2\) fold change due to factor \(l\) above control, where the control is modeled with just the intercept term. In the example above, sample 1 (\(x_{11} = x_{21} = 0\)) would be considered the control with \(\log_2 q_{i1} = \beta_{i0}\).

Then, if we consider sample 2 (\(x_{12} = 1, x_{22}= 0\)), for example, we see that \(\beta_{i1}\) gives the \(\log_2\) fold change in expression of gene \(i\) for a patient treated with the drug versus a control patient.

Similarly, if we consider sample 4 (\(x_{14} = 1, x_{24}= 1\)), for example, we see that \(\beta_{i12}\) is the \(\log_2\) fold change due to the interaction of the drug with disease.

Normalization

Raw alignment counts are often not directly comparable across samples or across genes due to the following biases:

• library size (sequencing depth): larger library sizes result in higher counts for the entire sample.
• uniquely high counts of a few genes can take up a large proportion of the sequencing “real estate,” leading to relatively lower counts for the remaining genes. [Robinson & Oshlack]
• gene length: longer genes are more likely to be “sampled” (sequenced) by RNA-seq.
  • This is an important concern when comparing the expression of different genes within the same sample (e.g., whether the expression of GLUT2 (glucose transporter 2) or GLUT4 (glucose transporter 4) higher in a particular sample).
  • This is not a significant concern for differential gene expression analysis, since gene length should exert the same relative bias across samples.

This article presents a few popular methods. The first three methods (naive normalization, median-of-ratios, and TMM) estimate a sample-specific factor \(s_j\) yielding normalized counts \(q_{ij} = \frac{\mathbb{E}(K_{ij})}{s_j}\) that can be compared across samples. There also exist methods for estimating sample- and gene-specific factors \(s_{ij}\) such that \(q_{ij} = \frac{\mathbb{E}(K_{ij})}{s_{ij}}\), but they are not discussed here. Additionally, spike-in controls (known quantities of artificial fragments) can be used to specifically address technical variability in library preparation and sequencing and are especially important when the assumption that the majority of the genes are not differentially expressed is invalid. However, spike-in controls may suffer from their own reliability issues and rely on the strict assumption that technical biases affect spike-ins in the same way as they do genes. [Risso et al.] The gene-length dependent metrics (RPKM, FPKM, and TPM) estimate normalized proportions best used for comparisons within samples.

For further readings, consider the following review articles:

• Dillies, M.-A. et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Briefings in Bioinformatics 14, 671-83 (2013). doi:10.1093/bib/bbs046.
• Liu X. et al. Normalization Methods for the Analysis of Unbalanced Transcriptome Data: A Review. Frontiers in Bioengineering and Biotechnology 7 (2019). doi:10.3389/fbioe.2019.00358.

Naive normalization

Assumptions / limitations

• Only corrects for library size
• Each gene in a sample is normalized by the same factor: \(s_{ij} = s_j\)

Method: Divide by the total count of (mapped) sequencing reads of a sample, then multiply by the mean total count across all samples. The multiplication step ensures that the normalized values are counts (as assumed by models like DESeq2) rather than proportions.

Problem: This potentially reduces the power of the method to pick out differentially expressed genes, and in extreme cases, might lead to false positives. Consider the following dataset where gene C is knocked-out in sample 3. Using the normalized counts \(q_{ij}\), one might erroneously conclude that genes A, B, D, and E are upregulated in sample 3 relative to samples 1 and 2.

More robust alternatives to normalizing to the total count include normalizing to the upper quartile (75th percentile) or median (50th percentile) of nonzero counts.

Median-of-ratios

This is the default normalization method in DESeq and DESeq2 (DESeq2::estimateSizeFactors()).

Assumptions / limitations

• Corrects for library size and uniquely high counts
• Most genes are not differentially expressed

Method: Consider a pseudo-reference sample \(K^R\) whose counts for each gene are obtained by taking the geometric mean across samples. The size factor \(\hat{s}_j\) for sample \(j\) is computed as the median of the ratios of the \(j\)-th sample’s counts to those of the pseudo-reference:

Example: Using the same example as above, the normalized counts of genes A, B, D, and E are the same across samples 1, 2, and 3, as desired.

Trimmed mean of M values (TMM)

This is the default normalization method in the edgeR R package (edgeR::calcNormFactors()).

A trimmed mean is the average after removing the upper and lower x% of the data. The TMM procedure is doubly trimmed, by log-fold-changes \(M_{ij}^r\) (sample \(j\) relative to sample \(r\) for gene \(i\)) and by absolute intensity (\(A_i\)). By default, we trim the \(M_i\) values by 30% and the \(A_i\) values by 5%. [Adapted from Robinson & Oshlack]

Symbols

• \(M_{ij}^r = \log_2 \frac{K_{ij} / N_j}{K_{ir} / N_r}\): gene-wise log-fold changes (sample \(j\) relative to sample \(r\) for gene \(i\))
• \(A_{ij}^r = \log_2 \sqrt{\frac{K_{ij} \cdot K_{ir}}{N_j \cdot N_r}}\): absolute expression (log-geometric mean of expression in samples \(j\) and \(r\))

Let \(M^*\) and \(A^*\) be the percentile thresholds for trimming genes by log-fold change and by absolute expression, respectively. For RNA-seq, Robinson & Oshlack propose default values of \(M^* = 30\) and \(A^* = 5\). Let \(\mathcal{Q}\) be the quantile function, i.e., \(\mathcal{Q}(v, p)\) returns the \(p\)th quantile of vector \(v\). Then, the set of untrimmed genes is

Finally, the size factors are calculated as

where the weights are the inverse of estimated asymptotic variances

In the edgeR source code, the size factor produced by the median-of-ratios method is further divided by the library size, but the TMM size factor is not. Why?

Gene length-dependent metrics

Symbols

• \(X_t\): number of reads or fragments mapping to transcript \(t\)
• \(N = \sum_{t \in T} X_t\): total number of mapped reads
• \(\tilde{l}_t = l_t - m + 1\): effective length of a transcript \(t\), i.e., the number of positions in a transcript in which a read of length \(m\) can start
• \(p_t = \frac{m_t}{\sum_{t \in T} m_t}\): relative abundance of transcript \(t\), where \(m_t\) is the copy number of transcript \(t\) in the sample
  • \(M = \sum_{t \in T} m_t\) generally cannot be inferred from sequencing data but can be measured using qPCR?

*RPKM and FPKM values are not comparable across samples because they do not sum (over all unique transcripts in a sample) to a fixed value: \(\sum_{t \in T} \text{RPKM}_t\) is not constant. By normalizing to total effective (with respect to gene length) number of mapped reads instead of total number of mapped reads \(N\), TPM always sums to \(10^6\) for a given sample. Note that an equivalent estimate of relative abundance is \(\hat{p}_t = \frac{\text{RPKM}_t}{\sum_{r \in T} \text{RPKM}_r}\).

• Since RPKM and FPKM are not comparable across samples, there is not a defined value corresponding to 1 transcript per sample. However, as a very rough estimate, values in the range of 1-10 have been suggested. [Pachter’s blog post] By definition, a relative abundance value \(p_t = \frac{1}{M}\) (or a TPM of \(\frac{ {10}^6 }{M}\)) corresponds to 1 transcript per sample.

However, even the relative abundance values \(p_t\) (and consequently TPM values) are not necessarily comparable across experiments due to the potential bias introduced by high counts of a few genes, as explored in the example in the naive normalization section. In short, consider samples 1 and 2 that are identical, except that one gene (say, gene C) is completely knocked-out in sample 2. Because the relative abundance values \(\hat{p}_t\) sum to one (by construction), the relative abundance values in sample 2 will be higher than those in sample 1 for all genes except gene C, for which \(p_C = 0\) in sample 2.

For a brief comparison of the three metrics, consider

1. Harold Pimentel’s blog post
2. Lior Pachter’s 2013 CSHL Keynote

Model

DESeq2

Generalized linear model (GLM) with logarithmic link

Hierarchical construction of negative binomial

• \(R_{ij} \sim \text{Gamma}(a_{ij} = \frac{q_{ij}^2}{v_{ij}}, \theta_{ij} = \frac{v_{ij}}{q_{ij}})\): normalized count of transcripts from sample \(j\) belonging to gene \(i\); the distribution is over biological replicates from experimental condition \(\rho(j)\)
  • \(\mathbb{E}(R_{ij}) = a_{ij} \theta_{ij} = q_{ij}\)
  • \(\text{Var}(R_{ij}) = a_{ij} \theta_{ij}^2 = v_{ij}\)
• \(K_{ij} \mid R_{ij} \sim \text{Poi}(s_j R_{ij}) \rightarrow K_{ij} \sim \text{NB}(\mu_{ij}, \alpha_i)\)
  • \(\mu_{ij} = \mathbb{E}(K_{ij}) = a_{ij} s_j \theta_{ij} = s_j q_{ij}\) [Eq. (2), DESeq paper]
  • \(\sigma^2_{ij} = \text{Var}(K_{ij}) = a_{ij} s_j \theta_{ij} + a_{ij} s_j^2 \theta_{ij}^2 = s_j q_{ij} + s_j^2 v_{ij}\) [Eq. (3), DESeq paper]
  • \(\alpha_i = \frac{1}{a_{ij}} = \frac{v_{ij}}{q_{ij}^2}\)

Parameters (in order of estimation) [DESeq2 vignette]

• \(s_{ij}\): gene- and sample-specific normalization factor
  • By default, DESeq2 uses only a sample-specific normalization factor by assuming \(s_{ij} = s_j\) for all \(i = 1, ..., n\). This will account for differences in total reads (sequencing depth). See the median-of-ratios method described under the Normalization section.
  • To account for effects of GC content, gene length, and other gene-specific properties during the sample preparation (e.g., PCR amplification) and sequencing processes, consider using packages/methods like cqn or EDASeq to calculate gene-specific normalization factors.
• \(q \in \mathbb{R}^{n \times m}\): \(q_{ij}\) is a normalized count of transcripts from sample \(j\) belonging to gene \(i\).
  • For a given gene \(i\), DESeq2 uses the same estimate \(\hat{q}_{i, \rho(j)}\) for all samples from the same experimental condition. In other words, for all samples \(\tilde{j}\) such that \(\rho(\tilde{j}) = \rho(j)\), where \(\rho(j)\) is the experimental condition of sample \(j\), we have \(q_{i\tilde{j}} = q_{i, \rho(j)}\).
  • \(\hat{q}_{i, \rho(j)} = \frac{1}{m_{\rho(j)}} \sum_{\tilde{j}: \rho(\tilde{j}) = \rho(j)} \frac{k_{ij}}{\hat{s}_j}\). Estimate as the average of normalized counts from samples in the same experimental condition.
• \(\alpha_i\): dispersion values for each gene
  • See next section for estimation method.
  • Note that \(v_{ij}\), \(a_{ij}\), and \(\theta_{ij}\) are not necessarily explicitly estimated but can be computed from \(\alpha_i\), \(q_{ij}\), and \(s_{ij}\).
• \(\beta \in \mathbb{R}^{n \times c}\): \(\beta_{il}\) is the \(\log_2\) fold change due to experimental condition \(l\).
  • DESeq2 fits a zero-centered normal distribution to the observed distribution of MLE (maximum-likelihood estimate) \(\beta\) values and uses that as a prior for final MAP (maximum a posteriori) estimation of \(\beta\). The normal prior leads to “shrunken” logarithmic fold change (LFC) estimates, which are more accurate when the amount of (Fisher) information for a gene is low (e.g., when the gene has very low counts and/or high dispersion.)

Sharing dispersion information across genes

When sample sizes are small, dispersion estimates \(\alpha_i\) are highly variable (noisy). DESeq2 addresses this problem by assuming that “genes of similar average expression strength have similar dispersion.” The counts \(K_{ij}\) are first fit to the negative binomial model (as a generalized linear model) using MLE. Next, DESeq2 fits a smooth curve regressing the MLE dispersion estimate against the mean of normalized counts. Finally, the MLE estimates are shrunk towards the smooth curve by treating the smooth curve as a prior and generating MAP estimates of the dispersion values. [DESeq2 paper]

Knowledge base-driven pathway analysis

Let \(G\) denote the set of all genes considered, with \(\abs{G} = N\). For example, if we consider all protein-coding human genes, \(N \approx 20438\) (according to statistics from Ensembl).

Input

• Expression matrix
  • Rows: genes \(G\)
  • Columns: samples labeled by phenotype (e.g., control1, control2, cancer1, cancer2, …)
• Pathways (gene sets)
  • In the sections below, I use the symbol \(S\) to denote a gene set of interest.
  • Following [Khatri et al.], the sections below focus “on methods that exploit pathway knowledge in public repositories such as GO … rather than on methods [such as WGCNA] that infer pathways from molecular measurements.” Consequently, the analyses are only as good as the pathways annotated in existing databases; see [Tomczak et al.] and [Haynes et al.].

Over-representation analysis (ORA)

Idea: “statistically evaluates the fraction of genes in a particular pathway found among the set of genes showing changes in expression” [Khatri et al.]

Assumptions

1. Each gene is independent of other genes.
2. Each pathway is independent of other pathways.

Limitations

1. Ignores real values associated with each gene, since it only considers whether a gene is part of some gene set or not.
2. Only uses inforation from the most significantly differentially expressed genes.

Method

1. Choose set \(\delta\) of differentially expressed genes (DEGs). Common threshold: fold-change > 2 and \(q\)-value < 0.05.
2. For each pathway gene set \(S\) (of size \(\abs{S} = N_S\)), test for enrichment of DEGs. The null hypothesis is that the set of DEGs is not enriched for genes in the pathway, i.e., the proportion \(\frac{\abs{\delta \cap S}}{\abs{\delta}}\) of pathway genes in the set of DEGs is the same as the proportion \(\frac{\abs{S}}{N}\) of pathway genes in the set of all genes.
   • Khatri et al. (Text S2.3) argues that ORA approaches test competitive null hypotheses because they compare “the proportion of differentially expressed genes in a pathway with the proportion of differentially expressed genes not in the pathway.” However, the second “proportion” should be the proportion of pathway genes in the set of all genes.
   • Test sample: \(x = \abs{\delta \cap S}\). \(\abs{\delta}\) genes are drawn from the set \(\delta\) of all DEGs.
   • Null distribution: \(\abs{\delta}\) genes are drawn from the set \(G\) of all genes.

     • Binomial: assumes genes are drawn with replacement. \(X \sim \mathrm{Binom}(n = \abs{\delta}, p = \frac{N_S}{N})\).

       \[p = P(X \geq x) = 1 - \sum_{k=0}^{x-1} P(X = k) = 1 - \sum_{k=0}^{x-1} {\abs{\delta} \choose k} \left(\frac{N_S}{N}\right)^k \left(1 - \frac{N_S}{N} \right)^{\abs{\delta} - k}\]

     • Hypergeometric (Fisher’s exact test): assumes genes are drawn without replacement. \(X \sim \mathrm{Hyper}(N = N, n = \abs{\delta}, K = N_S)\), where the notation follows Wikipedia.

       \[p = P(X \geq x) = 1 - \sum_{k=0}^{x-1} P(X = k) = 1 - \sum_{k=0}^{x-1} \frac{ {N_S \choose k}{N - N_S \choose \abs{\delta} - k}}{ {N \choose \abs{\delta}}}\]

Functional class scoring (FCS)

Idea: “weaker but coordinated changes in sets of functionally related genes (i.e., pathways) can also have significant effects” [Khatri et al.]

Method

1. Compute gene-level statistics. Examples:
   • correlation with phenotype
   • t-statistic
   • fold-change
2. Aggregate gene-level statistics into a single pathway-level statistic. Examples:
   • Kolmogorov-Smirnov statistic
   • sum, mean, or median of gene-level statistics
   • Wilcoxon rank sum
3. Assess the statistical significance of the pathway-level statistic.
   • Competitive null hypothesis: “genes in a pathway are at most as often differentially expressed as the genes not in the pathway” [Khatri et al., Text S2.3]
     • Null distribution: permute gene labels
   • Self-contained null hypothesis: “no genes in a given pathway are differentially expressed” [Khatri et al., Text S2.3]
     • Null distribution: permute class (phenotype) labels for each sample. This null incorporates the correlation structure of genes. This is reasonable if we assume that most pathways are not affected by different phenotypic conditions.

Assumptions

1. Each pathway is independent of other pathways.
2. [Only for the self-contained null hypothesis] The correlation structure of genes is preserved across phenotypic conditions.

GSEA

Acronym: Gene Set Enrichment Analysis

Algorithm

1. Compute gene-level statistic: fold-change of average expression between conditions, correlation between expression and phenotype, etc.
2. Rank genes from most-to-least differentially-expressed according to their gene-level statistic.
   • Let \(r_i\) denote the gene-level statistic of the \(i\)th ranked gene \(g_i\), where \(i \in [1, N]\).
3. Use a “random walk” to calculate enrichment score (ES).
   • Use labels \(y_i \in \{0, 1\}\) to denote whether ranked gene \(i\) is in the gene set \(S\) of interest.

   • Old method [GSEA 2003]: equal weights at every step; \(\text{ES}\) is a standard Kolmogorov-Smirnov statistic.

     \[\begin{gathered} p_\text{hit} = \sqrt{\frac{N - N_S}{N_S}} \\ p_\text{miss} = -\sqrt{\frac{N_S}{N - N_S}} \\ \text{ES} = \max_{n \in [1, N]} \left\lvert \sum_{i=1}^n y_i p_\text{hit} + (1 - y_i) p_\text{miss} \right\rvert \end{gathered}\]
   • New method [GSEA 2005]: steps are weighted according to each gene’s gene-level statistic; \(\text{ES}\) is a weighted Kolmogorov-Smirnov statistic.

     • Formulation to match the notation of the old method:

       \[\begin{aligned} p_\text{hit}(i) &= \begin{cases} \frac{\abs{r_i}^p}{N_R} & i \in S \\ 0 & \text{otherwise} \end{cases}, & N_R = \sum_{j=1}^N y_j \abs{r_j}^p \\ p_\text{miss} &= -\frac{1}{N - N_S} \\ \text{ES} &= \max_{n \in [1, N]} \left\lvert \sum_{i=1}^n y_i p_\text{hit}(i) + (1 - y_i) p_\text{miss} \right\rvert \end{aligned}\]

     • Formulation from the 2005 paper:

       \[\begin{aligned} P_\text{hit}(S, i) &= \sum_{g_j \in S,\, j \leq i} \frac{\abs{r_j}^p}{N_R}, & N_R = \sum_{g_j \in S}^N \abs{r_j}^p \\ P_\text{miss}(S, i) &= \sum_{g_j \notin S,\, j \leq i} \frac{1}{N - N_S} \\ \text{ES}(S) &= \max_{i \in [1, N]} \left\lvert P_\text{hit}(S, i) - P_\text{miss}(S, i) \right\rvert \end{aligned}\]
       • \(P_\text{hit}(S, i)\) and \(P_\text{miss}(S, i)\) are the fraction (with respect to genes \(i = 1, ..., i\)) of genes in / not in \(S\) weighted by their gene-level statistics.

     • \(p\) is a hyperparameter controlling the weight of each step in the random walk.

       • The authors recommend a default of \(p = 1\). [GSEA 2005 supplemental]
       • When \(p=0\), the new method becomes nearly identical to the old method with equal weights at every step: \(p_\text{hit}(i)\) becomes \(\frac{1}{N_S}\).
   • Note that the “random walk” always returns to 0:

     • Old method:

       \[\begin{aligned} \sum_{i=1}^{N} y_i p_\text{hit} + (1 - y_i) p_\text{miss} &= N_S p_\text{hit} + (N - N_S) p_\text{miss} \\ &= \sqrt{N_S (N - N_S)} - \sqrt{(N - N_S) N_S} \\ &= 0 \end{aligned}\]

     • New method:

       \[\begin{aligned} \sum_{i=1}^{N} y_i p_\text{hit}(i) + (1 - y_i) p_\text{miss} &= \left( \sum_{i=1}^N y_i \frac{\abs{r_i}^p}{N_R} \right) + (N - N_S) \left(- \frac{1}{N - N_S}\right) \\ &= \left( \frac{1}{N_R} \sum_{i=1}^N y_i \abs{r_i}^p \right) - 1 \\ &= \frac{1}{N_R} N_R - 1 \\ &= 0 \end{aligned}\]
4. Estimate significance with permutation test.
   1. Permute column labels of expression matrix.
   2. Perform steps 1-3 of the GSEA algorithm with the permuted expression matrix (including re-computing gene-level statistics and re-ranking genes) and compute an enrichment score.
      • [New method only] Discard this enrichment score sampled from the null distribution if its sign does not match the sign of the observed \(\text{ES}(S)\). Separately considering positively and negatively scoring gene sets is necessary because

        the use of weighted steps could cause the distribution of observed enrichment scores to be asymmetric in cases where many more genes are correlated with one of the two phenotypes. Why?

   3. Estimate nominal \(p\)-value for \(\text{ES}(S)\) from the null distribution.
   4. [New method only] Compute a normalized enrichment score \(\text{NES}(S)\) that adjusts for variation in gene set size.
5. Correct for multiple hypothesis testing (each gene set is a hypothesis).

Pathway-topology (PT)-based approaches

Idea: Take advantage of knowledge (e.g., KEGG, MetaCyc, Reactome, PantherDB) of how (subcellular localization; activation/inhibition; etc.) genes in a pathway interact, not just the knowledge that they are in the same pathway.

Method: Generally follows same three steps as functional class scoring but incorporates additional pathway information in computing gene-level and pathway-level statistics.

Limitations

• [T]rue pathway topology is dependent on the type of cell due to cell-specific gene expression profiles and [the] condition being studied. [Khatri et al.]

SPIA

Acronym: Signaling Pathway Impact Analysis

Idea

The impact analysis combines two types of evidence: (i) the over-representation of DE genes in a given pathway and (ii) the abnormal perturbation of that pathway, as measured by propagating measured expression changes across the pathway topology. These two aspects are captured by two independent probability values, \(P_{NDE}\) and \(P_{PERT}\). [SPIA paper]

DEG enrichment analysis: \(P_{NDE}\)

Any of the existing ORA or FCS approaches can be used to calculate \(P_{NDE}\), as long as this probability remains independent of the magnitudes of the fold-changes.

Pathway perturbation analysis: \(P_{PERT}\)

Define the following symbols:

• \(g_1, ..., g_{N_S}\): genes in the pathway \(S\)
• \(\Delta E(g_i)\): \(\log_2\) fold-change in expression of gene \(g_i\)
• \(N_{ds}(g_i)\): number of genes downstream of \(g_i\) in the pathway
• \[\beta_{i,j} = \begin{cases} -1 & \text{gene } j \text{ inhibits gene } i \\ 1 & \text{gene } j \text{ activates gene } i \\ 0 & \text{otherwise} \end{cases}\]
• \(PF(g_i)\): perturbation factor of gene \(i\) (this is the gene-level statistic)

Linear system

Observe that the perturbation factor of gene \(j\) is distributed among all its downstream interacting partners. (This is reminiscent of the PageRank system of equations.)

The net accumulations of the perturbations from other genes on gene \(i\) is \(\text{Acc}(g_i) = PF(g_i) - \Delta E(g_i)\). The total accumulated perturbation in the pathway \(S\) is then computed as the pathway-level statistic

Finally, \(P_{PERT} = P(T_A \geq t_A \mid H_0)\) is calculated as the \(p\)-value of the pathway-level statistic \(t_A\) following a bootstrap procedure described in the supplemental text of the [SPIA paper].

Combining DEG and pathway analyses: \(P_G\)

Let the observed / computed values of \(P_{NDE}\) and \(P_{PERT}\) be \(p_{nde}, p_{pert}\), respectively with \(p_{nde} \cdot p_{pert} = c\).

Since \(P_{NDE}\) and \(P_{PERT}\) are \(p\)-values, under the null hypothesis they are uniformly distributed on the interval \([0, 1]\). The support of their joint distribution is therefore the unit square. Since \(P_{NDE}\) and \(P_{PERT}\) are independent, points with the same probability lie as \((p_{nde}, p_{pert})\) lie along the hyperbola defined by \(P_{NDE} P_{PERT} = c\). Then, the probability of obtaining a set of \(p\)-values as extreme or more extreme than \((p_{nde}, p_{pert})\) is the area under and to the left of the hyperbola (see Figure S1 in the [SPIA paper]):

References

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

   • Depending on the experimental situation, Poisson seems appropriate for technical replicates and Negative Binomial may be appropriate for the additional variation observed from biological replicates.

   • Presents the “trimmed mean of M values” (TMM) library count normalization method.
2. Risso, D. et al. Normalization of RNA-seq data using factor analysis of control genes or samples. Nature Biotechnology 32, 896–902 (2014). doi:10.1038/nbt.2931
   • Explores the reliability of spike-in controls and presents a normalization strategy remove unwanted variation (RUV) (not explained here).
3. Lipp, J. Why sequencing data is modeled as negative binomial. Bioramble (2016). https://bioramble.wordpress.com/2016/01/30/why-sequencing-data-is-modeled-as-negative-binomial/.
4. Anders, S. & Huber, W. Differential expression analysis for sequence count data. Genome Biology 11, R106 (2010). doi:10.1186/gb-2010-11-10-r106.
   • DESeq paper.
5. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014). doi:10.1186/s13059-014-0550-8.
   • DESeq2 paper.
6. Love, M. I., Anders, S. & Huber, W. Analyzing RNA-seq data with DESeq2. https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html (2020).
   • DESeq2 vignette.
7. Pachter, L. Models for transcript quantification from RNA-Seq. arXiv:1104.3889 [q-bio, stat] (2011). http://arxiv.org/abs/1104.3889.
   • Derives general likelihood model for RNA-seq reads, likelihood estimation of \(p\) and \(\alpha\) values, and inference using the EM-algorithm.
8. Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L. & Wold, B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods 5, 621–628 (2008). doi:10.1038/nmeth.1226.
   • One of the original RNA-seq papers; introduces RPKM metric. See Wikipedia, Lior Pachter’s *Seq chronology, and this blog post.
9. Trapnell, C. et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology 28, 511–515 (2010). doi:10.1038/nbt.1621.
   • Cufflinks paper; introduces FPKM metric.
10. Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011). doi:10.1186/1471-2105-12-323.
    • RSEM paper; introduces TPM metric.
11. Pachter, L. Estimating number of transcripts from RNA-Seq measurements (and why I believe in paywall). Bits of DNA (2014). https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/.
    • Explains why the FPKM value of a single gene cannot be converted to “transcripts per cell” without knowing the FPKM values of all genes.
12. Holmes, S. & Huber, W. Modern Statistics for Modern Biology. (Cambridge University Press, 2018). https://web.stanford.edu/class/bios221/book/index.html.
    • Chapter 4 derives the negative binomial model as a hierarchical Gamma-Poisson model.
    • Chapter 8 describes the DESeq2 model.
13. Khatri, P., Sirota, M. & Butte, A. J. Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges. PLoS Computational Biology 8, e1002375 (2012). doi:10.1371/journal.pcbi.1002375.
    • Classifies existing pathway analysis methods (as of 2012) into one of three approaches: ORA, FCS, or PT; analyzes the assumptions and limitations of each approach.
14. Mootha, V. K. et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nature Genetics 34, 267–273 (2003). doi:10.1038/ng1180.
    • Preliminary GSEA paper.
15. Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. 102, 15545–15550 (2005). doi:10.1073/pnas.0506580102.
    • Robust GSEA paper. Describes software package GSEA-P and the creation of the Molecular Signature Database (MSigDB).
16. Tomczak, A. et al. Interpretation of biological experiments changes with evolution of the Gene Ontology and its annotations. Scientific Reports 8, 5115 (2018). doi:10.1038/s41598-018-23395-2.

    • Our analysis suggests that GO evolution may have affected the interpretation and possibly reproducibility of experiments over time.

17. Haynes, W. A., Tomczak, A. & Khatri, P. Gene annotation bias impedes biomedical research. Scientific Reports 8, 1362 (2018). doi:10.1038/s41598-018-19333-x.

    • Collectively, our results provide an evidence of a strong research bias in literature that focuses on well-annotated genes instead of those with the most significant disease relationship in terms of both expression and genetic variation.

18. Tarca, A. L. et al. A novel signaling pathway impact analysis. Bioinformatics 25, 75–82 (2009). doi:10.1093/bioinformatics/btn577.
    • SPIA paper.