史隨隨，羅越毅，唐夢丹，祁寶玉，周 劍*，巢國俊*

1.石河子大學醫(yī)學院，新疆 石河子 832000

2.中國中醫(yī)科學院眼科醫(yī)院，北京 100040

3.北京中醫(yī)藥大學東方醫(yī)院，北京 100078

年齡相關(guān)性黃斑變性（age-related macular degeneration，AMD）是與年齡相關(guān)的視網(wǎng)膜黃斑進行性退變疾病，造成患者不可逆性中心視力下降、視物變形，為世界衛(wèi)生組織三大重點致盲眼病之一[1]。近30 年中國AMD 疾病負擔居全球前列，2019 年因AMD 致中度視力障礙患者占4.56%，重度視力障礙者占4.14%[2]。早期AMD 與玻璃膜疣沉積有關(guān)，晚期分干性、濕性2 種，其發(fā)病機制復雜，治療方法有限，診療效果不佳[3-4]。越來越多的證據(jù)表明，AMD 與衰老和能量代謝密切相關(guān)，視網(wǎng)膜色素上皮（retinal pigment epithelial，RPE）細胞利用自噬過程釋放代謝能量，使成分降解為簡單分子，并產(chǎn)生維持能量穩(wěn)態(tài)的三磷酸腺苷（adenosine triphosphate，ATP），隨著RPE 細胞的衰老，自噬過程減少，代謝能量降低，代謝碎片和蛋白質(zhì)聚集體積累，從而導致AMD[5-6]。研究也證明改善能量代謝可以延緩衰老，因此其可能是治療AMD 的一種潛在方法[7]。

隨著快速、高通量轉(zhuǎn)錄組學測序方法的普及，對現(xiàn)有的高通量測序數(shù)據(jù)進行挖掘可為醫(yī)學研究提供有效信息。Bulk RNA 數(shù)據(jù)基于細胞群體水平獲得組織的平均特性，反映細胞群體的平均差異，而Single-cell RNA 通過單個細胞的mRNA 進行轉(zhuǎn)錄組分析，識別不同細胞亞群的細胞類型、狀態(tài)和譜系，繪制細胞圖譜、鑒定細胞異質(zhì)性，可以精準解析細胞亞類功能及對疾病的影響[8-9]。本研究結(jié)合Bulk RNA 測序和Single-cell RNA 測序數(shù)據(jù)探究AMD 發(fā)病機制中能量代謝關(guān)鍵基因并對潛在的調(diào)控中藥化合物進行預測，為AMD 藥物靶點和治療藥物的研究提供新視角[10]。

1 資料與方法

1.1 Bulk RNA 和Single-cell RNA 測序數(shù)據(jù)獲取、處理及差異基因篩選

借助GEO 數(shù)據(jù)庫（https://www.ncbi.nlm.nih.gov/gds/）以“age-related macular degeneration”“AMD”作為檢索詞，種屬選擇“human sapiens”，查找下載Bulk RNAseq數(shù)據(jù)集GSE135092[11]。該數(shù)據(jù)集共537個樣本，其中來自視網(wǎng)膜組織271 例；來自RPE 組織266例。選取數(shù)據(jù)集中266例RPE組織樣本［AMD組（AMD 患者）RPE 組織黃斑區(qū)26 例，非黃斑區(qū)23 例；對照組（正常RPE 組織）黃斑區(qū)105 例，非黃斑區(qū)112 例］納入本研究分析?；赗 軟件“DESeq2”包［R 軟件命令行輸入install.packages（“DESeq2”）后運行］對數(shù)據(jù)進行標準化和分析，以校正后P＜0.05 為閾值，且|log2FC|＞2［FC 表示差異倍數(shù)（fold change）］為初篩條件，對表達數(shù)據(jù)進行差異表達分析。

