Cnv檢測方法和裝置的制造方法
【專利摘要】本發明提供了一種檢測CNV的方法，包括：獲取目標個體的基因組測序結果；將所述測序結果與參考序列比對，獲得比對結果，所述參考序列包括多個窗口；基于所述比對結果，計算每個窗口的初始比對率，窗口的初始比對率＝比對上所述窗口的讀段數目/比對上所有窗口的讀段數目的平均值，所述比對上所有窗口的讀段數目的平均值＝比對上所有窗口的讀段總數/窗口個數；合并初始比對率無顯著差異的多個相鄰窗口，定義合并后的多個相鄰窗口為一次區域，剩余的每個獨立窗口也分別稱為一次區域；基于所述一次區域的比對率與預定比對率不相等，判定所述一次區域存在CNV。
【專利說明】
GNV檢測方法和裝置
技術領域
[0001] 本發明涉及生物信息領域，具體的，本發明涉及檢測CNV的方法和裝置。
【背景技術】
[0002] 單細胞測序技術是利用二代測序技術對單個細胞的微量核酸進行測序。該項技術 主要包括單細胞分離，單細胞核酸的提取與擴增和測序三部分。單細胞測序作為一項革命 性技術，近幾年來在科研及生物醫藥領域得以廣泛應用。例如，對腫瘤單細胞進行測序，揭 示腫瘤單細胞層面的異質性，推演腫瘤的演化過程；無創產前診斷；組裝不能培養的微生 物基因組；痕量細胞（可以運用于法醫學等）基因組的獲取；單細胞技術也被引入到胚胎 植入前診斷等等。單細胞測序技術解決了痕量細胞基因組獲取的難題，為疾病發病機制和 診斷學研究提供新的方法。
[0003] 在單細胞研究中，單細胞拷貝數變異（Copy Number Variants, CNV)扮演了很重要 的角色。CNV-般指染色體上大于1Kb的片段發生丟失或重復的現象。CNV是一種廣泛存 在動植物基因組中的遺傳多態，它的突變頻率遠高于SNP，基因組研究已證明CNV與部分人 類疾病相關，比如與腫瘤、肥胖、自閉癥、自身免疫性疾病和系統性紅斑狼瘡等多種復雜疾 病相關。在腫瘤的異質性和進化研究中，腫瘤單細胞的CNV的檢測，通過比較單細胞之間， 以及單細胞與對應組織的CNV的差異，揭示了腫瘤在單個細胞層面的異質性，為腫瘤演化 推演提供了依據；無創產前診斷，則需要對微量DNA進行檢測是否有染色體非整倍體變異 (CNV的一種）而導致的21-三體綜合征（47, +21)、18-三體綜合征（47, +18)、13-三體綜 合征（47,+13)等；胚胎植入前的診斷和篩查，需要對單個生殖細胞或胚胎細胞進行相關檢 測分析；法醫取證樣本（痕量的血液，精液等），需要進行痕量細胞的分析等等。綜合來看， 當前生物醫學領域對于痕量細胞，甚至單個細胞的大片段CNV的檢測提出了需求和挑戰。
[0004] 現有的CNV檢測方法大多是針對組織測序數據的，如CNV-seq，PenCNV，CNAseg和 Readdepth等。單細胞測序數據，尤其是低深度測序數據，具有低基因組覆蓋度和高擴增偏 向性，在基因組的不同區域短序列比對波動很大，這些CNV檢測方法并不太適合對單個細 胞拷貝數變異的檢測。

【發明內容】

