基于RNA测序的A亚型呼吸道合胞病毒毛细支气管炎转录组学特征分析及关键基因预测

李海燕，翁婷婷，王乐颖，梁春婵，董琳

温州医科大学附属第二医院育英儿童医院 儿童呼吸科，浙江 温州 325027

呼吸道合胞病毒（respiratory syncytial virus, RSV）是引起婴幼儿下呼吸道感染的重要病毒，在婴幼儿中有较高的发病率和致死率[1]。其中毛细支气管炎是RSV感染导致的最常见的下呼吸道感染，占所有毛细支气管炎住院患儿的50%～80%[2]。根据病毒G蛋白基因的不同，RSV可分为A、B亚型[3]。不同亚型RSV感染的婴幼儿毛细支气管炎的免疫学机制较为复杂，目前尚不明确。现仍无有效的疫苗和特效治疗药物，只能对症支持治疗。因此需对可能的发病机制进行深入探讨，为临床干预和治疗提供科学依据。

转录组学是研究病毒与宿主细胞相互作用的重要手段，为感染性疾病的预防和治疗以及揭示疾病的发病机制提供更多的研究方法[4-6]。目前国内尚少见不同亚型RSV感染转录组学的研究报道，鉴于不同国家人群的遗传背景不尽相同，有必要对我国不同亚型RSV感染儿童开展转录组学研究。本项目前期从收集的23个RSV毛细支气管炎患儿鼻咽分泌物中分离出11个RSV病毒样品，最终鉴定均是A亚型RSV，同时应用RNA高通量测序方法对A亚型RSV毛细支气管炎患儿全血基因转录组学进行分析，并与性别、年龄配对的健康儿童进行对比，筛选差异表达基因（differentially expressed genes, DEGs），进一步分析差异基因表达功能，以期从分子水平了解A亚型RSV感染的宿主遗传与免疫反应特点，阐明A亚型RSV毛细支气管炎的潜在发病机制。

1 对象和方法

1.1 对象 招募2020年1月1日至2020年12月31日RSV感染的毛细支气管炎患儿作为病例组。病例组纳入标准：①1～24月龄，汉族；②首次出现咳嗽、喘息[符合毛细支气管炎诊断、治疗与预防专家共识（2014版）诊断标准]；③鼻咽分泌物免疫荧光检测RSV阳性。病例组排除标准：①早产，合并慢性肺部疾病、先天性心脏病、免疫缺陷病、神经肌肉疾病等其他慢性疾病；②14 d内应用过全身糖皮质激素治疗；③合并呼吸系统其他病毒或细菌感染；④合并其他系统病毒或细菌感染。健康对照组纳入标准：①1～24月龄，汉族；②同期在儿童保健科体检的健康儿童；③既往无喘息病史，1个月内无呼吸道感染史。本研究经温州医科大学附属第二医院育英儿童医院伦理委员会批准（编号：LCKY2019-148），所有研究对象监护人知情同意。

1.2 标本采集、RSV病毒亚型鉴定、RNA提取 标本采集：入院当天采集毛细支气管炎患儿鼻咽分泌物1 mL置于2 mL 0.9%氯化钠溶液中，应用直接免疫荧光法筛查RSV抗原；健康儿童于就诊时采血，RSV阳性患儿在入院48 h内采集外周静脉血1 mL置于EDTA抗凝管，均置于-80 ℃冰箱保存，用于全血白细胞分离和RNA提取。

RSV病毒亚型鉴定：提取病毒总RNA并反转录，利用Nest-PCR法扩增目的基因，凝胶电泳初步鉴定RSV亚型；同时把剩余RNA送至生工生物工程（上海）股份有限公司进行反转录、扩增目的基因、1%琼脂糖凝胶电泳检测，通过条带回收、测序分析、序列比对等流程进一步明确鉴定RSV为A或B亚型。

外周血白细胞RNA提取：红细胞裂解液裂解全血中的红细胞，加入适量TRIzol裂解细胞，按照操作说明提取总RNA；取2 μL样品于核酸蛋白检测仪上测OD260和OD280值；剩余RNA用于建库测序。

1.3 RNA测序、测序数据前期处理 将检测合格的RNA样品送至上海浦东解码生命科学研究院利用磁珠富集mRNA，并将mRNA片段化，加入随机引物和反转录酶合成cDNA，最后利用核酸外切酶、核酸聚合酶进行末端修复，末端添加PolyA，加测序接头后进行PCR扩增，完成cDNA文库制备。将质量检测合格的cDNA文库用Illumina高通量测序平台进行测序，最终获得原始测序数据。

通过对原始下机数据进行质量评估、去除低质量序列和接头污染等质控和预处理过程，得到高质量测序数据用于下游分析。同时使用DEseq2软件包对测序数据进行标准化，并进行二维主成分分析（principal component analysis, PCA）和样本间的层次聚类分析，观察病例组和对照组之间的相似性和差异。

