微生物細胞工廠的智能設計進展

張震，曾雪城，秦磊，李春，2

（1清華大學化學工程系生物化工研究所/工業生物催化教育部重點實驗室，北京 100084；2清華大學合成與系統生物學研究中心，北京 100084）

引 言

“碳達峰、碳中和”作為我國的戰略發展規劃，對社會各行各業的發展將產生深遠影響，其中對于制造業而言，生物制造基于生物催化過程可生產各種高附加值產品，且以可再生的生物質為原料，是未來制造業可持續發展的重要方向。微生物細胞工廠通過在微生物底盤細胞中表達異源生物合成途徑來生產目標化合物，是綠色生物制造的重要實現形式。當下結構已報道的天然產物超過了30萬種，目前利用微生物細胞工廠實現了生物合成途徑已知的青蒿酸[1]、1,3-丙二醇[2]、法尼烯[3]、甘草次酸[4]、PHA[5]、1,3-丁二醇[6]等多種高附加值化合物的生產，產品廣泛應用于化工、醫藥、能源、食品等領域。然而，生物合成途徑完全解析的天然產物只有不到3萬種，尚有大量天然產物的生物合成途徑未解析，嚴重阻礙了微生物細胞工廠的發展[7]。傳統微生物細胞工廠的設計和構建方式通量小、效率低，亟需自動化的設計工具代替傳統依賴經驗和試錯的設計構建方式，以加速未知途徑化合物生物合成途徑的解析及其微生物細胞工廠的合成設計。

設計-構建-測試-學習(design-build-testlearn,DBTL)循環是開發微生物細胞工廠的基本研究思路，使用遞歸循環的方式使設計的微生物細胞工廠逐漸提高得率、速率及產量等生產目標。設計(design)是DBTL循環的第一步，對細胞工廠能否成功構建產生重要影響。基于有機化學和生物化學的專業知識、文獻報道及自身實踐經驗繪制潛在的生物合成途徑的設計方法通常僅適用于化學結構簡單的目標化合物。此外，當前微生物細胞工廠中所設計途徑與合成調控元件的組合較為盲目或隨機，只能逐次對元件進行優化，而且設計出的途徑不能很好與底盤細胞適配。這些因素導致傳統微生物細胞工廠的設計和構建通量小、效率低，且驗證過程周期長，例如，阿米瑞斯生物技術公司（Amyris Biotechnologies）構建微生物細胞工廠生產抗瘧疾前體青蒿素花費了十年時間、150人年的工作量，實現微生物生產法尼烯花費了四年時間、130～575人年；杜邦公司（DuPont）和杰能科公司（Genencor）分別花費約15和575來實現微生物細胞工廠生產1,3-丙二醇[8]。然而，組學時代的到來和生物數據資源的爆炸式增長已經改變了生物制造領域的研究模式，為生物制造和微生物細胞工廠設計提供了新的發展機遇。因此，在大數據基礎上的智能化微生物細胞工廠設計方法成為加速微生物細胞工廠設計構建的關鍵。

借鑒電子設計自動化（electronic design automation，EDA，指利用計算機輔助設計軟件，來完成超大規模集成電路芯片的功能設計、驗證、物理設計等流程的設計方式）的概念，微生物細胞工廠的智能設計可采用生物設計自動化（biological design of automation，BDA）的形式，立足于生物數據庫中的豐富資源，通過一系列算法完成細胞工廠的自動化設計，實現生物合成途徑的預測與篩選、調控元件的設計、途徑與元件組裝設計、設計途徑與底盤代謝網絡適配等功能。通過生物合成途徑的預測與篩選，對途徑未知化合物提供可靠性較強、效率較高的候選途徑方案，輔助途徑設計，此外，從全局出發預測潛在的分支途徑，以實現對實驗中可能出現的副產物進行快速定位，進而指導細胞工廠的構建。通過調控元件的設計及元件與途徑的組合設計，進一步提高預測途徑與底盤菌株的適配性及可靠性，為細胞工廠后續的構建、測試和學習指明方向。

1 逆合成算法預測生物合成途徑

化合物生物合成途徑的設計主要采用生物逆合成算法，其思想來源于化學逆合成，分為生物逆合成途徑的預測和途徑篩選兩個環節。在逆合成途徑預測階段，通過使用一組在原子水平上描述底物和產物分子之間化學轉化模式的生化反應規則，推測合成目標化合物的反應及催化該步反應的酶，實現將輸入化合物（即目標化合物）轉化為一系列中間化合物，并最終轉化為前體化合物的過程。生物逆合成算法按預測中間化合物的方式不同可分為兩大類：一是在數據庫中檢索已知代謝反應并預測反應和中間化合物，所預測的中間化合物種類受到化合物數據庫規模的限制，如FMM、DESHARKY和Metabolic tinker等工具；二是基于泛化的生化反應規則來預測新反應且可產生數據庫中不存在的新化合物，如XTMS[9]、RetroPath[10]、RetroPath2.0[11]、RetroPath RL[12]、novoPathFinder[13]等工具（表1），這一類基于反應規則的算法在生物逆合成預測中具有更大的應用潛力，因此下述生物逆合成預測流程主要介紹基于反應規則的預測算法。

表1 生物逆合成工具Table1 Retrobiosynthesis tools

1.1 生物逆合成途徑預測流程

