iMetaOmics | 魏来/贾慧珏/何明光-多组学助力揭示塑造转录组的基因型-微生物组相互作用

多组学整合揭示塑造结膜转录组的基因型-微生物组相互作用

●  宿主基因型和微生物组对宿主分子表型（例如基因表达）的影响可能存在相互依赖关系，但这些交互效应在很大程度上仍然未知；

●  我们通过分析来自120对双生子的多组学数据，探究了基因型和微生物组与结膜转录组的关联；

●  将基因型-微生物组交互项纳入基因表达模型，可提高大量基因的表达水平预测准确性，其中包括许多编码细胞粘附分子的基因;

●  我们开发了R语言工具包MicroGenix，用于在多组学数据中筛选可能参与塑造分子表型的基因型-微生物组相互作用。

引  言

基因型-环境交互作用（genotype-environment interaction, GEI）是指基因型和环境对表型的影响因彼此而异的现象。GEI被认为是表型变异（包括对复杂疾病的易感性）的关键机制。微生物组可以被视为宿主的环境暴露，并可能与宿主基因型相互作用以塑造宿主表型。为简洁起见，我们将这种交互效应称为基因型-微生物组交互作用（genotype-microbiome interaction, GMI）。基因表达是一种分子表型，可能由宿主基因型和微生物组共同调节。例如，一项基于原代人结肠上皮细胞和肠道微生物组提取物共培养的研究提供了宿主基因型和微生物组在调节宿主转录应答中相互作用的证据。此外，微生物组应答表达数量性状位点（microbiome response expression quantitative trait locus）的概念已被提出，它是指宿主基因表达对特定微生物组暴露的响应的eQTL。然而，宿主基因表达调控中具体的GMI在很大程度上仍然未知。

研究GMI的实验方法仅限于针对少数基因和微生物，主要基于体外或动物模型研究，而多组学数据整合分析在生成GMI假设方面具有较大潜力，并且在人类研究中是可行的，特别是对于基因表达等分子表型。在GMI研究中，考虑分子表型至关重要，因为它们在复杂性状的出现中起着基础作用，并可以帮助阐明潜在的分子机制。现有的计算方法评估基因组-转录组或微生物组-转录组的关联，而没有明确建模两者加性效应或GMI。在多元回归中，交互项捕捉两个变量之间的乘法效应，使得可以建模一个变量对结果的影响可能取决于另一个变量的值。乘法交互项已被用于在模型中表示基因型和环境在表型塑造中的交互作用。然而，这种方法尚未被用于在多组学数据中研究GMI对分子表型的影响。

在这项研究中，我们对广州双生子研究中的120对双生子的基因组、结膜转录组和结膜微生物组数据进行了整合分析。我们首先确定了宿主基因型与结膜转录组之间的关联，以及结膜微生物组与结膜转录组之间的关联。我们进一步开发了MicroGenix方法，并将其应用于结膜数据集，以筛选潜在的GMI。MicroGenix方法将基因表达建模为宿主基因型、微生物组及其交互作用的函数。这种方法还允许识别与疾病等复杂性状相关的GMI。

结  果

宿主遗传变异与结膜转录组的关联

迄今为止，大多数粘膜组织（包括结膜）的表达数量性状位点（eQTL）尚未确定。我们通过进行全基因组eQTL定位来表征遗传因素对结膜转录组的影响。分析仅纳入每对双生子中的一个个体，并平均左右眼的基因表达水平，因此分析数据集由不相关个体组成。为便于与GTEx数据集比较，我们采用与GTEx项目相同的分析方法和显著性阈值。我们鉴定了属于5946个基因的187,241个顺式eQTL（FDR < 0.05，表S1）。对于每个eGene（即具有eQTL的基因），检测到的eVariant（相应的显著变异）的数量的中位数为五个，73.4%的eGene有一个以上的eQTL（图1A）。随后，我们将结膜eQTL与GTEx项目中其他组织的eQTL进行比较。共695个eGene为结膜独有，例如IL5和FCRL4（图1B；表S1）。甲状腺与结膜共有的eGene数量最多，其次是日晒皮肤、胫神经、非日晒皮肤和食管粘膜（图1C）。对于结膜与其他组织共有的eVariant-eGene对，不同组织与结膜中作用方向相同的eVariant的百分比在90.77%到98.66%之间（图S1；表S1）。我们鉴定了编码细胞因子（例如IL5、TNF和VEGFC）和趋化因子（例如CCL22、CCL28、CXCL14、CXCL17和XCL1）及其受体（例如IL6R、IL12RB2、IL18R1、IL20RB、TNFRSF1A、CCR1和CCR2）的基因的顺式eQTL。这些结果还为研究结膜免疫性疾病的遗传基础提供了资源。

