glm模型参数量

2024-05-31 23:39:15

glm模型参数量

我非常愿意为大家解答关于glm模型参数量的问题。这个问题集合包含了一些复杂而有趣的问题，我将尽力给出简明扼要的答案，并提供进一步的阅读材料供大家深入研究。

1.广义线性模型有什么缺点？

2.转录组入门（7）：差异表达分析

glm模型参数量

广义线性模型有什么缺点？

广义线性模型（Generalized Linear Models，GLMs）是一种用于建立响应变量与预测变量之间关系的广泛应用的统计模型。虽然GLMs在许多领域都表现出色，但它们仍然存在以下一些缺点：

假设限制：GLMs的假设限制较多，如线性关系、等方差性、独立性等。如果实际数据不符合这些假设，则GLMs可能会失效或给出不准确的结果。

数据分布假设：GLMs通常假设响应变量服从特定的概率分布，如正态分布、泊松分布等。如果响应变量的实际分布与这些假设不符，则GLMs的拟合结果可能不准确。

局限性：GLMs只能处理连续型和离散型的响应变量，而不能处理其他类型的响应变量，如序数型、分类型和时间序列型的数据。

非线性关系：如果预测变量与响应变量之间存在非线性关系，则GLMs可能无法捕捉到这种关系，因为它们只能建立线性关系。

参数解释：GLMs的参数解释可能比较困难，因为不同的预测变量可能具有不同的比例因子、基线值和缩放因子，这使得参数解释变得更加复杂。

需要注意的是，这些缺点并不意味着GLMs是一种不好的模型，而是需要针对具体应用场景和数据特征来选择合适的模型。

转录组入门（7）：差异表达分析

最近看了不少统计基因组方面的资料，以为懂了，其实懵逼。实在是因为统计学基础太菜，似懂非懂，似是而非。记录下自己的理解，求轻喷。

GWAS是表型和基因型之间的相关性分析，然而这个相关性用什么统计方法？

GWAS模型的发展：

实际比例是否符合预期分离比例，若不符合则认为基因与表型相关。计算量大。

首先计算基因和表型之间的相关性（如pearson），再经过r换算为t值，进行t检验。计算量小。

直接将基因型x和表型y做回归拟合，即y=xb+e。

其中y是表型，x是SNP、Q矩阵和总均值的一个综合变量，e是残差。

GLM模型中，如果两个表型差异很大，但群体本身还含有其他的遗传差异（如地域等），则那些与该表型无关的遗传差异也会影响到相关性。MLM模型可以把群体结构的影响设为协方差，把这种位点校正掉。此外，材料间的公共祖先关系也会导致非连锁相关，可加入亲缘关系矩阵作为随机效应来矫正，即y=Xb+Zu+e，Z为亲缘关系矩阵，u为个体育种值，其他同GLM。

GLM是固定效应模型。MLM除了固定效应，还加入了随机效应，即亲缘关系矩阵。所谓的固定效应是有限水平的，易控制的，关心的是水平上的影响；而随机效应有很多水平，关心的是水平背后的群体（如均值，方差等）。

有种只可意会，不可言传的感觉。我百度了下，认为它解释的还可以：随机效应模型

MLM的矫正过于严格，会把一些真实相关的SNP标记也过滤掉，因此CMLM模型目的是重新检测到那些假阴性SNP标记。

方法是根据个体遗传关系的相似性将其分组，然后将压缩后的组当做协变量，替换原来的个体。组内个体的亲缘关系都是一样的。

问题来了，CMLM应该选择哪些SNP来计算亲缘关系矩阵，答案是使用所有跟表型相关的SNP（且排除了检测到的那个SNP）来构建亲缘关系矩阵的效果最好，这就是SUPER (Settlement of Kinship Under Progressively Exclusive Relationship, 逐步排他性亲缘关系解决方案)。

QTN（数量性状SNP），即控制QTL的SNP。