生物合成途徑由多步生化反應級聯組合而成，后一步反應的前體化合物作為前一步反應的目標產物，相當于多個特殊的單步反應的組合，因此生物合成途徑的預測關鍵在于對某一指定化合物的單步生物合成反應的預測。單步生物合成反應的預測以普遍酶具有底物雜泛性的假設條件為前提，即認為酶可以催化與底物具有相似化學結構的化合物，因此可借鑒已報道生化反應的轉化模式及相應的酶序列來預測與底物具有相似結構化合物發生的新生化反應。同時存在一些計算工具可對酶的底物雜泛性進行預測，如EPP-HMCNF[23]可根據BRENDA數據庫中酶-底物的互作預測可催化給定查詢化合物的不同種類的酶。基于上述假設單步反應的預測流程可以概括為，對于目標化合物Q(query)，在化合物數據庫中檢索與之結構相似的化合物M1、M2(match)，在生化反應數據庫中檢索M1、M2的相關反應R1、R2(reaction)，提取相應的反應規則及酶序列，將其應用在目標化合物Q上，從而得到目標化合物Q的生成反應、催化該步反應的酶序列及直接前體化合物Q11、Q21（圖1，步驟1～5）。

通過對預測前體迭代應用單步預測反應，可逐步延伸預測途徑的長度，直至達到規定的終止條件。通常以預測途徑達到規定步長、預測的前體化合物為指定化合物或指定底盤宿主內源化合物等易于計算的指標作為終止條件。由于化合物相關合成反應并不單一，因此預測出的合成反應通常具有許多分支，從而使得最終的逆合成途徑通常以樹的形式呈現，被稱為逆合成網絡（圖1，步驟6～7），其中目標化合物為根節點、中間化合物為子節點、指定起始化合物或底盤宿主可內源合成的化合物為葉子結點，生化反應及催化該步反應的酶序列為連接節點的邊。基于逆合成網絡的樹形特點，一些同樣呈現為樹形的人工智能算法被應用到逆合成網絡的生成中，如RetroPath RL[12]利用蒙特卡洛樹搜索算法（Monte Carlo tree search）的選擇-擴展-隨機模擬-反向傳播四個過程對預測途徑進行延伸，最終生成逆合成網絡。

圖1 生物逆合成途徑預測流程圖Fig.1 Workflow of retrobiosynthesis

1.2 生化數據庫