图1. 结膜eQTL定位分析

（A）具有一个或多个顺式eQTL的结膜基因数量。（B）维恩图表示结膜eGene（具有eQTL的基因）与GTEx项目所包括的所有组织的eGene的交集。（C）与GTEx项目中每个组织共有的eGene和eVariant-eGene对（eGene和相应的显著变异）的数量。

微生物组与结膜转录组的关联

图2. 结膜基因表达与微生物丰度在双生子内的相关性

（A）所有同卵双生子基因表达水平的双生子内倍数变化（FC）的分布。条形图的颜色表示log₂转换的FC范围。图中标注了log₂转换的FC > 2或 < 1的基因数量。（B）同卵双生子中表达水平平均log₂转换FC最高的前50个基因。数据表示为平均值±标准误差。误差条表示标准误差。（C, D）同卵双生子内微生物Bray-Curtis距离与LCN1（C）和LACRT（D）表达的log₂转换FC之间的Spearman秩相关性。（E）表达受环境影响的基因（中位数和平均FC > 2）的双生子内FC与微生物丰度的双生子内FC之间的Spearman秩相关性。标有星号的相关性满足Spearman ρ > 0.2和FDR < 0.1。（F）曼哈顿图展示了结膜基因表达与微生物丰度在双生子内的全转录组关联。线性混合模型的p值使用Satterthwaite方法计算。

将GMI纳入模型可提高基因表达预测的准确性

为了鉴定可能参与塑造结膜转录组的GMI，我们将宿主基因表达建模为基因型和微生物组的主要效应及其交互效应的组合。对于主要效应，我们包括了位于基因起始和结束位置± 1 Mb的单核苷酸多态性（SNP）的剂量，以及所有微生物物种的丰度。应用了两阶段拟合策略来确定要包括作为预测因子的SNP-微生物交互项。具体来说，首先拟合一个没有交互作用的模型，以选择重要的SNP和微生物预测因子，定义为在正则化模型中有非零系数或在标准线性回归模型中有名义显著性。然后通过将每个选定的SNP的剂量和微生物的丰度相乘来构建交互项。使用所有主要效应预测因子（不进行重要性筛选）加上交互项来拟合包含交互的模型。预测因子的高维性和多重共线性对标准线性回归构成挑战。对于包含和不包含交互的两种模型，通过决定系数R²进行评估（配对t检验，p < 2.2×10^-16；图S3），弹性网络在预测准确性方面优于多元线性回归和线性混合模型，线性混合模型与多元线性回归相比改进不显著（配对t检验，p > 0.05）。因此，我们在后续分析中使用了弹性网络。

在分析中包括的所有蛋白质编码基因（n = 14,197）中，18.9%的基因（n = 2681）在没有交互的模型中至少有一个SNP和一个微生物物种具有非零系数。在这些基因中，92.3%（n = 2475）至少有一个交互项具有非零系数。在包含交互的模型中预测良好的基因包括编码HLA类II分子的基因，如HLA-DQB1、HLA-DQA2和HLA-DOA（R² > 0.2，图S4）。值得注意的是，包含交互的模型的预测准确性高于不包含交互的模型（配对t检验，p < 2.2×10^-16；图3A），并且2158个基因在包含交互的模型中显示出准确性提高（表S4）。准确性提高的基因在细胞-基质粘附这一生物过程富集（图3B）。在细胞粘附分子中，ANTXR1（ANTXR细胞粘附分子1）的准确性提升最多，R²从0.012提升到0.194（图3C、D；表S4）。准确性提高的细胞粘附分子编码基因（表S4）的SNP-微生物交互形成了一个网络，包括217个顺式eQTL和35个与基因表达相关的微生物物种（后者在双生子对内分析中所识别）。例如，ANTXR1的表达水平与rs6760455和Corynebacterium ureicelerivorans之间的交互作用相关（图3E）。我们检查了2158个基因的最重要的交互项（定义为系数的绝对值最大的项）涉及的物种（图3F；表S4）。Porphyromonas gingivalis参与最多基因的最重要交互项，其中LRFN2是R²提升最多的基因。

