GSSG论文阅读
有源码kkdey/GSSG: Gene Set + S2G strategy annotations analyzed for disease architecture
简介
GSSG : Combine SNP-gene linking strategy with Gene Sets A pipeline to perform the following steps.
- Generate gene programs for analysis
- Link SNPs to genes in a gene program to generate a genome-wide annotation
- Perform disease heritability enrichment analysis using S-LDSC
这里我们将主要关注其第二篇文章的工作：将GWAS数据与scRNA结合以识别细胞种类和细胞过程(Identifying disease-critical cell types and cellular processes by integrating single-cell RNA sequencing and human genetics)。
作者开发了sc-linker，一个框架将单细胞数据、表型图谱(epigenomic maps)和GWAS数据集成用以识别潜在的细胞类型和细胞过程(cell process)。该框架再次捕捉了已知的生物学规律并突显重要的细胞-疾病相关性，包括重度抑郁症中的 GABA 能神经元、溃疡性结肠炎中依赖疾病状态的 M 细胞程序，以及多发性硬化症中疾病特异性的补体级联过程。在自身免疫疾病中，健康状态和疾病依赖的免疫细胞类型程序均表现出相关性，而只有疾病依赖的上皮细胞程序显著，这表明它们更多参与疾病的反应过程而非起始阶段。我们的框架提供了一种强有力的方法，用于识别遗传变异影响疾病的细胞类型及其细胞过程。
结论
Overview
作者构建了一个框架来将从单细胞中获得基因过程(gene programs)与疾病和复杂的特征相关联。首先，作者使用scRNA-seq构建了基因过程，该基因过程定义未一个具有连续值的基因集合，该基因集合特征化了: (1) 个体的细胞类型 (2) 疾病的依赖性(健康人与患者的同种细胞的比较) (3) 细胞过程(cellular processes)。基因过程中的连续值为0-1之间的数值，但并不表示概率。
然后，作者将这些基因过程中潜在的基因和调控它的SNP相关联，关联方法依赖于结合两种组织特异性的增强子基因关联策略：1. Roadmap增强子，2. Activity-by-Contact(ABC)模型。
最后，作者通过使用S-LDSC，在考虑基线LD模型提供的一系列编码区、保守曲、调控区及连锁不平衡相关注释的条件下，评估所得SNP注释对于疾病信息的贡献。
总之，作者的这些操作再次捕捉了已知的生物学规律并突显重要的细胞-疾病相关性。对于这个流程看的不是很明白，后续结合具体数学公式与代码分析。
评价标准(benchmark)
为了证明自己的模型优于别人的模型，作者自己做了一个benchmark。作者分析了5种学习的特征，这5个特征对应于特定的免疫细胞类型。作者构建6种免疫细胞类型程序（来自4个数据集：2个PBMC，1个骨髓，1个脐带血），作者从中识别出：红系细胞的富集和红细胞数量有关、巨核细胞和血小板数量有关，单核细胞和单核细胞有关，B、T细胞和淋巴细胞比例有关。这些富集是一种”预期中的富集“，因为这是符合既往的研究和生物学相关性的。
同时还定义了敏感性与特异性的index，用以量化现有富集方法和研发出的富集方法的表现效果。
Method
方法上都是常规思路，但是该论文中对于scRNA的处理策略可以作为后续进行scRNA处理流程的参考
scRNA-seq
正常处理，然后把基因表达矩阵均转化为$\log_{2}(\text{TP10K}+1)$
降维/去批次/聚类/标注
作者首先使用Scanpy’s highly_variable_genes函数识别出所有数据集中表达差异最大的2000个基因，然后使用PCA对着2000个基因提取40个主成分，使用Harmony包进行去批次操作，最后再用KNN对于前40个主成分进行聚类分析，聚类为10类。然后对于聚类结果采用Leiden图聚类方法进行社区检测(Community Dectection)，细粒度(resolution)为1。对于每个数据集，都采用UMAP进行可视化操作，如果之前的标注可以获得，那就可以作为参考对每个细胞进行标注，反之则获取每个细胞的独特表达印迹(signatures)与基因markers来进行标注。
Cell type gene programs
作者对于每个细胞类别进行一次非参数Wilcoxon值和检验比较该细胞类型与其他细胞来获得每个基因在该类细胞中的显著性。然后使用之前提出的策略将p值化为$X=-2\log(p)$，这个$X$服从卡方分布，然后再用最大最小归一化方法将X归一化到0-1之间。从健康人哪里获得的细胞类型基因程序被称为健康细胞类型程序，反之则成为疾病细胞类型程序。
Disease-dependent gene programs
作者进行健康组织与异常相同组织的非参数检验，然后采用相同的方法转化为0-1之间。作者同时还在饭呢西COVID-19数据集时构建了一个病毒进展程序（比较未感染和已感染的细胞类型）。作者观察到健康细胞类型基因程序与疾病依赖基因程序之间相关性弱。
Comment: 博主认为，这里的programs其实质上是采用一种特殊的方式将细胞进行一种独特的编码，显著性较高的部分标为1，而不高的部分标为0，这样就形成了一种专属于每种细胞的独特 one-hot encoding。这种独热编码描述了每个细胞的特征，且采用的两种方式的编码具有较弱的关联性，也体现出了：两种不同的编码所描述的细胞特征为不同维度，即正交的
Cellular process gene programs
使用从NMF（非负矩阵分解，non-negative matrix factorization）中获取的潜在因子。作者定义了一种细胞过程程序（cellular process program），该过程基于在不同细胞之间表达量和在每个因（Latent factors）中的贡献性呈现高相关性的基因。这个相关性然后被归一化到0-1（负相关被分配为0），然后作者对于每个因子采用这个因子中“驱动力最强”（most driving）的基因所富集出的通路进行标注。并根据该通路是只与一个细胞种类程序相关还是与多个细胞种类程序相关分为“内细胞类型（intra-cell type）”和”间细胞类型（inter-cell type）”因子。
- Comment: 这个NMF主要用于图像识别、文本处理，主要思路在于将一个大的矩阵分解为两个非负矩阵的乘积，该方法将一个大的矩阵中提取出重要信息。该思路也在深度学习中有所应用，LoRA作为大模型微调的低廉且快速的策略，其全名便是”Low-Rank Adaptation“。LoRA中，作者认为参数更新中存在一个“内在秩”，也就是说，我们更新的参数，可以采用两个低秩的向量相乘所得，于是我们仅仅需要训练这两个向量，而不需要训练巨大的权重矩阵。
我借用Chatgpt的解释：这两种方法都假设“变化在低维子空间里“。NMF：高维观测其实由少数 latent factors 线性组合而成；LoRA：微调任务只需要在权重空间里走一个低维方向。可将其类比为一种 ”基 \times 系数“的类比，也就是说我们将一个庞大的空间分为多个小的基向量，然后与系数相乘来描述原空间。这就保证了我们的预测/解释不会在整个参数空间内随意优化，即仅在rank-r子空间内学习这也就是为什么LoRA不易过拟合、NMF可解释性强。
再进一步考虑，这种思路其实是类似于谱分析思路，将高维度数据转移到一个低维度/低频率的空间进行考虑.
这里作者获得细胞过程基因程序采用以下优化函数进行：
$$
\begin{aligned}
\operatorname{argmin}\frac{1}{2}\left|X_{n, m}- \sum_p W_{n, p} \times H_{p, m} \right|_{F}^{2}+
\frac{1-\alpha}{2}\left|W_{n, p}\right|+ \frac{1-\alpha}{2}\left|H_{p, m}\right|+ \alpha\left|\operatorname{vec}(W_{n, p})\right|_{1}+ \alpha\left|\operatorname{vec}(H_{p, m})\right|_{1}
\end{aligned}
$$
上图中$X_{n,m}$代表着基因表达矩阵中的元素：基因m在样本n中的表达情况，$W_{n,p}$代表了样本n中的潜在因子p，$H_{p,m}$代表因子p在基因m中的权重。这里第一项表述采用了非矩阵的表示，NMF如果采用矩阵则可以直接如下表示：
$$
\operatorname{argmin}\left|X - WH\right|^2 _{F}
$$
其中 $X\in \mathbb{R}^{N\times M}, W\in \mathbb{R}^{N\times P},H\in \mathbb{R}^{P\times M}$ 最外层表示2范数，也被称为Frobenius范数(F范数)，也就是平方和的平方根。后续其还有多项，我们可以认为是某种正则项，保证优化得到的矩阵不会过拟合。
*
Comments: 这里的数学公式后面几项是错误的，同时argmin函数的书写不全，可能导致误解，因为W_{n,p}等应为标量(scalar)，而标量不能进行vec操作，取范数操作也不严谨。正确书写形式应该如下：
*
$$
\begin{aligned}
\min_{W, H} & \frac{1}{2}|X-W H|_{F}^{2} +(1-\alpha) \frac{1}{2}|W|_{F}^{2}+(1-\alpha) \frac{1}{2}|H|_{F}^{2} +\alpha|\operatorname{vec}(W)|_{1}+\alpha | \operatorname{vec}(H)|_{1}
\end{aligned}
$$
如上式子展开便可以发现，这个优化函数其实是类似于LASSO + 岭回归的正则化方法（对于系数分别进行1阶和2阶范数的惩罚，目前类似于这样的正则化方法被称为Elastic-Net）。
作者通过如下方法定义一个细胞过程基因程序：识别基因中 在每个细胞表达量和factors(也就是W中的一列)呈现强相关。潜在因子中，如果相关性大于0.8（潜在因子之间的，也就是说W中如果出现某2个列相关性大于0.8，那就只考虑其中之一），那么只考虑一个潜在因子。我们采用相关性最强的（这里的相关性是 表达水平和factors之间的相关性）且factors所对应的H中权重更高的基因所富集得出的通路来对于每个细胞过程程序进行标注。并将其标注为”intra-cell type”和”inter-cell type”，前者代表这个细胞过程基因程序只和单个细胞类型程序强相关，后者代表多个细胞类型程序都与之有关。
Cellular process gene programs constructed from healthy and disease tissues
作者这部分则是将目标放在提取每个数据集中健康人以及患疾病人中的共性与特性。它提出了一种NMF方法方法：
$$
H_{P\times N_{1}}=[L_{P\times K_{C}}^{CH}, L_{P\times K_{H}}^{UH}]F^H_{(K_{C}+K_{H})\times N_{1}}, L^{CH},L^{UH},F^H_{0}>0
$$
$$
H_{P\times N_{2}}=[L_{P\times K_{C}}^{CD}, L_{P\times K_{D}}^{UD}]F^D_{(K_{C}+K_{D})\times N_{2}}, L^{CD},L^{UD},F^D_{0}>0
$$
上文的公式是博主进行了一些“修正”得出的结果，原作者的公式有明显错误。$H_{P\times N_{1}}$代表健康个体中特定组织观察到的基因表达矩阵；$D_{P\times N_{2}}$为相同组织在患病个体中的基因表达矩阵；P代表基因的数目(非常坏的标注习惯，上面不都用过P了，为什么不能用N_Gene_H，这种无歧义的变量名)。这里原作者考虑是将H分解为：$L^{CH}$与$L^{UH}$其中C代表共性、U代表个性，H代表健康，D代表疾病，所以解释就是：$L^{CH}$代表健康人中与患病者的共性部分，$L^{UH}$代表健康人中个性的部分；$L^{CD}$代表患病者中与健康人的共性部分，$L^{UD}$代表患病者中的个性部分。所以作者希望共性部分能尽量接近，这也在后面的优化函数中有所体现。
$$
\begin{array}{l}\quad \underset{L^{H}, L^{D}, F^{H}, F^{D}}{\operatorname{argmin}} \frac{1}{2}\left|H-L^{H} F^{H}\right|_{F}^{2}+\frac{1}{2}\left|D-L^{D} F^{D}\right|_{F}^{2}+\frac{\mu}{2}\left(\left|L^{H}\right|_{F}^{2}+\left|L^{D}\right|_{F}^{2}\right) \ \ +\frac{\gamma}{2}\left(\left|{L}^{C H}-L^{C D}\right|_{F}^{2}\right)\end{array}
$$
其中，$L^H = [L^{CH}_{P\times K_{C}}, L^{UH}_{P\times K_{H}}]$，$L^D=[L^{CD}_{P\times K_{C}}, L^{UD}_{P\times K_{D}}]$，$\gamma$为可调参数，用以控制共性部分的相似性；$\mu$依旧控制正则化部分。
对于以上所有待优化参数进行求导得到倒数:
$$
\begin{array}{c}
\nabla _{L^{H}}=-H F^{H^{T}}+L^{H} F^{H} F^{H^{T}}+\mu L^{H}-\gamma\left[L^{CH} - L^{C D} ;0\right] \\
\nabla _{L^{D}}=-D F^{D^{T}}+L^{D} F^{D} F^{D^{T}}+\mu L^{D}-\gamma\left[L^{C D} -L^{C H}; 0\right] \\
\nabla_{F^{H}}=-L^{H^{T}} H+L^{H^{T}} L^{H} F^{H} \\
\nabla_{F^{D}}=-L^{D^{T}} D+L^{D^{T}} L^{D} F^{D}
\end{array}
$$
Comments: 作者后面使用了NMF自带的参数更新策略，我在验证后发现，最终的优化模型大致一致，但是模型中对于前项求导错误
作者应采用以下公式进行参数更新，而非原公式，直到100个周期或者误差小于一个特定值。
$$
\begin{align}
L_{i j}^{H} &\leftarrow L_{i j}^{H} \frac{\left(H F^{H^{T}}\right)_{i j}}{\left(L^{H} F^{H} F^{H^{T}}+\mu L^{H} + \gamma\left[L^{C H} - L^{C D}; 0\right]\right)_{i j}} \\
F_{i j}^{H} &\leftarrow F_{i j}^{H} \frac{\left(L^{H^{T}} H\right)_{i j}}{\left(L^{H^{T}} L^{H} F^{H}\right)_{i j}}
\end{align}
$$
Enhancer-gene linking stategies
作者使用一种增强子-基因连接策略，将0或多个基因分配给在1000 Genomes Project 欧洲人群参考中次等位基因计数大于5的每一个SNP。这里，我们主要考虑了Roadmap + ABC的策略，这两种方法一方面广泛可用，另一方面是优于其他方法。对于与相关组织对应的任何生物样本（细胞类型或组织），我们均采用组织特异性的 Roadmap 和 ABC 增强子–基因链接策略来构建相应的基因程序（gene programs）。基于对于免疫细胞的类型的分析，在单细胞测序中的87%的基因可以被观察到有增强子-基因连接。同时，我们也考虑了非组织特异的Roadmap和ABC策略，除了这个，作者还考虑标准的100kb窗宽的策略
这里的增强子-基因连接策略用于分析SNP与特定基因的相关性，是连接GWAS与scRNA的重要桥梁。
Genomic annotations and the baseline-LD models
作者将会给每个SNP赋予一个特定的值，这个过程便被称为标注，这个标注如果对于了某种功能，则被称为功能标注。基线-LD模型包括了86种标注。
Stratified LD score regression
分层连锁不平衡分数回归（S-LDSC）评估了一个基因标注对于疾病或者其他表征的贡献性。
作者使用了两个度量来评估标注c的信息量，分别为E-score和标准化效应量$\tau^*$
MAGMA 基因富集
MAGMA评估了疾病中基因富集的情况
评估敏感性与特异性
作者定义了一个敏感性与特异性的benchmark。1. sc-linker与MAGMA基因富集分析的比较；2. 不同版本中采用不同方法来定义细胞类型程序和SNP-to-gene连接策略的sc-linkers之间的比较。
为了和MAGMA比较，作者比较了1. 对于期待通路的sc-linker中的富集评分和MAGMA中的相关性之间的差值(敏感性，也就是该找的都该找到)；2. 对于其他通路在基因集等级上的sc-linker中的富集评分和MAGMA中的相关性之间的差值（特异性，也就是不该找的有没有没找到）。作者承认有未被观测到的未知富集，所以作者也会考虑在寻找目标富集的敏感性作为代替。
Identifying genes driving heritability enrichment
对于每个基因程序，我们首先从完整的基因列表中进行筛选，仅考虑那些在该基因程序中具有超过 80% 成员概率等级（probability grade of membership） 的基因。接着用MAGMA对于留下的基因进行排序，并考虑最高的50个基因来进行下一步的下游分析。而baseline则是取出前200个基因进行针对细胞类型的疾病通路富集。
