如何提高循環腫瘤DNA（ctDNA）信號的獲取能力以及判定低頻信號的真實性是腫瘤微小殘留病灶（MRD）檢測的難點問題。本研究開發了多變異聯合置信概率分析的MRD生物信息學算法（MinerVa），并在ctDNA標準品和早期非小細胞肺癌患者血漿DNA樣本中對算法效能進行評估。結果顯示，MinerVa算法對多變異共追蹤的特異性穩定在99.62%～99.70%，當追蹤30個變異時，可檢測低至6.3 × 10?5變異豐度的變異信號；進一步對27例非小細胞肺癌患者的數據進行分析，ctDNA-MRD動態監測復發的特異性為100%，靈敏度為78.6%。數據表明MinerVa算法可高效地捕獲血液ctDNA信號，在MRD檢測中具有較好的準確性。
引用本文： 楊飄, 張亞晰, 夏粱, 梅建東, 范銳, 黃宇, 劉倫旭, 陳維之. MinerVa——一種用于實體瘤微小殘留病灶檢測的高性能檢測算法. 生物醫學工程學雜志, 2023, 40(2): 313-319. doi: 10.7507/1001-5515.202303039 復制
0 引言
循環腫瘤DNA（circulating tumor DNA，ctDNA）攜帶了癌組織特定的基因組突變、染色體重排、拷貝數異常和表觀遺傳異常等信息[1-2]。越來越多的實體瘤臨床研究證實，ctDNA可作為檢測微小殘留病灶（minimal residual disease，MRD）的理想生物標志物；ctDNA-MRD監測比傳統影像學能更及時地提示復發，為提前治療干預提供參考，提高患者治愈可能[3-9]。相較傳統的用藥檢測，MRD對ctDNA檢測技術的靈敏度提出了更高的要求[10]。由于單變異的檢出上限具有明顯的“天花板”效應，選擇覆蓋更多變異的大panel通過多個變異共追蹤可顯著提升檢測靈敏度[11]。MRD的主要技術難點是準確判定血漿低頻變異的真實性，由于基于二代測序技術的血漿ctDNA檢測受文庫制備、上機測序以及基因組比對等多個環節的影響，不同堿基、不同變異呈現出不同的技術噪音[12]。精準劃定每個變異的判定閾值是提升檢測性能的關鍵，但當前依然缺乏在單堿基層面上對變異的噪音水平進行有效評價的方法。另外，由于患者腫瘤突變負荷存在差異，在共追蹤不同數量的變異時，怎樣合理設定樣本的MRD陽性判斷標準也是目前ctDNA-MRD檢測中的一個難點。為了解決上述問題，我們采用腫瘤先驗技術路線，基于固定大panel開發了一種嶄新的生物信息學算法——多變異聯合置信概率分析（MinerVa）。MinerVa采用了目前已知的最大規模的陰性血漿數據庫實現了大panel覆蓋區間內單堿基分辨率的技術噪音基線模型搭建，通過共追蹤變異的整合分析對樣本的MRD狀態進行判定，并在小規模臨床隊列中探索MinerVa算法在ctDNA-MRD檢測中的價值。
1 資料與方法
1.1 MinerVa算法應用方案
MinerVa算法采用腫瘤先驗技術路線，使用固定大panel進行血漿ctDNA-MRD檢測。首先獲取患者腫瘤組織和配對白細胞的測序數據，構建患者的個性化腫瘤變異圖譜；然后在血漿游離DNA的測序數據中，根據腫瘤變異圖譜提取相應的變異信號；對提取到的變異信號，使用噪音基線分布模型進行單變異顯著性分析；在樣本層面上，綜合單變異顯著性結果進行多變異聯合置信概率分析，判定樣本的MRD狀態。
1.2 生信分析算法及模型搭建
1.2.1 腫瘤組織變異檢測
本研究通過配對的腫瘤組織和白細胞測序樣本識別體細胞突變。首先使用BWA（0.7.17）將原始的讀段比對到參考基因組hg19上，然后采用變異檢測軟件VarDict（1.5.1）識別突變（包括SNV、Indel）信息，并得到相應的覆蓋深度和突變支持數目。要求位點深度不低于200 X，SNV類型突變頻率不低于5%，Indels類型突變頻率不低于4%。此外，對于repeat/low complexity等區間易重復出現的假陽性信號，基于30個正常樣本構建了PoN模型（Panel of Normal），過濾檢出頻率低于閾值的突變。
1.2.2 基于MinerVa算法的血漿ctDNA-MRD判定
（1）單堿基分辨率的技術噪音基線模型的搭建：噪音基線模型的建立使用了由1 132例陰性血漿數據構建的數據庫。陰性血漿數據是處理后的腫瘤患者及健康志愿者的血漿數據，腫瘤患者的實測血漿數據扣除了配對組織及白細胞變異，健康志愿者的實測血漿數據則扣除了配對的白細胞變異。針對panel覆蓋范圍內的每個堿基都進行了三個變異方向的SNV基線模型的搭建，對每個SNV，提取基線血漿數據庫中所有未扣除該指定SNV的個體數據進行模型搭建，同時為了排除腫瘤異質性信號以及其他未知來源的變異信號，搭建前對指定SNV的離群值做處理。每個變異的噪音分布模型為混合模型，由無變異人群占比模型以及變異模型兩部分組成，第一部分是無指定SNV變異的人群占比（Pzero），另一部分是對含有指定SNV變異人群的等位基因變異頻率（variant allele frequency，VAF）分布擬合的韋伯分布模型或逆伽馬分布模型。Indel噪音模型的構建不區分具體Indel，只定位堿基坐標，將所有關聯指定堿基坐標的Indel變異合并計算Pzero以及擬合分布模型，方法同SNV模型。
（2）單變異信號相對噪音基線的顯著性水平評估：根據血漿變異的堿基坐標信息，調取該變異的噪音基線模型，通過蒙特卡洛抽樣N=10 000，生成N*Pzero個0，同時由變異模型部分生成N*（1–Pzero）個隨機VAF，分別以這N個VAF作為先驗的噪音頻率，根據二項分布計算患者血漿變異信號計算來自于噪音信號的可能性：
 |