图3. 基于基因型、微生物组及其相互作用的基因表达预测

（A）比较包含和不包含SNP-微生物交互项的弹性网络模型的预测准确性。散点图显示了两种模型类型的决定系数（R²），每个点代表一个基因。呈现了每次交叉验证中的平均R²。p值是使用配对t检验计算的。密度图说明了基因准确性提高的分布，定义为包含交互的模型中的R²减去不包含交互的模型中的R²。（B）显示在包含交互的模型中预测准确性提高的基因的富集通路。点的大小表示基因数量，点的颜色表示通路中预测准确性提高的基因比例。（C, D）使用包含和不包含交互的模型预测的ANTXR1表达水平与观测值的对比。点的颜色表示交叉验证中每个折叠的测试集。图中标注多个折叠的平均R²。线条表示使用线性模型的最佳拟合，阴影区域表示95%置信区间。（E）网络图展示编码细胞粘附分子的基因（在包含交互的模型中准确性提高）与模型中交互项所涉及的微生物。（F）参与构成最多基因的最重要交互项（弹性网络中系数绝对值最大的交互项）的前30种微生物物种。

识别与复杂性状相关的GMI

为了便于使用我们的多组学数据筛选GMI的方法，我们开发了一个名为MicroGenix的R包，它由三个模块组成：使用基因型、微生物组及其交互作用对参考转录组进行训练；使用基因型和微生物组数据预测基因表达水平；以及将预测的表达与复杂性状相关联（图4A）。预测的表达水平可以被视为由基因型、微生物组及其交互作用解释的转录变异的估计。将解释的基因表达组分与复杂性状相关联有助于优先考虑GMI，这些GMI在关联基因的建模中发挥作用。基于基因的关联测试使用逻辑和线性回归分析定性和定量性状，使用单个基因的预测表达水平作为变量。

我们通过识别可能在过敏性结膜炎中起作用的GMI来举例说明MicroGenix的应用。为此，我们预测了一个包含健康个体和过敏性结膜炎或干眼病患者样本的宏基因组数据集中的结膜基因表达水平。我们使用根据纳入分析的总基因数量的Bonferroni校正来确定显著关联阈值（p < 2.4×10^-5）。我们鉴定了148个与过敏性结膜炎相关的基因，其中39个基因显示出疾病特异性，即不与干眼病相关联。此外，我们使用不包含交互的模型、仅用SNP或微生物预测因子拟合的模型进行了相同的分析。39个基因中有22个是仅通过包含交互的模型鉴定的，并且所有这些基因在包含交互的模型中的R²都高于其他三种模型（表S5）。值得注意的是，22个基因中有11个与过敏性结膜炎呈正相关（图4B），其中5个位于先前已与过敏相关联的遗传区域，包括PTTG1IP、SCIN、TMEM18、CFLAR和LAMC2，其中LAMC2是一种细胞粘附分子。其他参与免疫反应的基因包括CD276和CXCR2。CD276在过敏性哮喘儿童中升高，并调节Th2/Th17细胞分化。基因表达预测模型将这些基因与涉及24种微生物物种和31个SNP的SNP-微生物交互作用联系起来（图4C；表S6）。这些发现为进一步研究眼表疾病中的GMI提供了基础。

图4. MicroGenix用于筛选可能与疾病发病机制有关的SNP-微生物交互作用