1.4 DEGs筛选 通过三种R软件包（DESeq2、edgeR、limma）进行病例组和对照组之间基因表达差异分析，并将FDR值≤0.05且差异倍数（fold change, FC）绝对值≥2作为DEGs的筛选条件。最终将三种R包同时筛选出的DEGs用于后续分析，并使用R包ggplot2将最终筛选的DEGs绘制成可视化火山图。

1.5 功能富集分析 利用R包clusterProfiler中的enrichGO模块基于基因本体论（gene ontology,GO）对差异基因进行生物学功能富集分析，以校正后P值≤0.05作为阈值筛选出显著富集的通路。同时利用Metascape在线工具基于京都基因与基因组百科全书（Kyoto Encyclo-pedia of Genes and Genomes, KEGG）数据库对DEGs进行信号通路富集分析，以校正后P值≤0.05作为阈值筛选出显著富集的通路。

1.6 蛋白相互作用（p r o t e i n-p r o t e i n interaction，PPI）网络分析、关键模块鉴定、高风险关键基因（Hub基因）筛选 通过STRING数据库（https∶//cn.string-db.org/）预测DEGs的蛋白相互作用关系，在Cytoscape中可视化PPI网络，利用通过MCODE插件对整个PPI网络进行聚类分析，筛选出差异显著的核心模块。同时分别使用Cytoscape中的插件“Cytohubba”和“CytoNCA”获取高风险Hub基因。

1.7 统计学处理方法 本研究所涉及的统计学检验均由R软件（v4.1.0）各种分析包及在线分析工具完成。P＜0.05为差异有统计学意义。

2 结果

2.1 RSV亚型鉴定 通过对23个RSV毛细支气管炎患儿鼻咽分泌物中分离出的RSV病毒样品进行RNA提取、反转录、扩增目的基因、1%琼脂糖凝胶电泳检测，发现有11个样本扩增出RSV病毒特异性片段，约499 bp。通过条带回收、测序分析、序列比对等流程，11个RSV病毒样品全为A亚型（见图1）。

图1 1%琼脂糖凝胶电泳检测RSV亚型

2.2 RNA测序数据质控及样本聚类分析 对11个鉴定出RSV亚型的患儿全血样品以及10个健康对照提取RNA进行RNA测序，通过对原始数据质控分析，发现有1例编号为3的患儿测序数据质量较差，因此排除该病例，最终得到10例患儿和10例健康对照的RNA测序数据。同时利用PCA的方法对测序数据分析，结果显示A亚型RSV感染的病例组与健康对照组的表达数据具有明显的差异性，两个主成分对数据差异的贡献度达到51.6%（见图2A）。样本的聚类分析结果也显示出病例组与健康对照组形成了两个明显不同的簇，进一步反映了这两组样品的差异性（见图2B）。

图2 样本测序数据质量控制分析

2.3 筛选DEGs 将FDR≤0.05且|FC|≥2作为DEGs的筛选阈值，利用三种差异表达分析的R包（DESeq2、edgeR、limma）对病例组和对照组之间进行基因表达差异分析。结果如图3A所示，DESeq2包筛选出648个DEGs，edgeR包筛选出610个DEGs，limma包筛选出318个DEGs。为了进一步筛选出可靠的DEGs，我们从三种方法中同时鉴定出286个DEGs，火山图显示有166个基因上调，120个基因下调，并用于后续的功能分析（见图3B）。

图3 DEGs筛选

2.4 差异基因的GO和KEGG富集分析 对286个DEGs进行功能注释和通路富集分析，GO功能富集显示DEGs主要参与过氧化氢代谢过程、细胞对生物刺激的应激反应、维生素D受体信号通路、氧气转运、受体信号活性的调控等生物学过程（见图4A）。KEGG通路富集分析显示DEGs显著富集在免疫系统中的细胞因子信号、适应性免疫系统、中性粒细胞脱粒、淋巴细胞与非淋巴细胞免疫调控互作、细胞循环等与免疫反应有关的信号通路（见图4B）。

图4 DEGs功能富集分析

2.5 PPI网络构建及筛选核心子网络 将286个DEGs导入STRING数据库构建PPI网络，并利用Cytoscape软件可视化与分析，得到239个节点（Node）和662条线（Edge）的互作网络（见图5A）。然后使用MCODE对整个网络进行聚类分析，以K-core值为3作为阈值，最终筛选到4个关键的PPI模块，这些模块内的基因主要涉及细胞循环、病毒蛋白与细胞因子及受体的相互作用、淋巴细胞与非淋巴细胞的免疫调控互作、激酶级联反应等生物学免疫相关过程（见图5B）。

图5 PPI网络筛选

2.6 与A亚型RSV感染相关的Hub基因筛选 利用Cytoscape软件的CytoHubba插件中MCC算法从上述的PPI网络中筛选出排名前15位的节点基因（见图6A）。同样采用CytoNCA插件中的Degree算法从PPI网络中筛选出排名前15位的节点基因（见图6B）。最终将两种算法共同筛选到的11个基因作为潜在的Hub基因，即RRM2、BUB1B、BUB1、MCM10、CDC45、MKI67、ASPM、NCAPG、NUSAP1、ESPL1、CDT1基因（见图6C）。

