[發明專利]用于計算癌癥樣本純度和染色體倍性的方法和裝置

有效

• 申請號: 201710312237.7
• 申請日: 2017-05-05
• 公開（公告）號: CN108804876B
• 公開（公告）日: 2022-03-15
• 發明（設計）人: 黃宇;羅志輝;蘇瑤;范新平
• 申請（專利權）人: 中國科學院上海藥物研究所
• 主分類號: G16B20/10
• 分類號: G16B20/10;G16B30/00;G16B40/00;G16B50/30
• 代理公司: 北京金信知識產權代理有限公司 11225
• 代理人: 張皓;徐琳

1.一種用于計算癌癥樣本中癌癥細胞純度和染色體倍性的方法，所述方法包括以下步驟：

步驟A：

獲取配對的癌癥組織樣本和正常組織樣本的全基因組測序數據，并將測序數據比對到參考基因組；

步驟B：

從步驟A得到的比對結果文件中，提取read位置和長度信息，HGSNV位點和覆蓋該位點的read數量信息，計算所有HGSNV的MAF，其中，read是指高通量測序平臺產生的測序序列；HGSNV是指雜合生殖系細胞單堿基變異；MAF表示主要等位基因分數，指兩種等位基因分數中的較大值，等位基因分數的計算方法為覆蓋某一等位基因的read數量除以覆蓋該位點總read數量的比值；MAF計算公式如(1.1)所示：

公式(1.1)中，n^r為包含與參考基因組相同等位基因的read數量，n^a為包含另一種等位基因的read的數量，n^t表示覆蓋該HGSNV位點的總read數量，C為該HGSNV的MAF值；

步驟C：

根據步驟B得到的read位置和長度信息，以window為單位統計各window內包含的read數量，使用基因組GC含量校正所有window內read數量，其中，window是指按照一定長度劃分的基因組片段，該長度代表window大小；

步驟D：

使用步驟C校正后的read數量，使用公式(1)計算每一個window的TRE，然后運用TRE，通過BIC-seq軟件對基因組進行片段化，獲得以拷貝數劃分的基因組片段，其中，TRE表示癌癥片段讀長富集程度e_s，指癌癥樣本中某一片段s內read數量與相應正常樣本中對應片段read數量的比值，定義公式如下：

公式(1)中，和分別表示在癌癥樣本中覆蓋片段s的read數量和在正常樣本中覆蓋片段s的read數量，N_t表示癌癥樣本總read數量，N_n表示相應正常樣本總read數量，e_s為TRE值；

步驟E：

以步驟D中BIC-seq處理后的基因組片段為單位，統計片段內所有window的TRE的均值、方差和該片段內window數量，根據均值和方差對基因組每個片段的window數量進行平滑化處理，使TRE的分布更均勻，然后將平滑化處理后所有片段的window分布匯總，得到基因組上window隨TRE變化的分布結果；同時以片段為單位，計算片段中所有HGSNV的MAF的均值和標準差；

步驟F：

使用公式(12)、(13)所示的類自回歸模型，計算兩個相鄰拷貝數片段對應的peak的TRE的差值即P，其中，peak是指基因組所有window的TRE分布中，聚集在一起的TRE簇；以及其中，遍歷一定范圍的P，計算Y(P)，在Y(P)的分布中，選擇第二高峰內Y(P)的最大值對應的P作為P的計算結果：

公式(12)和(13)中，X_t表示0到M_t之間的TRE值；t表示擴大了1000倍的TRE值；M_t表示TRE的最大值；變量P表示兩個相鄰拷貝數片段對應的peak的TRE的差值；C(X_t)表示在TRE為X_t的位點，對應的window數量；C(X_t+1000×P)表示在TRE為X_t+1000×P的位點，對應的window數量；Y(P)表示在變量P下，類自回歸模型的函數值；

步驟G：

根據步驟F得到的P，計算TRE分布中第一個實際觀測peak的TRE均值，然后計算在第一個實際觀測peak之前最多可能存在理論peak的數量N，最后當第一個實際觀測peak之前存在的理論peak數量為n時，計算Q的值，以Q_n表示，其中步驟G包括：

G1：

根據步驟F計算的P，使用公式(13.1)，選取使公式(13.1)取最大值的X_f作為第一個實際觀測peak的TRE均值：

公式(13.1)中，i表示第i個peak，C(X_f+P×i)表示在TRE為X_f+P×i的位點，對應的window數量，n_m表示M_t以內peak的最大數量，M_t表示TRE的最大值；

G2：

使用公式(13.2)，根據步驟F計算的P和步驟G1計算的X_f，計算在X_f之前最多可能存在的peak數量N：

公式(13.2)中，X_f表示第一個實際觀測peak的TRE均值，P表示兩個相鄰拷貝數片段對應的peak的TRE的差值，floor表示向下取整數；

G3：