[0005] 本發明旨在至少解決上述問題至少之一或者提出至少一種商業選擇。
[0006] 依據本發明的一方面，本發明提供一種檢測CNV的方法，所述方法包括以下步驟： 獲取目標個體的基因組測序結果，所述測序結果包括多個讀段；將所述測序結果與參考序 列比對，獲得比對結果，所述參考序列包括多個窗口，所述比對結果包括比對上每個所述窗 口的讀段的數目；基于所述比對結果，計算每個窗口的初始比對率，窗口的初始比對率=比 對上所述窗口的讀段數目/比對上所有窗口的讀段數目的平均值，所述比對上所有窗口的 讀段數目的平均值=比對上所有窗口的讀段總數/窗口個數；合并初始比對率無顯著差異 的多個相鄰窗口，定義合并后的多個相鄰窗口為一次區域，剩余的每個獨立窗口也稱為一 次區域；基于所述一次區域的比對率與預定比對率不相等，判定所述一次區域存在CNV，所 述一次區域的比對率為所述一次區域包含的窗口的初始比對率的均值，所述預定比對率為 所有窗口的比對率中頻率最高的窗口的比對率，所述窗口的比對率為其所在的一次區域的 比對率。在本發明的一個實施例中，所述基因組獲自所述目標個體的單個細胞。通過構建 單細胞的基因組測序文庫，并對所述測序文庫進行序列測定獲得所述測序結果。任選的， 構建所述測序文庫包括對所述基因組進行簡并寡核苷酸引物PCR，多重置換擴增和/或多 次退火環狀循環擴增，以獲得足夠建庫的核酸量和/或足夠上機測序的核酸量。序列測定 可以利用現有測序平臺，包括但不限于CG (Complete Genomics)、Illumina/Solexa、Life Technologies ABI SOLiD和Roche 454測序平臺，可根據所選用的測序平臺進行相應的測 序文庫制備，可選擇單端或雙端測序，由此獲得的測序結果由多個短序列組成，將各個短序 列稱為讀段。所說的比對可以利用已知比對軟件進行，比如利用Bowtie、SOAP、BWA和/或 TeraMap等進行。在本發明的一個實施例中，只利用所述比對結果中的比對到所述參考序列 唯一位置的讀段進行比對率的計算，以提高數據準確性以利于提高CNV檢測的準確性。
[0007] 所稱的窗口可以預先確定，也可以在進行目標個體檢測時同時確定。在本發明的 一個實施例中，窗口是預先確定的。所述窗口的確定包括：將短序列集與參考序列比對， 確定比對上所述參考序列的短序列的起始位置，所述短序列集包括多個短序列；在所述參 考序列上劃定窗口，使每個所述窗口包含相同數目的所述起始位置，任選的，所述窗口之 間沒有重疊。在比對過程中，根據比對參數的設置，一條短序列最多允許有m個堿基錯配 (mismatch)，m優選為1或2,若一條短序列中有超過m個堿基發生錯配，則視為該短序列無 法比對到參考序列。所稱的起始位置為比對上參考序列的各個短序列的起始堿基與參考序 列的匹配位置，當有多個短序列的起始堿基比對到參考序列同一位置時，只記錄一次，即記 錄所說的起始位置為一個。這里，所稱的短序列的起始堿基，即涉及的短序列的方向，是以 參考序列的方向為參照的，例如，將一條短序列上的匹配到參考序列最前位置（位置編號 最小）的堿基稱為該短序列的起始堿基。使每個所述窗口包含相同的起始位置數目，而不 限制其包含的沒有讀段匹配的位點數目，所以一般的各個窗口的大小不一樣，這樣，有利于 減少單細胞基因組擴增帶來的偏向性。據此，在這一構思下，也可以使每個窗口包含相同數 目的特定位置來進行窗口劃定，所說的特定位置為比對上參考序列的各個短序列的相同位 置堿基與參考序列的匹配位置，例如，使所說的特定位置為比對上參考序列的各個短序列 的末端堿基與參考序列的匹配位置。
[0008] 所稱短序列集可來自模擬序列集和/或測序結果，這里所說的測序結果可以是自 己測定的人核酸的測序數據，也可以是他人公開的核酸樣本的測序結果，核酸可以是基因 組DNA也可以是游離DNA。較佳的，使所說的模擬序列集中的模擬序列在比對到參考基因組 上能有相對均勻的分布，在本發明的一個實施例中，模擬序列可以這樣獲得：從所述參考序 列的長度為Q的染色體的一端的堿基開始，拷貝所述染色體的P個堿基，以獲得第一條模擬 序列，沿所述染色體的另一端方向移動一個堿基拷貝所述染色體的P個堿基，以獲得第二 條模擬序列，沿所述染色體的另一端方向移動兩個堿基拷貝所述染色體的P個堿基，以獲 得第二條模擬序列，依此獲得第Q _P+1條模擬序列，所述第Q_P+1條模擬序列的末端堿基與 所述染色體的另一端的堿基重合，其中，P為模擬序列的長度，較佳的，P多10。在本發明的 一個實施例中，所述窗口之間沒有重疊且所述窗口總數不大于100, 〇〇〇。窗口的大小的設置 可以基于CNV檢測精度調整，在人參考基因組大小確定的情況下，窗口的大小與窗口數目 成反比。在該實施例中，窗口的總數不少于10, 〇〇〇且不大于100, 〇〇〇,而且之間沒有重疊， 利于準確檢測出一般定義的不小于1K的CNV。
[0009] 在本發明的一個實施例中，所述目標個體為人類，人類為二倍體生物，其染色體組 數為2,優選對應的所述參考序列為人參考基因組的至少一部分，例如為HG19，HG19可以從 NCBI數據庫獲取，或者為所有窗口構成的參考序列。在本發明的另一個實施例中，以N替 代所述人參考基因組的Y染色體的擬常染色體區的每個堿基，N表示A、T、C和G中的任一 種，這樣，有利于避免性染色體的擬常染色體區域CNV檢測的假陽性。
[0010] 在本發明的一個實施例中，在合并初始比對率無顯著差異的多個相鄰窗口之前， 利用比對率-GC含量的關系對每個所述窗口的初始比對率進行GC校正，獲得各個窗口的校 正比對率，以消除或減少GC含量對測序結果、比對率的影響，并且以窗口的校正比對率替 代所述窗口的初始比對率進行后續檢測，例如，所述一次區域的比對率變為其包含的所有 窗口的校正比對率的均值，而在確定所稱的預定比對率時，將一次區域的比對率賦給其所 包含的窗口，即處于同一一次區域的各個窗口的比對率均相等，為其所在的一次區域的比 對率，這樣，對所有窗口的比對率進行計數，確定各個比對率出現的次數，將出現次數最多 即頻率最高的窗口比對率定為所稱的預定比對率。所述比對率-GC含量的關系可以預先利 用對照樣本的測序數據來建立、保存，用以校正各個待測樣本測序結果，優選的對照樣本為 與目標個體同物種的組織樣本，也可以在檢測目標樣本時同時利用目標樣本基因組的測序 結果來建立。在本發明的一個實施例中，直接利用檢測所需的目標樣本的測序結果來建立 比對率-GC含量的關系，所說的比對率-GC含量的關系的建立如下：獲得至少一個樣本的核 酸的測序數據，所述測序數據由多個讀段組成；將所述測序數據與參考序列進行比對，獲得 比對結果，所述參考序列包括多個窗口，所述比對結果包括比對上每個所述窗口的讀段的 數目；計算每個所述窗口的初始比對率，窗口的初始比對率=比對上所述窗口的讀段數目 /比對上所有窗口的讀段數目的平均值，所述比對上所有窗口的讀段數目的平均值=比對 上所有窗口的讀段總數/窗口個數；基于多組的窗口的初始比對率和該窗口的GC含量的數 值，利用二維回歸分析法建立所述比對率-GC含量的關系。在本發明的一個實施例中，利用 的二維回歸分析法為局部加權回歸散點平滑法（Lowess)。
[0011] 所稱的初始比對率無顯著差異的多個相鄰窗口指初始比對率無實質差異的相鄰 窗口，例如，由于初始比對率或者校正比對率是圍繞著"1"波動的一組數值，可以以1或 者以1± 1*10%為有無實質差異的界限，如相鄰的、校正比對率都在0.9以下，或0.90~ 1. 10,或1. 10以上的窗口為無實質差異的窗口。在本發明的一個實施例中，所述合并初 始比對率無顯著差異的多個相鄰窗口為合并符合以下描述的相鄰窗口，多個相鄰窗口 的校正比對率都大于1或者都小于1。進一步的，為確定所檢測的CNV的大小和確切發 生位置（斷點），該方法還包括：確定所述一次區域中的二次區域，包括，（1)基于公式
算所述一次區域中的子區域Μ與該一次區域中的所有其它窗口的比對 率的差異，獲得所有Zi j，取4= max i i j~| |，（2)將Ζ。與第一臨界值比較，當Ζ。超過 第一臨界值時，相應的子區域Μ為所述二次區域，所述二次區域為CNV區域，所述二次區域 的邊界即為CNV的發生位置，（3)去除所述一次區域中的二次區域，更新i、j和n，進行步 驟（1)和（2)，直至無 Z。超過第一臨界值；其中，i和j為所述一次區域中的窗口的編號，η 為所述一次區域中的窗口的數目，所述子區域Μ為所述一次區域中的第i+1個窗口到第j 個窗口之間的區域，氏為所述一次區域中第i個窗口的校正比對率，所述第一臨界值是Z ^ 分布中的第一預定概率的概率密度，所述第一預定概率多95%，1 < i < j < mSii L+… +民，S,= R片…+R,，Sn= R彳…+Rn。假設子區域Μ為正常非變異區域，Zy分布指Z u服從標 準正態分布，第一預定概率和第一臨界值一一對應，一般統計書籍都包含第一預定概率和 第一臨界值對應的表格供查閱。在本發明的一個實施例中，當Z。落入拒絕域，即Z。超過第 一預定概率例如為99. 9%對應的第一臨界值，可知發生了小概率事件，否定原假設，即子區 域Μ為變異區域。上述過程，依據窗口的校正比對率的同向性，即都大于1或者都小于1， 對窗口進行合并，獲得大的一次區域，再在各個一次區域中進行循環判斷以確定其中的CNV 的發生邊界，即從其中確定二次區域，這樣同時在多個一次區域中并行確定二次區域，利于 快速檢測CNV。在本發明的一個實施例中，所述二次區域的比對率為所述二次區域包含的所 有窗口的校正比對率的均值。在本發明的一個實施例中，該方法還包括，基于比較所述二次 區域的比對率與所述預定比對率的大小，判定CNV的類型，其中包括，當所述二次區域的比 對率大于所述預定比對率時，判定所述二次區域為拷貝數增加區域，當所述二次區域的比 對率小于所述預定比對率時，判定所述二次區域為拷貝數減少區域。在本發明的另一個實 施例中，利用以下公式計算所述二次區域的拷貝數，二次區域的拷貝數=該二次區域的比 對率/預定比對率*目標個體的染色體組數，所述二次區域的比對率為其包含的所有窗口 的校正比對率的均值。
[0012] 無顯著差異，也可以指統計學上的對數據差異性的評價一一差異無顯著性，例如 設定預定概率，通常預定概率可以設為不小于95 %，對相鄰多個窗口的校正比對率進行統 計檢驗，例如可以利用ζ檢驗或t檢驗，多個校正比對率之間的差異無顯著性（ρ > 0. 05)， 即認為達到所說的無顯著差異。在本發明的一個實施例中，所述合并初始比對率無顯著差 異的多個相鄰窗口為合并滿足如下描述的相鄰窗口一一校正比對率的差異無統計意義，使 合并得的一次區域為CNV區域。合并初始比對率無顯著差異的多個相鄰窗口具體包括：(a)
計算區域N與其它所有窗口的比對率的差異，獲得所有Zxy， 取Zb= max i < x < y J Zxy |，（b)將Zb與臨界值比較，當Z b超過所述臨界值時，相應的區域N 為所述一次區域，（c)將所述一次區域去除，更新x、y和w，進行步驟（a)和（b)，直至無 Zb 超過所述臨界值，其中，X和y為窗口的編號，w為窗口總數，所述區域N為第x+1個窗口到 第y個窗口之間的區域，Rx為第X個窗口的校正比對率，所述臨界值是Z xy分布中的預定概 率的概率密度，所述預定概率彡95%，1彡X < y彡w，Sx=心+…+艮，Sy= RJ…+Ry，Sw= ^+???+1。所說的Zxy分布為Z xy服從標準正態分布，預定概率與臨界值一一對應。在本發 明的一個實施例中，假設區域N為正常非變異區域，當Z b落入拒絕域，即Z b超過預定概率例 如為99. 9%對應的臨界值，可知發生了小概率事件，否定原假設，即區域N為變異區域。上 述過程，基于對所有窗口進行循環判斷確定CNV的發生邊界，確定出的一次區域即為CNV區 域。在本發明的一個實施例中，該方法還包括：基于比較所述一次區域的比對率與所述預 定比對率的大小，判定所述CNV的類型，其中包括，當所述一次區域的比對率大于所述預定 比對率時，判定所述一次區域為拷貝數增加區域；當所述一次區域的比對率小于所述預定 比對率時，判定所述一次區域為拷貝數減少區域。在本發明的另一個實施例中，該方法還包 括：利用以下公式計算所述一次區域的拷貝數，一次區域的拷貝數=該一次區域的比對率 /預定比對率*目標個體的染色體組數，所述一次區域的比對率為其包含的所有窗口的校 正比對率的均值。
[0013] 利用上述的本發明一方面的或者任一【具體實施方式】中的CNV檢測方法，能夠解決 上述現有的CNV檢測流程中，存在一些不足，例如現有方法中的采用固定長度的窗口，不能 很好地解決單細胞測序中全基因組擴增所帶來的偏性問題以及重復序列問題，不能很好的 用于二倍體單細胞的CNV的檢測等。上述的本發明一方面的或者任一【具體實施方式】中的 CNV檢測方法，非常適用于基于單細胞測序數據的CNV檢測，特別是基于單細胞低深度測序 的CNV檢測，對于不同測序平臺采用不同擴增方法進行單細胞測序或組織測序的數據都是 有效的，適用性廣泛。在不同測序平臺使用不同全基因組擴增方法進行單細胞測序時，本發 明的方法在檢測CNV的敏感性和特異性都很好，尤其是基于成環循環擴增技術（MALBAC)的 Proton平臺的測序數據。而且，利用本發明的方法的檢測結果具有高重復性，結果可信。與 現有的CNV檢測方法比較，本發明的方法采用長度變化的窗口，有利于保持所有窗口比對 上的短序列數的平均值的穩定性，還可以避免重復序列區域帶來的影響，使得CNV檢測更 準確。
[0014] 依據本發明的另一方面，本發明提供一種檢測CNV的裝置，所述裝置能夠用以執 行或完成上述本發明一方面或者任一【具體實施方式】中的CNV檢測方法，所述裝置包括：數 據輸入單元，用以接收數據；數據輸出單元，用以輸出數據；處理器，用以執行計算機可執 行程序，執行所述計算機可執行程序包括實現上述本發明一方面或者任一【具體實施方式】中 的CNV檢測方法；以及，存儲單元，用以存儲數據，其中包括所述計算機可執行程序。所說的 計算機可執行程序可以保存在存儲介質中，所稱存儲介質可以包括：只讀存儲器、隨機存儲 器、磁盤或光盤等。本發明還提供一種計算機可讀存儲介質，其用于存儲供計算機執行的程 序，所述程序的執行包括完成前述本發明一方面的或者其任一【具體實施方式】中的CNV檢測 方法。前述對本發明的CNV檢測方法的優點和技術特征的描述也適用于該CNV檢測裝置和 計算機可讀存儲介質，在此不再贅述。
【附圖說明】
[0015] 本發明的上述和/或附加的方面和優點從結合下面附圖對實施方式的描述中將 變得明顯和容易理解，其中：
[0016] 圖1是本發明的一個【具體實施方式】中的窗口合并后的每個窗口的ratio的密度分 布圖；
[0017] 圖2是本發明的一個【具體實施方式】中的基于MDA擴增的CG平臺單細胞測序數據 檢測CNV的結果不意圖；
[0018] 圖3是本發明的一個【具體實施方式】中的基于MDA擴增的Proton平臺單細胞測序 數據的CNV檢測的結果示意圖；
[0019] 圖4是本發明的一個【具體實施方式】中的基于MALBAC擴增的Proton平臺單細胞測 序數據的CNV檢測的結果示意圖。
【具體實施方式】
[0020] 以下對本發明方法的一般步驟或者相關信息的獲取方式進行介紹。
[0021] 1.首先從 UCSC 的網站（//hgdownload. cse. ucsc. edu/goldenPath/hgl9/ bigZips/)下載hgl9參考基因組的序列文件chromFa. tar. gz。在此，CNV檢測方法將只使 用那些恰好比對到參考基因組一個位置的短序列，我們將Y染色體上擬常染色體區的序列 用N代替。對于Y染色體，后續過程使用的都是這個經修改的版本。擬常染色體區是X染 色體和Y染色體間唯一可發生互換的位置，這也是它的名稱由來，由于會發生互換的是常 染色體，X染色體和Y染色體一般是沒有互換的現象，只有在擬常染色體區反常地出現了互 換，導致男性和女性都帶有兩個該區域基因的復本。這使擬常染色體區的基因表達類似常 染色體，而非性染色體的伴性遺傳模式，因而得名。
[0022] 2.確定每一個檢測窗口（window)的大小。
[0023] a)針對Proton測序平臺的下機數據，可采用單端模擬數據來劃分。以hgl9參考 基因組為基準，模擬單端測序短序列，從基因組染色體的第一個堿基開始，每50個堿基為 一條讀段（reads)，并為其生成ID和質量值生成fastaq格式。然后依次往后移一個堿基， 直到短序列的末端為染色體的最后一個堿基。使用bowtie把模擬數據對到參考基因組上， 結果只保留那些唯一比對的短序列（即去除可重復比對上的短序列），使用samtools把比 對結果轉換為BAM格式。
[0024] 比對參數可設置為：bowtie-S-t-n 2_e 7〇-m 1 -best - strata，后續單細胞測序 數據的比對參數也可以一樣，參數意義為：-n 2表示高保真區域內錯配數不能超過2個，-e 70表示錯配位點質量值不能超過70, 一best報告文件中，每個短序列的匹配結果將按匹配 質量從高到低排序，一srtata與一best -起使用報告質量最高的那部分，-m 1表示報告所 有比對的短序列。
[0025] b)對于CG平臺的下機數據，通過對正常人的細胞系的一團細胞DNA的測序數據來 劃分。CG的下機數據，進行流程信息分析，例如利用Teramap軟件比對到參考基因組，然后 將比對結果的格式轉換成BAM格式結果。
[0026] 不同平臺的數據最后都可使用軟件samtools去除重復的短序列。記錄參考基因 組上每一個被短序列起始堿基覆蓋的位置，并把這些位置劃分到10, 000至100, 000個窗 口，每一個窗口內包含的位置數目完全相同，但其區間長度是變化的。然后分別計算每個窗 口所包含參考基因組序列的GC含量。
[0027] 3.提取單個細胞的DNA，進行全基因組擴增，然后建庫上機測序，得到下機數據， 并進行相應分析處理得到bam格式比對的結果。
[0028] a)Proton平臺，其下機的數據（BAM格式），我們使用BEDTools轉換為FASTQ格 式數據，然后使用Trimmomatic軟件對長于50bp的短序列從3'端截取有效長度（50bp加 上與全基因組擴增方法的primer等長的短序列序列），如多重退火和成環循環擴增技術 (MALBAC)的引物為35bp，有效長度為85bp，同時過濾掉長度小于有效長度的短序列。使用 bowtie把截取后的短序列比對到參考基因組，并用samtools view轉換成bam文件排序后 去除重復短序列。
[0029] b)CG平臺，使用其平臺研發的分析流程把下機短序列數據與參考基因組hgl9進 行比對。然后把比對結果轉成BAM格式并進行排序，samtools的單端模式去除重復短序列。
[0030] 4.對每一個窗口中比對上的短序列進行統計計數并進行標準化處理，即計算每個 窗口的比對率（ratio)=比對上的短序列數目/所有窗口比對上的短序列數的平均值。
[0031] 5.使用L0WESS算法確定的ratio-GC含量關系對每個窗口中標準化處理后得到的 ratio進行GC校正，獲得校正ratio。
[0032] 6.每個樣本根據所有窗口校正后的ratio值，可使用CBS segment軟件對窗口進 行合并形成無重疊的區域（segment)并計算其ratio值，把此ratio值賦給區域（segment) 內的每個窗口。具體的包括，（a)
計算區域N與其它所有窗口 的比對率的差異，獲得所有Zxy，其中，區域N為第x+1個窗口到第y個窗口之間的區域，Zxy 呈標準正態分布，取Zb= maXl<x<y<w|Zxy|，（b)將Z b與臨界值比較，當Zb超過臨界值時， 相應的區域N為預計的窗口合并區域，即區域N為發生CNV的區域（c)將窗口合并區域去 除，更新X、y和《，進行上述兩步驟（a)和（b)，直至無 Zb超過臨界值，即循環劃分窗口，直 到窗口不能再進行合并；其中，X和y為窗口的編號，w為窗口總數，所述R x為第X個窗口的 校正比對率，所述臨界值是Zxy分布中的預定概率的概率密度，所述預定概率多95%，1 < X < y彡w，Sx= R片…+RX，Sy= R彳…+Ry，Sw= R片…+RW。所說的Zxy分布為Z xy服從標準正 態分布，預定概率與臨界值一一對應。上述過程可理解為，假設區域N為正常非變異區域， 當Z b落入拒絕域，即Z b超過99. 9%對應的臨界值，可知發生了小概率事件，否定原假設，即 區域N為變異區域。各合并區域的ratio為其所包括窗口的校正ratio的均值，然后把該 合并區域的ratio值賦值給其包括的所有窗口，即為窗口的比對率。
[0033] 接著，對所有窗口的ratio畫密度曲線分布圖，如圖1所示。對于近二倍體細胞或 者二倍體是所有倍型的眾數的細胞，密度分布圖中最大峰值對應的ratio值則為該細胞拷 貝數為2的ratio值。
[0034] 7.把每個區域的ratio除以拷貝數為2的ratio,再乘以2,則得到每個區域的拷 貝數。
[0035] 8.計算CNV檢測的敏感性和特異性。敏感性=LT/LC，特異性=LT/L，其中，L :指 單細胞測序找到的CNV(彡1Mb)的總長度，LC :表示組織測序中找到的CNV(彡1Mb)的總長 度，LT :表示單細胞測序和組織測序共同找到的CNV(彡1Mb)的總長度。
[0036] 以下結合具體個體樣本對依據本發明的檢測方法及檢測結果進行詳細的描述。下 面示例，僅用于解釋本發明，而不能理解為對本發明的限制。在本發明的描述中，"一次"、 "二次"等為指代或描述方便，不能理解為有順序關系或者相對重要性指示，除非另有說明， "多個"的含義是兩個或兩個以上。
[0037] 除另有交待，以下實施例中涉及的未特別交待的試劑、序列（接頭、標簽和引物）、 軟件及儀器，都是常規市售產品或者公開的，比如購自Illumina公司的hise q2000測序平 臺建庫相關試劑盒等。
[0038] 實施例一：基于MDA擴增的CG平臺低深度測序數據的CNV檢測方法測試
[0039] 隨著高通量測序技術的蓬勃發展，以Complete Genomics (CG)、IlluminaSolexa 和Roche454為代表的二代測序，以及三代測序技術（即單分子測序技術）所包括的 HelicosGenetic Analysis System、單分子實時測序技術（SMRT)和納米孔單分子測序技 術等各種測序技術已成為單細胞組學研究的重要工具。CG平臺作為一種專注于人類基因 組的二代測序技術，能夠完整并精確地測序人類全基因組，其測序通量大，在業內被高度認 可。CG平臺其主要包括測序平臺，高通量過程自動化技術和完整的數據管理解決方案三個 部分，其測序平臺包括DNA納米陣列（DNANanoball arrays，DNB? arrays)和組合探針錨定 連接測序法（combinatorial probe-anchor ligation，cPAL?)，這兩項技術的應用大大減 少試劑的消耗和縮短成像的時間。所以我們首先在CG平臺、利用CG平臺的下機數據對本 發明的CNV檢測方法進行試驗驗證。
[0040] 從病人的膠質母細胞瘤組織中分離出了 3個單細胞，組織樣本來自北京天壇醫院 提供，提取每個單細胞的DNA并利用MDA全基因組擴增技術進行擴增，再進行文庫構建，然 后在CG平臺進行單細胞低深度測序。最后按本發明的方法進行單細胞CNV的檢測分析。為 了驗證本發明方法的CNV檢測效率，我們提取了組織的DNA進行文庫構建，然后在CG平臺 進行全基因組測序，并使用CG的標準分析流程檢測得到組織的CNV結果。3個單細胞樣本 (P1-T2-SC#)和組織樣本（P1-T2)檢測到的CNV如圖2所示，平行X軸的粗黑線表示各區域 的拷貝數，其大于2則說明該區域內發生拷貝數增加，小于2則說明此區域內發生拷貝數減 少，等于2則表示拷貝數正常，每個窗口的ratio值用散點表示。
[0041] 進一步對CNV檢測方法的敏感性和特異性進行估計，敏感性=LT/LC，特異性= LT/L。估算5個樣本的平均敏感性為91. 01%，特異性為74. 47%，結果如表1所示。
[0042] 表 1
[0043]
[0044] 然后，對基于CG測序平臺MDA擴增的CNV檢測方法的重復性進行統計，發現樣本 間的重復性高于0. 7,結果詳見表2。
[0045] 表 2
[0046]
[0047] 從敏感性、特異性以及重復性統計計算結果可以得出，本發明的CNV分析檢測流 程的檢測結果的有效性，其在CG測序平臺上是可行的。
[0048] 實施例二：基于MDA擴增的Proton平臺的單細胞低深度測序的CNV檢測方法測試
[0049] 現有單細胞測序數據多由Illumina測序平臺產出。盡管Illumina測序儀的測序 通量大，但是其上機測序時間周期長，測序成本高，這些會限制單細胞CNV檢測分析的快速 發展。而一些研究往往對測序通量沒有需求，相對的，對測序的時間和成本具有更高的需 求，此時Proton測序平臺將是更好的選擇。Proton測序平臺運行速度快，測序周期只需幾 個小時，測序成本低，更適合部署到醫院或第三方檢測機構，縮短檢測時間，降低成本，從而 提尚檢測效率。而對于Proton的單細胞測序CNV檢測鮮有報道。
[0050] 提取來自北京大學腫瘤醫院的人類胃腺癌細胞系（BGC823)5個細胞（MDA-2_ BGC#)，并用多重置換擴增（MDA)技術進行文庫構建后在Proton平臺進行單細胞低深度測 序。同時提取人類胃腺癌細胞系（BGC823) -團細胞（BGC)的DNA，進行常規文庫構建后在 Proton平臺進行測序。然后按我們的方案進行CNV的檢測分析，在BGC823的5個單細胞樣 本和組織樣本（BGC)檢測到的CNV如圖3,每個區域的拷貝數（平行X軸的粗黑線）大于2 則為說明區域發生拷貝數增加，小于2為拷貝數減少，等于2則為拷貝數正常，三種拷貝數 變化區域內窗口的ratio值分別用不同灰度深度的散點表示。
[0051] 在整個基因組上多個染色體上檢測到大片段的CNV，與細胞團的CNV檢測結果保 持一致，驗證了本發明方法檢測CNV的有效性。
[0052] 然后，根據五個單細胞和一團細胞CNV檢測結果，我們進一步對檢測CNV方法的敏 感性和特異性進行估計，敏感性=LT/LC，特異性=LT/L。估算5個單細胞樣本的平均敏感 性為85. 86%，特異性為81. 18%，結果如表3所示。
[0053] 表 3
[0054]
[0055] 對基于Proton測序平臺MDA擴增的CNV檢測方法的重復性進行統計，發現樣本間 的重復性高于0. 7,結果詳見表4。
[0056] 表 4
[0057]
[0058] 從敏感性、特異性以及重復性統計計算結果可以得出，本發明的CNV檢測流程分 析結果的有效性，其對于Proton測序平臺MDA擴增方法的測序數據是可行的。
[0059] 實施例三：基于MALBAC擴增的Proton平臺的單細胞低深度測序的CNV檢測方法 測試
[0060] 提取5個人類胃腺癌細胞系細胞（BGC823)，用MALBAC全基因組擴增方法進行常規 文庫構建后在Proton平臺進行單細胞低深度測序；同時提取人類胃腺癌細胞系（BGC823) 一團細胞（BGC)的DNA在Proton平臺進行測序。得到的下機數據按我們的方案進行CNV 的檢測分析，在BGC823五個樣本找到CNV，結果如圖4所示，其中，橫坐標表示染色體；右側 縱坐標為5個單細胞樣本及團細胞樣本，左側縱坐標為拷貝數，圖上的黑色粗線代表劃分 區域的ratio值，其值大于2的說明此區域拷貝數增加，小于2則說明區域內拷貝數減少， 等于2則說明區域內拷貝數正常。三種拷貝數變化區域分別用不同灰色深度的散點表示窗 口的 ratio 值。
[0061 ] 進一步對本發明的檢測CNV方法的敏感性和特異性進行估計，敏感性=LT/LC，特 異性=LT/L。估算5個樣本的平均敏感性為84. 72%，特異性為85. 18%，結果如表5。
[0062] 表 5
[0063]
[0064] 對基于Proton測序平臺MALBAC擴增的CNV檢測方法的重復性進行統計，發現樣 本間的重復性高于〇. 92,詳見表6。
[0065] 表 6
[0066]
[0067] 從敏感性、特異性以及重復性統計計算結果可以得出，本發明的CNV檢測流程分 析結果的有效性，其對于Proton測序平臺MALBAC擴增方法的測序數據是可行的。
【主權項】
1. 一種檢測CNV的方法，其特征在于，包括W下步驟： 獲取目標個體的基因組測序結果，所述測序結果包括多個讀段； 將所述測序結果與參考序列比對，獲得比對結果，所述參考序列包括多個窗口，所述比 對結果包括比對上每個窗口的讀段的數目； 基于所述比對結果，計算每個窗口的初始比對率，窗口的初始比對率=比對上所述窗 口的讀段數目/比對上所有窗口的讀段數目的平均值，所述比對上所有窗口的讀段數目的 平均值=比對上所有窗口的讀段總數/窗口個數； 合并初始比對率無顯著差異的多個相鄰窗口，定義合并后的多個相鄰窗口為一次區 域，剩余的每個獨立窗口也分別稱為一次區域； 基于所述一次區域的比對率與預定比對率不相等，判定所述一次區域存在CNV， 所述一次區域的比對率為所述一次區域包含的窗口的初始比對率的均值， 所述預定比對率為所有窗口的比對率中頻率最高的窗口的比對率，所述窗口的比對率 為其所在的一次區域的比對率。2. 權利要求1的方法，其特征在于，所述基因組獲自所述目標個體的單個細胞； 任選的，通過構建所述細胞的基因組測序文庫，并對所述測序文庫進行序列測定獲得 所述測序結果； 任選的，構建所述測序文庫包括對所述基因組進行簡并寡核巧酸引物PCR，多重置換擴 增和/或多次退火環狀循環擴增。3. 權利要求1的方法，其特征在于，所述參考序列為人參考基因組； 任選的，W N替代所述人參考基因組的Y染色體的擬常染色體區的每個堿基，N表示A、 T、C和G中的任一種。4. 權利要求1的方法，其特征在于，所述窗口的確定，包括， 將短序列集與參考序列比對，確定比對上所述參考序列的短序列的起始位置，所述短 序列集包括多個短序列，所述短序列集來自模擬序列集和/或測序結果； 在所述參考序列上劃定窗口，使每個所述窗口包含相同數目的所述起始位置； 任選的，所述模擬序列集中的模擬序列的獲取包括， 從所述參考序列的長度為Q的染色體的一端的堿基開始，拷貝所述染色體的P個堿基， W獲得第一條模擬序列， 沿所述染色體的另一端方向移動一個堿基拷貝所述染色體的P個堿基，W獲得第二條 模擬序列， 沿所述染色體的另一端方向移動兩個堿基拷貝所述染色體的P個堿基，W獲得第=條 模擬序列， 依此獲得第Q-P+1條模擬序列，所述第Q-P+1條模擬序列的末端堿基與所述染色體的 另一端的堿基重合，其中， P為模擬序列的長度，P > 10 ; 任選的，所述窗口之間沒有重疊； 任選的，所述窗口總數不大于100,000。5. 權利要求1的方法，其特征在于，在所述合并初始比對率無顯著差異的多個相鄰窗 口之前，利用比對率-GC含量的關系對每個所述窗口的初始比對率進行GC校正，獲得各個 窗口的校正比對率， W所述窗口的校正比對率替代該窗口的初始比對率。6. 權利要求5的方法，其特征在于，建立所述比對率-GC含量的關系，包括， 獲得至少一個樣本的核酸的測序數據，所述測序數據由多個讀段組成； 將所述測序數據與參考序列進行比對，獲得比對結果，所述參考序列包括多個窗口，所 述比對結果包括比對上每個所述窗口的讀段的數目； 計算每個所述窗口的初始比對率，窗口的初始比對率=比對上所述窗口的讀段數目/ 比對上所有窗口的讀段數目的平均值，所述比對上所有窗口的讀段數目的平均值=比對上 所有窗口的讀段總數/窗口個數； 基于多組的窗口的初始比對率和該窗口的GC含量，利用二維回歸分析法建立所述比 對率-GC含量的關系； 任選的，所述二維回歸分析法為局部加權回歸散點平滑法。7. 權利要求5的方法，其特征在于，所述合并初始比對率無顯著差異的多個相鄰窗口 是指，合并滿足W下的相鄰窗口， 校正比對率的差異無統計意義。8. 權利要求7的方法，其特征在于，所述合并初始比對率無顯著差異的多個相鄰窗口， 包括， (a) 基于公^5 計算區域N與其它所有窗口的比對率的差異，獲得 所有 Zxy,取 Zb= max Jz巧 I， (b) 將Zb與臨界值比較，當Z b超過所述臨界值時，相應的區域N為所述一次區域， (C)將所述一次區域去除，更新x、y和W，進行步驟（a)和化），直至無 Zb超過所述臨界 值，其中， X和y為窗口的編號， W為窗口總數， 所述區域N為第X+1個窗口到第y個窗口之間的區域， Rx為第X個窗口的校正比對率， 所述臨界值是Zyy分布中的預定概率的概率密度，所述預定概率> 95%， 1《X < y《W， Sx= R 1+. .. +Rx, Sy= R 1+. .. +Ry, S"= R 1+. . . +斬。9. 權利要求8的方法，其特征在于，還包括， 基于比較所述一次區域的比對率與所述預定比對率的大小，判定所述CNV的類型，其 中包括， 當所述一次區域的比對率大于所述預定比對率時，判定所述一次區域為拷貝數增加區 域， 當所述一次區域的比對率小于所述預定比對率時，判定所述一次區域為拷貝數減少區 域。10. 權利要求7-9任一方法，其特征在于，還包括， 利用W下公式計算所述一次區域的拷貝數， 一次區域的拷貝數=該一次區域的比對率/預定比對率*目標個體的染色體組數， 所述一次區域的比對率為其包含的所有窗口的校正比對率的均值。11. 一種檢測CNV的裝置，其特征在于，包括， 數據輸入單元，用W接收數據； 數據輸出單元，用W輸出數據； 處理器，用W執行可執行程序，執行所述可執行程序包括完成權利要求1-10任一方 法；化S， 存儲單元，用W存儲數據，其中包括所述可執行程序。
【文檔編號】C12Q1/68GK105986008SQ201510039685
【公開日】2016年10月5日
【申請日】2015年1月27日
【發明人】李甫強, 史旭蓮, 謝國云, 魯娜, 趙至坤, 蔣潤澤, 梁瀚, 侯勇, 吳逵
【申請人】深圳華大基因科技有限公司