3 讨论

毛细支气管炎是婴幼儿常见的下呼吸道感染，常发生于2岁以内[7]。其中RSV是毛细支气管炎最常见病原。近年来，转录组测序已初步应用于RSV

毛细支气管炎的发病机制及抗感染药物的筛选研究，并预测与发病密切相关的Hub基因，从而有助于RSV感染的临床干预和治疗[4-6]。BUCASAS等[5]应用RNA测序方法检测了21例RSV毛细支气管炎患儿和37例健康儿童全血转录组，结果显示RSV毛细支气管炎患儿的固有免疫反应基因，尤其是干扰素信号通路广泛激活，而与抗原递呈相关的转录表达下调。翁婷婷等[4]利用22例RSV毛细支气管炎患儿的全血RNA测序数据，通过构建加权基因共表达网络最终识别与RSV毛细支气管炎密切相关的Hub基因模块及Hub基因。另外不同亚型RSV毛细支气管炎的转录组学研究在国外也有被报道。有研究应用微阵列芯片分析方法对不同RSV亚型和基因型毛细支气管炎患儿全血转录组测序，发现与B亚型感染者相比，A亚型感染患儿干扰素、炎症及中性粒细胞相关基因表达明显上调[6]。然而，目前国内尚未见不同亚型RSV感染儿童全血转录组学的研究报道，鉴于不同国家和种族人群的遗传背景不尽相同，因此，本研究进行了A亚型RSV毛细支气管炎儿童转录组学分析。

本项目通过高通量RNA测序技术，发现A亚型RSV毛细支气管炎的患儿和健康对照组存在286个DEGs，其中IL10、CCL2、MUC1等细胞因子已被报道参与RSV感染后的固有免疫反应[8-9]，另外，SIGLEC1、HMOX-1等DEGs也不同程度参与细胞抗病毒的免疫反应[10-11]。基于GO富集分析和KEGG通路分析，我们发现这些DEGs主要参与宿主免疫应答、炎症反应、代谢调控、信号转导等生物学过程或信号通路。其中DEGs也参与维生素D受体信号通路和过氧化氢代谢过程等生物学通路，已有研究发现维生素D具有RSV感染后固有免疫应答的调控作用[12-13]，RSV感染呼吸道上皮细胞后，过氧化氢代谢过程产生的氧化应激反应具有调控抗炎基因TLR3表达的作用[14]。另外，这些DEGs也参与到中性粒细胞脱粒的信号通路中，有研究显示中性粒细胞的激活和脱粒在RSV感染过程具有提高激活标志物基因表达、延长细胞凋亡、病毒清除的作用[15-16]。

多个基因之间的相互作用常常参与了疾病的发生发展过程，本研究进一步深入分析所有DEGs的PPI网络，并从中鉴定出4个包含Hub基因的蛋白质相互作用模块。这些核心子网络的基因主要与细胞循环、细胞因子-细胞因子受体信号通路、淋巴细胞-非淋巴细胞免疫调控作用、RAF/MAP激酶级联反应等通路有关，其中细胞因子-细胞因子受体信号通路已经被大量研究证实与各种免疫细胞相互作用，在病毒感染后发挥抗感染和保护机体的作用[17-18]。另外有研究发现Raf/MEK/ERK级联反应及MAP激酶的激活可能参与到RSV感染后的炎症反应、宿主细胞信号的转导等过程[19-21]。因此以上通络均在RSV感染过程中发挥重要的作用。通过进一步对PPI模块的网络拓扑分析，我们最终确定了11个Hub基因，包括RRM2、BUB1B、BUB1、MCM10、CDC45、MKI67、ASPM、NCAPG、NUSAP1、ESPL1、CDT1。其中RRM2、BUB1B、BUB1、MCM10、MKI67、ASPM、NUSAP1等基因已被报道参与呼吸系统病毒感染或其他系统病毒感染过程，并扮演一定的调控作用[22-27]，不过其参与RSV感染进程的具体机制需要进一步的探索。

本研究的不足之处是从收集的23个RSV毛细支气管炎患儿鼻咽分泌物中只分离出11个RSV病毒样品，而且鉴定都为A亚型RSV，鉴于研究发现不同亚型RSV感染患儿在临床表型上的差异，有必要进一步收集B亚型RSV毛细支气管炎患儿血液样本，进行转录组学层面的研究，并与A亚型RSV毛细支气管炎患儿全血转录组学进行比较，分析其异同。

综上所述，本研究通过RNA测序及生物信息学分析手段比较了A亚型RSV毛细支气管炎与健康儿童之间的基因表达差异，观察到这些DEGs涉及到多个有意义的生物学过程和通路，并筛选出一批潜在的与A亚型RSV毛细支气管炎相关的Hub基因，这些信号通路及Hub基因未来可成为A亚型RSV毛细支气管炎的干预和治疗靶点。