同樣查找并下載AMD Single-cell RNA 數(shù)據(jù)集GSE188280[12]，該數(shù)據(jù)集包含2 例AMD 患者和1例正常人（對照）的視網(wǎng)膜和RPE 組織的黃斑區(qū)和非黃斑區(qū)樣本，選取6 例RPE 組織樣本納入本研究。使用R 軟件“Seurat”包［R 軟件命令行依次輸入if (!requireNamespace ("BiocManager", quietly＝TRUE))和install.packages ("BiocManager") 和Bioc Manager::install ("Seurat") 后運行］將GSE188280 的表達矩陣創(chuàng)建為Seurat 對象。過濾線粒體基因含量＞20%的細胞以及feature＜500 或＞6 000 的細胞，僅選取RPE 組織的細胞，最終獲得32 714 個細胞的轉(zhuǎn)錄組數(shù)據(jù)。通過“NormalizeData”函數(shù)對GSE188280 的測序深度進行標準化，調(diào)用“FindVariableFeatures”函數(shù)使用“vst”方法檢測數(shù)據(jù)集的2 000 個可變特征，并使用“ScaleData”對數(shù)據(jù)進行縮放以排除測序深度的影響。

1.2 Bulk RNA 測序中AMD 相關(guān)的能量代謝相關(guān)基因（energy metabolism-related genes，EMRGs）的蛋白相互作用（protein-protein interaction networks，PPI）、轉(zhuǎn)錄因子（transfer factor，TF）、微小RNA（microRNA，miRNA）網(wǎng)絡(luò)構(gòu)建

EMRGs 來自Gene Set Enrichment Analysis 數(shù)據(jù)庫（GSEA，https://www.gsea-msigdb.org/），最終有48 個AMD 相關(guān)EMRGs 從GSE135092 的轉(zhuǎn)錄組數(shù)據(jù)集中提取出來。使用 Wilcoxcon 檢驗比較GSE135092 RPE 組織黃斑區(qū)和非黃斑區(qū)來源的EMRGs 的表達，P＜0.05 被認為具有統(tǒng)計學差異，使用R 軟件“pheatmap”包繪制熱圖，并用Pearson分析不同EMRGs 之間的相關(guān)性。

借助STRING（http://www.string-db.org/）在線工具和Cytoscape 軟件將差異表達的EMRGs 上傳構(gòu)建PPI 網(wǎng)絡(luò)。進一步分析其TF 和miRNA 作用網(wǎng)絡(luò)，將差異表達EMRGs 上傳至Network Analyst（https://www.networkanalyst.ca/），參考TF databases ENCODE（http://cistrome.org/BETA/）和miRTarBase v8.0（https://mirtarbase.cuhk.edu.cn） 數(shù)據(jù)庫，分別使用TF-gene Interactions 和Gene-miRNA Interactions模塊進行分析。

1.3 Logistic 回歸分析評估AMD 患者EMRGs 診斷能力

Logistic 回歸模型正是作為被解釋變量與其他自變量相互關(guān)系的一種線性模型，廣泛用于醫(yī)學領(lǐng)域[13]。受試者工作特征曲線（receiver operating charateristic，ROC）作為疾病診斷的常用統(tǒng)計學分析方法，是一種全面、準確評價診斷指標的工具[14]。將有差異的EMRGs 作為自變量，將AMD 發(fā)生情況作為因變量，使用Logistic 回歸分析差異表達EMRGs與是否患有AMD的關(guān)聯(lián)性，并根據(jù)Logistic模型擬合結(jié)果利用R 軟件的“pROC”包繪制ROC曲線，并計算曲線下面積（area under curve，AUC），評估AMD 患者EMRGs 診斷能力。AUC＝0.5 時，說明關(guān)鍵靶基因無診斷價值；AUC＞0.5 時，面積越接近于1，說明關(guān)鍵靶基因診斷的準確性越高；0.5＜AUC≤0.7，說明關(guān)鍵靶基因準確性較低；0.7＜AUC≤0.9 時，說明關(guān)鍵靶基因有準確性中等；AUC＞0.9 時，說明關(guān)鍵靶基因準確性高[15]。

1.4 最小絕對值收斂和選擇算子（least absolute shrinkage and selection operator，LASSO）篩選AMD 診斷標志物及風險評分構(gòu)建

LASSO 算法是一種既能對變量進行選擇和壓縮來進行特征選擇，又能得出參數(shù)估計值的基于線性回歸模型的方法，可顯著降低特征空間維度，有效地避免過擬合和多重共線性問題，被廣泛應用和推廣[16]。使用LASSO 回歸在差異表達基因中篩選診斷標志物。使用“glmnet”包的“glmnet”函數(shù)進行變量篩選，采用“cv.glmnet”函數(shù)進行交叉驗證，然后篩選出使得變異系數(shù)最小的預后標志物組合。每個樣本的風險評分的計算公式如下。用Wilcoxcon檢驗比較AMD 和對照的風險評分，以箱式圖展示。

