生物统计常见概念，记录 P 值、方差分析和检验方法的学习理解。

Understand P-value

  • P value: a ==sum of probability== of those events ==rare and rarer==.

  • P-value Explained Clearly — Regression, PDF, Discrete

  • statQuest: p-value

  • P-values, False Discovery Rate (FDR) and q-values

    • Thus, if the difference is large but the variance is also large, then the difference may not be significant. On the other hand, a small difference coupled with a very small variance could be significant.
    • We use Anova tests to formalise this calculation.
    • The tests return a p-value that takes into account the mean difference and the variance and also the sample size. The p-value is a measure of how likely you are to get this compound data if no real difference existed.

  • p-value 不是“这个结果为真的概率”，而是在 null hypothesis 成立的情况下，观察到当前结果以及更极端结果的概率。

  • rare and rarer 这个说法很重要：它不是单点概率，而是一侧或两侧尾部概率的累积。

  • p-value 同时受 mean difference、variance、sample size 影响，所以不能直接等同于 effect size 或 biological importance。

  • 在 DEG 分析里，log2FC 大但方差/标准误也大时，p-value 可能仍然不显著。

FDR / q-value

  • https://stats.stackexchange.com/questions/252937/benjamini-hochberg-choosing-the-false-discovery-rate-q-value

    • FDR的值可以不是0.05，但有可能会受到审稿人的批评
    • Moreover, the formula in the original paper shows that you just select an upper bound for FDR, but that actual FDR is below the product between the chosen Q and the proportion of null hypotheses that are true, in the paper’s terminology): https://rss.onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1995.tb02031.x
    • https://onlinelibrary.wiley.com/doi/10.1111/psyp.12616) with a higher coverage than 95% and not including the null. Instead, with a Q above 0.05, you could find significance in comparisons with a p-value above 0.05 (and with a 95% unadjusted CI including the null). Moreover, referees might not accept a Q above 0.05, for analogy between FDR and type-I error.

  • RNA-seq/omics 数据通常会同时检验成千上万个 gene。单独使用 p < 0.05 会产生大量由多重检验带来的 false positives。

  • FDR/q-value 关注的是：在被判定为 significant 的结果中，预期 false positives 所占的比例。

  • FDR threshold 可以根据研究目的调整。0.05 更常见；0.1 有时用于 exploratory analysis，但需要说明理由。

  • DESeq2 结果里通常看 padj，也就是 multiple testing correction 后的 adjusted p-value。

1
2
res <- results(dds, alpha = 0.05)
sig <- res[which(res$padj < 0.05), ]

Anova test

  • https://www.youtube.com/watch?v=R7xd624pR1A

  • ANOVA 的直觉仍然是比较组间差异和组内变异。

  • 组间差异大，但组内变异也大，结果不一定 significant；组间差异小，但组内很稳定，也可能 significant。

  • 这一点和前面 p-value 的解释一致：variance 会直接影响统计检验结果。

PCA

  • What Is Principal Component Analysis (PCA) & How It Works? (analyticsvidhya.com)

    • Key Takeaways
    • Principal Component Analysis (PCA) is used to overcome feature redundancy in a data set. These features are low dimensional in nature. These features, a.k.a components, are a result of normalized linear combinations of original predictor variables.
    • The first component has the highest variance, followed by second, third, and so on. The components must be uncorrelated (remember orthogonal direction ?). See above.
    • Normalizing data becomes extremely important when the predictors are measured in different units.PCA works best on data sets having 3 or higher dimensions. Because, with higher dimensions, it becomes increasingly difficult to make interpretations from the resultant data cloud.

  • PCA 不用于判断某一个 gene 是否 differential expression，而是用于观察样本整体结构。

  • 在 RNA-seq 里常用于检查：

    • biological replicates 是否聚在一起；
    • treatment/control 是否分开；
    • 是否存在 outlier；
    • 是否存在 batch effect。

  • PC1 是解释最多 variance 的方向，PC2 是第二多的方向，并且和 PC1 orthogonal。

  • RNA-seq 的 PCA 一般不直接用 raw count。raw count 受 library size 影响，需要先 normalization，再使用 vst/rlog 等转换后的表达矩阵。