（A）MicroGenix方法包括使用参考数据集训练基因表达预测模型，这些数据集包括宿主基因型、微生物组和转录组数据。这一过程生成了一个包含预测模型和SNP-微生物交互的数据库。这些模型将被应用在包括宿主基因型和微生物组数据的新数据集上，以预测基因表达。随后进行关联测试，将预测的基因表达水平与特定目标性状相关联。这有助于优先考虑基因，并在数据库中识别关键的SNP-微生物交互作用。G×M表示基因型（G）和微生物组（M）之间的交互作用，捕捉两者之间的乘法效应。宿主基因表达（Y）被建模为截距（β₀）、基因型（β₁）和微生物组（β₂）的主要效应以及交互效应（β₃）的线性组合。（B）与过敏性结膜炎正向且特异性相关的11个基因，在健康个体和过敏性结膜炎患者中的预测基因表达水平。（C）过敏性结膜炎相关基因的弹性网络模型中，拥有非零系数的SNP-微生物交互项的系数绝对值。对于具有多个交互项的基因-微生物对，系数绝对值最大的SNP在图中展示。

讨  论

结  论

方  法

多组学数据

结膜多组学数据来自广州双生子研究（GTES）。从2006年起，居住在广州的双生子被招募参加GTES。所有参与者在入组时提供了书面知情同意。研究遵循赫尔辛基宣言的人体研究原则，并获得了中山大学中山眼科中心伦理委员会的批准。本研究纳入110对同卵双生子和10对异卵双生子的数据，包括117名男性和123名女性，平均年龄为18.9岁（18.9 ± 2.9，平均值±标准差）。本研究分析了480个结膜样本（来自240名个体的双眼）和240个血液样本。所有双生子都是健康受试者，没有系统性疾病或眼病的证据。至少有一名个体在过去6个月内有抗生素治疗或隐形眼镜佩戴史的双生子被排除在外。230名健康受试者、89名过敏性结膜炎患者和166名干眼病患者的宏基因组数据在我们的既往研究中收集。

样本收集和处理

按照我们之前建立的方法收集结膜印迹细胞学样本。简而言之，在采样前用1-2滴Alcaine眼药水（Alcon, Fort Worth, TX, USA）对眼睛进行局部麻醉。将一个无菌的半圆形膜滤器（REF: HAWP01300, 0.45 µm; Merck Millipore, Burlington, MA, USA）放置在下球结膜上10秒钟。然后将膜立即放入含有300 µl组织和细胞裂解液（Epicentre, Ambleside, UK）的无菌管中，并储存在-80°C。未用于采样并暴露于与用于采样的膜相同环境的膜被收集作为阴性空白对照。使用MasterPure Complete DNA和RNA纯化试剂盒（Epicentre, Ambleside, UK）按照制造商的说明从结膜样本中提取RNA。阴性空白对照与结膜样本一起使用相同的程序进行处理，包括储存、RNA提取、文库制备和测序。血液样本被收集并储存在-80°C。使用TIANamp Blood DNA试剂盒（DP348, TianGen Biotech, Beijing, China）按照制造商的说明提取DNA，并进行SNP分型。

测序

使用VAHTS Total RNA-seq文库制备试剂盒（Vazyme, Nanjing, China）按照制造商的说明为RNA-seq准备文库。这个试剂盒可以耗尽核糖体RNA和DNA，留下mRNA和其他非编码RNA。使用KAPA SYBR FAST qPCR试剂盒（Kapa Biosystems, Wilmington, MA, USA）通过qPCR定量文库。在NovaSeq 6000仪器（Illumina, San Diego, CA, USA）上进行2×150-bp的双端测序。

SNP阵列数据的插补

