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

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

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

〔關(guān)鍵詞〕 陽虛;全基因組表達;基因集;通路;生物信息學(xué)

〔中圖分類號〕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

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

1 資料與方法

1.1? 數(shù)據(jù)獲取

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

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

1.2? 基因表達譜數(shù)據(jù)

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

1.3? 數(shù)據(jù)集預(yù)處理

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

1.4? 差異表達基因的篩選

數(shù)據(jù)集的差異基因需同時滿足|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數(shù)據(jù)集進行GSEA分析。設(shè)定顯著性基因富集的臨界值為P<0.05。

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

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

2 結(jié)果

2.1? 多組數(shù)據(jù)的獲得和數(shù)據(jù)集差異基因篩選

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

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

其中GSE56116數(shù)據(jù)集中差異最大的30個基因里包括4個上調(diào)基因和26個下調(diào)基因。GSE87474數(shù)據(jù)集中差異最大的30個基因里包括2個上調(diào)基因和28個下調(diào)基因。

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

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

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

2.3? GO功能富集

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

2.4? KEGG通路富集

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

2.5? GSEA基因富集分析

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

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

2.6? 多組數(shù)據(jù)集的共同富集基因集和通路的篩選

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

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

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

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

hsa05321;其功能分別為造血細胞系（免疫功能相關(guān)）、移植物抗宿主?。庖吖δ芟嚓P(guān)）、利什曼?。庖吲c炎癥反應(yīng)相關(guān)）、I型糖尿病、炎癥性腸?。ㄑ装Y反應(yīng)相關(guān)）。

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

3 討論

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

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

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

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

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

參考文獻

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

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

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

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

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

[6] 譚從娥.腎陽虛證膝骨關(guān)節(jié)炎的差異基因表達譜研究[D].成都：成都中醫(yī)藥大學(xué)，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] 鄭? 玫，連? 方，孫振高.補腎中藥對腎虛不孕患者卵巢顆粒細胞基因表達干預(yù)的研究[J].遼寧中醫(yī)雜志，2017，44（1）：2-15，221.

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

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

[16] 賈文睿.針刺對應(yīng)激性高血壓前期大鼠心臟基因表達譜的影響[D].北京：北京中醫(yī)藥大學(xué)，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.