基于多組全基因組表達譜的陽虛人群關鍵候選基因集和通路篩選及生物信息學分析

林景峰　高強　劉甘露　馬華萍　胡文悅　韓振蘊

〔摘要〕 目的 基于多組全基因組表達譜篩選陽虛人群關鍵候選基因集和通路，探討不同陽虛人群與對照組差異基因的異同，提出陽虛與基因表達關系的可能結論。方法 查找GEO基因表達數據庫及PubMed、Embase、CNKI、萬方、維普等中英文文獻數據庫，篩選出其中存在陽虛人群及其對照組的基因表達譜數據或基因表達譜分析結果。應用R語言和生物信息學方法篩選差異表達基因，并進行GO、KEGG和GSEA富集分析，利用韋恩圖展現不同陽虛人群的差異表達基因結果的關系，并分析差異基因、富集結果與陽虛之間的關系。結果 共得到2個數據集（GSE87474、GSE56116）和4個基因表達譜分析結果，數據集GSE56116中存在350個差異表達基因，數據集GSE87474中存在138個差異表達基因。4個基因表達譜進行去重與合并后，形成3個差異基因集，分別報道了190個、66個和21個差異表達基因。對差異表達基因結果取交集，未發現重合的基因。對其中差異表達基因相關的基因集和通路取交集分析，尋找到8個共同的基因集和通路。這些通路的功能集中在免疫功能、細胞質囊泡等方面。結論 陽虛人群與非陽虛人群的差異基因集和通路主要與能量代謝、細胞質囊泡及免疫調節相關。不同陽虛人群間關鍵候選基因集和通路存在一定的相似性，但其作用的具體基因存在差異。中醫的陽虛概念更可能與一系列基因集和通路存在緊密的聯系，而不是單一的基因。

〔關鍵詞〕 陽虛;全基因組表達;基因集;通路;生物信息學

〔中圖分類號〕R2-0? ? ? ? 〔文獻標志碼〕A? ? ? ?〔文章編號〕doi：10.3969/j.issn.1674-070X.2021.03.020

〔Abstract〕 Objective To screen the key candidate gene sets and pathways in the Yang deficiency population based on multiple genome-wide expression profile sets. To explore the differences and similarities between the different genes in the Yang deficiency syndrome population and the control group， and to propose a possible hypothesis on the relationship between Yang deficiency and gene expression. Methods Searching for GEO gene expression database and PubMed， Embase， CNKI， Wanfang， Viper databases in language of Chinese and English， and screened out gene expression profile data sets or gene expression profile analysis results of Yang deficiency group and its control group. Differentially expressed genes were screened by R software and bioinformatics methods， and gene ontology （GO） analysis， kyoto encyclopedia of genes and genomes （KEGG） and gene set enrichment analysis （GSEA） were carried out. The Venn diagram was used to show the relationship between the results of differentially expressed genes in different people with Yang deficiency， and the relationship between the results of differentially expressed genes， enrichment and Yang deficiency was analyzed. Results Two data sets （GSE87474， GSE56116） and four gene expression profiles were obtained. There were 350 differentially expressed genes in the data set GSE56116 and 138 differentially expressed genes in the data set GSE87474. Four gene expression profiles were selected and combined to form three differentially expressed gene sets. 190， 66 and 21 differentially expressed genes were reported respectively. The results of differentially expressed genes were overlapped and no coincidence genes were found. The gene sets and pathways related to differentially expressed genes were analyzed and 8 common gene sets and pathways were found. The functions of these pathways were concentrated in immune function， cytosolic vesicles and so on. Conclusion The different gene sets and pathways between Yang deficiency and non-Yang deficiency populations are mainly related to energy metabolism， cytoplasmic vesicles and immune regulation. The key candidate gene sets and pathways have some similarities among different Yang deficiency populations， but the expressions of specific genes are different. The concept of Yang deficiency in Chinese medicine is more likely to be related to a series of gene sets and pathways rather than a single gene.

〔Keywords〕 Yang deficiency; genome-wide expression; gene set; pathway; bioinformatics