Coefficient 為LASSO 回歸系數(shù)，mRNA Expressiopn 為該基因的表達水平（log2轉(zhuǎn)換）

1.5 Single-cell RNA 測序分析及細胞類型注釋

主成分分析法（principal component analysis，PCA）是一種將高維空間復雜數(shù)據(jù)轉(zhuǎn)化為低維空間的常用線性降維方法[17]。均勻流形逼近和投影（uniform manifold approximation and projection，UMAP）是基于流形學習技術(shù)和拓撲分析思想的降維算法，耗時短、對嵌入維數(shù)沒有限制，具有可視化的優(yōu)點，可用于醫(yī)學、生物學數(shù)據(jù)分析[18]。研究應用PCA 識別顯著主成分（principal component，PC），用“Elbowplot”函數(shù)可視化P值分布。然后選22 個PC進行UMAP 分析以進行降維，用“FindNeighbors”默認參數(shù)以及22 個PC 維度參數(shù)來構(gòu)建PCA 空間中歐幾里得距離的 K-最近鄰域。然后通過調(diào)用“FindClusters”和“clustree”函數(shù)尋找到最佳分辨率0.5，將細胞分為19 個不同的簇。最后用“RunUMAP”函數(shù)進行降維以便探索數(shù)據(jù)集和可視化。

通過文獻中已報道在RPE 和各種脈絡(luò)膜細胞類型中表達的基因列表以及人工注釋進行本研究的細胞類型注釋，例如，T 細胞的標記基因有C 型蛋白酪氨酸磷酸酶受體（protein tyrosine phosphatase receptor type C，PTPRC）、CD3D（CD3 delta）、CD2（CD2 molecule），B 細胞的標記基因有跨膜4 域亞家族 A1（membrane spanning 4-domains A1，MS4A1）、CD79A，NK 細胞的標記基因有自然殺傷細胞顆粒蛋白 7（natural killer cell granule protein 7，NKG7）、顆粒酶A（granzyme A，GZMA），骨髓細胞的標記基因有溶菌酶（lysozyme，LYZ）、CD68、Fc 區(qū)III 受體γ（Fc gamma RIIIA，F(xiàn)CGR3A）、CD14，內(nèi)皮細胞的標記基因有血小板內(nèi)皮細胞粘附分子1（ platelet endothelial cell adhesion molecule 1，PECAM1）、血管性血友病因子（von willebrand factor，VWF），成纖維細胞的標記基因有α1-1 型膠原（collagen, type I, alpha 1，COL1A1）、COL1A2、核心蛋白聚糖（decorin，DCN）、胰島素樣生長因子2（insulin-like growth factor 2，IGF2）、纖維蛋白1（fibulin 1，F(xiàn)BLN1），施萬細胞的標記基因有髓鞘質(zhì)蛋白1（proteolipid protein 1，PLP1），黑色素細胞的標記基因有黑色素A（melan-A，MLANA）、黑素細胞蛋白（premelanosome protein，PMEL）、酪氨酸酶相關(guān)蛋白1（tyrosinase related protein 1，TYRP1），巨噬細胞的標記基因有同種異體移植炎癥因子1（allograft inflammatory factor 1，AIF1），肥大細胞的標記基因有原癌基因KIT（KIT proto-oncogene，KIT），平滑肌細胞的標記基因有蛋白信號調(diào)節(jié)因子5（regulators of G-protein signaling 5，RGS5）、α 肌動蛋白2（actin alpha 2，ACTA2）。

1.6 細胞亞群中EMRGs 的表達

通過“DoHeatmap”函數(shù)展現(xiàn) EMRGs 在GSE188280 數(shù)據(jù)集不同細胞類型的表達情況，用熱圖展示。然后使用“DotPlot”函數(shù)繪制6 個診斷標志物在不同細胞類型間的表達水平。最后使用“VlnPlot”函數(shù)繪制診斷標志物在不同細胞類型的表達情況的小提琴圖。

1.7 調(diào)控中藥預測及計算機分子對接、分子動力學模擬