Modern statistics for modern biology

DEG/DE test

  • https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8

    • DESeq2 论文

  • https://hbctraining.github.io/DGE_workshop/lessons/05_DGE_DESeq2_analysis2.html

    • Lessons: https://github.com/hbctraining/DGE_workshop_salmon/blob/master/lessons/05_DGE_DESeq2_analysis2.md 提到的是01和05的部分。 这里01有讲解负二项分布, 05讲解DE test。
    • ==拟合模型后，使用公式估算每个样本组的系数及其标准误差==
    • 系数是各样本组折合变化对数的估计值。然而，这些估计值并没有考虑到我们观察到的低读数的巨大离散性。为避免出现这种情况，需要对模型计算出的对数折合变化进行调整
      • 我们首先是用负二项分布对数据建模，也就是确定其系数，但不理解这个是什么意思，好像是和GLM的建模有关，对于GLM来说，其检验的东西就是模型的系数。
    • Shrinking the log2 fold changes will not change the total number of genes that are identified as significantly differentially expressed.
      • The shrinkage of fold change is to help with downstream assessment of results. For example, if you wanted to subset your significant genes based on fold change for further evaluation, you may want to use shruken values. Additionally, for functional analysis tools such as GSEA which require fold change values as input you would want to provide shrunken values.
    • In our case is, the null hypothesis is that there is no differential expression across the two sample groups (LFC == 0) ^5b406c
      • 这里假设的是LFC服从正态分布，计算了LFC的均值和标准误，从而计算Wald Test的z值，来得到p值
    • The Wald test can also be used with continuous variables. If the variable of interest provided in the design formula is continuous-valued, then the reported log2 fold change is per unit of change of that variable. Wald 检验也可用于连续变量。如果设计公式中提供的相关变量是连续值，那么报告的 log2 折合变化就是该变量的单位变化 (per unit of change)
    • Thus, if the difference is large but the variance is also large, then the difference may not be significant.

  • DESeq2 对每个 gene 拟合 negative binomial GLM。condition 是模型里的变量，模型系数表示 condition 对表达量的影响。

  • 因为模型使用 log link，condition 的系数可以和 log fold change 对应起来：

1
2
fold change = treatment / control
log(fold change) = log(treatment) - log(control)
  • 因此，log normalized gene count 和 log fold change 的关系在于：模型是在 log scale 上估计组间差异。
  • Wald test 可以理解为检验 coefficient 是否显著偏离 0：
1
Wald statistic = estimated coefficient / standard error
  • DESeq2 结果中的关系大致是：
1
stat = log2FoldChange / lfcSE
  • LFC 大但 lfcSE 也大时，stat 不一定大，p-value 也不一定小。
  • continuous variable 的情况：如果 design formula 中的变量是连续值，reported log2 fold change 表示该变量每增加一个单位对应的 log2 fold change，而不是两组之间的 fold change。
  • LFC shrinkage 不改变 significant genes 的数量。它主要用于让 fold change 的估计更稳定，特别是 low count genes；后续排序、可视化、GSEA 等分析更适合使用 shrunken values。

如何理解正态分布，二项分布，泊松分布，负二项分布

  • 前面提及的lesson01 的这一篇在后面就提到了负二项分布。

  • Modern Statistics for Modern Biology - 1 Generative Models for Discrete Data (stanford.edu)

    • 这本数学书里在统计建模的章节详细讲了这一问题。

  • 正态分布：常用于连续变量和近似连续的误差项。

  • 二项分布：固定试验次数下的成功次数。

  • 泊松分布：常用于 count data，但要求 mean 和 variance 相等。

  • 负二项分布：可以处理 overdispersion，也就是 variance > mean。

  • RNA-seq read counts 常见 overdispersion，所以 DESeq2 使用 negative binomial distribution，而不是简单的 Poisson distribution。