陽虛人群包括陽虛證患者和陽虛體質人群。陽虛證、陽虛體質是中醫常見的一種證候和體質，通常都有畏寒怕冷、手足不溫、喜熱飲食等一系列表現[1]。近年來，中醫證候和體質的基因組學、轉錄組學和蛋白組學研究逐漸成為分子機制研究的熱點，而基因組學整體性及穩定性特點與中醫證候特性有異曲同工之處[2]，利用全基因組表達譜對具有陽虛證候群的人群進行生物信息學分析具有可行性。既往亦有不同陽虛人群與對照組基因表達狀況的相關研究[3-7]。但此類研究有樣本量較少等不足，同時，不同陽虛人群之間的基因表達差異尚不明確，而其共同差異基因、共同基因集和通路可能與陽虛證候群有較為確切的聯系。本研究團隊運用生物信息學相關方法對陽虛人群相關基因芯片數據進行分析，結合既往陽虛人群的基因表達研究結果，探討不同陽虛人群與對照組差異基因、基因集和通路的異同，提出陽虛與基因表達關系的可能結論。

1 資料與方法

1.1? 數據獲取

查找基因表達數據庫（gene expression omnibus， GEO）及PubMed、Embase、CNKI、萬方、維普等中英文文獻數據庫，篩選出其中存在陽虛證候或體質及其對照組的基因表達譜數據（包括基因芯片與轉錄組測序數據）或基因表達譜分析結果（差異基因集合）。檢索時限均為建庫至 2020年2月1日。檢索采用主題詞和自由詞相結合的方式，并追溯納入文獻的參考文獻，以補充獲取相關文獻。中文檢索詞包括：陽虛、基因表達、轉錄組測序、基因芯片等，英文檢索詞包括：RNA expression， gene expression， expression profiling， Yang deficiency， TCM syndrome，

humans。未進行灰色文獻的人工檢索。納入標準：（1）研究對象為人類。（2）若為分析結果，則報道的差異基因數量需>10個。（3）試驗組為“陽虛體質”“陽虛證”或陽虛相關證型，對照組無陽虛或中醫其他證型。得到基因表達數據或基因表達譜分析結果后，對異質性較小的文獻結果進行合并，以增加其差異基因數量。

1.2? 基因表達譜數據

陽虛人群基因表達數據集來源于GEO數據庫，分別為GSE87474和GSE56116。GSE87474包含32個中國漢族人的外周血單核細胞樣本，分別為12個陽虛體質、12個陰虛體質和8個平和體質樣本。選取其中陽虛體質和平和體質樣本進行生物信息學分析。GSE56116包含13份中國絕經后骨質疏松癥患者外周血樣本，其中辨證為腎陰虛者4份、腎陽虛3份、無腎虛患者3份、健康對照組3份。選取其中腎陽虛和無腎虛患者人群進行生物信息學分析。另有數據集GSE67090為陽虛四逆證人群與健康人群的比較，但其原始文件與表達矩陣數據均非基因表達值，無法進行二次分析，故剔除。使用R語言（版本 3.6.2）進行數據集的生物信息學分析。

1.3? 數據集預處理

下載GSE56116和GSE87474的探針表達矩陣，對其進行Log2轉化和質量控制分析。然后根據soft注釋文件將表達矩陣與基因名和entrez ID進行對應。多個探針對應同一個基因的，取最大表達量探針。

1.4? 差異表達基因的篩選

數據集的差異基因需同時滿足|Log2 fold change （Log2FC）|>1且P<0.05。利用pheatmap包里對前30顯著差異基因（以|Log2FC|為標準）繪制熱圖，直觀地展示每個差異基因在每個樣本中的表達情況。利用Enhanced Volcano包對所有基因繪制火山圖，直觀地展示每個基因在每個樣本中的表達情況。同時提取文獻中報道的差異基因集合。

1.5? 差異表達基因的基因本體論（gene ontology， GO）、通路富集分析（kyoto encyclopedia of genes and

genomes， KEGG）和基因富集分析（gene set enrichment analysis， GSEA）

對多組差異基因集合取并集得到差異基因并集。利用R語言clusterProfiler 包對多組差異基因集合及差異基因并集進行GO分析、KEGG分析。對GSE56116和GSE87474數據集進行GSEA分析。設定顯著性基因富集的臨界值為P<0.05。

1.6? 多組數據集的共同差異基因、富集基因集和通路的篩選

利用Venn Diagram包繪制多組共同差異基因、基因集和通路的交集。利用intersect函數篩選得到共同基因集和通路。

2 結果

2.1? 多組數據的獲得和數據集差異基因篩選