Coremine Medical 數(shù)據(jù)庫（https://www.pubgene.com/coremine-medical/）是國際上最先進的醫(yī)學信息檢索平臺，可查詢基因所影響的生物學過程[19]。通過Coremine Medical 數(shù)據(jù)庫，以P＜0.05 為閾值，下載對所篩選的關(guān)鍵靶基因具有統(tǒng)計學意義的相關(guān)中藥。中藥系統(tǒng)藥理學與分析平臺 TCMSP（Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform ） 數(shù) 據(jù) 庫（https://www.tcmsp-e.com/）能篩選中藥化學成分，查找499 味中藥口服生物利用度、藥物相似度等藥動學特性[20]，通過TCMSP 查找關(guān)鍵中藥主要成分，從化學分子數(shù)據(jù)庫PubChem（https://pubchem.ncbi.nlm.nih.gov/）獲取中藥活性小分子成分配體的3D結(jié)構(gòu)。從蛋白質(zhì)數(shù)據(jù)庫UniPro（thttps://www.uniprot.org/）查詢基因蛋白Entry 號，RCSB PDB（Protein Data Bank）數(shù)據(jù)庫（https://www.rcsb.org/）中下載關(guān)鍵基因?qū)鞍资荏w3D 結(jié)構(gòu)。將中藥小分子配體及關(guān)鍵基因?qū)鞍资荏w運用AutoDock Tools-1.5.6 軟件進行分子對接及分子動力學模擬。

2 結(jié)果

2.1 Bulk RNA 測序篩選差異表達AMD 相關(guān)EMRGs 基因

研究比較GSE135092 數(shù)據(jù)集中AMD 組和對照組的黃斑區(qū)和非黃斑區(qū)組織EMRGs 的表達水平。組間差異采用Wilcoxcon 檢驗，如圖1-A 所示，在黃斑區(qū)發(fā)現(xiàn)了10 個EMRGs 的表達都具有顯著差異（P＜0.05）；如圖1-B 所示，在非黃斑區(qū)發(fā)現(xiàn)了9 個EMRGs 的表達具有顯著差異（P＜0.05）。如圖1-C所示，兩類組織差異表達基因中有6 個交集基因，其中叉頭蛋白O3（forkhead box O3，F(xiàn)OXO3）、絲氨酸蘇氨酸激酶3β（glycogen synthase kinase 3 beta，GSK3B）、肌細胞增強因子2A（myocyte enhancer factor 2A，MEF2A）、核受體共激活因子1（nuclear receptor coactivator 1，NCOA1）在AMD 組高表達，組蛋白脫乙酰酶1（histone deacetylase 1，HDAC1）、過氧化物酶體增生激活受體 γ（peroxisome proliferator-activated receptor gamma，PPARG）在AMD 組低表達。如圖1-D 所示，基因表達水平顏色從淺到深表示表達水平升高，藍色表示呈負相關(guān)，紅色表示呈正相關(guān)。EMRGs 表達水平的相關(guān)性矩陣如圖1-E 所示，F(xiàn)OXO3與NCOA1、GSK3B、MEF2A呈顯著正相關(guān)，MEF2A與NCOA1、GSK3B呈顯著正相關(guān)，GSK3B和NCOA1以及PPARG和HDAC1分別呈顯著的正相關(guān)；而FOXO3和HDAC1呈顯著的負相關(guān)；紅色表示正相關(guān)，藍色表示負相關(guān)，顏色越深表示相關(guān)程度越高。對差異表達EMRGs 之間的相關(guān)性進行分析，相關(guān)系數(shù)大于0.65 的相關(guān)性結(jié)果的散點圖如圖2-A～E 所示。

圖1 GSE135092 數(shù)據(jù)集EMRGs 的表達情況及相關(guān)性分析Fig.1 Expression and correlation analysis of EMRGs in GSE135092 dataset

圖2 差異表達EMRGs 相關(guān)性分析Fig.2 Correlation analysis of differentially expressed EMRGs

2.2 差異表達AMD 相關(guān)EMRGs 的PPI、TF、miRNA 網(wǎng)絡(luò)

將6 個差異表達AMD 相關(guān)的EMRGs 進行PPI網(wǎng)絡(luò)分析，PPI 網(wǎng)絡(luò)如圖3-A 所示，紅色表示上調(diào)，藍色表示下調(diào)，節(jié)點越大表示和其他基因聯(lián)系越多。mRNA-TF 互作結(jié)果顯示，F(xiàn)OXO3與85 種TF、GSK3B與31 種TF、HDAC1與33 種TF、MEF2A與2 種TF、NCOA1與3 種TF、PPARG與4 種TF存在互作關(guān)系，如圖3-B 所示。mRNA-miRNA 分析發(fā)現(xiàn)，F(xiàn)OXO3與64 種miRNA、GSK3B與68 種miRNA、HDAC1與10 種miRNA、MEF2A與8 種miRNA、NCOA1與17 種miRNA、PPARG與11 種miRNA 存在互作關(guān)系，如圖3-C 所示。