使用Illumina Infinium Asian Screening Array进行基因分型。我们使用PLINK在SNP和样本水平上对基因型数据进行质量控制。在SNP水平上，我们排除了次要等位基因频率（MAF）< 1%、基因型缺失率> 5%和偏离哈代-温伯格平衡（HWE）的SNP（p < 10^-5）。在样本水平上，我们移除了基因型缺失率> 5%和杂合度值超出平均值±3标准差的样本。我们进行了同系身份分析，并在存在一级或二级关系的情况下保留了基因型缺失率较低的样本。我们进行了主成分分析（PCA），并在前五个主成分（PC）中偏离平均值±6标准差的样本被排除。使用SHAPEIT4和IMPUTE5以及1000 Genomes Project数据集（第三阶段）的参考面板进行质量控制的基因型数据的阶段和插补。排除了MAF < 1%或INFO值< 0.8的SNP。插补后概率< 0.9的基因型被视为缺失，基因型缺失率> 5%的SNP被移除。

宏基因组测序数据的插补

使用内部管道将宏基因组数据作为浅全基因组测序数据进行基因分型。简而言之，质量控制后的读取与人类参考基因组（GRCh37）使用BWA进行比对。消除PCR重复。使用GLIMPSE进行基因型插补。用于插补的单倍型参考面板来自1000 Genomes Project数据集。排除了MAF < 1%和INFO分数< 0.8的SNP。基因型概率< 0.9的位点被视为缺失。

人类基因表达分析

RNA-seq数据的比对和定量是遵循GTEx项目的RNA-seq管道进行的。简而言之，RNA-seq数据与人类参考基因组GRCh37使用STAR进行比对。使用Picard套件标记重复读取。使用RNA-SeQC进行基因水平表达定量。定量基于GENCODE Release 41注释。

微生物组成分析

在质量过滤后，使用KneadData去除来自人类基因组和核糖体RNA的读取。使用Kraken 2将读取映射到自定义数据库。总共从国家生物技术信息中心（NCBI）下载了29,943个完整的微生物基因组，其中19,362个是细菌，368个是古菌，9346个是病毒，867个是真菌。使用置信度分数0.20过滤分类结果。只有在至少一个样本中读取数超过10个且在> 10%的样本中存在的物种才被保留。使用decontam R包的isContaminant函数的基于频率和基于流行度的方法进行污染物识别。在空白样本中检测到的物种被进一步移除。

eQTL定位

我们遵循GTEx项目的流程进行顺式eQTL定位。读取计数使用修整均值的M值（TMM）进行标准化。将在≥ 20%的样本中表达量≥ 0.1 TPM且在≥20%的样本中≥ 6个读取的基因纳入分析。每个基因在样本中进行逆正态转换。我们使用在peer R包中实现的PEER方法估计基因表达水平的潜在协变量。在顺式eQTL定位中，每个家庭只包括双生子的其中一个个体。对于每对双生子，选择基因分型率较高的双生子。个体的左眼和右眼的基因表达水平被平均。使用修改版的FastQTL进行顺式eQTL定位，参数与GTEx项目相同：定位窗口设置在转录起始位点（TSS）的±1 Mb，应用FDR阈值≤ 0.05。FDR由FastQTL内部使用置换方法估计。我们包括了前15个PEER因子、前五个基因型PC、年龄和性别作为协变量。GTEx项目的eQTL数据从GTEx门户网站下载（https://www.gtexportal.org）。

转录组和微生物组的双生子对内相关性

我们计算了同卵双生子内基因表达和微生物丰度的倍数变化之间的Spearman秩相关。通过计算双生子中一个个体的左眼和右眼的基因表达或微生物丰度的平均值，将一个双生子的平均值除以共生双生子的平均值，并进行log₂转换来获得倍数变化。对于全转录组关联分析，我们使用lme4 R包拟合线性混合模型来分析所有宿主基因表达与微生物丰度的关联。在每个模型中，我们将宿主基因对微生物物种进行回归，将家庭作为随机截距校正。这个测试可以被视为配对设计（双生子对中的两个个体）的配对t检验的泛化，线性混合模型的p值使用lmerTest R包计算。与之前的研究一致，使用Benjamini-Hochberg FDR < 0.1来确定显著的基因-微生物对。

通路富集分析

我们使用clusterProfiler R包基于京都基因与基因组百科全书（KEGG）数据库和基因本体（GO）数据库进行通路富集。