GWAS的瓶颈一是计算速度，二是统计准确性。FarmCPU能提升速度和准确性，首先把随机效应的亲缘关系矩阵（Kinship）转换为固定效应的关联SNP矩阵（S矩阵/QTNs矩阵），使计算速度大大加快；其次利用QTN矩阵当做协变量，重新做关联分析，提升准确率。

SUPER和FarmCPU都是把bin（一段区域，比如10kb）当做SNP单位，而不是单个SNP。

Blink是进阶版GWAS，也是为提高速度和准确率。如下图所示：先用上方的GLM模型获得QTNs，然后用右侧的GLM以QTNs当做协变量进行SNP检测，得到的SNP根据LD信息确定QTNs的信息（根据染色体实际位置来选择对应的bin大小），进而利用左侧的GLM以BIC（Bayesian information criterion）策略进行QTNs准确性检测，排除假设错误的部分，保留真实的QTNs，不断循环这一过程，直到检测到所有关联SNP（即QTNs）。

模型总结：

这里是以张志武老师《统计基因组学》课程整理，所以重点介绍的是他们课题组开发的模型。

GWAS常用软件：

GS模型的发展：

当控制表型的基因数量比较少，同时遗传力比较高的性状，可以用少量标记对表型进行很好的预测。

利用个体亲缘关系构建协变量矩阵，然后根据个体育种值对表型进行预测。

如果把GBLUP中构建协变量的个体亲缘关系矩阵换成SNP标记构成的关系矩阵，构建模型，然后对个体进行预测，就是rrBLUP的思路。

也就是把模型y=Xb+Zu+e变为y=Xb+Ms+e：

其中M是SNP构建的矩阵（替换个体亲缘关系矩阵Z），s就是标记（替换个体u）。

rrBLUP假设所有标记效应符合一个正态分布，同时标记效应的方差（σ^2）相等，这可能与基因的实际效应值不相符。

这时引入了贝叶斯的分析方法：预期控制表型的基因数目未知，基因效应值的分布未知。根据预先假定基因的数量和基因效应值分布的不同，建立了不同的贝叶斯模型，如BayesianA,B,C,Cπ,LASSO等。

其差别主要在于：标记效应是否符合相同分布；是否所有标记都有效应值；标记效应方差服从什么分布。

不同贝叶斯模型的假设和分布可参考之前的推文： # GS文献基因组选择技术在农业动物育种中的应用

用图形直观表示不同的贝叶斯方法，下图中π表示没有标记效应的比例。

岭回归的所有标记效应方差都相等；贝叶斯A是所有标记都有效应方差，但不同标记有不同的效应方差；贝叶斯B是部分标记有效应方差，同时具有差异；贝叶斯Cπ是部分标记效应有方差，同时所有方差都相等。

GS常用免费开源R包：

如果GWAS使用的是MLM模型及其进阶版，那么它和GS用的其实是同一模型。

只是GWAS的重点在挖掘显著性位点，关注固定效应Xb，加入随机效应只是为了控制与表型显著相关位点的假阳性；而GS重点在计算育种值，关注随机效应Zu，加入固定效应是为了控制不同个体相同的部分。

两种方法虽然模型一样，但应用方向不同，所以视为不同技术，实际上是一回事儿。我们可以将二者结合起来，比如先用GWAS找出显著位点，再加入GS的固定效应中，可能会使预测的准确性更高。或者是如果标记数目太多，我们不想过滤掉一些无关标记，就可以结合GWAS的结果来过滤。

在GS和GWAS的结合方面，也已经有了一些模型和算法，如SSGBLUP，GS + de novo GWAS，GS + historical GWAS等。

致谢：

感谢张志武老师提供的最新教程以及开发的工具，感谢知乎张帆的课程总结。教程PPT参考： ## Statistical Genomics: 2020

原先三个样本的HTSeq-count计数的数据可以在我的GitHub中找到，但是前面已经说过Jimmy失误让我们分析的人类就只有3个样本，另外一个样本需要从另一批数据获取（请注意batch effect)，所以不能保证每一组都有两个重复。