圖3 差異表達EMRGs 互作網(wǎng)絡(luò)分析Fig.3 Interaction network analysis of differentially expressed EMRGs

2.3 差異表達AMD 相關(guān)EMRGs 診斷AMD 能力評估

使用Logistic 回歸分析差異表達AMD 相關(guān)EMRGs 與是否患有AMD 的關(guān)聯(lián)性，發(fā)現(xiàn)NCOA1、MEF2A、FOXO3、GSK3B表達分別與AMD 發(fā)生呈正相關(guān)，HDAC1表達與AMD 發(fā)生呈負相關(guān)，如圖4-A 所示。根據(jù)Logistic 模型擬合結(jié)果分別繪制了ROC 曲線，AUC 分別為FOXO30.720、GSK3B0.654、HDAC10.695、MEF2A0.647、NCOA10.655、PPARG0.577，如圖4-B～G。

圖4 差異表達EMRGs Logistic 回歸分析與ROC 曲線Fig.4 Logistic regression analysis and ROC curve of differentially expressed EMRGs

2.4 LASSO 分析差異表達EMRGs 篩選AMD 診斷標志物及風險評分

對獲得的差異表達EMRGs 進行LASSO 回歸分析，并利用具有10 倍交叉驗證的部分似然偏差來計算最佳Lambda 值（λ），其中紅點代表每個λ的目標參量，虛線與橫坐標相交處為最優(yōu)λ值。如圖5-A 所示，當均方誤差取最小值時，預測模型描述實驗數(shù)據(jù)精確度越高，對應的λ值即可確定LASSO 回歸復雜度調(diào)整的程度。如圖5-B 所示，6 個差異EMRGs 都被篩選出來作為診斷標志物。根據(jù)LASSO 回歸模型的分析結(jié)果確定候選預后標志物的系數(shù)，計算風險評分，風險評分＝?1.821HDAC1－0.144NCOA1＋1.185MEF2A＋2.243FOXO3－1.662GSK3B＋0.229PPARG。分別計算AMD 和對照組樣本的風險評分，組間差異采用Wilcoxcon 檢驗，結(jié)果表明AMD 組的風險評分顯著高于對照組，如圖5-C 所示。采用Logistic 回歸模型構(gòu)建診斷模型，模型擬合結(jié)果顯示ROC AUC 為0.725，如圖5-D 所示。

圖5 LASSO 回歸篩選診斷AMD 相關(guān)EMRGsFig.5 LASSO regression screening for diagnostic AMD-related EMRGs

2.5 Single-cell RNA 測序分析RPE 組織的細胞異質(zhì)性

本研究對單細胞GSE188280 數(shù)據(jù)集進行質(zhì)控，過濾線粒體基因含量＞20%，feature＜500 或＞6 000的細胞后獲得了105 424 個細胞，提取來源RPE 組織的細胞共計32 714個，以供下游分析。應用UMAP降維進行可視化，將32 714 個細胞成功分類為19個獨立的簇，如圖6-A、B 所示，顏色越深表示表達水平越高，圓圈越大表示基因在細胞群內(nèi)的表達比例越高。使用marker 基因?qū)⒓毎罔b定為共9 種細胞類型，如圖6-C 所示，cluster 1、5 注釋為T 細胞（7 481，22.92%）；cluster 8、16 注釋為B 細胞（1 884，5.77%）；cluster 11 注釋為內(nèi)皮細胞（593，1.81%）；cluster 0、9、10 注釋為成纖維細胞（8 773，26.88%）；cluster 7、12、15 注釋為巨噬細胞（2 433，7.45%）；cluster 6 注釋為黑色素細胞（2 049，6.27%）；cluster 3、13 注釋為NK 細胞（3 016，9.24%）；cluster 2、14 注釋為施萬細胞（3 766，11.54%）；cluster 4注釋為平滑肌細胞（2 634，8.07%）。每個樣本間的各細胞比例如圖6-D 所示。

圖6 單細胞RNA 測序數(shù)據(jù)集GSE188280 鑒定RPE 組織細胞亞群Fig.6 Identification of RPE cell subsets using single-cell RNA sequencing dataset GSE188280