共獲得2個基因組表達譜數據集（GSE56116和GSE87474）。兩組數據集分別有32 696個、28 267個有效基因探針結果，對其進行Log2轉化。清洗未標注基因名的探針、多個探針對應同一個基因取最大表達量探針后，分別得到擁有19 749個、17 447個基因的表達矩陣。R語言進行差異基因篩選，GSE56116共篩選出350個差異基因，其中92個上調基因、258個下調基因。GSE87474共篩選出138個差異基因，其中96個上調基因、42個下調基因。所有基因表達的火山圖見圖1。

選取兩個數據集中差異最大的30個基因繪制熱圖，見圖2。

其中GSE56116數據集中差異最大的30個基因里包括4個上調基因和26個下調基因。GSE87474數據集中差異最大的30個基因里包括2個上調基因和28個下調基因。

在文獻的基因表達譜分析結果方面，提取到4個分析結果。其中，李艷艷[4]報道了差異最大的20個基因，其余3個分析結果[3，5，7]報道了其研究中的所有差異基因。GSE87474在其注釋文件中未對納入人群的性別和年齡進行說明。Cheng[7]在論文中未說明納入受試者的性別和年齡。其中湯朝暉[3]和李艷艷[4]的研究為同一機構、相近時間的研究，且納入試驗的受試者信息相近，為增加差異基因總數，便于更好地進行富集分析，將兩者差異基因取并集進行合并。6個差異基因集合進行去重與合并后，形成5個差異基因集，分別為GSE87474 138個，GSE56116 350個，Cheng[7]190個，湯朝暉[3]+李艷艷[4]66個和楊嘉慧[5]21個差異表達基因。

2.2? 多組差異基因的交集分析及差異基因并集的建立

對獲取的5組差異基因取交集進行分析。結果顯示，沒有基因同時存在于5組、4組及3組的差異基因集交集中。

2.3? GO功能富集

將5組差異基因集取并集后共得到742個差異基因。GO功能富集分析包括生物過程（biological process， BP）、細胞組成（cellular component， CC）和分子功能（molecular function， MF）3個方面。利用clusterProfiler包對差異基因并集進行GO功能富集分析，選取前10位的功能富集類別，見圖3。富集功能主要涉及中性粒細胞介導免疫、粒細胞趨化性、趨化因子介導的信號通路、細胞質囊泡腔等。

2.4? KEGG通路富集

對差異基因并集分別進行KEGG通路富集分析，選取前10位的通路，見圖4。基因通路主要涉及造血細胞譜系、移植物抗宿主病、抗原處理及呈遞等。

2.5? GSEA基因富集分析

對GSE87474和GSE56116進行GSEA富集分析。選取P<0.05的富集通路，其中GSE87474存在1個富集通路，GSE56116存在2個富集通路。

在GSEA分析中，從2個數據集中篩選出了3個基因富集功能，無共有富集功能。相比于GO基因富集分析和KEGG基因通路分析來說，其結果數量太少，未具有代表性。

2.6? 多組數據集的共同富集基因集和通路的篩選

對各組的差異基因分別進行GO功能3個方面的富集分析。由于楊嘉慧[5]發表的差異基因數量過少，不適于進行基因富集和通路分析，對其進行剔除。選取P<0.05的剩余4組GO功能富集結果，并繪制韋恩圖，見圖5。

在有關細胞組成的GO功能中可見有3個GO功能集合為4個組所共有。這3個功能集合分別為GO：0060205、GO：0031983、GO：0034774;其功能分別為細胞質囊泡腔、囊泡腔、分泌顆粒腔。

對剩余4個組的差異基因分別進行KEGG通路富集分析。選取P<0.05的剩余4組KEGG通路富集分析結果，并繪制韋恩圖，見圖6。

圖中可見有5個通路為4組所共有。這5個通路分別為hsa04640、hsa05332、hsa05140、hsa04940、

hsa05321;其功能分別為造血細胞系（免疫功能相關）、移植物抗宿主病（免疫功能相關）、利什曼病（免疫與炎癥反應相關）、I型糖尿病、炎癥性腸病（炎癥反應相關）。

對超過3個組所共有的差異基因、GO富集基因集、KEGG富集通路進行統計。共有的差異基因數目為0，共有的GO富集集合（細胞組成）數為22個，共有的GO富集集合（生物過程）數為39個，共有的GO富集集合（分子功能）數為1個，共有的KEGG通路為23條。利用χ2檢驗比較其間的差異，發現差異有統計學意義（P<0.05）。結果表明：不同陽虛人群之間更有可能也更容易存在共有的基因集和通路，而非單個差異基因。

