結合Bulk RNA 和Single-cell RNA 測序探索年齡相關性黃斑變性能量代謝關鍵基因及調控中藥化合物預測

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

1.石河子大學醫學院，新疆 石河子 832000

2.中國中醫科學院眼科醫院，北京 100040

3.北京中醫藥大學東方醫院，北京 100078

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

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

1 資料與方法

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

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

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

1.2 Bulk RNA 測序中AMD 相關的能量代謝相關基因（energy metabolism-related genes，EMRGs）的蛋白相互作用（protein-protein interaction networks，PPI）、轉錄因子（transfer factor，TF）、微小RNA（microRNA，miRNA）網絡構建

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

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

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

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

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

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

Coefficient 為LASSO 回歸系數，mRNA Expressiopn 為該基因的表達水平（log2轉換）

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

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

通過文獻中已報道在RPE 和各種脈絡膜細胞類型中表達的基因列表以及人工注釋進行本研究的細胞類型注釋，例如，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 區III 受體γ（Fc gamma RIIIA，FCGR3A）、CD14，內皮細胞的標記基因有血小板內皮細胞粘附分子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，FBLN1），施萬細胞的標記基因有髓鞘質蛋白1（proteolipid protein 1，PLP1），黑色素細胞的標記基因有黑色素A（melan-A，MLANA）、黑素細胞蛋白（premelanosome protein，PMEL）、酪氨酸酶相關蛋白1（tyrosinase related protein 1，TYRP1），巨噬細胞的標記基因有同種異體移植炎癥因子1（allograft inflammatory factor 1，AIF1），肥大細胞的標記基因有原癌基因KIT（KIT proto-oncogene，KIT），平滑肌細胞的標記基因有蛋白信號調節因子5（regulators of G-protein signaling 5，RGS5）、α 肌動蛋白2（actin alpha 2，ACTA2）。

1.6 細胞亞群中EMRGs 的表達

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

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

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

2 結果

2.1 Bulk RNA 測序篩選差異表達AMD 相關EMRGs 基因

研究比較GSE135092 數據集中AMD 組和對照組的黃斑區和非黃斑區組織EMRGs 的表達水平。組間差異采用Wilcoxcon 檢驗，如圖1-A 所示，在黃斑區發現了10 個EMRGs 的表達都具有顯著差異（P＜0.05）；如圖1-B 所示，在非黃斑區發現了9 個EMRGs 的表達具有顯著差異（P＜0.05）。如圖1-C所示，兩類組織差異表達基因中有6 個交集基因，其中叉頭蛋白O3（forkhead box O3，FOXO3）、絲氨酸蘇氨酸激酶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 所示，基因表達水平顏色從淺到深表示表達水平升高，藍色表示呈負相關，紅色表示呈正相關。EMRGs 表達水平的相關性矩陣如圖1-E 所示，FOXO3與NCOA1、GSK3B、MEF2A呈顯著正相關，MEF2A與NCOA1、GSK3B呈顯著正相關，GSK3B和NCOA1以及PPARG和HDAC1分別呈顯著的正相關；而FOXO3和HDAC1呈顯著的負相關；紅色表示正相關，藍色表示負相關，顏色越深表示相關程度越高。對差異表達EMRGs 之間的相關性進行分析，相關系數大于0.65 的相關性結果的散點圖如圖2-A～E 所示。

圖1 GSE135092 數據集EMRGs 的表達情況及相關性分析Fig.1 Expression and correlation analysis of EMRGs in GSE135092 dataset

圖2 差異表達EMRGs 相關性分析Fig.2 Correlation analysis of differentially expressed EMRGs

2.2 差異表達AMD 相關EMRGs 的PPI、TF、miRNA 網絡

將6 個差異表達AMD 相關的EMRGs 進行PPI網絡分析，PPI 網絡如圖3-A 所示，紅色表示上調，藍色表示下調，節點越大表示和其他基因聯系越多。mRNA-TF 互作結果顯示，FOXO3與85 種TF、GSK3B與31 種TF、HDAC1與33 種TF、MEF2A與2 種TF、NCOA1與3 種TF、PPARG與4 種TF存在互作關系，如圖3-B 所示。mRNA-miRNA 分析發現，FOXO3與64 種miRNA、GSK3B與68 種miRNA、HDAC1與10 種miRNA、MEF2A與8 種miRNA、NCOA1與17 種miRNA、PPARG與11 種miRNA 存在互作關系，如圖3-C 所示。

圖3 差異表達EMRGs 互作網絡分析Fig.3 Interaction network analysis of differentially expressed EMRGs

2.3 差異表達AMD 相關EMRGs 診斷AMD 能力評估

使用Logistic 回歸分析差異表達AMD 相關EMRGs 與是否患有AMD 的關聯性，發現NCOA1、MEF2A、FOXO3、GSK3B表達分別與AMD 發生呈正相關，HDAC1表達與AMD 發生呈負相關，如圖4-A 所示。根據Logistic 模型擬合結果分別繪制了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 值（λ），其中紅點代表每個λ的目標參量，虛線與橫坐標相交處為最優λ值。如圖5-A 所示，當均方誤差取最小值時，預測模型描述實驗數據精確度越高，對應的λ值即可確定LASSO 回歸復雜度調整的程度。如圖5-B 所示，6 個差異EMRGs 都被篩選出來作為診斷標志物。根據LASSO 回歸模型的分析結果確定候選預后標志物的系數，計算風險評分，風險評分＝?1.821HDAC1－0.144NCOA1＋1.185MEF2A＋2.243FOXO3－1.662GSK3B＋0.229PPARG。分別計算AMD 和對照組樣本的風險評分，組間差異采用Wilcoxcon 檢驗，結果表明AMD 組的風險評分顯著高于對照組，如圖5-C 所示。采用Logistic 回歸模型構建診斷模型，模型擬合結果顯示ROC AUC 為0.725，如圖5-D 所示。

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

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

本研究對單細胞GSE188280 數據集進行質控，過濾線粒體基因含量＞20%，feature＜500 或＞6 000的細胞后獲得了105 424 個細胞，提取來源RPE 組織的細胞共計32 714個，以供下游分析。應用UMAP降維進行可視化，將32 714 個細胞成功分類為19個獨立的簇，如圖6-A、B 所示，顏色越深表示表達水平越高，圓圈越大表示基因在細胞群內的表達比例越高。使用marker 基因將細胞簇鑒定為共9 種細胞類型，如圖6-C 所示，cluster 1、5 注釋為T 細胞（7 481，22.92%）；cluster 8、16 注釋為B 細胞（1 884，5.77%）；cluster 11 注釋為內皮細胞（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 測序數據集GSE188280 鑒定RPE 組織細胞亞群Fig.6 Identification of RPE cell subsets using single-cell RNA sequencing dataset GSE188280

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

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

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

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

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

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

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

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

3 討論

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

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

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

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

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

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