我一直坚信”你并不孤独“这几个字，遇到这种情况的人肯定不止我一个，于是我找到了几种解决方法

以上方法都会在后续进行介绍，但是我们DESeq2必须得要有重复的问题亟待解决，没办法我只能自己瞎编了。虽然是编，我们也要有模有样，不能直接复制一份，要考虑到高通量测序的read是默认符合泊松分布的。我是这样编的。

这仅仅是一种填坑的方法而已，更好模拟数据的方法需要参阅更加专业的文献， 有生之年 我希望能补上这一个部分。

这部分内容最先在 RNA-Seq Data Analysis 的8.5.3节看到，刚开始一点都不理解，但是学完生物统计之后，我认为这是理解所有差异基因表达分析R包的关键。

基本上，统计课都会介绍如何使用t检验 用来比较两个样本之间的差异，然后在样本比较多的时候使用方差分析 确定样本间是否有差异。当然前是样本来自于正态分布的群体，或者随机独立大量抽样。

对于基因芯片的差异表达分析而言，由于普遍认为其数据是服从正态分布，因此差异表达分析无非就是用t检验和或者方差分析应用到每一个基因上。高通量一次性找的基因多，于是就需要对多重试验进行矫正，控制假阳性。目前在基因芯片的分析用的最多的就是limma 。

但是，高通量测序(HTS)的read count普遍认为是服从泊松分布（当然有其他不同意见），不可能直接用正态分布的t检验 和方差分析 。当然我们可以简单粗暴的使用对于的非参数检验 的方法，但是统计力不够，结果的p值矫正之估计一个差异基因都找不到。老板花了一大笔钱，结果却说没有差异基因，是个负结果，于是好几千经费打了水漂，他肯定是不乐意的。因此，还是得要用参数检验的方法，于是就要说到方差分析和线性模型之间的关系了。

线性回归和方差分析是同一时期发展出的两套方法。在我本科阶段的田间统计学课程中就介绍用方差分析 （ANOVA）分析不同肥料处理后的产量差异，实验设计如下

这是最简单的单因素方差分析，每一个结果都可以看成 yij = ai + u + eij，其中u是总体均值，ai是每一个处理的差异，eij是随机误差。

注：方差分析(Analysis of Variance, ANAOVA)名字听起来好像是检验方差，但其实是为了判断样本之间的差异是否真实存在，为此需要证明不同处理内的方差显著性大于不同处理间的方差。

线性回归 一般是用于量化的预测变量来预测量化的响应变量。比如说体重与身高的关系建模：

当然线性回归也可用处理名义型或有序型因子（也就是离散变量）作为预测变量，如果要画图的话，就是下面这个情况。

如果我们需要通过一个实验找到不同处理后对照组和控制组的基因变化，那么基因表达可以简单写成， y = a + b · treament + e。和之前的 yij = ai + u + eij 相比，你会发现公式是如此的一致。这是因为线性模型和方差分析都是广义线性模型 (generalizing linear models, GLM)在正态分布的预测变量的特殊形式。而GLM本身只要采用合适的连接函数 是可以处理对任意类型的变量进行建模的。

目前认为read count之间的差异是符合负二项分布，也叫gamma-Possion分布。那么问题来了，如何用GLM或者LM分析两个处理件的差异呢？其实可以简单的用上图的拟合直线的斜率来解释，如果不同处理之间存在差异，那么这个拟合线的斜率必定不为零，也就是与X轴平行。但是这是一种便于理解的方式（虽然你也未必能理解），实际更加复杂，考虑因素更多。

注1 负二向分布有两个参数，均值(mean)和离散值（dispersion). 离散值描述方差偏离均值的程度。泊松分布可以认为是负二向分布的离散值为1，也就是均值等于方差（mean=variance）的情况。

注2 这部分涉及大量的统计学知识，不懂就用维基百科一个个查清楚。