2.6 RPE 組織細胞亞群中EMRGs 的表達

利用熱圖展現(xiàn)48 個EMRGs 在各類型細胞中的表達情況，如圖7-A 所示。其中6 個差異EMRGs 在各類型細胞中的表達如圖7-B 所示，顏色越深表示表達水平越高，圓圈越大表示基因在細胞群內(nèi)的表達比例越高；6 個差異EMRGs 在各類型細胞中的表達小提琴圖如圖7-C 所示。HDAC1、MEF2A在內(nèi)皮細胞和B 細胞中存在較高表達，F(xiàn)OXO3在內(nèi)皮細胞、成纖維細胞、施萬細胞、平滑肌細胞、巨噬細胞中都存在較高表達，NCOA1在內(nèi)皮細胞中存在較高表達。

圖7 差異表達EMRGs 在單細胞RNA 測序數(shù)據(jù)集GSE188280 中的表達水平Fig.7 Expression level of differentially expressed EMRGs in single-cell RNA sequencing dataset GSE188280

2.7 潛在中藥化合物篩選與計算機分子對接、分子動力學模擬

從Coremine Medical 數(shù)據(jù)庫，以P＜0.05 為閾值，篩選出與關(guān)鍵靶基因具有統(tǒng)計學意義的相關(guān)中藥黃芪、人參、山藥、山茱萸、桑白皮、桑葉、淫羊藿、西紅花和丹參。從TCMSP 數(shù)據(jù)庫根據(jù)中藥化學物質(zhì)口服生物利用度（oral bioavailability，OB）、類藥性（drug-likeness，DL）篩選關(guān)鍵中藥主要活性物質(zhì)285種，選擇生物活性最高的9 個活性物質(zhì)（丹參醇A、丹參醇B、薯蕷皂素、富馬堿、馬卡因、白樺脂酸、馬爾肯久納醇苯酯、正庚醛、淫羊藿A）及對應的關(guān)鍵基因靶點蛋白（FOXO3、HDAC1、MEF2A）進行分子對接。所有的對接實驗運用AutoDock Tools-1.5.6軟件，并根據(jù)結(jié)合自由能評估最終對接結(jié)構(gòu)，最佳前4 組利用Pymol 2.1 軟件進行可視化，如圖8-A～D 所示。其中，丹參醇B 能夠與蛋白HDAC1 的GLY-149氨基酸形成氫鍵距離短（0.28 nm），結(jié)合能力強；還能與蛋白口袋的PHE-205、HIS-141、PHE-150、PRO-206、LEU-271、TYR-303 形成很好的疏水作用，特別是化合物的苯環(huán)能夠與PHE-205、HIS-141 形成p-p 共軛作用，對穩(wěn)定小分子有著重要貢獻，如圖8-D 所示。

圖8 分子對接可視化結(jié)果Fig.8 Visualization results of molecular docking

借助Gromacs 2020 軟件包對分子對接最佳的丹參醇B 與HDAC1 蛋白進行分子動力學模擬，使用Gromacs 2020 程序的g_MMPBSA 方法對結(jié)合自由能進行計算。結(jié)果顯示，復合物的平均方根偏差（RMSD）均小于0.2 nm，并且基本在60 ns 左右達到動態(tài)平衡，表明小分子與靶點蛋白匹配較好能形成穩(wěn)定的復合物，如圖9-A 所示。由平均方根波動（RMSF）結(jié)果可知，蛋白與小分子作用形成的復合物中少部分氨基酸構(gòu)象變化較大，大部分變化較小，說明小分子與蛋白相互作用形成穩(wěn)定的復合物，如圖9-B 所示。HDAC1 蛋白的回轉(zhuǎn)半徑（Rg）在前20 ns 略有升高，之后緩慢下降，蛋白與化合物的結(jié)合促使蛋白維持了更多的疏水接觸，內(nèi)部形成更多有效相互作用，而此時蛋白空腔擴大以更好地匹配化合物，從而促進復合物的穩(wěn)定性，如圖9-C 所示。統(tǒng)計蛋白與小分子的氫鍵數(shù)量發(fā)現(xiàn)，化合物均能夠與蛋白口袋氨基酸形成一個氫鍵相互作用，這些氫鍵在穩(wěn)定的小分子與蛋白結(jié)合方面起到了重要作用，如圖9-D 所示。丹參醇B 與HDAC1 蛋白的結(jié)合自由能為（?75.780±11.918）kJ/mol，范德華力發(fā)揮主要作用［（?149.094±9.569）kJ/mol］。氫鍵使靜電相互作用在穩(wěn)定小分子方面也存在貢獻［（?23.388±6.477）kJ/mol］。綜上，丹參醇B 能夠與HDAC1 蛋白有著很強的親和力，可以促使小分子與蛋白形成為穩(wěn)定的復合物，從而發(fā)揮活性作用。