3 討論

中醫陽虛概念與全基因組表達的關系目前尚不明確。本研究選取了GSE87474、GSE56116兩組GEO數據集中陽虛體質與正常體質，絕經后骨質疏松癥腎陽虛患者與無腎陽虛患者樣本進行分析，將陽虛證候或體質作為變量，對其進行生物信息學分析，以探討陽虛證候或體質與其基因表達的關系。本研究發現，陽虛人群與非陽虛人群的差異基因集和通路主要與能量代謝、細胞質囊泡及免疫調節相關，且不同陽虛人群間關鍵候選基因集和通路存在一定的相似性。各個納入分析的差異基因集的通路分析結果具有較好的一致性。與相關研究相比較，研究結果亦具有一致性。在GO基因富集分析中，研究團隊選出了3個共同GO富集集合：“GO：0060205”“GO：0031983”“GO：0034774”，主要涉及細胞質囊泡。在KEGG基因通路分析中選出了5個共同基因通路：“hsa04640”“hsa04940”“hsa05240”“hsa05321” “hsa05332”，主要與免疫調節等相關。既往研究中，GSE87474與GSE56116數據集下尚無數據集上傳者的論文，李艷艷[4]報道了“hsa04640”“hsa04940”，湯朝暉[3]報道了“hsa04940”，楊嘉慧[5]報道了囊泡運輸途徑和“hsa04640”。其他相關研究中，TANG等[8]利用四逆湯治療腎陽虛證模型大鼠腎上腺基因表達變化尋找相關基因集與通路，發現四逆湯治療腎陽虛證涉及能量代謝與線粒體功能相關通路。XU等[9]利用 GSE57273和GSE56116數據集探討六味地黃丸治療絕經后骨質疏松癥腎陰虛證的靶點與基因通路，得到的差異基因、基因集和通路主要涉及免疫功能等方面。LIU等[10]對氣虛血瘀型和陰虛血瘀型缺血性腦卒中雄性大鼠模型的基因表達譜進行分析，發現兩組間存在一定數量的相同差異基因，這些差異基因可能與氣虛、血瘀或陰虛有關，功能主要集中在能量代謝方面，并利用免疫印跡法驗證了富集基因和通路的產物。但其實驗對象全部為大鼠，組間差異小，難以體現出人類人群間的多樣性對實驗結果的影響。GUAN等[11]對腎陰虛大鼠基因表達情況進行分析，發現其差異基因功能亦主要集中于能量代謝等方面。此外，多項研究[12-19]亦從中醫方藥與針灸治療角度分析了治療前后基因表達變化，但對證型間的差異尚無涉及。盡管不同證候、不同物種的基因表達譜應有差異[6]，但本研究團隊總結既往類似研究結果，綜合本研究結果，認為陽虛人群與非陽虛人群間的基因表達差異可能多與能量代謝、免疫調節相關基因集和通路有關。

從差異基因的篩選結果來看，GSE56116共篩選出350個差異基因，GSE87474共篩選出138個差異基因，疾病狀態下的差異基因數量亦大于健康人群中的差異基因數量。本研究團隊認為部分差異基因的表達差異可能由疾病或腎虛造成，從而導致篩選出差異基因數量增多。

對多組差異基因取交集后發現并無共同差異基因。而進行GO功能富集和KEGG基因通路分析后，發現其存在多個共同基因集和通路。見圖5-6。對共同差異基因或基因集合數量進行χ2檢驗，發現差異有統計學意義（P<0.05）。本研究團隊認為陽虛人群在基因組表達層面更有可能與一些基因集和基因通路相關，而非單一的基因。

本研究的不足之處：（1）本數據集數據量偏少，且各個數據集使用的平臺并不完全一致，納入受試者的年齡、性別亦存在一定的差異。本研究的研究結果需要更大、更多元、相同平臺的中醫陽虛人群基因表達數據來進行證實。（2）納入的兩個GEO數據集存在一定的異質性，此異質性可能會對分析結果形成影響。

本研究表明：陽虛人群與非陽虛人群的差異基因集和通路主要與能量代謝、細胞質囊泡及免疫調節相關，不同陽虛人群間關鍵候選基因集和通路存在一定的相似性，且本研究的結果與既往相關研究具有較強一致性。疾病狀態下的陽虛證候、體質人群的基因表達與健康人的陽虛體質基因表達存在一定的差異，今后開展類似研究可能需考慮此方面的差異。同時，不同陽虛人群的基因表達與基因集和通路的關系可能比單個基因更為密切，此觀點亦需大數據層面的驗證。最后，如何更好地從全基因組差異表達層面闡釋陽虛概念乃至中醫證候、體質、證型的形成機制，需在大數據支持下進一步探索。