利用步驟G2計算的N值，當n取0到N之間的整數時，使用公式(13.3)計算Q_n的值：

Q_n＝X_f-n×P+2×P＝X_f+(2-n)×P，n∈[0，N] (13.3)

公式(13.3)中，n表示第一個實際觀測peak之前存在的理論peak數量，取值范圍是0到N之間的整數，P表示兩個相鄰拷貝數片段對應的peak的TRE的差值，X_f表示第一個實際觀測peak的TRE均值，Q_n表示在X_f之前理論上存在n個peak時的Q值；

步驟H：

使用步驟F計算的P與步驟G計算的Q_n，使用公式(10)、(11)計算癌癥樣本純度γ和染色體倍性κ：

公式(10)、(11)中，γ表示樣本純度，κ表示染色體倍性，由此對(P，Q_N)得到對應的(γ，κ)；

步驟I：

當n取[0，N]之間的某個整數值時，使用公式(13.4)計算第i個peak的TRE均值：

T_i＝X_f-n×P+i×P＝X_f+(i-n)×P，n∈[0，N] (13.4)

公式(13.4)中，n表示第一個實際觀測peak之前存在的理論peak數量，取值范圍是0到N之間的整數，P表示兩個相鄰拷貝數片段對應的peak的TRE的差值，X_f表示第一個實際觀測peak的TRE均值，T_i表示第i個peak的TRE均值，

對于落在T_i附近的片段，認為該片段具有拷貝數i；對于沒有落在T_i附近的片段，將其歸類為亞克隆片段，在后續分析中剔除所有亞克隆片段；然后根據步驟H計算的癌癥樣本純度γ和peak對應的拷貝數，計算peak的MAF的期望f_b，不同peak的MAF期望不同，對基因組上的所有peak，最終得到MAF期望的集合{f_b}；同時計算各個peak的TRE均值和標準差；

步驟J：

根據步驟F計算的P和步驟I計算的{f_b}構建如公式(19)所示的用“貝葉斯信息準則”校正后的最終的混合高斯模型，然后對最終的混合高斯模型極大似然估計；其中，步驟J包括如下幾步：

J1：

以步驟F計算的P構建如公式(17)所示的高斯分布模型：

公式(17)中，L(e_s；γ，κ)表示基因組片段TRE的似然函數，N_w表示基因組上的所有window的數量，I表示基因組中所有片段的最大的拷貝數，σ_i表示拷貝數為i的所有片段的TRE的標準差由步驟I得到，e_s為第s個window的TRE觀測值，Sⁱ表示第i個peak的TRE均值即步驟I中的T_i，p_i表示第s個window的拷貝數為i的權重，對所有的i，p_i均取值為1；

J2：

以步驟I計算的f_b構建如公式(18)所示的第一混合高斯模型：

公式(18)中，L(f_s；γ，κ)表示HGSNV的似然函數，M表示基因組中所有HGSNV數量，S表示第S個HGSNV，I表示基因組中所有片段的最大的拷貝數；F^i，j表示拷貝數為i，主要等位基因的拷貝數為j的片段內HGSNV的MAF期望值，由步驟I得到；f_s表示該片段內所有HGSNV的MAF的均值，由步驟E得到；σ_i，j表示該片段內所有HGSNV的MAF的標準差，由步驟E得到；p_i，j表示在主要等位基因的拷貝數為j時，高斯分布的權重，對所有的i和j，p_i，j取值均為1，p_i表示第S個HGSNV所在片段的拷貝數為i的權重，對所有的i，p_i取值均為1；

J3：

將(17)與(18)相加得到第二混合高斯模型，然后對第二混合高斯模型進行BIC(Bayesian Information Criterion)校正得到最終的混合高斯模型如公式(19)：

BIC(e_s，f_s；γ，κ)＝-2×logL(f_s；γ，κ)-2×logL(e_s；γ，κ)+I×log(N_w)+J×log(M) (19)

公式(19)中，BIC(e_s，f_s；γ，κ)表示最終的混合高斯模型的似然函數，I表示基因組中所有片段的最大的拷貝數，J是公式(18)中j的取值個數，N_w是基因組中window的數量，M是基因組中HGSNV的個數，

對[0，N]范圍內的每一個整數值n，通過步驟G得到Q_n，或者通過步驟I得到所有peak的MAF期望的集合{f_b}，由一對(P，{f_b})構建一個公式(19)所示的模型；

步驟K：

以0.001為分辨率，對[P-m，P+m]區間的所有P值，重復步驟G～J，得到一系列不同的(P，Q_n)與對應的似然函數值，取最大的似然函數值對應的(P，Q_n)作為最合適的P和Q值，m是0到0.5之間的一個值；

步驟L：

查詢步驟H的結果，找到在步驟K得到的(P，Q)下，對應的癌癥樣本純度和染色體倍性。