圖9 分子動力學模擬可視化結(jié)果Fig.9 Visualization results of molecular dynamics simulation

3 討論

RPE 是視網(wǎng)膜的基本組成部分，RPE 細胞具有高度的新陳代謝活性，其依賴能量代謝來完成運輸營養(yǎng)、吞噬光感受器脫落的外節(jié)盤、消除活性氧、分泌因子等功能，對維持視網(wǎng)膜的健康至關(guān)重要，其功能異常將引起AMD 等視網(wǎng)膜退行性疾病[21-22]。1 項研究比較從流產(chǎn)胎兒或老年供體眼球中分離出來的RPE 細胞的能量代謝，發(fā)現(xiàn)線粒體耗氧量（ATP 生產(chǎn)）在老化細胞中完全喪失儲備能力，而糖酵解在老化細胞中也被削弱，說明RPE細胞的能量代謝隨年齡增加而減弱[23]。煙酰胺腺嘌呤二核苷酸（nicotinamide adenine dinucleotide，NAD+）是一種中心代謝輔助因子，在調(diào)節(jié)細胞代謝和能量平衡方面起關(guān)鍵作用，有研究評估小鼠RPE中關(guān)鍵的NAD+生物合成基因的時間表達和相關(guān)的NAD+水平，發(fā)現(xiàn)NAD+水平隨著年齡的增長而下降，也說明能量對RPE 健康的重要性[24]。與健康對照組的RPE 細胞相比，從AMD 患者中分離出來的RPE 細胞的能量代謝有障礙，該研究指出衰老加上外部壓力擾亂RPE 細胞的代謝是AMD 重要病理特征[25]。總之，這些研究支持增強能量代謝對延緩衰老和保持RPE 活力有潛在影響，或許是治療AMD的有效方式。

高通量技術(shù)有低成本、測序快、通量高等優(yōu)勢，是新一代測序技術(shù)，其中轉(zhuǎn)錄組測序?qū)膊》肿訖C制以及藥物靶點的預測有一定幫助。本研究借助轉(zhuǎn)錄組測序數(shù)據(jù)，比較了AMD 患者和非AMD 患者的mRNA 的表達譜，并使用相關(guān)計算方法篩選，并進行分子對接發(fā)現(xiàn)FOXO3和HDAC1是與AMD 有較大關(guān)聯(lián)的關(guān)鍵基因。研究發(fā)現(xiàn)FOXOs 轉(zhuǎn)錄因子在細胞能量代謝、線粒體自噬等方面發(fā)揮著重要的調(diào)控作用，F(xiàn)OXO3 可調(diào)節(jié)葡萄糖代謝中的炎癥標志物，F(xiàn)OXO3a 影響糖酵解、三羧酸循環(huán)中關(guān)鍵酶的表達[26-27]。FOXOs 主要通過能量穩(wěn)態(tài)的中樞調(diào)節(jié)因子（adenosine monophosphate，AMP）依賴的蛋白激酶（adenosine 5′-monophosphate-activated protein kinase，AMPK）參與能量代謝，當能量缺失后細胞內(nèi)的高AMP/ATP 比率會導致AMPK 激活，而活化的AMPK 可以磷酸化FOXO3 的6 個不同氨基酸殘基，導致抗氧化反應和能量利用通路的相關(guān)基因表達上調(diào)[28]。查閱FOXO3 相關(guān)文獻，未見與AMD 相關(guān)報道，因此FOXO3 或許是作用于AMD 的新靶點。

