Science|中国人的千万级免疫细胞图谱|亚群|位点|免疫性疾病|免疫细胞|染色质|特异性

摘要：超大超全的PBMC的免疫细胞图谱！整合scRNA-seq、scATAC-seq和WGS测序数据，结合脂质组学、蛋白组学和血液生化检测，共同解析遗传变异对基因表达和表型的影响。

总而言之，不用翻来覆去的找数据了！！！

研究背景

人类免疫系统是一个复杂的细胞网络。不同个体之间免疫系统的差异源于遗传和非遗传因素（年龄和生活习惯等）的共同作用，其中遗传因素约占个体间差异的20%至40%，而这些差异和个体的生理状态的差异息息相关。因此解析这种差异对于免疫相关疾病的机制的阐明是十分重要的。

之前对于这种差异的的遗传因素的研究主要基于xQTL研究（影响某个细胞类型的基因表达的DNA位点的研究）、细胞分选和批量RNA测序（RNA-seq），但是我们知道RNA-seq不能区分精细的细胞亚群和状态，所以不能解析具体的某种遗传变异是否影响某个细胞亚型的变化。

因此，整合单细胞转录组（scRNA-seq）和单细胞基因组（scATAC-seq和WGS）对于联系

遗传变异

细胞类型和特定背景下的基因调控

，揭示这些变异在健康和疾病状态下以及在不同人群中的作用机制是十分必要的，这也是此工作的研究重点。另外，作者开发了一个细胞语言模型（CIMA-CLM），该模型整合了染色质序列特征和单细胞转录组数据，以预测染色质可及性并评估非编码变异的功能效应。

研究结果

1 一份涵盖428名个体的综合的中国免疫多组学图谱

作者对来自428名中国人的外周血进行了scRNA-seq、scATAC-seq、WGS、脂质组学、代谢组学和生理生化的检测（图1A）。这些人的年龄在20-77之间，包括189名男性和239名女性（图1B）。PBMCs经过质控和去双细胞后，保留了来自scRNA-seq的6,484,974个细胞和来自scATAC-seq的3,762,242个细胞。

通过逐层注释（先注释成6类细胞-L1，然后分别筛选各类细胞进一步细分注释-L2，再分别筛选细分——L3，L4），最终界定了73种免疫细胞类型（L4，第四层）。

作者使用scGlue深度学习模型，基于基因组邻近性构建先验网络，实现了scRNA-seq与scATAC-seq数据的整合，使得scRNA-seq与scATAC-seq的注释一致（图1C，D）。并且发现了和年龄/性别相关的亚群（图1E：最外层是和年龄的相关性；向内一层是和性别的相关性。*代表显著性，颜色代表和年龄的政府相关性以及性别的倾向性）。

2 血液中的生物分子和免疫特征表现出性别和年龄相关的差异为研究个体间差异，作者对5类数据（脂质组学、代谢组学和血液生化取其指标，scRNA-seq和scATAC-seq取细胞比例）分别进行了PCA降维（图2A）。细胞比例前两个主成分（PC1和PC2）能够按年龄区分个体，而代谢和生化分子则按性别分离（图2B），也就是说细胞比例和年龄相关，而代谢和生化指标和性别相关。随后使用MOFA以428名的5类数据（包括19个血液生化特征、1228个脂质物种、321种代谢物，以及来自scRNA-seq的73种细胞类型和来自scATAC-seq的65种细胞类型的比例）作为输入进行分析（高级PCA），发现一种和性别相关的因子（什么是因子？可以简单理解为MOFA算出来的一个和每个个体的所有指标相关的一个总指标），这个因子的特征是男性中的2，3-二磷酸甘油酸水平升高，CD8 CTL-GZMB比例升高（图2C，D）。调整性别影响后，发现因子2和年龄相关，例如ncMono-FCGR3A和CD4 CTL-GZMH的比例随年龄增加（图2E）（？有个问题，文章说因子2也和脂质组相关，但是没看到图，而且figs4特别糊，看不清）。作者也使用scRNA-seq的基因表达数据和scATAC-seq的染色质可及性数据进行MOFA分析，发现因子7和年龄高度相关（图2F）。其中对于DC2-CD1C这个细胞亚群，因子7的方差解释率较高，且这个亚群的marker中CXC3R1&CCR2随年龄表达上调，SERTAD1&IER3表达下调。为解释基因表达变化是由于细胞类型还是个体因素，作者使用variancePartition工具进行了变化的贡献度计算，发现个体原因构成的变化主要集中于核心生理过程（性别相关，氧气转运/二氧化碳转运），细胞类型构成的变化主要是免疫相关功能（图2I，J）。年龄相关基因PER1和PPP1R15A和个体因素更相关，CXC3R1和细胞类型更相关，主要在NK细胞表达（图2K）。