其中，VSM為血漿中支持變異的分子數，TSM為覆蓋該堿基的總分子數。
為第i次模特卡洛抽樣獲取的先驗噪音頻率。
綜合N次計算結果，進一步計算 
的加和平均值 
以衡量患者血漿單堿基變異的顯著性程度。P ≤ 0.05作為血漿單變異的陽性判定值。
（3）多變異聯合置信概率分析：該分析檢驗所有P ≤ 0.05的單變異同時是噪音來源的可能性。聯合多個突變，經過多變異聯合置信概率分析計算其聯合置信概率combine_P，combine_P ≤ 0.01則認為樣本MRD陽性。計算公式如下所示，其中，n為追蹤的突變總數，k為陽性變異數目，
為陽性位點的p_value：
 |
1.3 驗證使用樣本
本次驗證包括2套陽性標準品（GW_OGTM006、多變異共檢測標準品）和1套陰性標準品（GM12878），在初步的臨床應用驗證中，我們回顧性地對27例非小細胞肺癌患者的樣本和臨床復發結果進行一致性評估。
1.3.1 GW_OGTM006
GW_OGTM006（菁良科技）是商業化標準品（見附件1）。用正常二倍體細胞系GM12878（科佰生物）梯度稀釋，用于熱點靈敏度和特異性分析。在本次性能驗證中使用的GM12878，實驗數據顯示KRAS p.G12C在多次特異性重復中有低頻信號出現，且頻率低于GM12878出廠質控陰性頻率值（0.1%）。通過定制化Panel進行高深度驗證后證實確實存在低頻變異，因此后續數據分析未納入KRAS p.G12C位點。另外，本方法主要評估的是常見的熱點插入缺失變異，一般在20 nt以內，而MET Exon 14 Skipping是一個長度140 nt的缺失變異，并不在本方法的分析范圍內。因此，后續在熱點變異靈敏度及特異性數據分析中使用的是去除了KRAS p.G12C及MET Exon 14 Skipping之外的12個熱點變異。
1.3.2 多變異共檢測標準品
由于常規標準品的變異數量較低，不能夠滿足多變異共追蹤的MRD分析性能驗證，因此我們使用差異化的單核苷酸多態性（single nucleotide polymorphism，SNP）信號模擬ctDNA的變異，以提高可追蹤的變異的數量。方案為通過將已知SNP信息的隨機一個個體（以個體S指代）的白細胞DNA用GM12878進行梯度稀釋，得到一系列可用于整體分析性能驗證的樣品。將個體S特有的（不同于hg19以及GM12878）的100個SNP變異納入陽性變異集（見附件2），將個體S的白細胞DNA和細胞系GM12878 DNA中與參考基因組hg19具有相同基因型的547個SNP位點納入陰性變異集（見附件3）。
1.3.3 臨床樣本
1 132例基線樣本包含以非小細胞肺癌為主的泛癌種患者的腫瘤組織、血液白膜層和血漿樣本，以及健康人的血漿白膜層和血漿樣本，來自四川大學華西醫院、四川省人民醫院等多個中心的臨床隊列，以及臻和科技公司的健康志愿者。
27例接受根治性手術治療的Ⅰ～Ⅲ期非小細胞肺癌患者在四川大學華西醫院納入，其中Ⅰ期7例，Ⅱ期14例，Ⅲ期6例（詳見附件4），采集術中腫瘤組織樣本及多個時間點的血液樣本，對患者進行為期30個月的隨訪信息收集。臨床樣本收集和使用已通過四川大學華西醫院倫理委員會批準，患者提供書面知情同意。
1.4 驗證實驗流程
1.4.1 GW_OGTM006標準品實驗流程
首先使用正常二倍體細胞系GM12878基因組對陽性標準品GW_OGTM006進行6個梯度的稀釋，按照熱點變異的理論變異頻率均值由高到低分別為1%、0.5%、0.3%、0.1%、0.05%和0.03%，這6個梯度的樣本被分別命名為PC-1P、PC-05P、PC-03P、PC-01P、PC-005P、PC-003P。
對稀釋好的6個梯度的DNA樣本使用Covaris超聲波DNA破碎儀（M220，Covaris，美國）進行片段化，取33 ng片段化的DNA使用 KAPA Hyper Preparation Kit（KK2621，羅氏，瑞士）試劑盒進行文庫構建，構建好的文庫使用xGEN lockdown（1080584，IDT，美國）探針捕獲法進行目標區域富集，使用MGISEQ-T7（華大科技，中國）進行PE100雙端測序，下機后平均測序深度 ≥ 4 500 X，熱點區域深度 ≥ 30 000 X。每個樣本進行了4次技術重復。
1.4.2 多變異共檢測標準品實驗流程
將個體S的白細胞DNA用GM12878進行9個梯度的稀釋，得到一系列可用于整體分析性能驗證的多變異共檢測標準品。預期的理論變異頻率均值由高到低分別為1%、0.5%、0.3%、0.1%、0.05%、0.03%、0.01%、0.005%和0.003%。建庫富集方法同GW_OGTM006標準品，同樣每個樣本進行了4次技術重復。
1.4.3 27例臨床樣本及1 132例基線樣本實驗流程
術中組織標本及白膜層的提取使用天根血液/組織/細胞基因組提取試劑盒（DP304-02，天根，中國），血漿樣本的提取使用MagMAX?游離DNA分離試劑盒（A29319，Applied Biosystems，美國），建庫富集方法同前述標準品。
2 結果
2.1 熱點變異檢測性能
2.1.1 熱點變異檢測靈敏度及檢出限評估
標準品GW_OGTM006統計了4次技術重復中各個變異的檢出情況（見表1），綜合所有的熱點變異（包括SNV和Indel）進行分析。在33 ng起始建庫條件下，預期稀釋變異頻率為0.1%的標準樣本PC-01P中，熱點變異的檢測靈敏度為98%；即使在預期稀釋變異頻率為0.05%的標準樣本PC-005P中，熱點變異的檢測靈敏度仍可達85%。通過probit回歸方法進行檢出限（limit of detection，LoD）評估（見圖1），對于肺癌常見的熱點驅動變異EGFR L858R和19DEL，LoD可達0.055%；綜合標準品中所有的EGFR變異，LoD為0.074%；考慮標準品中12個熱點變異，LoD為0.178%。
表1
GW_OGTM006各熱點變異的檢出率
Table1.
Detection rate of hotspots in GW_OGTM006
圖1
GW_OGTM006熱點變異LoD
Figure1.
Limit of detection of hotspots in GW_OGTM006
2.1.2 熱點變異檢測特異性
使用正常細胞來源的細胞系GM12878基因組DNA進行熱點檢測特異性評價。將GM12878的基因組DNA打斷作為熱點變異特異性檢測樣本，進行8次獨立的建庫捕獲和測序分析，統計在所有樣本中的12個熱點變異的檢出情況，用以衡量熱點變異的檢出特異性。
在驗證中共進行了8*12 = 96次變異檢測，96次檢出結果均為陰性，對熱點變異的檢出特異性為100%（95% CI：96.15%～100%）。
2.2 非熱點變異的檢測性能
分析多變異共檢測標準品的陽性變異集的100個變異作為非熱點變異，統計其在各稀釋梯度的檢出率，通過probit回歸分析獲得最低檢出限。
2.2.1 非熱點變異檢測的靈敏度及檢出限
在33 ng起始建庫的條件下，不同稀釋梯度的多變異共檢測標準品變異的檢出率見附件5，probit回歸分析表明，非熱點單變異的最低檢出限為0.353%（見圖2）。
圖2
非熱點單變異LoD
Figure2.
Limit of detection of non-hotspots
2.2.2 非熱點變異檢測的特異性
陰性變異集中共納入了547個SNP位點，這些多態性位點在個體S的白細胞DNA和細胞系GM12878 DNA中與參考基因組hg19具有相同基因型。對9個梯度的標準品分別進行4次重復獨立建庫和測序分析，在對SNP位點的547*4*9 = 19 692次的變異檢測中，僅有45個假陽性判定結果，非熱點單變異的檢測特異性為99.77%（95% CI：99.69%～99.83%）。
2.3 多變異追蹤時樣本的檢出性能
在多變異共檢測標準品系列樣本中，我們利用差異化的SNP位點分別設置陽性變異集和陰性變異集，用以評估MinerVa算法在多變異追蹤時的檢出靈敏度及特異性。
2.3.1 多變異共追蹤時樣本的檢出靈敏度
模擬先驗的腫瘤變異圖譜，設定變異追蹤數目為1、2、3、6、10、20、30、40和50。對每個指定的變異數目，隨機抽取陽性變異集中的變異進行組合，每個梯度1 000次隨機組合，即形成了1 000組獨立的腫瘤先驗變異圖譜，分別按照每一組圖譜進行MRD檢測，最終統計每個稀釋梯度樣本的這1 000個組合的陽性檢出率（見表2）。結果顯示隨著建庫起始量的升高，變異追蹤數目的增多，檢出靈敏度顯著提高。
表2
追蹤不同數目變異時的檢出靈敏度（%）
Table2.
Sensitivity (%) for tracking different number of variants
進一步經過probit回歸，統計出在33 ng起始建庫的條件下，不同追蹤變異數目時的95檢出限（見圖3）。結果顯示隨著追蹤位點數增加，LoD逐漸降低，當追蹤30個變異時，可檢測低至6.3 × 10?5變異豐度的變異信號。
圖3
追蹤不同數目變異時的LoD
Figure3.
Limit of detection for tracking different number of variants
2.3.2 多變異共追蹤時樣本的檢出特異性
同樣設定變異追蹤數目為1、2、3、6、10、20、30、40和50，對每個指定的變異數目，隨機抽取陰性變異集中的變異進行組合，每個梯度1 000次隨機組合，綜合1 000次隨機組合的MRD判定結果，評估不同追蹤變異個數下的檢測特異性（見表3）。當追蹤多個變異時（變異數目：2、3、6、10、20、30、40、50），特異性穩定在99.62%～99.70%。
表3
追蹤不同數目變異時的特異性
Table3.
Specificity for tracking different number of variants
2.4 肺癌隊列ctDNA-MRD檢測性能
27例非小細胞肺癌患者中，14例患者在隨訪中出現了復發，復發患者的中位無病生存期（disease-free survival，DFS）為337天（166～632天），13例患者未復發。13例未復發的患者在術后動態ctDNA監測中，ctDNA檢測結果均為陰性，特異性為100%（95% CI：77.19%～100%）（見附件6）。14例復發的患者在術后一個月ctDNA檢測陽性的比例為35.7%（5/14），在動態ctDNA監測中，11例患者ctDNA陽性，靈敏度達78.6%（95% CI：52.41%～92.43%）。該隊列研究表明使用MinerVa算法分析的ctDNA-MRD狀態與臨床復發具有較高的一致性。
3 討論
使用固定大panel的MRD檢測方法在臨床研究中已顯示了其可行性[13]。固定panel通過充分的前期驗證，保證了操作的穩定性和臨床的成功率，同時也保證了檢測時效性，可滿足在關鍵時間節點為制定治療決策提供參考的需要。
靈敏度不足是MRD檢測的主要技術挑戰。用藥熱點突變檢測的常規方法是以支持reads數或變異頻率為固定閾值，但MRD檢測中包括了大量個性化的變異，簡單套用這種“一刀切”的方法忽略了不同變異噪音水平的差異，無法平衡不同變異的靈敏度和特異性，從而影響MRD判定的準確性。細化不同變異的技術噪音是合理設定閾值、提升檢測性能的關鍵。既往的研究也嘗試構建技術噪音模型，但由于只使用了少量健康人的數據，差異分層只能集中在變異亞型（subtype，如C > G、C > A等）或者三堿基組合（context）層面[14-15]。Capp-Seq研究嘗試建立單堿基分辨率的噪音模型，但由于只使用了12例健康人血漿，模型構建并不充分，只作為簡單過濾工具使用[16-17]。MinerVa算法創新性地將腫瘤患者及健康志愿者的血漿數據扣除配對組織和白細胞的變異，構建了陰性血漿數據，同時利用固定panel累積數據的優勢，使用超千例的陰性血漿數據，對panel覆蓋區域的每個單堿基都構建了技術噪音分布模型；基于模型，只有高于噪音基線的變異才會被判定為陽性。分析性能數據顯示，EGFR L858R和19DEL的LOD低至0.055%，顯著低于常規熱點變異的檢出限，而非熱點單變異的檢測LoD也達到了0.353%。
樣本水平MRD陽性判斷值目前并沒有統一的標準。常用的方法要求滿足大于等于1個或2個ctDNA單變異檢出[18-19]，但該方法在追蹤變異數量較多時，檢測特異性會下降。假設每個單變異的檢測特異性是99.6%，在追蹤50個變異、使用 ≥ 1個變異的陽性判定條件時，樣本水平特異性低于85%，而使用比較嚴苛的 ≥ 2的陽性變異作為陽性判斷條件，特異性也跌至99%（見圖4），且損傷了檢測的靈敏度。MinerVa算法開發了多變異聯合置信概率分析，將每個變異源于噪音的概率進行整合，通過控制全部血漿變異都來源于技術噪音的概率，保證了樣本檢測的特異性。分析性能驗證顯示，無論追蹤2個變異還是追蹤50個變異，樣本檢測特異性穩定在99.62%～99.70%之間。
圖4
不同MRD陽性判定標準的檢測特異性
Figure4.
Specificity of different MRD positive determination thresholds
本研究也具有一定的局限性。首先，基于固定panel的MRD檢測由于覆蓋的區間依然有限，對一些低突變負荷癌種或患者的變異數量覆蓋不足，可能會影響MRD檢測靈敏度。此外，本研究僅采用小的臨床隊列對MinerVa的臨床價值進行了評估，仍需通過更大樣本的臨床實踐探索本算法在不同癌種及不同臨床場景（如術后復發監測、新輔助/輔助治療療效評估等）中的適用性。
重要聲明
利益沖突聲明：楊飄、張亞晰、范銳、陳維之均是MinerVa相關專利的共同發明人，其他作者均聲明不存在利益沖突。
作者貢獻聲明：楊飄負責分析性能實驗設計、實驗開展，以及論文撰寫；夏粱、梅建東、劉倫旭負責臨床隊列設計、臨床信息收集和臨床數據分析；張亞晰、范銳、黃宇、陳維之負責算法開發和生信數據分析；陳維之負責整體實驗設計。
倫理聲明：本研究通過了四川大學華西醫院生物醫學倫理分委會的審查審批（批文編號：2017年 審（34）號）。
本文附件見本刊網站的電子版本（biomedeng.cn）。
0 引言
循環腫瘤DNA（circulating tumor DNA，ctDNA）攜帶了癌組織特定的基因組突變、染色體重排、拷貝數異常和表觀遺傳異常等信息[1-2]。越來越多的實體瘤臨床研究證實，ctDNA可作為檢測微小殘留病灶（minimal residual disease，MRD）的理想生物標志物；ctDNA-MRD監測比傳統影像學能更及時地提示復發，為提前治療干預提供參考，提高患者治愈可能[3-9]。相較傳統的用藥檢測，MRD對ctDNA檢測技術的靈敏度提出了更高的要求[10]。由于單變異的檢出上限具有明顯的“天花板”效應，選擇覆蓋更多變異的大panel通過多個變異共追蹤可顯著提升檢測靈敏度[11]。MRD的主要技術難點是準確判定血漿低頻變異的真實性，由于基于二代測序技術的血漿ctDNA檢測受文庫制備、上機測序以及基因組比對等多個環節的影響，不同堿基、不同變異呈現出不同的技術噪音[12]。精準劃定每個變異的判定閾值是提升檢測性能的關鍵，但當前依然缺乏在單堿基層面上對變異的噪音水平進行有效評價的方法。另外，由于患者腫瘤突變負荷存在差異，在共追蹤不同數量的變異時，怎樣合理設定樣本的MRD陽性判斷標準也是目前ctDNA-MRD檢測中的一個難點。為了解決上述問題，我們采用腫瘤先驗技術路線，基于固定大panel開發了一種嶄新的生物信息學算法——多變異聯合置信概率分析（MinerVa）。MinerVa采用了目前已知的最大規模的陰性血漿數據庫實現了大panel覆蓋區間內單堿基分辨率的技術噪音基線模型搭建，通過共追蹤變異的整合分析對樣本的MRD狀態進行判定，并在小規模臨床隊列中探索MinerVa算法在ctDNA-MRD檢測中的價值。
1 資料與方法
1.1 MinerVa算法應用方案
MinerVa算法采用腫瘤先驗技術路線，使用固定大panel進行血漿ctDNA-MRD檢測。首先獲取患者腫瘤組織和配對白細胞的測序數據，構建患者的個性化腫瘤變異圖譜；然后在血漿游離DNA的測序數據中，根據腫瘤變異圖譜提取相應的變異信號；對提取到的變異信號，使用噪音基線分布模型進行單變異顯著性分析；在樣本層面上，綜合單變異顯著性結果進行多變異聯合置信概率分析，判定樣本的MRD狀態。
1.2 生信分析算法及模型搭建
1.2.1 腫瘤組織變異檢測
本研究通過配對的腫瘤組織和白細胞測序樣本識別體細胞突變。首先使用BWA（0.7.17）將原始的讀段比對到參考基因組hg19上，然后采用變異檢測軟件VarDict（1.5.1）識別突變（包括SNV、Indel）信息，并得到相應的覆蓋深度和突變支持數目。要求位點深度不低于200 X，SNV類型突變頻率不低于5%，Indels類型突變頻率不低于4%。此外，對于repeat/low complexity等區間易重復出現的假陽性信號，基于30個正常樣本構建了PoN模型（Panel of Normal），過濾檢出頻率低于閾值的突變。
1.2.2 基于MinerVa算法的血漿ctDNA-MRD判定
（1）單堿基分辨率的技術噪音基線模型的搭建：噪音基線模型的建立使用了由1 132例陰性血漿數據構建的數據庫。陰性血漿數據是處理后的腫瘤患者及健康志愿者的血漿數據，腫瘤患者的實測血漿數據扣除了配對組織及白細胞變異，健康志愿者的實測血漿數據則扣除了配對的白細胞變異。針對panel覆蓋范圍內的每個堿基都進行了三個變異方向的SNV基線模型的搭建，對每個SNV，提取基線血漿數據庫中所有未扣除該指定SNV的個體數據進行模型搭建，同時為了排除腫瘤異質性信號以及其他未知來源的變異信號，搭建前對指定SNV的離群值做處理。每個變異的噪音分布模型為混合模型，由無變異人群占比模型以及變異模型兩部分組成，第一部分是無指定SNV變異的人群占比（Pzero），另一部分是對含有指定SNV變異人群的等位基因變異頻率（variant allele frequency，VAF）分布擬合的韋伯分布模型或逆伽馬分布模型。Indel噪音模型的構建不區分具體Indel，只定位堿基坐標，將所有關聯指定堿基坐標的Indel變異合并計算Pzero以及擬合分布模型，方法同SNV模型。
（2）單變異信號相對噪音基線的顯著性水平評估：根據血漿變異的堿基坐標信息，調取該變異的噪音基線模型，通過蒙特卡洛抽樣N=10 000，生成N*Pzero個0，同時由變異模型部分生成N*（1–Pzero）個隨機VAF，分別以這N個VAF作為先驗的噪音頻率，根據二項分布計算患者血漿變異信號計算來自于噪音信號的可能性：
|
其中，VSM為血漿中支持變異的分子數，TSM為覆蓋該堿基的總分子數。為第i次模特卡洛抽樣獲取的先驗噪音頻率。
綜合N次計算結果，進一步計算 的加和平均值 以衡量患者血漿單堿基變異的顯著性程度。P ≤ 0.05作為血漿單變異的陽性判定值。
（3）多變異聯合置信概率分析：該分析檢驗所有P ≤ 0.05的單變異同時是噪音來源的可能性。聯合多個突變，經過多變異聯合置信概率分析計算其聯合置信概率combine_P，combine_P ≤ 0.01則認為樣本MRD陽性。計算公式如下所示，其中，n為追蹤的突變總數，k為陽性變異數目， 為陽性位點的p_value：
|
1.3 驗證使用樣本
本次驗證包括2套陽性標準品（GW_OGTM006、多變異共檢測標準品）和1套陰性標準品（GM12878），在初步的臨床應用驗證中，我們回顧性地對27例非小細胞肺癌患者的樣本和臨床復發結果進行一致性評估。
1.3.1 GW_OGTM006
GW_OGTM006（菁良科技）是商業化標準品（見附件1）。用正常二倍體細胞系GM12878（科佰生物）梯度稀釋，用于熱點靈敏度和特異性分析。在本次性能驗證中使用的GM12878，實驗數據顯示KRAS p.G12C在多次特異性重復中有低頻信號出現，且頻率低于GM12878出廠質控陰性頻率值（0.1%）。通過定制化Panel進行高深度驗證后證實確實存在低頻變異，因此后續數據分析未納入KRAS p.G12C位點。另外，本方法主要評估的是常見的熱點插入缺失變異，一般在20 nt以內，而MET Exon 14 Skipping是一個長度140 nt的缺失變異，并不在本方法的分析范圍內。因此，后續在熱點變異靈敏度及特異性數據分析中使用的是去除了KRAS p.G12C及MET Exon 14 Skipping之外的12個熱點變異。
1.3.2 多變異共檢測標準品
由于常規標準品的變異數量較低，不能夠滿足多變異共追蹤的MRD分析性能驗證，因此我們使用差異化的單核苷酸多態性（single nucleotide polymorphism，SNP）信號模擬ctDNA的變異，以提高可追蹤的變異的數量。方案為通過將已知SNP信息的隨機一個個體（以個體S指代）的白細胞DNA用GM12878進行梯度稀釋，得到一系列可用于整體分析性能驗證的樣品。將個體S特有的（不同于hg19以及GM12878）的100個SNP變異納入陽性變異集（見附件2），將個體S的白細胞DNA和細胞系GM12878 DNA中與參考基因組hg19具有相同基因型的547個SNP位點納入陰性變異集（見附件3）。
1.3.3 臨床樣本
1 132例基線樣本包含以非小細胞肺癌為主的泛癌種患者的腫瘤組織、血液白膜層和血漿樣本，以及健康人的血漿白膜層和血漿樣本，來自四川大學華西醫院、四川省人民醫院等多個中心的臨床隊列，以及臻和科技公司的健康志愿者。
27例接受根治性手術治療的Ⅰ～Ⅲ期非小細胞肺癌患者在四川大學華西醫院納入，其中Ⅰ期7例，Ⅱ期14例，Ⅲ期6例（詳見附件4），采集術中腫瘤組織樣本及多個時間點的血液樣本，對患者進行為期30個月的隨訪信息收集。臨床樣本收集和使用已通過四川大學華西醫院倫理委員會批準，患者提供書面知情同意。
1.4 驗證實驗流程
1.4.1 GW_OGTM006標準品實驗流程
首先使用正常二倍體細胞系GM12878基因組對陽性標準品GW_OGTM006進行6個梯度的稀釋，按照熱點變異的理論變異頻率均值由高到低分別為1%、0.5%、0.3%、0.1%、0.05%和0.03%，這6個梯度的樣本被分別命名為PC-1P、PC-05P、PC-03P、PC-01P、PC-005P、PC-003P。
對稀釋好的6個梯度的DNA樣本使用Covaris超聲波DNA破碎儀（M220，Covaris，美國）進行片段化，取33 ng片段化的DNA使用 KAPA Hyper Preparation Kit（KK2621，羅氏，瑞士）試劑盒進行文庫構建，構建好的文庫使用xGEN lockdown（1080584，IDT，美國）探針捕獲法進行目標區域富集，使用MGISEQ-T7（華大科技，中國）進行PE100雙端測序，下機后平均測序深度 ≥ 4 500 X，熱點區域深度 ≥ 30 000 X。每個樣本進行了4次技術重復。
1.4.2 多變異共檢測標準品實驗流程
將個體S的白細胞DNA用GM12878進行9個梯度的稀釋，得到一系列可用于整體分析性能驗證的多變異共檢測標準品。預期的理論變異頻率均值由高到低分別為1%、0.5%、0.3%、0.1%、0.05%、0.03%、0.01%、0.005%和0.003%。建庫富集方法同GW_OGTM006標準品，同樣每個樣本進行了4次技術重復。
1.4.3 27例臨床樣本及1 132例基線樣本實驗流程
術中組織標本及白膜層的提取使用天根血液/組織/細胞基因組提取試劑盒（DP304-02，天根，中國），血漿樣本的提取使用MagMAX?游離DNA分離試劑盒（A29319，Applied Biosystems，美國），建庫富集方法同前述標準品。
2 結果
2.1 熱點變異檢測性能
2.1.1 熱點變異檢測靈敏度及檢出限評估
標準品GW_OGTM006統計了4次技術重復中各個變異的檢出情況（見表1），綜合所有的熱點變異（包括SNV和Indel）進行分析。在33 ng起始建庫條件下，預期稀釋變異頻率為0.1%的標準樣本PC-01P中，熱點變異的檢測靈敏度為98%；即使在預期稀釋變異頻率為0.05%的標準樣本PC-005P中，熱點變異的檢測靈敏度仍可達85%。通過probit回歸方法進行檢出限（limit of detection，LoD）評估（見圖1），對于肺癌常見的熱點驅動變異EGFR L858R和19DEL，LoD可達0.055%；綜合標準品中所有的EGFR變異，LoD為0.074%；考慮標準品中12個熱點變異，LoD為0.178%。
表1
GW_OGTM006各熱點變異的檢出率
Table1.
Detection rate of hotspots in GW_OGTM006
圖1
GW_OGTM006熱點變異LoD
Figure1.
Limit of detection of hotspots in GW_OGTM006
2.1.2 熱點變異檢測特異性
使用正常細胞來源的細胞系GM12878基因組DNA進行熱點檢測特異性評價。將GM12878的基因組DNA打斷作為熱點變異特異性檢測樣本，進行8次獨立的建庫捕獲和測序分析，統計在所有樣本中的12個熱點變異的檢出情況，用以衡量熱點變異的檢出特異性。
在驗證中共進行了8*12 = 96次變異檢測，96次檢出結果均為陰性，對熱點變異的檢出特異性為100%（95% CI：96.15%～100%）。
2.2 非熱點變異的檢測性能
分析多變異共檢測標準品的陽性變異集的100個變異作為非熱點變異，統計其在各稀釋梯度的檢出率，通過probit回歸分析獲得最低檢出限。
2.2.1 非熱點變異檢測的靈敏度及檢出限
在33 ng起始建庫的條件下，不同稀釋梯度的多變異共檢測標準品變異的檢出率見附件5，probit回歸分析表明，非熱點單變異的最低檢出限為0.353%（見圖2）。
圖2
非熱點單變異LoD
Figure2.
Limit of detection of non-hotspots
2.2.2 非熱點變異檢測的特異性
陰性變異集中共納入了547個SNP位點，這些多態性位點在個體S的白細胞DNA和細胞系GM12878 DNA中與參考基因組hg19具有相同基因型。對9個梯度的標準品分別進行4次重復獨立建庫和測序分析，在對SNP位點的547*4*9 = 19 692次的變異檢測中，僅有45個假陽性判定結果，非熱點單變異的檢測特異性為99.77%（95% CI：99.69%～99.83%）。
2.3 多變異追蹤時樣本的檢出性能
在多變異共檢測標準品系列樣本中，我們利用差異化的SNP位點分別設置陽性變異集和陰性變異集，用以評估MinerVa算法在多變異追蹤時的檢出靈敏度及特異性。
2.3.1 多變異共追蹤時樣本的檢出靈敏度
模擬先驗的腫瘤變異圖譜，設定變異追蹤數目為1、2、3、6、10、20、30、40和50。對每個指定的變異數目，隨機抽取陽性變異集中的變異進行組合，每個梯度1 000次隨機組合，即形成了1 000組獨立的腫瘤先驗變異圖譜，分別按照每一組圖譜進行MRD檢測，最終統計每個稀釋梯度樣本的這1 000個組合的陽性檢出率（見表2）。結果顯示隨著建庫起始量的升高，變異追蹤數目的增多，檢出靈敏度顯著提高。
表2
追蹤不同數目變異時的檢出靈敏度（%）
Table2.
Sensitivity (%) for tracking different number of variants
進一步經過probit回歸，統計出在33 ng起始建庫的條件下，不同追蹤變異數目時的95檢出限（見圖3）。結果顯示隨著追蹤位點數增加，LoD逐漸降低，當追蹤30個變異時，可檢測低至6.3 × 10?5變異豐度的變異信號。
圖3
追蹤不同數目變異時的LoD
Figure3.
Limit of detection for tracking different number of variants
2.3.2 多變異共追蹤時樣本的檢出特異性
同樣設定變異追蹤數目為1、2、3、6、10、20、30、40和50，對每個指定的變異數目，隨機抽取陰性變異集中的變異進行組合，每個梯度1 000次隨機組合，綜合1 000次隨機組合的MRD判定結果，評估不同追蹤變異個數下的檢測特異性（見表3）。當追蹤多個變異時（變異數目：2、3、6、10、20、30、40、50），特異性穩定在99.62%～99.70%。
表3
追蹤不同數目變異時的特異性
Table3.
Specificity for tracking different number of variants
2.4 肺癌隊列ctDNA-MRD檢測性能
27例非小細胞肺癌患者中，14例患者在隨訪中出現了復發，復發患者的中位無病生存期（disease-free survival，DFS）為337天（166～632天），13例患者未復發。13例未復發的患者在術后動態ctDNA監測中，ctDNA檢測結果均為陰性，特異性為100%（95% CI：77.19%～100%）（見附件6）。14例復發的患者在術后一個月ctDNA檢測陽性的比例為35.7%（5/14），在動態ctDNA監測中，11例患者ctDNA陽性，靈敏度達78.6%（95% CI：52.41%～92.43%）。該隊列研究表明使用MinerVa算法分析的ctDNA-MRD狀態與臨床復發具有較高的一致性。
3 討論
使用固定大panel的MRD檢測方法在臨床研究中已顯示了其可行性[13]。固定panel通過充分的前期驗證，保證了操作的穩定性和臨床的成功率，同時也保證了檢測時效性，可滿足在關鍵時間節點為制定治療決策提供參考的需要。
靈敏度不足是MRD檢測的主要技術挑戰。用藥熱點突變檢測的常規方法是以支持reads數或變異頻率為固定閾值，但MRD檢測中包括了大量個性化的變異，簡單套用這種“一刀切”的方法忽略了不同變異噪音水平的差異，無法平衡不同變異的靈敏度和特異性，從而影響MRD判定的準確性。細化不同變異的技術噪音是合理設定閾值、提升檢測性能的關鍵。既往的研究也嘗試構建技術噪音模型，但由于只使用了少量健康人的數據，差異分層只能集中在變異亞型（subtype，如C > G、C > A等）或者三堿基組合（context）層面[14-15]。Capp-Seq研究嘗試建立單堿基分辨率的噪音模型，但由于只使用了12例健康人血漿，模型構建并不充分，只作為簡單過濾工具使用[16-17]。MinerVa算法創新性地將腫瘤患者及健康志愿者的血漿數據扣除配對組織和白細胞的變異，構建了陰性血漿數據，同時利用固定panel累積數據的優勢，使用超千例的陰性血漿數據，對panel覆蓋區域的每個單堿基都構建了技術噪音分布模型；基于模型，只有高于噪音基線的變異才會被判定為陽性。分析性能數據顯示，EGFR L858R和19DEL的LOD低至0.055%，顯著低于常規熱點變異的檢出限，而非熱點單變異的檢測LoD也達到了0.353%。
樣本水平MRD陽性判斷值目前并沒有統一的標準。常用的方法要求滿足大于等于1個或2個ctDNA單變異檢出[18-19]，但該方法在追蹤變異數量較多時，檢測特異性會下降。假設每個單變異的檢測特異性是99.6%，在追蹤50個變異、使用 ≥ 1個變異的陽性判定條件時，樣本水平特異性低于85%，而使用比較嚴苛的 ≥ 2的陽性變異作為陽性判斷條件，特異性也跌至99%（見圖4），且損傷了檢測的靈敏度。MinerVa算法開發了多變異聯合置信概率分析，將每個變異源于噪音的概率進行整合，通過控制全部血漿變異都來源于技術噪音的概率，保證了樣本檢測的特異性。分析性能驗證顯示，無論追蹤2個變異還是追蹤50個變異，樣本檢測特異性穩定在99.62%～99.70%之間。
圖4
不同MRD陽性判定標準的檢測特異性
Figure4.
Specificity of different MRD positive determination thresholds
本研究也具有一定的局限性。首先，基于固定panel的MRD檢測由于覆蓋的區間依然有限，對一些低突變負荷癌種或患者的變異數量覆蓋不足，可能會影響MRD檢測靈敏度。此外，本研究僅采用小的臨床隊列對MinerVa的臨床價值進行了評估，仍需通過更大樣本的臨床實踐探索本算法在不同癌種及不同臨床場景（如術后復發監測、新輔助/輔助治療療效評估等）中的適用性。
重要聲明
利益沖突聲明：楊飄、張亞晰、范銳、陳維之均是MinerVa相關專利的共同發明人，其他作者均聲明不存在利益沖突。
作者貢獻聲明：楊飄負責分析性能實驗設計、實驗開展，以及論文撰寫；夏粱、梅建東、劉倫旭負責臨床隊列設計、臨床信息收集和臨床數據分析；張亞晰、范銳、黃宇、陳維之負責算法開發和生信數據分析；陳維之負責整體實驗設計。
倫理聲明：本研究通過了四川大學華西醫院生物醫學倫理分委會的審查審批（批文編號：2017年 審（34）號）。
本文附件見本刊網站的電子版本（biomedeng.cn）。