基于基因型和微生物组的基因表达建模

基因表达预测模型以表达为响应变量，以SNP剂量、微生物物种丰度和交互项为预测因子进行拟合。模型定义为：

其中Y是基因的表达，w_k是SNP k对基因的效应大小，S_k是SNP k的等位基因剂量；w_j是微生物j对基因的效应大小，M_j是微生物j的log₂转换丰度；w_i是SNP-微生物交互i对基因的效应大小，X_i是SNP-微生物对i的交互项，ε是假设与遗传和微生物组分独立的其他因素的贡献。

基因表达使用TMM进行标准化，并使用线性模型调整年龄、性别、前五个基因型PC和前15个PPER因子，残差用于拟合模型。纳入从基因起始上游1 Mb到基因结束下游1 Mb的SNP和所有微生物物种。SNP-微生物交互项，即SNP剂量和微生物丰度的乘积，用于表示GMI。在拟合包含交互的模型之前，先拟合不包含交互的模型，以识别用于构建交互项的重要SNP及微生物：

弹性网络使用glmnet R包进行拟合。对于在不包含交互的模型中具有非零系数的每个SNP和微生物，构建了交互项。使用stats和lme4 R包分别拟合了多元线性回归和线性混合模型，其中线性混合模型将家庭作为随机截距校正。对于未正则化的模型，为在不包含交互的模型中具有名义显著性（p < 0.05）的每个SNP和微生物构建了交互项。

基因预测准确性的交叉验证

嵌套交叉验证的执行方式如下：我们首先从100对双生子（n = 400）中随机选择样本作为训练数据，剩余的样本作为保留的测试数据。通过在家庭水平选择样本，我们确保了同一对双生子的样本在相同的集合中（无论是作为训练还是测试数据），以避免数据泄露。该过程针对每个折（fold）重复，测试集中的决定系数（R²）的平均值用于衡量准确性。

MicroGenix方法已经实现为一个R包，包含三个模块：训练、预测和关联。训练模块接受SNP剂量、微生物丰度和分子表型（如基因表达）作为输入，并执行两阶段拟合以构建包含交互的模型。预测模块使用这些构建的模型从SNP剂量和微生物丰度预测分子表型。关联模块使用逻辑或线性回归测试预测的分子表型与复杂性状（如疾病）之间的关联，这取决于性状是定性的还是定量的。它还从与性状相关的分子表型模型中检索SNP-微生物交互项。

识别与过敏性结膜炎相关的GMI

利用在双生子数据集中拟合的模型，预测来自健康个体和过敏性结膜炎或干眼病患者的宏基因组数据集的基因表达。分析中包括了年龄和性别匹配的健康个体和患者的样本。使用逻辑回归分析表达量预测值与疾病的关联。分析纳入在包含交互的模型中R² > 0.01的2086个基因。采用了全基因显著性阈值p < 2.4×10^-5（0.05/2086），遵循全转录组关联研究的常见做法。

（▼ 点击跳转）

“iMeta” 是由威立、肠菌分会和本领域数百千华人科学家合作出版的开放获取期刊，主编由中科院微生物所刘双江研究员和荷兰格罗宁根大学傅静远教授担任。目的是发表所有领域高影响力的研究、方法和综述，重点关注微生物组、生物信息、大数据和多组学等。目标是发表前10%(IF > 20)的高影响力论文。期刊特色包括视频投稿、可重复分析、图片打磨、青年编委、前3年免出版费、50万用户的社交媒体宣传等。2022年2月正式创刊发行！发行后相继被Google Scholar、ESCI、PubMed、DOAJ、Scopus等数据库收录！2024年6月获得首个影响因子23.7，位列全球SCI期刊前千分之五(107/21848)，微生物学科2/161，仅低于Nature Reviews，同学科研究类期刊全球第一，中国大陆11/514！

“iMetaOmics” 是“iMeta” 子刊，主编由中国科学院北京生命科学研究院赵方庆研究员和香港中文大学于君教授担任，是定位IF>10的高水平综合期刊，欢迎投稿！