3 构建免疫细胞类型特异性的基因调控网络（GRNs）

作者使用MACS3对scATAC-seq进行峰值检测，识别到338,036个染色质开放区域（cCRE），仅6.9%的cCRE位于启动子区，而大多数位于内含子区（31.2%）和基因间区（20.7%），或与转座元件重叠（32.9%）（图3A）。CIMA的cCRE有70%以上和已发表的scATAC-seq的cCRE重叠，40%以上和ENCODE数据库重叠。重叠的cCRE在10种以上细胞类型都开放，非重叠的部分有50%以上在单一细胞类型开放（图3B）。通过整合scRNA-seq（基因表达）和scATAC-seq（染色质开放）数据，基于基因组邻近性表达相关性，将增强子区域与靶基因连接起来，并推断主导该调控的转录因子。作者一共识别到404个GRNs，包括84,625个调控区和13,645个目标基因。特定转录因子（TF）调控的基因数量与其相关增强子区域的数量呈正相关（图3C）。在所有GRNs中，97%的基因与1至15个预测增强子相关，95%由1至20个TF调控（图3D）。 GRNs为免疫细胞类型的细胞身份提供了重要见解。基于eRegulon富集评分的降维有效区分了不同免疫细胞类型。作者选取了237种eRegulons，并鉴定了针对不同免疫细胞类型的的特异转录因子（图3E），确认了参与谱系分化的转录因子和罕见细胞类型种活性较强的转录因子。另外作者探讨了单核细胞分化过程中转录因子活性的变化，在从经典单核细胞向非经典单核细胞过渡过程中，作者观察到TF活性的动态变化，包括MITF、FOS、MBD2和CUX1（图3F，G，H）。为评估个体间差异，作者评估了GRN与年龄之间的相关性。热图展示了CD8 T cells 和unconventional T cells（图3I）以及单核细胞中的转录因子和年龄的相关性（图3J），比如CD8 T cells 和unconventional T cells中EOMES、RUNX3和TBX21等TFs与年龄呈显著正相关。在CD8 CTL-GZMK中，NFATC3表现出高调控活性，靶向PRF1和SLFN5等基因（图3K）；在intMono-GFRA2中，ZBTB7A调控了C3AR1和TLR1等基因（图3L）。在不同性别中有比例差异的CD4 T细胞中，女性的HIF1A、NFKB1和STAT5B活性更高，男性的BHLHE40更高。在ncMono-FCGR3A中，女性的STAT1活性升高（图3M）。

4 免疫细胞类型层面的eQTLs（影响基因表达的遗传变异位点）和caQTLs（影响染色质区域开放程度的遗传变异位点）的鉴定为探究遗传变异如何影响染色质可及性和基因表达，作者对部分样本进行了WGS测序，获得505万常见的遗传变异。为了绘制顺式表达定量性状位点（cis-eQTLs）和顺式染色质可及性质（cis-caQTLs），将每个个体中同一细胞类型的单细胞表达/染色质数据进行汇总，为每个细胞类型生成个体水平的表达谱和开放染色质谱。最终在69种细胞类型中寻找cis-eQTL，在42种细胞类型中寻找cis-caQTL，分别识别出9600个eGene（与主变异显著相关的受调控的基因）和52,361个caPeak（与主变异显著相关的染色质可及性峰）（图4A，B）。且检测到的eGene和caPeak的数量和细胞类型及细胞类型的占比呈正相关（图4C，D）研究还关注到仅在单一细胞类型识别到的eGene/caPeak和在一种细胞类型显著被影响但在另一种细胞类型中不被影响的eGene/caPeak（成对特异性）。通过计算π1（共享信号比例）rb（效应值相关性）统计量，量化了xQTL在不同细胞类型间的共享程度。结果表明，eQTL和caQTL在细胞类型间存在中等至高度的共享，但共享并不意味着调控变异相同（图4E，F）

CIMA队列的遗传变异以东亚人群常见变异为主。与全球人群（ALFA数据库）比较发现，约10.4%的lead xQTLs在欧洲或非洲人群中为罕见变异（MAF < 0.01），提示存在人群特异的调控效应。

