平質,張顥齡,陳世宏,倪鳴,徐訊,朱砂,沈玥,5
(1深圳華大生命科學研究院,廣東 深圳 518083;2深圳市合成生物學創新研究院,中國科學院深圳先進技術研究院,廣東 深圳 518055;3廣東省高通量基因組測序與合成編輯應用重點實驗室,深圳華大生命科學研究院,廣東 深圳518120;4(廣東省)華大基因合成基因組學院士工作站, 深圳華大基因科技有限公司,廣東 深圳 518120;5深圳國家基因庫,廣東 深圳 518120;6英國牛津大學大數據研究所,牛津 OX3 7LF;7英國TAICHI AI Ltd., 倫敦 N1 7GU)

作為一種備受關注的新型數據存儲介質,DNA分子與傳統存儲介質相比在信息容量、存儲時間及維護投入等方面都極具優勢。尤其是隨著DNA合成和測序技術的快速發展,基于DNA的數據存儲在編碼方法與全技術集成系統方面均有了相當的進展[1-2]。如圖1(a)所示,常規的DNA存儲流程包括從數字文件的編碼,生成編碼DNA序列的從頭合成,通過測序確定DNA文件的序列信息讀取信息,根據編碼方法進行解碼以恢復原始數字文件。其中,作為實現DNA數據存儲的關鍵步驟,編解碼方法的開發在過去十年間積累了大量的研究基礎。自2012年起,該研究方向聚焦于提升編碼密度、通過控制DNA序列的特定生化參數提升DNA存儲與現有技術的兼容性、提升編解碼方案的穩健性,George Church、Goldman、Grass、Erlich等研究團隊提出了不同策略的DNA存儲編解碼方案[3-7]。在兼容性方面,極端(極高或極低)的GC含量或單堿基長串重復對現有上下游銜接技術都非常不利,會造成DNA合成困難以及DNA測序錯誤。因此,Church轉碼采用兩種堿基信息對應一種二進制信息的方式,盡管編碼密度在所有算法中較低,但其利用隨機替換的方式均衡GC含量并避免較長的單堿基重復[3]。Goldman轉碼利用霍夫曼(Huffman)編碼將二進制信息首先轉化為三進制信息,再利用輪轉編碼的方式生成DNA序列以避免單堿基重復出現[4]。Grass轉碼則利用伽羅瓦域(Galois field,GF)進行編碼,使最長單堿基重復長度不超過3[5]。DNA Fountain轉碼方法首次引入了通信編碼中的噴泉碼,使得其凈信息密度盡可能逼近香農極限,并通過有效性篩選篩除了GC含量不均衡或最長單堿基重復長度超過4的DNA序列,使得生成的DNA序列與現有技術的兼容程度大幅提高[6]。Yin-Yang轉碼利用2條二進制序列通過不同規則組合成1條DNA序列的方式,通過有效性篩選,使輸出序列的GC含量、最長單堿基重復長度以及DNA二級結構自由能都在指定范圍內[7]。在編碼穩健性方面,為了應對存儲過程中可能出現的錯誤,Grass轉碼通過引入Reed-solomon(RS)糾錯碼,增加了編碼的穩健性,這種添加糾錯碼的策略被后續的許多轉碼方案用來確保數據的準確性[6-9]。在編碼學領域也有基于受限編碼(constrained coding)相關理論對DNA存儲的信道模型進行理論研究的相關報道[10-11]。當然,這些編解碼方案也存在局限性。在遇到某些特殊數據結構時,Goldman轉碼仍然會生成GC含量極端的DNA序列,同時,由于堿基信息間相互關聯,單個編碼位出現錯誤,存在后續編碼位也有發生連鎖錯誤的風險。對Grass編碼來說,GC含量不均衡仍然是其主要風險。DNA Fountain由于其編碼方法的特性,在追求極高編碼密度的同時,存在有些數字文件無法被編碼生成對應DNA序列的問題,或DNA文件無法解碼導致完全損壞的風險。因此,該領域研究的研究重點也從以往追求高信息存儲密度和技術兼容性的實現,逐漸在編碼方法對不同文件的通用性以及高穩健性等方面進行了拓展。
然而,目前已開發的多種轉碼方法所采用的編程語言、編碼所設定的技術參數和標準各不相同,不利于基于已有研究基礎的后續開發優化,從應用端看針對不同類型的數據文件最適配的方案選擇也缺乏相應的評價與選擇標準,從而對促進該領域的交流與發展帶來了阻礙。因此,對于已開發的不同轉碼方法的系統評價標準應當形成明確的共識,并建立相應標準化的體系。此外,現有的轉碼方法在應對不同數據結構的文件或在不同應用場景下時,表現亦有差異,通過建立系統集成評價平臺,針對不同類存儲需求以靈活調用的方式提供相應的存儲方案也將促進更多DNA存儲應用的普及。同時一個開放性、可拓展的集成系統也將有利于該領域的研究學者基于前人研究基礎持續提升優化,逐漸形成領域共識,促進這一新興領域的規范發展。
為解決前文所述的問題與需求,本研究開發了一個用于評價與輸出推薦編碼方案的系統集成軟件平臺Chamaeleo[圖1(a)],針對不同類型的輸入文件對模塊化集成在該平臺中的已報道經典轉碼方法提供了相應的評價方案[圖1(b)]。通過選擇此前經典轉碼方案關注的編碼密度、GC含量以及最大單堿基重復長度,對輸出序列進行定性定量的分析統計。進一步地,通過建立存儲過程中的出錯模型和序列丟失模型,對存儲過程進行模擬,并據此分析轉碼方法的穩健性。針對未來某些特殊應用場景,Chamaeleo也提供了不同轉碼方法所對應數據加密性的理論分析。此外,轉碼方案的編解碼計算效率也將作為額外的參數進行統計,為實際應用提供指導。考慮DNA存儲領域還處于快速發展初期,Chamaeleo采用開源、可擴展的設計思路,以支持未來新轉碼方法的持續嵌入,以及新評估參數的加載。針對特定文件類型,Chamaeleo可通過多維度的分析,對擬采用的轉碼方法進行較全面的評估,從而為基于特定需求的最優轉碼方案的選擇提供依據。
Chamaeleo作為一個DNA存儲轉碼方法的集成與評估平臺,主要功能定位于基于同一文件量化分析不同轉碼方法之間的特征參數系統性評價并提供方案選擇指導。如圖1(c)所示,Chamaeleo平臺分為三個主要的模塊:轉碼模塊、糾錯模塊和流程模塊。