HDAC 家族成員參與血管生成、炎癥等生物過程，在AMD 的病理生理學中發(fā)揮著重要作用[29]。在內(nèi)皮細胞中，HDAC1 通過調(diào)節(jié)主要的內(nèi)皮細胞功能，如血管生成、炎癥信號、氧化還原平衡和一氧化氮信號，介導外部和環(huán)境刺激的影響。血管生成最常被內(nèi)皮細胞HDAC1 抑制，HDAC1 通過組蛋白和非組蛋白的去乙?；瘜?nèi)皮功能的依賴性調(diào)節(jié)發(fā)揮作用[30]。HDAC1 表達減少血管新生增多，脈絡(luò)膜新生血管是濕性AMD 主要病理特征。相比于正常人，HDAC1、2、5 和6 在AMD 患者的視網(wǎng)膜中的表達減少[31]。結(jié)合mRNA-TF、mRNA-miRNA調(diào)控網(wǎng)絡(luò)可進一步去探討AMD 發(fā)病機制。

本研究針對關(guān)鍵基因篩選治療AMD 的潛在中藥化合物，得到黃芪、人參、山藥、山茱萸、桑白皮、桑葉、淫羊藿、西紅花和丹參共9 味中藥，主要歸脾、肝、腎經(jīng)，有補氣健脾、清潤明目、補陽活血之功。中藥調(diào)控能量代謝主要通過葡萄糖代謝信號通路、線粒體能量代謝等途徑來緩解病變，并且發(fā)現(xiàn)多數(shù)中藥與脾、肝、腎關(guān)系密切[32]。脾為氣血化生之源，氣血是人體活動的物質(zhì)基礎(chǔ)，黃芪、人參、山藥益氣健脾，有助于水谷精微物質(zhì)被輸布散播至全身及目竅，發(fā)揮濡養(yǎng)作用?，F(xiàn)代代謝組學研究指出，黃芪、人參可能通過轉(zhuǎn)化生長因子β1（transforming growth factor-β1，TGF-β1）/Notch 信號，改善線粒體能量信號，對線粒體動力學相關(guān)蛋白發(fā)揮作用[33]。肝開竅于目，肝氣調(diào)達氣血和順，則目竅得養(yǎng)，肝氣不舒，氣郁化火，肝火上炎，耗傷精血，則目竅失養(yǎng)。桑葉清肝明目，研究發(fā)現(xiàn)桑葉參與單磷酸腺苷激活的AMPK 能量代謝信號通路、胰島素調(diào)控相關(guān)的信號通路等[34]。腎為命門之火，腎之陽又稱元陽、真陽，溫養(yǎng)五臟六腑官竅，腎陽虛衰，則臟腑官竅失于溫養(yǎng)。淫羊藿溫陽補腎，其發(fā)揮溫補效應機制與淫羊藿苷調(diào)控PTEN 誘導激酶1（PTEN induced kinase 1，PINK1）/parkin 信號通路降低細胞線粒體自噬水平、改善能量代謝有關(guān)[35]。丹參活血通經(jīng)，可以提高線粒體膜的完整性，降低其通透性，抑制大量鈣離子涌入線粒體，保護參與代謝的酶類，在缺血缺氧時，促進三羧酸循環(huán)，維護呼吸鏈功能，增加線粒體ATP 酶活性[36]。以上研究均支持所預測的健脾益氣、清肝明目、溫陽補腎、活血通經(jīng)類中藥可改善機體能量代謝，從而防治疾病。結(jié)合分子動力學模擬證明丹參有效成分丹參酮B 與其對應基因轉(zhuǎn)錄蛋白HDAC1 有著很強的親和力，可促使小分子與蛋白形成為穩(wěn)定的復合物而發(fā)揮活性作用，故丹參通過能量代謝過程對改善AMD 具有潛在作用。因此，對這些中藥化合物的關(guān)注可能會給AMD 的防治帶來新的研究方向。

綜上所述，本研究通過Bulk RNA 和Single-cell RNA 測序數(shù)據(jù)探究AMD 的能量代謝相關(guān)關(guān)鍵基因，并通過關(guān)鍵基因分析得到潛在治療中藥，借助分子對接技術(shù)及分子動力學模擬技術(shù)選出丹參酮B是最有潛在治療價值的中藥化合物，為AMD 后續(xù)藥物開發(fā)和臨床應用提供一定的研究基礎(chǔ)。但本研究仍然存在不足，單細胞測序較貴，尚未普及，AMD相關(guān)基因芯片數(shù)量較少，除此之外本研究還需進一步體內(nèi)或體外實驗驗證關(guān)鍵基因及相關(guān)中藥的治療效果。

利益沖突所有作者均聲明不存在利益沖突