聊完了线性模型和方差分析，下面的设计矩阵（design matrix）就很好理解了，其实就是用来告诉不同的差异分析函数应该如何对待变量。比如说我们要研究的KD和control之间变化，设计矩阵就是

那么比较矩阵（contrast matrix)就是告诉差异分析函数应该如何对哪个因素进行比较，这里就是比较不同处理下表达量的变化。

其实read count如何标准化的方法有很多，最常用的是FPKM和RPKM，虽然它们其实是错的-- FPKM/RPKM是错的。

我推荐阅读 Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data , 了解不同标准化方法之间的差异。

有一些方法是要求原始数据，有一些则要求经过某类标准化后的数据，记得区分。

关于DESeq2分析差异表达基因，其实在 position bias)标准化。edgeR的 calcNormFactors 函数使用TMM算法 对DGEList标准化

注大部分的mRNA-Seq数据分析用TMM标准化就行了，但是也有例外，比如说single-cell RNA-Seq(Lun, Bach, and Marioni 2016), 还有就是global differential expression，基因组一半以上的基因都是差异表达的，请尽力避免，(D. Wu et al. 2013)，不然就需要用到内参进行标准化了(Risso et al. 2014).

第四步：实验设计矩阵(Design matrix)，类似于DESeq2中的design参数。 edgeR的线性模型和差异表达分析需要定义一个实验设计矩阵。很直白的就能发现是1vs0

第五步：估计离散值（Dispersion）。前面已经提到负二项分布（negative binomial，NB)需要均值和离散值两个参数。edgeR对每个基因都估测一个经验贝叶斯稳健离散值（mpirical Bayes moderated dispersion），还有一个公共离散值（common dispersion，所有基因的经验贝叶斯稳健离散值的均值）以及一个趋势离散值

还可以进一步通过quasi-likelihood (QL)拟合NB模型，用于解释生物学和技术性导致的基因特异性变异 (Lund et al. 2012; Lun, Chen, and Smyth 2016).

注1 估计离散值这个步骤其实有许多 estimate*Disp 函数。当不存在实验设计矩阵(design matrix)的时候， estimateDisp 等价于 estimateCommonDisp 和 estimateTagwiseDisp 。而当给定实验设计矩阵(design matrix)时， estimateDisp 等价于 estimateGLMCommonDisp , estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp 。其中tag与gene同义。

注2 其实这里的第三，四，五步对应的就是DESeq2的 DESeq 包含的2步，标准化和离散值估测。

第六步: 差异表达检验（1）。这一步主要构建比较矩阵，类似于DESeq2中的 results 函数的 contrast 参数。

这里用的是 glmQLFTest 而不是 glmLRT 是因为前面用了glmQLTFit进行拟合，所以需要用QL F-test进行检验。如果前面用的是 glmFit ，那么对应的就是 glmLRT . 作者称QL F-test更加严格。多重试验矫正用的也是BH方法。

后续就是提取显著性差异的基因用作下游分析，做一些图看看

第六步：差异表达检验（2）。上面找到的显著性差异的基因，没有考虑效应值，也就是具体变化了多少倍。我们也可用找表达量变化比较大的基因，对应的函数是 glmTreat 。

经过上面两个方法的洗礼，基本上套路你也就知道了，我先简单小结一下，然后继续介绍limma包的 voom 。

Limma原先用于处理基因表达芯片数据，可是说是这个领域的老大 :sunglasses: 。如果你仔细看edgeR导入界面，你就会发现，edgeR有一部分功能依赖于limma包。Limma采用经验贝叶斯模型（ Empirical Bayesian model）让结果更稳健。

在处理RNA-Seq数据时，raw read count先被转成log2-counts-per-million (logCPM)，然后对mean-variance关系建模。建模有两种方法：

数据预处理 ： Limma使用edgeR的DGEList对象，并且过滤方法都是一致的，对应edgeR的第一步,第二步，第三步

差异表达分析 : 使用”limma-trend“

差异表达分析 : 使用”limma-voom“

如果分析基因芯片数据，必须好好读懂LIMMA包。