在生物合成途徑的預測中生化數據庫為預測算法的實現奠定了基礎，所需要的生化數據庫包括化合物數據庫、生化反應數據庫、代謝數據庫、酶數據庫及細胞模型。KEGG[24]、KNApSAcK[25]等化合物數據庫及天然產物詞典(Dictionary of Natural Products,DNP)(https://dnp.chemnetbase.com)和Super NaturalⅡ[26]等天然產物數據庫可被用于相關反應已有文獻報道的相似化合物的搜索。生化反應數據庫用于提取反應規則，如ATLAS[27]數據庫中收錄了超過14萬條反應，ATLASx[22]數據庫中包含超過500萬個預測反應。基于代謝途徑中反應間的級聯關系可將反應規則構建為反應規則網絡，為長途徑的預測提供可能。常用的代謝途徑數據庫包括：KEGG PATHWAY、MetaCyc、Reactome和UM-BBD等。酶數據庫將提供催化生化反應可用的酶序列，如BRENDA是一個綜合酶信息數據庫，包含了酶促反應和相關的代謝通路。細胞模型為基因組規模代謝網絡模型（genome-scale metabolic model，GEM），GEM模型可提供底盤宿主中存在的化合物、生化反應等信息，并可對代謝途徑產量進行模擬計算、為預測途徑的優化提供手段。大腸桿菌、釀酒酵母等常用模式生物均發展了系列GEM模型。目前也創建了一些綜合數據庫，如biochem4j圖數據庫中包含36765條反應信息、19735條化合物信息、245704條酶序列和8431個細胞模型。

1.3 化合物描述符

生物合成途徑由一系列基于可進行酶催化轉化的化合物組成，化學結構式是這些化合物的基礎表示方式，可表示所有原子通過化學鍵與其相鄰原子連接的原子鍵合環境信息。在生物逆合成途徑預測中通常采用化學模式語言和分子指紋兩種方式對化合物的結構式進行數字化編碼，以滿足后續提取反應規則及搜索結構相似性化合物的需求。

1.3.1 化學模式語言 在化學信息學中，SMILES（simplified molecular-input line-entry system）和SMARTS（smiles arbitrary target specification）是兩種已有明確定義的化學模式語言，其中SMILES可以將化合物的二維結構式表示為ASCⅡ字符串，主要由原子和化學鍵兩種基本符號組成；SMARTS是SMILES的延伸，允許使用通配符表示原子和化學鍵，如SMARTS中符號[C,N]表示該原子是碳（C）或氮（N），符號～可匹配任何化學鍵[28-30][圖2(a)]，這兩種化學模式語言在生物逆合成算法中常被用于表示化學反應中反應物和產物的結構變化[9，11，13，15]。

圖2 化學模式語言與分子指紋Fig.2 Chemical model language and molecular fingerprints

1.3.2 分子指紋 通過預定義KEGG原子類型[31]、官能團等子結構將化合物結構分解為一個個子結構的累積，并將化合物中子結構的數量和各種物理化學特性編碼為二進制變量的位串，形成化合物的分子指紋。可以使用RDKit和PaDEL-descriptor等軟件包來生成分子指紋，RDkit[32]是一個用于化學信息學的開源工具包，可以生成RDKit指紋、Morgan指紋、Avalon指紋和MACCS指紋四種分子指紋，例如MACCS指紋可表示為長度為166位的向量，每一位對應一個分子特征，當化合物中存在此類特征時該位數值為1，否則為0[圖2(b)]。PaDEL-descriptor[33]可 生 成PubChem指 紋、CDK指紋、CDKextend指紋、子結構指紋和GraphOnly指紋五種分子指紋。這些指紋一方面可基于Tonimoto等相似性算法用于從化合物數據庫中快速搜索與查詢化合物具有相似結構的化合物，另一方面可用作相似性搜索、分類和回歸等各種機器學習任務的輸入。

1.4 反應規則

反應規則（也稱為反應描述符）為產物描述符與底物描述符之間的凈差，通過底物結構式與產物結構式間的原子-原子映射比對而得到，描述了底物向產物轉化時反應中心原子的鍵合環境變化（圖3）。目前對反應規則的研究多以反應中某一個主要的反應物向相應產物的一對一轉化模式為主。反應規則可作為一種模塊化操作，適用于與底物結構相似的化合物上，可預測目標產物的合成反應及相應前體化合物。反應規則可以從已知反應的數據庫中自動提取[34-35]，也可以通過手動輸入生產精簡的專家反應規則集[36]。

1.4.1 基于生化反應數據庫的反應規則 反應規則的自動提取需要經歷（1）反應物-產物對的識別；（2）原子-原子映射；（3）反應中心原子、反應基團及保守基團的識別和（4）提取反應規則四個過程。PathPred中將每個反應中反應對（反應物-產物對）中的匹配區和非匹配區之間的邊界原子分別定義為反應中心R原子、差異區域D原子和匹配區域M原子，引入了R-D-M原子模式。KEGG為反應對提供了基于原子映射的結構對齊信息，并構建了KEGG RPAIR數據庫。上述反應規則表示了反應物-產物間的最大結構差異，而XTMS[9]基于SMILES化學模式語言描述化學反應，同時將原子映射編號附加到反應物和產物側的相應原子上，以反映原子身份并跟蹤反應中原子的轉移；通過調整反應中心原子周圍環境的大小（即直徑d）可獲得不同泛化水平的反應規則，當d=0時僅包括反應中心原子，d=1時包括反應中心原子及與之直接相鄰的化學鍵及原子，如此類推，可見隨著直徑d的增加，反應規則變得更加具體（圖3）。該方法中反應規則的泛化水平是逆合成途徑預測的關鍵技術之一，過于具體的規則會限制預測新路線的潛力，而過于籠統的規則可能會使預測偏離實際[37]。RetroRules[34]反應規則數據庫按照SMARTS標準格式收錄了超過40萬條包含立體化學信息的反應規則，且每條反應規則均可以不同雜泛水平呈現。

圖3 反應規則的提取Fig.3 The extraction of reaction rules

1.4.2 專家反應規則集 自動提取的反應規則通常存在大量冗余，將降低途徑延伸時的預測計算效率，且使產生的化合物和反應的數量呈現指數增長，造成組合爆炸問題。此時可通過人工精簡產生規模較小但更精確的專家反應規則集，以限制途徑延伸時反應的數量且減少網絡規模，有利于提高途徑延伸的計算效率及反應規則覆蓋所有可能轉換的全面性。RetroBioCat[37]人工構建了由83個反應組成的反應規則集，使用107個反應的SMARTS編碼進行描述。Broadbelt[38]從基于原子映射的MetaCyc反應中自動提取反應規則，并通過人工精簡獲得最小但全面的1224條通用反應規則，經驗證可唯一地覆蓋所有常見的酶促轉化，且能夠重現KEGG和BRENDA數據庫中超過85%的所有反應，有利于探索已知酶促轉化的更大空間，加速生物合成途徑的設計。

1.4.3 反應規則網絡 代謝途徑是生物體內的級聯反應，基于代謝途徑數據庫可以將從數據庫中學習到的所有反應規則按照反應規則網絡(reaction rule network,RRN)的形式進行整合[16]。將各個反應規則均視為節點，若兩個規則在已知途徑上呈級聯狀態或具有形成級聯反應的潛力，則在兩個規則間添加邊來連接，最終形成反應規則網絡并將其應用在長途徑化合物的預測上。例如，ReactionMiner[16]基于反應規則網絡對衣康酸酯、柚皮素、1,3-丙二醇、木糖醇等高附加值化合物的生物逆合成途徑進行了預測，發現可復原這些化合物的已知途徑或預測出更短且生物學上更合理的逆合成途徑。

1.5 酶序列的獲取

化學催化和生物催化的主要區別在于生物催化采用酶作為催化劑，同樣的，生物逆合成區別于化學逆合成的關鍵在于生物逆合成需要為預測的反應提供可能的催化該步反應的酶序列。然而，許多預測的新反應沒有相關酶報道，預測與反應規則相關聯、可催化新反應且與底盤宿主進化親緣較近的酶的計算方法可大大加快生物合成途徑的開發。

EC號是國際生物化學與分子生物學聯盟(International Union of Biochemistry and Molecular Biology，IUBMB)中酶學委員會（Enzyme Commission）根據每種酶所催化的化學反應為分類基礎制定的一套由四個級別組成的編號分類法，而反應規則同樣體現了酶催化的功能，因此反應規則與EC號間存在關聯，并通過查詢EC號對應的酶序列可將反應規則與酶序列相關聯，為預測出的新反應提供酶催化數據。

反應同源性是預測催化目標反應候選酶序列的基礎，且不同酶序列預測工具在反應相似性矩陣的計算及反應的輸入格式上有所區別。Yamanishi等[39]采用KEGG化合物ID作為輸入內容，通過比較化合物間的轉化模式，構建了E-zyme工具，可對化學反應分配EC號前三位數字。Goto等[40]對E-zyme工具進行了拓展并設計了E-zyme2工具，可基于RDM模式對底物-產物對進行全結構比對，當已知的相似反應與至少一個基因序列相鏈接時，可為輸入的查詢反應分配催化相似反應的酶序列。此外，Thornton等[41]開發了EC-BLAST工具，采用KEGG反應ID、SMIRKS反應規則和EC號作為輸入內容，利用最大公共子圖算法（maximal common subgraph,MCS）進行原子-原子比對（atom-atom mapping,AAM），并根據鍵的變化、反應中心及結構相似性從KEGG REACTION數據庫中查詢相似反應，以此為新反應分配前三位EC號，但EC-BLAST不能輸出酶序列。Faulon等[42]開發了Selenzyme工具，采用SMIRKS化學模式語言作為輸入內容，基于biochem4j圖數據庫，通過綜合考慮序列相似性、酶序列的物理化學性質、酶來源物種與底盤宿主間的種群發生距離等方面來選擇最優的候選酶序列。

BridgIT[43]考慮了輔因子在反應機制中的關鍵作用，是目前唯一可以區分不同酶催化反應機制且可預測從頭設計反應的關聯酶序列的方法。BridgIT采用SMILES表示的化學反應及BNICE.ch反應規則作為輸入內容，可基于Daylight分子指紋計算目標反應與天然反應的Tanimoto相似性，并從KEGG中獲取相似天然反應及其對應的酶EC號，進而提供一個候選酶EC號名單。對于KEGG（2011）中未注釋酶序列且在KEGG（2018）中補充了注釋的反應，BridgIT對其中90%的反應能夠正確預測其相應的酶三級EC號。

基于以上工具，一些報道將酶數據庫與反應數據庫相結合，構建出反應規則-酶相互關聯的綜合數據庫，方便生物逆合成算法的計算，如EnzyMine[34]將酶催化反應的特征與酶序列和結構注釋相聯系，構建了包含7767個EC號、267345條蛋白序列和9831個反應的綜合酶-反應規則關聯數據庫。Fenner等[44]將從化學污染物生物轉化反應數據庫Eawag-BBD和KEGG數據庫中獲取的泛化生物轉化規則與前三位EC號相關聯，建立了具有316條反應規則-酶分類鏈接的enviLink數據庫。

2 途徑篩選指標與綜合評價

逆合成網絡中包含了大量預測途徑，但并非所有預測途徑都能夠實現目標催化功能，目前已開發的工具所預測的途徑假陽性過高，使得生物逆合成工具的應用仍不夠普及，因此提高生物逆合成途徑預測的準確率和可靠性是促進逆合成算法廣泛應用的關鍵。從大量預測途徑中推薦最佳候選途徑，需要基于一些評價指標對預測途徑進行評價、排序及篩選，目前已經報道的一些途徑評價指標，主要從途徑的理論可行性和與底盤宿主的適配性兩方面對途徑進行評價（圖1，步驟8）。

2.1 途徑的可行性

在途徑篩選中，首先需要采用如底物相似性、熱力學可行性、酶序列、途徑長度等定量指標，對一些理論上不可行的途徑進行排除。

2.1.1 底物相似性 考慮到生物逆合成算法采用酶-底物雜泛性假設，選擇與已知底物結構相似性較高的輸入化合物將更有可能被相應的已知酶催化；相反，若輸入的查詢化合物與化合物數據庫中已知化合物的結構相似度較低，則相關聯的已知酶對該查詢化合物成功催化的可能性同樣較低。基于化合物描述符，可以通過Tanimoto相似性等算法對預測的中間化合物與數據庫中已知化合物間的結構相似性系數進行計算，從而在化合物數據庫中檢索結構相似化合物[7，12]。

2.1.2 熱力學可行性 吉布斯自由能變化ΔG可表示反應的熱力學勢能變化，決定了酶促反應的方向性和效率，是檢測和選擇預測途徑熱力學可行性及評估生物合成途徑熱力學驅動力的重要手段。一些生物逆合成預測工具基于數據庫中的反應吉布斯自由能數據或熱力學計算工具檢測途徑的熱力學可行性并對途徑進行篩選，如Metabolic tinker基于CHEBi和RHEA數據庫，使用此前報道的基于基團貢獻（group contribution,GC）的熱力學計算工具計算并評估途徑的可行性[45]；XTMS[9]基于MetaCyc database數據庫提供的反應吉布斯自由能數據來評價途徑可行性。此外，OptMDFpathway[46]利用基于約束的模型以途徑的最大最小驅動力（max-min driving force）為優化目標，通過混合整數線性規劃來識別具有最高熱力學驅動力的途徑，這類途徑具有較高的代謝通量且對酶的表達強度要求較小。eQuilibrator3.0[47]可利用組分貢獻（component contribution,CC）工具計算反應的生化平衡常數和ΔG，同時給出衡量預測不確定性的協方差矩陣用于基于約束的熱力學模型計算。dGPredictor[48]基于KEGG數據庫可利用基團貢獻工具計算不同pH和離子強度下的代謝途徑中酶催化反應的ΔG，且考慮了化合物結構中的立體化學信息，增加了熱力學預測的精準度。

2.1.3 酶序列可行性 有無催化目標反應的酶序列對于預測途徑的實現十分重要，尤其是對于基于反應規則預測的新反應。在延伸預測途徑時可能一些反應規則不存在相關聯的酶序列，此時預測出的新反應需要依靠人工查詢文獻以尋找酶序列，降低了設計的效率，因此在選擇反應規則時應當增加與酶序列相關聯的反應規則的權重，提高酶催化反應的可行性。

2.1.4 途徑長度 途徑長度是最直接的篩選指標，長途徑意味著在底盤宿主中引入了更多的酶，從而使代謝負擔增加，而結構復雜化合物需要經歷多種后修飾過程，需要較長的生物合成途徑才能完成修飾作用，因此需要對途徑長度進行合理篩選。目前一些算法直接利用途徑長度對預測途徑進行打分和排序，而基于圖論的生物逆合成預測工具通常采用混合整數線性規劃等基于約束計算方法尋找底物到目標產物的前k條最短基元模式分析（elementary flux mode，EFM）[49]、最 短 碳 流 量 途 徑（carbon flux path，CFP）[50]及最短活性途徑（active pathway）[51]。例如，NICEpath[52]將反應物-產物對中保守的原子數量作為反應權重結合k-最短圖搜索（k-shortest graph search）算法可實現KEGG中途徑的篩選。PATHcre8[53]構建了包含可逆反應的雙向圖，采用Yen算法與PathLinker算法相結合的前K條無環最短路徑算法篩選目標途徑。

2.2 途徑與底盤宿主的適配性

由于預測途徑在底盤宿主中的實現會受到內源化合物及調控網絡影響，從而呈現出偏離預測的現象，因此為了合理設計一個高效的異源生物合成細胞工廠，必須考慮外源反應在底盤宿主中特定內源性代謝網絡影響下的穩定性，需要對預測途徑與底盤宿主的適配性進行評價，以增加預測途徑的可行性。OptStrain[54]、DESHARKY[55]、FMM[49]、Metabolic tinker[45]、GEM-Path[56]、XTMS[9]、MRE[57]、RetroPath2.0[11]等生物逆合成預測工具均對途徑與大腸桿菌、酵母、藍細菌等底盤宿主的適配性進行了探討。

2.2.1 化合物毒性 中間化合物對細胞的毒性將妨礙途徑中的酶在底盤宿主中正常表達，因此需要對化合物毒性進行預測，避免預測途徑中包含高毒性中間化合物。通常采用化合物的半數抑制濃度(the half inhibitory concentration，IC50)作為化合物毒性的評價指標，表示一半細胞種群的生長受到抑制時的化合物濃度[58]。目前有一些收錄了化合物毒性的 數 據 庫，如TOXNET[59]、DSSTox[60]、T3DB[61]以 及RTECS[62]等，但其中的毒性數據多以動物細胞為對象，缺少對微生物細胞的毒性數據。目前已報道可預測化合物對微生物底盤宿主毒性的軟件較少，其中EcoliTox[63]可基于化學結構與活性的定量關系預測中間化合物在大腸桿菌中的毒性，在多種逆合成工具中均有應用[9,11-12,64-65]。此外Toxicity Estimation Software Tool(TEST)工具也可實現對化合物毒性的預測。

2.2.2 代謝負擔 代謝途徑中引入的異源途徑將增加細胞的代謝負擔，從而妨礙細胞的生長和生產。DESHARKY[55]使用蒙特卡洛啟發式算法對生物合成途徑進行預測，同時對大腸桿菌中細胞資源和內源代謝情況進行建模，基于對核糖體和RNA聚合酶的消耗量計算異源代謝途徑對底盤宿主的負擔，從而可選取對宿主產生較小代謝負擔的途徑。

2.2.3 理論產量 通過構建包含底盤菌株內源化合物及預測途徑相關化合物相對應的化學計量矩陣來構建GEM模型，并利用流量平衡分析（flux balance analysis,FBA）及基于約束的混合整數線性規劃算法可計算途徑在目標底盤宿主中化合物的理論產量，且常以產量最大化為優化目標來選擇途徑，這類生物逆合成預測工具包括OptStrain[54]、DESHARKY[55]、FMM[49]、GEM-Path[56]及XTMS[9]等。需要注意的是基于FBA的工具在評估途徑時需要提供豐富的信息來為給定底盤細胞模型設定嚴格的反應通量邊界，因而僅適用于大腸桿菌、釀酒酵母等經過充分研究的模式微生物。

2.2.4 內源起始化合物 底盤宿主內結構更簡單的起始化合物通常具有更大的代謝通量，有利于增加目標產物的產量。SCScore[66]是一項衡量分子復雜性的指標，基于大量合成化學反應訓練的神經網絡對化合物分子的結構復雜性進行評分，從而有利于指導生物逆合成途徑選擇更簡單的起始化合物。此外化合物官能團與反應中心碳原子之間的相對位置，將影響官能團的電子構型變化以及酶活性位點與底物的結合，進而影響化學反應的可行性。Lee等[67]基于ChemAxon Reactor工具通過比較所選氨基酸前體中官能團（氨基和羧基）與L-纈氨酸、L-亮氨酸和L-異亮氨酸的相對位置，為短鏈伯胺生物合成途徑選擇了最佳的氨基酸前體。

2.2.5 內源競爭途徑 底盤宿主內源代謝網絡可能會對異源生物合成途徑產生競爭作用，進而影響目標化合物的生產。MRE[57]將預測出的異源生物合成途徑整合到底盤宿主內源代謝網絡中，并基于熱力學可行性確定代謝網絡中反應的方向，形成有向圖，基于標準化的Boltzmann因子計算可對內源前體化合物（節點）進行轉化的內源和異源競爭反應（邊）的概率分布并由此進行賦權，以考慮特定底盤宿主中內源代謝反應對異源途徑的競爭作用。

2.3 途徑的綜合評價

在實際情況中，常使用不同指標的加權組合來對途徑進行綜合評價，且一些基于綜合評價獲得的預測途徑已被實驗驗證并用于微生物細胞工廠的構建。

XTMS[9]對食品工業中使用的昂貴風味成分樹莓酮的生物合成途徑進行了預測，并利用通量平衡分析(FBA)根據熱力學可行性（吉布斯自由能）、酶性能（基因評分）、途徑可行性（反應步驟的數量）、中間化合物的毒性及目標化合物產量等指標對預測途徑進行了綜合評價和排序，恢復了以香豆酰輔酶A為底物合成樹莓酮的兩步天然合成途徑，同時考慮到作為底盤宿主的大腸桿菌中不存在香豆酰輔酶A，XTMS給出了從內源性化合物到樹莓酮的生物合成途徑。

PATHcre8可選擇藍藻為底盤宿主[53]，預測從乙酰乙酰輔酶A到IPP、從甲羥戊酸到異戊二烯的生物合成途徑及從可卡因到偽雌二醇輔酶A（pseudoecgonyl-CoA）的生物降解途徑，并根據反應熱力學可行性、途徑中的潛在毒性產物（化合物毒性）、競爭反應消耗的途徑中的中間產物（產物消耗）以及拷貝數等指標對途徑進行綜合評價，結果顯示從乙酰乙酰輔酶A到IPP的天然途徑在預測的前15條候選途徑之中，經實驗驗證的基于磷酸異戊烯酯的(R)-甲羥戊酸到異戊二烯途徑在候選途徑中排名靠前，所預測的可卡因生物降解途徑為尚未經過實驗測試的潛在異源降解途徑。

Ahsanul Islam等[68]利用ReactPRED和RetroPath2.0兩種工具對苯、苯酚和1，2-丙二醇的生物途徑進行預測，結合底物可行性和熱力學可行性兩個篩選指標，最終共獲得49條生產苯、苯酚和1，2-丙二醇的預測途徑，包含了從乙酸鹽、葡萄糖和丙酮酸鹽起始到苯、苯酚和1，2-丙二醇的106個反應，且25條預測途徑完全由新反應組成，表明生物逆合成預測加速了潛在新反應的發現。

BioNavi-NP[7]對倍半萜類衍生物Sterhirsutin J和戊二酸的生物合成途徑進行了預測，并結合底物相似性和途徑長度兩個指標對候選途徑進行篩選，所得戊二酸的新生物合成途徑已被實驗驗證[69]。Lee等[67]利用Park等[70]報道的工具對短鏈伯胺的生物合成途徑進行了預測，結合底物相似性、反應位點相似性、熱力學可行性、路徑距離及酶與底盤的適配性五個評價指標對預測途徑進行排序，所預測異丁胺生物合成途徑在大腸桿菌中產量最高達到10.6 7g/L。Smolke等[71]利用BNICE.ch對那可丁衍生物的生物合成途徑進行了預測，通過化合物引用次數、相關專利數對候選衍生物進行篩選，并結合熱力學可行性、底物相似性、候選衍生物的生理功能對候選產物的預測途徑進行評價，最終得到了(S)-四氫巴馬汀、(S)-armepavine、(S)-laudanine和(S)-nandinine四種衍生物的候選途徑，并在釀酒酵母中成功地進行了途徑構建。

3 調控元件的設計和優化

微生物細胞工廠的設計不僅包括生物合成途徑的設計，途徑中所需的編碼基因在底盤宿主中表達時還需要一系列必要的轉錄和翻譯調節元件。這些元件將在一定程度上決定途徑中酶的表達活性，并進一步影響菌株的生長和目標化合物的產量。因此需要對調節元件進行設計和優化，以精確控制酶活性且提供途徑與底盤宿主的適配性，而人工智能的出現加速了元件從頭設計的研究，并使人工定制遺傳元件成為可能。由于真核生物轉錄和翻譯調控十分復雜，對調控元件的研究集中在原核生物，尤其是大腸桿菌表達體系。

3.1 轉錄水平元件

啟動子是在轉錄水平調控基因表達的關鍵元件，可驅動對基因表達的調控。先前尋找新啟動子的研究主要集中在通過誘變或調控元件組合對已知啟動子進行改造并形成啟動子文庫，結合人工智能手段對啟動子的強度進行預測，以實現為細胞工廠提供不同轉錄強度的啟動子元件（圖4）。

圖4 啟動子的設計Fig.4 The design of the promoter

在大腸桿菌的啟動子設計研究中，SelProm[72]將擁有120個質粒的BglBrick文庫中的誘導型啟動子替換為組成型啟動子，構建了10種不同表達強度水平的組成型表達質粒，覆蓋的表達強度水平中最弱為未誘導的PlacUV5的1/5.6 ，最強比誘導的Ptrc高4.3 倍，最強與最弱表達強度水平之差為156倍，以良好的分辨率提供了廣泛的表達水平。基于該數據，SelProm利用偏最小二乘回歸(partial least squares regression)算法建立了預測選擇模型，可對不同的質粒成分參數（啟動子、抗性基因）下誘導型和組成型質粒的表達強度水平進行預測，并推薦目標表達強度水平相應的誘導型和組成型啟動子，實驗結果驗證了啟動子推薦工具的有效性。此外，Deng等[73]基于易錯PCR技術對pTrc99a質粒上的Ptrc啟動子進行誘變，產生了由3665個突變體組成的大腸桿菌人工啟動子文庫，所跨越的表達強度水平超過兩個數量級，最強的啟動子是1mmol/L IPTG誘導的PT7強度的1.52 倍。使用該合成啟動子庫作為輸入數據集，構建并優化了基于XGBOOST機器學習算法的啟動子強度預測模型，可對所設計的人工啟動子的轉錄強度水平進行預測，且經比較發現，理性設計的一百個人工啟動子的預測強度和實際強度十分接近（R2=0.88 ），從而驗證了XgBoost模型在啟動子表達強度的預測上的可靠性。

隨著啟動子突變體文庫的增多，越來越多的啟動子序列及其表達強度數據被公開報道，同時結合合成生物學及生物信息學領域的快速發展，啟動子的從頭設計成為可能。Wang等[74]基于生成對抗網絡（generative adversarial networks，GANs）從大腸桿菌天然啟動子中學習關鍵特征（k-mer頻率、-10和-35基序及其間距限制），以捕獲不同位置的核苷酸之間的相互作用，從而建立了大腸桿菌啟動子的從頭設計方法。人工啟動子可基于大腸桿菌中的啟動子活性和預測模型進行優化，兩輪優化后高達70.8%人工啟動子被實驗驗證了其調控轉錄水平的功能，且多數人工啟動子與大腸桿菌基因組在序列上具有正交性。此外，其中一些人工啟動子顯示出與大多數天然啟動子及其最強突變體相當甚至更高的活性，表明深度學習的方法可以為細胞工廠的設計提供更廣泛的遺傳元件來源。

真核生物的啟動子設計研究報道較少，且集中于釀酒酵母表達體系，主要通過從釀酒酵母內源啟動子中獲取保守的模體（motif），并對模體之間的間隔序列進行設計，得到大型人工啟動子文庫，并結合人工智能算法建立啟動子的預測模型。如Smolke等[75]以酵母內源TDH3啟動子為研究對象，獲取了TDH3啟動子中轉錄因子結合位點、TATA框、轉錄起始位點等保守序列，對保守序列之間的間隔序列進行隨機設計，得到超過675000條基于TDH3啟動子的酵母人工啟動子文庫，測量了其中327000條序列的基因表達活性，并利用卷積神經網絡算法建立了具有較高預測準確性的人工啟動子表達強度預測模型。

3.2 翻譯水平元件

對于原核生物而言，在翻譯水平上的調控主要通過核糖體結合位點(ribosome binding sit,RBS)的設計來實現。RBS Calculator[76]、RBSDesigner[77]、RedLibs[78]、PartsGenie[79]等工具可預測RBS序列的翻譯起始速率以估計給定mRNA序列的蛋白質表達水平，并被用于設計符合所需翻譯起始速率的RBS序列。

例如，RBS Calculator[76]基于翻譯啟動階段關鍵分子相互作用的吉布斯自由能建立了平衡統計熱力學模型，通過將熱力學模型與隨機優化方法相結合，可設計具有特定的翻譯起始速率或使得編碼序列具有盡可能高的翻譯起始速率的RBS序列，此外RBS Calculator還可以通過手動設計或復制強大的自然序列來設計比以前更強大的人工RBS序列。RedLibs[78]基于RBS Calculator可生成全局優化的簡并RBS文庫，以減少RBS文庫中的冗余元件，且在文庫規模盡量小的情況下包括更多中等或高強度的RBS序列，所得的兼并RedLibs文庫中的核糖體結合位點樣本能夠以線性方式均勻覆蓋整個翻譯起始速率(rates for translation initiation,TIR)空間，充分滿足不同強度RBS序列的選擇需求。

此外，Ding等[80]將RBS的設計與生物傳感器相結合，利用DNA微陣列構建了包含12000個RBS的葡萄糖酸生物傳感器，通過熒光激活細胞分選（FACS）檢測了生物傳感器中綠色熒光蛋白的熒光強度，利用卷積神經網絡（convolutional neural network,CNN）對其中7053個RBS的七個特征（RBS的GC頻率，堿基A、T、C、G的頻率，SDn的GC頻率及SDm的GC頻率）進行訓練，建立了可預測RBS序列表達強度的神經網絡模型CLM-RDR，能夠快速確定與RBS序列對應的生物傳感器的平均動態表達范圍及其序列特征（圖5）。

圖5 核糖體結合位點的設計Fig.5 The design of the RBS

4 生物合成途徑與調控元件的組合設計

在微生物細胞工廠的設計階段需要確定所構建基因線路中啟動子、終止子、標簽、整合位點、載體等元件的組合方式，使目標途徑能夠發揮正常功能。在傳統實驗設計中由于不同特性的元件使得元件-途徑的組合空間變得十分巨大，組合優化將使得實驗量呈現指數級增加，在缺乏高通量構建平臺的情況下難以實現，因此通常只能一次對某一類元件進行優化（one factor/one variable at a time,OFAT or OVAT），但元件間的內在聯系和約束使所獲取的組合方式只能達到局部最優。實驗設計（design of experiments,DoE）是在生物工藝工程中廣泛應用的高效探索大規模設計空間的系統方法，可以對組合空間進行優化、壓縮，得到精簡的組合空間，進而為獲取全局組合優化提供了可能。在代謝工程及生物制造領域，DoE被用于優化實驗條件以提高目標化合物產量（圖6）。

圖6 實驗設計Fig.6 Design of experiments

4.1 DoE的原理

DoE可以對有不同取值（levels）的因子（factors）的組合空間進行評價，并通過一些代表性實驗來高效探索設計空間。DoE包括識別相關因子和對相關因子的取值組合優化兩個過程。在微生物細胞工廠構建過程中，啟動子、拷貝數、抗生素抗性等遺傳相關變量以及碳源、氮源、添加物等培養基相關變量均可以作為具有不同取值的因子。而與預測途徑及底盤宿主相關因子的識別和選擇依賴于先驗知識，其中在DBTL循環的早期階段，先驗知識較為匱乏，需要考慮較多的因子，并通過盡量少的DBTL的循環對設計空間進行精簡。不同因子組合要實現的典型優化目標包括目標化合物產量的最大化（如增加得率、速率和產量）、使中間化合物毒性最小化等。

4.2 DoE優化途徑與元件組合空間

Carbonell等[81-82]基于R程序包planor和DoE.base先后采用了以啟動子強度、途徑中基因位置次序及質粒拷貝數為因子和以底盤菌株和培養基組成為因子的方式對大腸桿菌中黃酮化合物(2S)-喬松素的產量進行優化，將(2S)-喬松素的產量提高了500倍。Singleton等[83]結 合DoE軟件JMP Pro（SAS Institute Inc.USA）和人工神經網絡（artificial neural networks,ANN），在搖瓶和96孔板兩種生長環境條件下，將碳源、氮源、磷酸、維生素、氨基酸濃度等培養基相關變量作為因子，對熱葡糖苷土芽孢桿菌（Geobacillus thermoglucosidans）的發酵培養基組分進行優化，實現熱葡糖苷土芽孢桿菌利用己糖單糖和二糖進行生長，并可產生乳酸或乙酸鹽。Radivojevi?等[84]結合機器學習和貝葉斯概率模型建立了自動推薦工具（Automated Recommendation Tool,ART），輸入組學數據或啟動子的組合可預測最終產量各種可能值的概率分布情況，同時基于逆向設計可提供使得目標化合物產量最高的候選組學數據或啟動子組合，從而指導下一輪DBTL循環的實驗設計，但ART目前僅支持單目標優化。Jensen等[85]將上述ART算法與釀酒酵母GEM、機器學習算法相結合，以啟動子為因子對色氨酸（tryptophan）的生產進行優化，最終滴度和產量分別提高74%和43%。

5 結論與展望

生物制造作為制造業可持續發展的重要方向受到人們越來越多的關注，而微生物細胞工廠是生物制造的有力手段。當前持續積累的生物大數據極大地促進了計算機輔助設計工具的發展，對微生物細胞工廠的智能化設計將起到革命性的作用。本文依據細胞工廠在實際構建中的先后次序對細胞工廠中生物逆合成途徑的預測與篩選、轉錄水平和翻譯水平上遺傳調控元件的設計、途徑與元件的組合優化三個環節相關的智能設計工具進行了綜述。

在生物合成途徑的預測與篩選環節，生物逆合成算法基于泛化的反應規則擴展逆合成網絡，并利用多種指標對預測途徑的途徑可行性和底盤適配性進行綜合評價，給出最具實際可行性的推薦途徑，幫助人們進行已知途徑化合物的途徑優化及提產和未知途徑化合物的途徑預測及設計。此外，蒙特卡洛樹等人工智能算法的引入為生物逆合成算法的發展提供了新的思路。但值得注意的是，當前生物逆合成工具預測途徑的假陽性率仍然較高，其主要原因在于途徑評價算法不能充分模擬底盤宿主對途徑的選擇，如不能精確計算化合物對指定微生物宿主的毒性、未考慮化合物及酶的區室化對途徑表達的影響以及缺乏酶催化底物雜泛性數據等。目前AlphaFold2[86]和RoseTTAFold[87]的出現可實現基于序列預測蛋白質晶體結構，為酶與底物的適配性問題提供了解決思路。

在遺傳調控元件設計環節，主要研究對象是原核生物的啟動子和核糖體結合位點設計。對原核啟動子的設計主要從對已有啟動子進行改造（誘變、易錯PCR等）或從頭設計兩種手段獲取新啟動子，并基于啟動子調控的熒光蛋白的表達強度數據建立預測模型及人工啟動子庫，為所需表達強度提供推薦的啟動子序列。由于大腸桿菌等原核生物中啟動子的長度較短，一般小于150bp，而酵母等真核生物的啟動子較長，通常為數百個核苷酸，且調控機制更為復雜，因此開發適用于真核生物啟動子設計和表達強度精準預測工具是遺傳調控元件設計的重要挑戰。此外，原核生物核糖體結合位點的設計主要基于熱力學模型，可計算RBS序列的翻譯起始效率，并提供推薦的RBS序列。值得注意的是，對這兩種遺傳調控元件進行設計的工具均很少考慮目標序列的實際表達環境參數（菌株、pH、質粒類型、標簽等），減少了推薦的遺傳元件可靠性。

在途徑與元件組合優化環節，DoE方法簡化了遺傳元件與途徑的組合空間，可以在多輪迭代設計中優化途徑表達及目標產物的產量，但DoE方法難以充分考慮途徑設計及構建中的相關因子，如底盤宿主中可能影響途徑表達的干擾因子，使得DoE推薦的組合方案具有一定的局限性。

此外，目前微生物細胞工廠相關的設計工具只能相對獨立地進行特定環節的設計，不能實現微生物細胞工廠的一站式自動化設計“流水線”，需要靠人力來完成各部分設計工具間的連接，因而不能有效地提高細胞工廠的設計效率，通過整合現有工具資源或創制新的工具、統一接口、建立標準化的自動化設計工作站，將是微生物細胞工廠智能設計的重要發展方向。