研究不仅限于基因附近的顺式调控，还探索了遗传变异对远处染色体上基因的调控。鉴定出84个反式eGenes，其中69%具有细胞类型特异性，且顺式基因和反式基因数量成正比（图4G，H）。例如在初始CD4 T细胞中，rs11886530变异同时调控其邻近的生物钟基因NPAS2（顺式效应）和位于另一条染色体上的另一个核心生物钟基因NR1D1（反式效应）（图4I，J），且这个变异在中国人和非洲人群中比较常见（图4K）。

5 细胞类型特异性和共享的多效性关联揭示了炎症和疾病易感性的调控机制为探究遗传调控与免疫特征之间的潜在联系，作者利用xQTL结果，通过基于汇总数据的孟德尔随机化（SMR）方法，系统探索了遗传变异、eGene/caPeak与免疫相关性状（血浆脂质/代谢物、循环炎症蛋白、疾病风险等）之间的多效性关联。研究共评估了154个性状，在68种免疫细胞类型中识别出1196个显著的SMR关联，涉及833个变异、138个eGene、280个caPeak和68个性状（图5A, B）。其中73.2%的关联仅在单一细胞类型中显著，凸显了细胞类型特异性遗传调控的普遍性（图5B）。典型案例包括：rs34415530变异通过特异性下调CD4 Treg-FOXP3细胞中的IKZF4基因表达，进而影响循环IL-12B蛋白水平并增加哮喘风险，阐明了从遗传变异到特定细胞功能失调再到疾病的完整通路（图5C-G）。此外，rs2235922变异通过特异性上调经典单核细胞（cMono-CD14）中的PADI2表达来增加类风湿关节炎（RA）风险（图5H-K）。研究也发现部分关联在多个细胞类型共享（如rs10946216在B细胞和树突状细胞中调控RA风险基因CCR6）（图5L, M），并揭示了人群特异性（如rs312457与T2D风险）和中国南北方纬度选择相关的调控信号。这些整合分析揭示了遗传变异如何通过细胞类型特异性的分子机制影响免疫反应和疾病易感性。

6CIMA-CLM预测染色质可及性并评估单细胞中非编码变异的影响作者开发了CIMA-CLM，这是一个细胞语言模型。它整合了我们scATAC-seq和scRNA-seq数据集中的501-bp染色质序列单细胞基因表达数据，用以预测染色质可及性并评估非编码变异效应。

该模型采用双分支编码器架构：一支通过预训练的HyenaDNA将501-bp染色质序列编码为DNA嵌入，另一支通过预训练的scGPT将单细胞基因表达数据编码为细胞嵌入；再通过基于Transformer的编码器优化和多头交叉注意力融合解码器，整合序列与表达信息以预测染色质可及性（图6A）。为评估模型，作者基于20名随机个体的数据构建了32种细胞类型的独立测试集，其余数据用于训练与验证。

评估显示，CIMA-CLM在所有32种细胞类型中均能准确预测染色质可及性，预测峰与实验观测的scATAC-seq图谱高度一致（图6B）。皮尔逊相关系数分析显示中位PCC值介于0.7661至0.9612之间（整体均值=0.8951），且PCC值随单细胞捕获深度下降而降低，表明准确预测需要足够数据量（图6C）。受试者工作特征曲线下面积分析进一步证实模型高性能，AUROC值介于0.9058至0.9927（整体均值=0.9560）（图6D）。与现有模型相比，CIMA-CLM在预测染色质可及性信号方面表现更优。模型成功捕捉到跨细胞类型的共性信号及细胞类型特异性模式，例如B细胞、T细胞、髓系细胞及细胞毒性细胞分别在特征基因位点显示特异峰信号，且对FADS2、TMEM258等多个基因的预测准确性高（图6E）。

利用模型的单核苷酸嵌入能力，作者对开放染色质区域进行计算机模拟诱变，以预测非编码变异对染色质可及性的调控效应。通过比较突变与参考序列的预测峰，评估可及性变化（图6F，G）。大多数侧翼突变影响微小，而靠近QTL位点的替换常显著改变可及性。例如，模型预测Bn-TCL1A细胞中VASH1基因内的非编码变异rs3783996会显著降低染色质可及性（图6F），与SMR结果一致；预测类风湿关节炎风险变异rs2069235的核苷酸替换会增强SYNGR1基因附近区域的染色质可及性（图6G），且该预测与caQTL及GWAS结果相符。这些发现彰显了CIMA-CLM在阐明非编码变异功能及解析疾病位点机制方面的潜力。