圖1 Chamaeleo平臺簡介[(a)DNA存儲的常規流程:該流程分為終端部分和生化部分,Chamaeleo針對DNA存儲中的轉碼步驟進行設計,串聯上下游DNA合成與測序;(b)Chamaeleo搭建的評價體系,從多維度對轉碼方案進行量化分析;(c)Chamaeleo的程序架構及模塊間的相互關系。Chamaeleo開源軟件平臺可通過https:∥github.com/ntpz870817/Chamaeleo進行訪問,或通過pypi:pip install Chamaeleo的方式在計算機終端進行安裝和使用]Fig.1 Brief introduction of Chamaeleo[(a)General procedure of DNA storage includes in-silico transcoding and experimental DNA synthesis and sequencing.Chamaeleo focuses on in-silico transcoding part and connect the related technologies.(b)The evaluation system created by Chamaeleo quantitatively analyzes coding schemes from multiple aspects.(c)Software architecture and relations between modules and classes.The source code is available at https://github.com/ntpz870817/Chamaeleo,which can bealso installed by thecommand of pip.exe,"pip install chamaeleo"]
轉碼模塊集成了6種經典的轉碼算法,包括Church轉碼算法[3]、Goldman轉碼算法[4]、Grass轉碼算法[5]、Blawat轉碼算法[12]、DNA Fountain轉碼算法[6]和Yin-Yang轉碼算法[7]。此外,為提供評估指標的基準參考,Chamaeleo提供Base coding作為基準轉碼算法(即使用A-00、C-01、G-10、T-11的無約束的映射關系進行轉碼的算法)參與到現有的轉碼算法的評估體系中。基準算法的評估可以表征無約束的和帶有約束的轉碼算法對特定類型或格式文件中編碼得到序列的兼容性差異。
為使轉碼方法易讀并保持后續新開發編碼方法可拓展集成,本研究對已開發原始轉碼算法進行了系統性重構[13]。所有經典轉碼算法都通過抽象轉碼算法(abstract coding algorithm)實現Chamaeleo平臺系統中的特定接口[14]。未來開發的轉碼算法也可通過該路徑實現對應接口的建立或通過繼承平臺上已實現轉碼算法進行特定接口的修改,以完成代碼層面的快速搭建(如表S1)。

表S1 使用代碼調用流程完成轉碼過程實例Tab.S1 Example commands to invoke transcoding process
糾錯模塊采用轉碼模塊相同的架構,目前包含了DNA存儲轉碼方案中最常用的兩種糾錯碼:Hamming碼[2]和RS碼[5,15]。 通 過 抽 象 糾 錯 碼(abstract error correction)的接口,Chamaeleo實現了糾錯碼與校正序列信息兩種功能的嵌入。同時,糾錯模塊在處理比特序列之上,實現了擴展接口,以滿足二維及以上比特矩陣的處理。
流程模塊則用于進行實際轉碼/評估任務的執行,轉碼模塊中的轉碼算法以及糾錯模塊中的糾錯碼都會通過實例化的方式被流程模塊中的具體流程所使用。通常,一個流程會包含至少一個轉碼算法以及根據用戶來確定是否采用及指定所采用的糾錯碼。除了最基本轉碼流程,流程模塊還包含三個統計評估流程:算法屬性分析流程、序列特征分析流程和轉碼穩健性分析流程(見下文“評估體系”),用于轉碼方法的系統性分析,為后續基于不同文件采用的最適配轉碼方法提供參考依據。
為了進一步增強實用性,Chamaeleo平臺亦編寫了如數據操作工具(data handler)和流程監控工具(monitor)等DNA存儲所常用的工具。數據操作工具支持讀/寫命令行字符串或文件、壓縮/解壓二進制信息、保存或加載轉碼算法/糾錯碼/流程模塊。流程監控工具用于訂制化地監控算法或流程執行過程中的完成進展,以及記錄特定的過程信息。
Chamaeleo平臺建立了轉碼方案評估體系,對已集成轉碼方案進行系統評估,主要涉及算法基本信息、兼容性和穩健性這三方面的系統性評估分析。為建立可比較的評估策略,該平臺提供了基準測試文件和基準轉碼算法,用于收集評估信息并進行統計分析和對比。在基準測試文件方面,考慮到不同的文件具有不同的數據特征,繼而會對不同轉碼方案的實際性能產生影響。因此,該平臺目前選擇了10個基準文件并盡可能覆蓋更多的實際數據特征。基于操作系統最常用的文件類型[16],基準測試文件被分為5類:圖片、音頻、視頻、文本和可執行文件(或壓縮包)。由于文件越大,寡核苷酸庫中的單條序列所需索引長度就越大,同時包含特殊數據特征的可能性也就越大。而文件比特頻率的分布會對DNA存儲的轉碼產生影響,例如,當特殊數據(如大量重復的“0”或“1”)的比例較高時,轉碼算法不易得到與現有技術兼容性較好的DNA序列。同一類型的文件的不同文件格式會使得信息的比特頻率分布出現差異(如圖S1圖片、音頻、視頻)。同種格式文件,因為內容有差異,也會導致比特頻率的分布出現差異(如圖S1文本類)。而如果采用了壓縮算法,比特頻率的分布將趨向平均[17](如圖S1二進制類),但需要承擔壓縮包損壞導致文件無法恢復的風險。綜合以上因素,本研究中所選的不同測試文件在文件格式、文件大小或比特頻率上存在顯著差異以確保評估體系的可參考性。

圖S1 10個典型測試文件的文件大小和字節頻率[每個分布中的4條橫線從下到上表示1%~4%的出現比例。部分比特在文件中的出現概率超過4%,最高可達42.67%(BMP代表文件中的比特“00101110”)]Fig.S1 Sizes and byte-frequency distributions of 10 typical test files[The 4 lines in each distribution indicates 1%to 4%from bottom to top.The probabilities of occurrence of few specific bytes exceed 4%,which are labeled with digits.The highest probability of occurrence is 42.67%("00101110"in BMPfile)]
Chamaeleo對轉碼方案的評估目前分為三個維度:①算法基本信息,通過算法屬性分析流程統計信息密度、轉碼所需時間以及轉碼成功率等;②技術兼容性,通過序列特征分析流程評估轉碼方案與上下游技術的兼容性,通過收集被轉碼方案編碼(一個或多個文件)獲得的DNA序列,統計并分析針對這些DNA序列與上下游銜接技術的兼容性相關參數,比如GC含量和單堿基重復等;③穩健性,由于針對編碼DNA序列的合成與測序實際操作過程中,不可避免地會因為所采用的技術本身的局限性而產生錯誤(堿基插入/替換/刪除)或序列丟失。通過轉碼穩健性分析流程,在轉碼獲得的DNA序列文庫中分別隨機引入錯誤和序列丟失,計算源文件的成功恢復率,進而表征該轉碼方案對錯誤和序列丟失的可容錯性。
為降低DNA存儲技術的成本,轉碼方法通常需要在保證數字信息完備的情況下盡可能減少DNA序列合成量,因此,編碼密度是轉碼方法最重要的評價指標之一。在DNA存儲中,除了攜帶源文件二進制信息的數據區,為了確定混合DNA分子中不同DNA序列對應源文件中信息的位置,在編碼DNA序列中會引入索引區。DNA擴增是目前DNA存儲流程中實現獲取部分數據所常用的技術方法,因此也需要設計引物區,結合特定引物對文庫進行擴增,以便數據備份和測序建庫。此外,為了使這些DNA序列在解碼過程中可以一定程度上糾正堿基突變帶來的問題,編碼DNA序列往往還會選擇引入糾錯區以提升數據可恢復的穩定性[圖2(a)]。基于特定編碼設定下的二進制-堿基間映射關系,可以計算出不同的轉碼方案對應的理論編碼密度。但如上所述,除數據區外,DNA序列還會包括索引區、引物區和糾錯區,因此實際的編碼密度會低于理論編碼密度。例如,當信息區的長度固定,伴隨著數字信息(文件)的總比特數的增加,其索引區占整條序列的比例會增加,從而導致編碼密度的降低。與其他轉碼方案不同,DNA Fountain轉碼采用Luby Transform(LT)碼[6,18],使用長度為16個堿基的索引區保存每條DNA序列所包含的隨機數種子。因此,在數字信息大小不超過500 MB的情況下[6],其實際信息密度不變。而當數字信息的大小大于這個閾值,DNA Fountain轉碼需要使用更長的索引區(超過16個堿基)去保存每條DNA序列包含的隨機數種子。

圖2 DNA序列結構設計與已收錄編碼方法的凈信息密度評估[(a)輸出DNA序列設計,糾錯碼為可選項,其中,Hamming和RS碼在結構上的分布略有不同;(b)測試文件通過不同轉碼算法獲得的凈信息密度,其中Goldman轉碼算法使用的報道中生成的哈夫曼樹[4],DNA Fountain轉碼算法使用偽隨機數生成器中默認度分布參數(δ=0.05,c_dist=0.1)以及默認的7%的冗余[6],Yin-Yang轉碼算法則選擇相關研究中采用的第888號規則[7]]Fig.2 Structure of output DNA sequences and their net information densities[(a)Basic design of output DNA sequences with/without optional error-correction codes(Hamming code and RS code).(b)Distribution of net information density using different coding schemes.The setting of parameters is identical to that of original report:For Goldman's coding scheme,the Huffman tree used in the evaluation process is set as the same in Goldman,et al[4].For DNA Fountain scheme,the degree distribution tuning parameter(δ=0.05 and c_dist=0.1)and redundancy(7%)used in the evaluation process is set as the same in Erlich,et al[6].For Yin-Yang coding scheme,rule No.888 wasused as reported in Ping,et al[7]]
此前有相關研究圍繞信息密度提出了一些相關計算方法與相應的概念用于評價轉碼方案的性能[3,6,19]。 如 堿 基 編 碼 率 (base coding density,bit/nt)描述的是在不考慮分子拷貝數的情況下,實際操作中一個堿基攜帶的有效比特數量,即:

而物理信息密度(physical information density,petabyte/gram,即PB/g)則是基于流程操作考量下提出的概念[20],其描述的是經過換算后,實驗中每克堿基攜帶的有效信息量(字節數),即:

此外,還有一些概念,如編碼潛力(coding potential)和實際容量(realized capacity,凈信息密度與信道的香農容量之比)[6],也相繼被提出,但由于它們涉及信息論相關的復雜理論推導,與實際應用層面數據差別較大,因此在DNA存儲領域未被廣泛認可。參考該領域中如上概念的普遍認可程度,參考引物區長度在不同存儲的應用需求下會出現差異的情況,并結合考慮目前一般采用長度為200 nt的寡核苷酸文庫作為DNA存儲的主要方式(數據區長度約120~140 nt),Chamaeleo平臺中選擇采用數據區編碼256 bit數據條件下的凈信息密度(net information density,即輸入信息的比特數除以轉碼完成后排除引物區序列的堿基數)[6]作為評估轉碼方法編碼密度的參數:

在編碼密度評估中,本研究基于無糾錯碼、含漢明碼和含RS碼這三種情況分別進行了評估。為和此前的報道統一,本研究在二進制信息層面進行糾錯碼的插入。程序默認原始比特序列長度為128 bit,采用經典設置的(7,4)Hamming碼,其校驗位長度為8 bit,且校驗位插入在原始比特數據的2的冪次位上,理論可以糾正至多1個堿基替換。目前,此前提出的轉碼算法皆采用的是RS碼,通常設置3個字節以糾正至多3個堿基替換錯誤[2,21]。此外,按照文獻的設計,RS碼通常會被放到原始比特序列的后面。通過針對目前已開發的6種編碼方法進行對應編碼密度的量化,我們發現由于測試文件的大小不同,受到索引區長度變化的影響,絕大多數編碼方法的編碼密度都在一定范圍內發生了波動,而DNA Fountain轉碼由于索引區長度固定,編碼密度最為穩定。另一方面,Yin-Yang轉碼表現出相對較高的編碼密度,接近于無拘束的基準轉碼方法。此外,總體趨勢上不同的編碼方法在引入糾錯功能后均顯示出較其理論編碼密度有5.19%~13.49%的下降。采用RS碼比Hamming碼的下降更為顯著,說明轉碼方法的糾錯能力越強,越需要支出更多的信息密度[22]。
DNA存儲中的兼容性體現在與現有技術(DNA合成技術[23]、DNA測序技術[24]、PCR擴增技術[25]等)的適配程度。極端(極高或極低)的GC含量或單堿基長串重復對現有上下游銜接技術都非常不利,會造成DNA合成困難以及DNA測序錯誤,進而導致數據無法恢復的問題。因此避免極端GC含量和單堿基長串重復是目前各轉碼方法重點關注的方向。其他可供選擇的評價參數還包括DNA序列的二級結構自由能、規律性重復序列與特殊序列(如酶切位點、毒性序列)的出現頻率等。
在兼容性評估實驗的過程中,由于添加糾錯碼在一定程度上能夠消除部分對轉碼算法兼容性不友好的數據特征[26],本研究中以不采用糾錯碼的方案作為兼容性下限的評估。
如圖3所示,針對測試文件所涉及到的數據特征,對比基準算法,6個轉碼算法都對兼容性做了相應的優化。依據目前的測試結果,已收錄的轉碼算法都回避了單條序列可能存在的極端GC含量。其中,DNA Fountain轉碼算法和Yin-Yang轉碼算法將GC含量控制在了40%~60%,在兼容性方面表現較為突出。而Church轉碼算法、Goldman轉碼算法、Grass轉碼算法以及Blawat轉碼算法轉碼所獲DNA序列都存在40%以下或者60%以上GC含量的情況,因此存在生成DNA序列與上下游技術不兼容的風險。Yin-Yang轉碼算法的平均GC含量分布在40%和60%附近相比DNA Fountain轉碼算法更低,因此相比DNA Fountain轉碼算法,它更容易進一步提高GC含量的限制。控制單堿基長串重復方面,無約束的基準算法所獲序列的單堿基重復長度存在于1~42之間,而已收錄轉碼算法對其都有嚴格的限制。由于Goldman轉碼方法的單堿基重復設置為1,因此其單堿基重復的峰值聚集在1的區間,而其他方法的單堿基重復并沒有較大的差異,都處于2~4之間,說明目前所有的轉碼方案都在控制單堿基長串重復方面具備良好的兼容性。

圖3 轉碼方案的兼容性評估[在無糾錯碼情況下,不同轉碼算法從10個測試文件中編碼所得DNA序列的GC含量分布(a)和最大單堿基重復長度分布的統計(b)]Fig.3 Compatibility evaluation of coding schemes[Distribution of GCcontent(a)and maximum homopolymer length(b)of DNA sequencesfrom transcoding of 10 test-filesby different coding schemes.For Basecoding in(b),themaximum homopolymer length is42 and not shown in thispanel for aclarity purpose]
由于成本與訪問速度的限制,DNA存儲目前仍然被認為適用于長期的冷數據(無需頻繁訪問的數據)存儲,同時它可以應用基于宿主菌的體內存儲或寡核苷酸庫的形式進行體外存儲。而在這些存儲過程中由于宿主菌自身可能發生的變異積累又或者合成過程引入的錯誤,DNA分子不可避免地會發生堿基插入、刪除、替換等變異以及自然降解。因此,如何從轉碼方案層面應對這些可能帶來信息丟失的風險尤為重要。在轉碼后得到的DNA序列文庫中隨機引入定量的堿基錯誤和序列丟失,再使用對應方案進行解碼,Chamaeleo通過收集和計算所得正確解碼信息對原始信息的覆蓋率作為穩健性評估的指標。
如圖4(a)所示,當引入1%的序列丟失后,大多數轉碼方案的數據恢復率都在98.98%~99%之間。DNA Fountain轉碼與其他轉碼算法相比,單個文件的數據恢復率(恢復出的正確數據量/源數據的總量)波動較大(11.75%~99.99%)。其原因在于:為了滿足篩選條件,DNA Fountain轉碼算法所獲的每條DNA序列(編碼數據包)通常會包含多條比特序列信息,同時這些數據包之間具有相互關聯的拓撲結構,當一條DNA序列丟失后,可能會造成更多DNA序列無法滿足解碼條件,從而丟失更多的比特序列信息。Yin-Yang轉碼算法的數據恢復率在98.95%~99%之間,其恢復率下限略低于除DNA Fountain外的其他算法。在Yin-Yang轉碼算法中,一條DNA序列包含2條比特序列信息,所以當單條DNA序列出現錯誤不可修正,其損失可能會是其他轉碼算法(Church轉碼算法、Goldman轉碼算法、Grass轉碼算法、Blawat轉碼算法)損失的2倍。對任一轉碼方案,利用大量物理冗余或邏輯冗余均可以很好地應對引入的錯誤。以DNA Fountain轉碼算法為代表的抹除碼,依據其編碼原理特征,只要在解碼階段接收到足夠的編碼數據包,源數據即能完全恢復。值得注意的是,由于噴泉碼中的度分布設置,單個編碼數據包與其他源數據分片存在相互關聯,因此當邏輯冗余不足的情況下,解碼過程中單個編碼包的錯誤或丟失可能會連鎖影響其他編碼包的解碼,而多個編碼數據包的丟失可能造成整個文件幾乎無法恢復。另外,直接對源數據進行編碼的轉碼算法,例如Church轉碼算法,在添加源數據量10%的邏輯冗余后,轉碼算法對于<10%的序列錯誤或丟失將具備較強的耐受性,然而不能保證源數據的完全恢復。

圖4 轉碼方案的穩健性評估[(a)在測試文件轉碼獲得的DNA序列文庫中引入1%的隨機序列丟失,對比解碼后文件與源文件,獲得該情況下文件的恢復率。右側圖為98.94%~99.02%區間放大。(b)在測試文件轉碼獲得的DNA序列文庫中引入3%的隨機堿基錯誤(插入、刪除、替換各1%),對比解碼后文件與源文件,獲得該情況下文件的恢復率。右側圖為97%~98.4%區間放大]Fig.4 Robustness evaluation of coding schemes[Distribution of file retrieval rate under condition of(a)1%random sequence loss and(b)3%of nucleotide error(1%for insertion/deletion/substitution,respectively).Figureson theright isthezoom in part for acloser view]
在DNA存儲的生化反應過程中,穩健性評估
需要考慮的因素為由于測序深度、PCR隨機性等生化操作造成的DNA分子的突變或者丟失。這些突變和丟失通常分為系統誤差和隨機誤差。在DNA存儲中,隨機誤差一般由測序產生,而測序過程的隨機錯誤通常可以用序列比對的方式進行相互校正,系統誤差一般由合成或分子生物學操作產生,無法通過常規測序數據處理方式進行校正。因此,在Chamaeleo評估體系中,這里的穩健性評估指的是系統誤差對轉碼算法造成的影響。根據此前的文獻報道[6,20,25],經過測序后,序列的丟失率一般在1%~2%左右,因此本文用1%序列丟失進行穩健性評估。通過常規DNA合成的錯誤率分析,一般認為錯誤率為0.3%左右,而大片段DNA組裝合成中錯誤率會更高,因此本文用各1%堿基錯誤(插入、刪除、替換)進行穩健性評估。當引入堿基插入、刪除、替換各1%的錯誤后,大多數轉碼方案的數據恢復率都在97.05%~98.62%之間。DNA Fountain轉碼算法在此情況下的文件恢復率較低(3.78%~27.52%),顯示其在應對堿基錯誤的穩健性不足。另外,糾錯碼的使用對穩健性的提升有較為明顯的作用。對所有轉碼方案而言,糾錯碼的糾錯能力越強,穩健性越好。值得注意的是,目前常規的糾錯碼,如本文使用的Hamming碼、RS碼等,僅對堿基替換有效,而由堿基插入或刪除導致的錯誤,常規的糾錯碼無法進行錯誤的發現和糾正。但在當前DNA合成流程中,每種DNA序列的合成拷貝數大約為107個,所有拷貝同時出現插入或刪除錯誤的可能性極低,因此,可以通過測序后的序列比對找出共有序列,從而糾正堿基插入或刪除錯誤。
針對現有編碼方法基于特定規則的轉碼過程,如果規則相對簡單易推導,則數據可加密性會較低,相應的有特殊存儲需求的加密數據的安全可控性也會較低。對于部分轉碼算法來說,其比特與堿基之間的映射關系并不是唯一的,如果解碼過程和編碼過程使用的映射關系不同,就無法獲得正確的原始比特信息。因此,這些映射關系可以被視為密碼學的秘鑰,在基本的轉碼任務以外,用于特殊信息的加密和解密。不同的轉碼算法包含的映射關系數量是不同的。通常,映射關系(秘鑰)的數量和信息被破譯的時間成正相關[27]。本研究統計了已集成轉碼算法的映射關系數量,用于評估轉碼方案的數據安全性。
Church轉碼算法和Blawat轉碼算法的映射關系是固定的且只有1種。Goldman轉碼算法的映射關系可以視作Huffman三叉樹[28](k+1=3)與字節(n=256)之間的映射關系,其映射關系的數量可以使用Huffman三叉樹的形態種類表示。通過弗斯-卡特蘭(Fuss-Catalan)數[29-30]進行換算,可得映射關系共有種。Grass轉碼算法的映射規則是在3個堿基(由于規定后兩個堿基不可相同,4×4×3=48個組合)中篩選47個組合作為伽羅瓦域,并將其余比特信息進行映射,所以會存在種映射關系。在DNA Fountain轉碼算法中,4種堿基可以通過不同的排列方式和2個比特組合(00,01,10,11)進行映射,因此該轉碼算法包含=24種映射關系。Yin-Yang轉碼算法由虛擬堿基、Yang規則和Yin規則共同構成。其中,虛擬堿基有=4個選擇;Yang規則是在4種堿基中選擇2種映射為0,剩下2種映射為1,共計有=6種組合;Yin規則是在Yang規則的基礎上進行0和1的分配,共有種組合,因此共有有4×6×256=6144種映射關系。
綜上所述,Chamaeleo平臺將通過計算轉碼方案的映射關系數量(秘鑰數量),作為其數據安全性的參考依據。
DNA存儲轉碼方案的理論層面研究仍存在空白,結合DNA的生物特性以及合成與測序的技術進行的理論研究也相對較少。因此,除了量化分析實際轉碼后的參數,本研究嘗試從理論層面對轉碼算法進行評估分析,通過轉碼算法的映射關系可視化(圖5),并使用圖論[31]進行理論層面的分析與評估。除了使用信息論[24]宏觀描述轉碼算法的信息密度外,圖論分析有助于更細微地剖析轉碼算法的特征與傾向性。

圖5 使用圖論對不同轉碼算法進行可視化[實心紅點表示起點或虛擬起點(即Goldman轉碼及Yin-Yang轉碼中為編碼第一位堿基所用的虛擬位)。在每一輪編碼過程中,已知當前的堿基(節點)和輸入的二進制(箭頭),跳轉至下一個節點并輸出其所指代的堿基,最后將下一個節點作為當前節點。在對應的每一輪解碼過程中,已知當前節點(堿基)和下一節點(堿基),獲得兩個節點之間的箭頭,并輸出箭頭指代的二進制信息,最后將下一節點作為當前節點]Fig.5 Visualization of different coding schemes using Graph-theory-based presentation[Red dot represents starting point or virtual starting point(for Goldman&Yin-Yang coding schemes).During each step of encoding,with known previous base(node)and input bit value(arrow),the current basewould be obtained from graph.During each step of decoding,with known previous and current nodes(base),the bit value(arrow)would beobtained from graph]
對于一個長度為n比特的序列來說,0和1的組合方式是有限的,即2n種。通過轉碼過程,可以獲得轉碼方案可生成的所有DNA序列。進一步將這些DNA序列作為圖的節點,并使用k-mer組裝[32]的方式進行連接作為圖的邊,獲得轉碼算法的圖表示。根據生成圖的各種屬性可以從理論層面分析轉碼算法的性能,而通過圖的出度情況[31],可以近似估計轉碼算法的凈信息密度。圖中每個節點的兼容性即為DNA序列的局部兼容性,基于貪心算法[33],局部特征可以反映全局特征,即全長DNA序列的兼容性。通過圖中的環結構[31],不僅可以統計單堿基重復的情況,還可以統計重復序列、回文重復等復雜情況。對于圖中具有周期性且具有統計意義的子圖或模式[34],可以通過計算,例如標準分數(z-score),分析算法對不同比特數據信息的傾向性。這些屬性可以結合生物特性進行進一步的分析。
目前基于圖論的轉碼算法映射關系可視化有助于理解轉碼過程,轉碼算法可生成的圖種類數量亦能在一定程度上反映出轉碼規則的復雜程度,從而與其在加密存儲的應用中數據安全性的性能形成對應關系。
DNA存儲集編碼學、密碼學、信息學、分子生物學、生物信息學、計算機科學等多學科高度交叉發展而來,其全流程的實現仍高度依賴DNA合成與測序技術的支撐,相關的理論基礎的研究還處于早期階段,是一個充滿想象力的新興研究領域。其中,DNA存儲中的編解碼方法是作為連接數字信息和DNA分子的關鍵步驟,也是過去十年該領域的主要研究方向。不同的編碼方法在存儲信息密度、技術兼容性與存儲穩健性方面各有千秋,但根據不同的存儲需求,目前該領域還缺乏方法間標準化比較評估體系的建立,不利于研究人員基于已有研究的再次開發以及DNA存儲應用端的普及。
因此,本研究開發了一個DNA存儲堿基編解碼算法的可拓展集成與系統評估平臺Chamaeleo,取變色龍可針對不同環境快速適應進行特征變換之義,旨在促進該領域的開發者進行協同開發,為應用端提供一個輔助的指導工具以實現不同存儲需求的應用。Chamaeleo集成了系列已開發的DNA存儲編解碼方法并對其算法在軟件工程層面進行了優化,提供了一個通用的DNA存儲編解碼應用工具。同時,聚焦DNA存儲與現有技術的兼容性以及長期存儲情況下的穩定性,通過選取領域內已有共識的關鍵參數,建立了一套DNA存儲編解碼方案的評估分析體系。利用不同類型、不同格式的文件對經典DNA存儲編解碼方案進行多維度評估,得到的數據可以進行后續的系統性分析。值得指出的是,本研究中采用的評估參數是相互關聯、相互影響的。例如,編碼密度會受到兼容性、穩健性、安全性等方面的影響。從編碼學中constrained coding的角度看,拘束條件必然會導致編碼密度的降低[11]。為提高DNA存儲整個流程的穩健性或保證數據安全性,轉碼方案一般會采用增加糾錯碼、多拷貝[4]或者對比特序列預先進行異或操作[35]的策略。但這些策略會導致冗余的增加,一定程度上也會降低轉碼方案的編碼密度。兼容性則意味著轉碼方法可以通過設置相關參數以降低與上下游銜接技術不兼容序列出現的可能性,以提升轉碼方案的整體可行性。因此,在選擇最優轉碼方案時,需要對所有指標進行綜合考量。我們希望Chamaeleo平臺的建立能促進領域內學者的交流以及新研究者的融入,有助于形成標準化的行業流程與評價指標,從而推動該領域規范有序的快速發展。同時,本研究中首次提出基于圖論的理論評估方法及“特征”“傾向性”等評價指標,旨在促進DNA存儲整體評價體系的發展。在未來,我們期待更多DNA存儲領域的研究者將其獨特的DNA存儲轉碼方法嵌入Chamaeleo開源工具平臺中,也希望能通過廣泛的交流與討論形成更多的有指導意義的評價指標和策略,幫助該領域編解碼方法的理論體系逐漸形成。
盡管DNA存儲目前已經有了諸多的文獻報道,數據存儲體量也逐漸逼近GB級別,但相比于傳統的光盤、硬盤存儲,還需要在很多方面進行突破,Chamaeleo平臺也需要相應的完善和優化。從理論層面,結合DNA合成與測序技術特點,通過數學、編碼學等對DNA存儲的最優轉碼方法進行研究。轉碼應用方面,在某些編解碼過程中,為解決序列內部的錯誤和序列丟失問題,設計糾錯機制時會設置不同的內碼和外碼[5],因此在已有糾錯模塊的基礎上,Chamaeleo也將提供內外碼設置的選項以提高糾錯機制的靈活性。評估體系方面,根據不同應用場景,可以對目前已有評估參數或者新評估參數進行標準化,并設置相應的權重,從而根據實際需求提供最優轉碼方案。從安全體系層面,需要進一步完善數據安全體系,針對個人隱私安全、數據安全、軍事應用等不同場景,建立相應的加密方式,保障數據安全。應用層面,建立高效低成本且功能化的全流程DNA存儲體系需要軟件平臺進一步集成上下游銜接技術參數來進一步提升兼容性方案的輸出,同時也需要考慮數據存儲架構、快速尋址訪問、低成本復寫策略等方面的需求。在DNA存儲之外,目前也有基于質譜技術,利用寡肽或代謝物進行信息存儲的策略被陸續提出[36-37],Chamaeleo平臺也可以提供相應轉碼方案的程序端加載接口,為硅基存儲與碳基存儲建立橋梁。
綜上,基于DNA介質的新型數據存儲作為一種具有劃時代意義的存儲方式,或將打開全球海量數據存儲的新紀元。DNA存儲技術的發展,既需要針對編解碼方法持續取得創新突破,建立顛覆性DNA存儲信息學理論基礎,也依賴于如DNA合成組裝與測序的上下游銜接技術的快速發展以形成高效低成本規模化存儲能力。我們希望Chamaeleo平臺的建立能作為向該目標發展的一個起點,助力生物技術與大數據信息技術的協同發展,在歷史文化傳承信息的永久保存、經濟大數據的長期保存、軍事部署等戰略信息的傳送和保存等領域發揮作用。
致謝:在程序設計和文章撰寫過程中,莊乾龍協助了部分程序模塊的填寫,哈佛大學Geor ge Chur ch教授、周廣宇博士和首都師范大學葛根年教授及蘭昭君參與了算法方面的討論,維也納工業大學的Nat al ie Fr eiber ger協助我們在IOS平臺進行了系統測試,在此一并致謝。
補充材料

圖S2 不同轉碼算法編碼和解碼速度的比較(測試環境:Windows 7系統;i7 CPU處理器;12GB內存;Python版本,3.7.3;IDE,PyCham 2019.1)Fig.S2 Encoding and decoding speed of different transcoding algorithms(Test environment:Windows 7 environment including an i7 CPU and 12 GB of RAM using Python 3.7.3 in the IDEPyCharm 2019.1)