參考文獻

[1] 黃瑞聰，李美紅，陳夢華，等.膏肓灸對陽虛質人人群質量表積分及紅外熱圖的影響研究[J].中醫臨床研究，2019，11（30）：70-72.

[2] 曾召瓊，易? 帆，李? 萍，等.基因組學技術在中醫藥研究中的應用[J].國際檢驗醫學雜志，2018，39（24）：3089-3092.

[3] 湯朝暉.老齡腎陽虛證候的診斷及其差異表達基因譜研究[D].成都：成都中醫藥大學，2008.

[4] 李艷艷.老年腎陽虛證的差異基因表達研究[D].成都：成都中醫藥大學，2009.

[5] 楊嘉慧.腎陽虛證排卵障礙性不孕的差異基因表達譜研究[D].成都：成都中醫藥大學，2012.

[6] 譚從娥.腎陽虛證膝骨關節炎的差異基因表達譜研究[D].成都：成都中醫藥大學，2006.

[7] CHENG H T， CHEN C R， LI C Y， et al. The classification of sini decoction pattern in traditional Chinese medicine by gene expression profiling[J]. Evidence-Based Complementary and Alternative Medicine： ECAM， 2016， 2016： 8239817.

[8] TANG N， LIU L H， QIU H， et al. Analysis of gene expression and functional changes of adrenal gland in a rat model of kidney Yang deficiency syndrome treated with Sini decoction[J]. Experimental and Therapeutic Medicine， 2018， 16（4）： 3107-3115.

[9] XU F， GAO F. Liuwei Dihuang pill cures postmenopausal osteoporosis with kidney-Yin deficiency[J]. Medicine， 2018， 97（31）： e11659.

[10] LIU T L， LIU M N， XU X L， et al. Differential gene expression profiles between two subtypes of ischemic stroke with blood stasis syndromes[J]. Oncotarget， 2017， 8（67）： 111608-111622.

[11] GUAN W， LIU Y， LI X M， et al. iTRAQ-based proteomics to reveal the mechanism of hypothalamus in kidney-Yin deficiency rats induced by levothyroxine[J]. Evidence-Based Complementary and Alternative Medicine， 2019， 2019： 3703596.

[12] 周? 嵐，李? 楊，汪? 典，等.利用基因芯片分析活血中藥、破血中藥對ApoE基因敲除小鼠動脈粥樣硬化模型的差異表達基因[J]. 疑難病雜志，2017，16（1）：18-22.

[13] 鄭? 玫，連? 方，孫振高.補腎中藥對腎虛不孕患者卵巢顆粒細胞基因表達干預的研究[J].遼寧中醫雜志，2017，44（1）：2-15，221.

[14] 魏歆然，魏高文，鄭雪娜，等.不同經穴組合針刺對失眠大鼠下丘腦生物鐘基因Clock和Bmal 1表達的影響[J].針刺研究，2017，42（5）：429-433.

[15] 王? 俊，虞彬艷，周斯斯，等.粗針神道穴平刺對周圍性面癱大鼠面神經基因表達譜的影響[J].中華中醫藥雜志，2016，31（1）：287-291.

[16] 賈文睿.針刺對應激性高血壓前期大鼠心臟基因表達譜的影響[D].北京：北京中醫藥大學，2017.

[17] WU P， HUANG R， XIONG Y L， et al. Protective effects of curcumin against liver fibrosis through modulating DNA methylation[J]. Chinese Journal of Natural Medicines， 2016， 14（4）： 255-264.

[18] GUO Y， XIE X J， GUO C Q， et al. Effect of electro-acupuncture on gene expression in heart of rats with stress-induced pre-hypertension based on gene chip technology[J]. Journal of Traditional Chinese Medicine， 2015， 35（3）： 285-294.

[19] HUANG Y L， WAN M Y， LIANG X S， et al. Effect of acupuncture along affected meridian on the MME gene expression of migraine patients without aura of Gan-Yang hyperactivity syndrome[J]. Chinese Journal of Integrated Traditional and Western Medicine， 2015， 35（3）： 294-298.