新一代測序的拷貝數變異檢測算法研究與設計

李 燕，李垚垚（哈爾濱醫科大學大慶校區，黑龍江大慶163319）

doi：10.3969/j.issn.1672－5565.2015.03.07

新一代測序的拷貝數變異檢測算法研究與設計

李 燕?，李垚垚
（哈爾濱醫科大學大慶校區，黑龍江大慶163319）

基于不同的測序技術，基因拷貝數變異的檢測方法有多種，但時間復雜度較高，而新一代測序技術的發展為基因拷貝數變異檢測的研究開辟了新領域。通過仿真實驗、置換檢驗設計出一種新的基于新一代測序的拷貝數變異檢測算法。不同于其它算法，本算法無需參考樣本，通過直接研究比對后的序列以及reads與拷貝數的關系，來研究檢測拷貝數變異，實驗結果表明在時間復雜度上能提高50％以上的運算速度，這對今后拷貝數與疾病的研究具有重要意義。

新一代測序；拷貝數變異；仿真；置換檢驗

LIYan?，LIYaoyao

（Harbin Medical University Daqing campus，Daqing Heilongjiang 163319，China）

新一代測序（New generation sequencing，NGS）技術的發展越來越成熟，各測序平臺層出不窮，基因序列的測序成本大幅度地下降，測序的速度越來越高，這使得測序產生的DNA序列數據非常龐大，怎樣理解數據成為當務之急。

伴隨著人類基因組計劃及 1 000 genomes project的實施與發展，蛋白質、DNA、RNA的序列數據的規模日趨增加，僅僅依靠生物實驗來研究生物基因變異及疾病產生早已不能滿足現實需要，因此必須借助計算機、數學等學科的理論及思想方法從海量數據中來研究和闡明生物學問題。拷貝數變異（Copy number alterations，CNAs）檢測是生物信息學中研究生物基因結構改變的有效方法之一。

迄今為止，在HapMap計劃的樣本研究基礎上，已經基本構建成人類第一代基因組CNV圖譜［1］。隨著測序技術的發展，新一代測序技術更成熟，從NGS數據出發，更多的拷貝數變異可能被檢測，這也為研究CNV檢測算法開辟了新領域。

新一代測序技術在對數據的處理過程中，會產生許多數據格式：FASTQ文件、SAM文件、VCF （Variant call format）文件、TXT文件和 BED文件等［2］。本文算法的重點研究對象是txt文件，當利用samtools工具中mpileup命令處理數據時，無“－g”或“－u”參數時會輸出類似“.txt”文本文件，此文本文件統計了參考序列上每一堿基位點的比對結果，每一行表示reference中某一堿基位點的比對情況［3，14］。

1 拷貝數變異概述

1.1 拷貝數變異含義

誘發基因變異的因素有多個方面，基因的遺傳變異的方式也多種多樣［4］。大部分研究都表明，CNV指大小從Kb到Mb范圍內的亞微觀（Submicroscopic，指的是在普通電子顯微鏡下能分辨的范圍）片段發生了拷貝數突變，這些拷貝數的復制、缺失、倒置等變異，統稱為拷貝數變異（Copy number alterations，CNAs），但不包括轉座子的插入和缺失引起的基因變異［5－7］（見圖1）。

圖1 基因組中的拷貝數變異Fig.1 Copy number variation in genome

1.2 目前檢測方法

目前拷貝數變異的檢測方法主要分為三大類：一是定量PCR技術；二是 基于芯片的 array?based comparative genomic hybridization和SNPs芯片；三是新一代測序技術。

對于目標基因CNV檢測常常采取基于定量PCR技術和雜交技術的方法。其中熒光定量PCR技術應用比較廣泛，它的一個反應只測得一個拷貝，通過將檢測樣本的目標基因與參照基因定量后的檢測值的比值相比較來估計此樣本基因的拷貝數［6］。

基于芯片技術的CNV檢測方法主要有：比較基因組雜交（Comparative genomic hybridization，CGH）技術、aCGH技術、oaCGH技術和SNPS芯片技術。其中，aCGH是基于微陣列的CGH技術，其芯片探針可以覆蓋整個基因組，因此這種高通量分析法的準確度、敏感度和分辨度更高，結果更加準確［8］。SNPs芯片技術不同于CGH技術，僅僅使用單雜交就可實現檢測。它是通過被測試的樣本信號強度跟其他樣本個體的強度作比較來確定每一位點對應的基因拷貝數［9］。這些方法都比較適合在全基因組范圍內尋找CNV。

目前基于新一代測序數據的CNV方法主要有：分解讀段（Split read）、讀段深度（Read depth，RD）、末端配對法（Pair?end mapping，PEM）和重組（Assembly）等。由于新一代測序技術具有高通量、門檻低、簡單等特點，因此基于NGS的CNV檢測方法克服了雜交固有的某些缺點，即不需要太多特別復雜的設計工作，可以直接處理比對后數據，無需參考樣本并可應用自身測序鑒定基因變化，而且費用相對低于aCGH技術。所以基于新一代測序的拷貝數檢測方法具有良好的發展前景，這也為本次論文的研究內容提供了方向。

2 算法介紹

本文算法的目的是檢測基于新一代測序的拷貝數變異，無需參考樣本，這既減少了實驗樣本數量，還降低實驗成本與時間。實驗分為兩大部分：（1）算法的設計及仿真實驗；（2）真實數據的應用。

2.1 仿真實驗

2.1.1 檢驗標準

在新一代測序時，高通量測序儀器一個反應得到的測序序列片段稱為reads［10］。不同的測序儀器產生的reads數長度也不同，reads數的長度大小在36～200 bp不等。正常在沒有發生拷貝數變異時，當測序depth和coverage一定時，同一測序儀器測序得到的一條染色體上的堿基序列上的reads數是基本相同的，若該序列上的reads數有一段區域不同于其他大部分區域，則可能說明這段reads數異常區域可能發生了拷貝數變化［11，15］。Reads數的異常主要表現在拷貝數的缺失、擴增等。因此本實驗選取reads數作為衡量是否發生拷貝數的標準［15］。為了產生模擬數據這里自行定義reads數S＝40bp，為測得正常序列的reads。若測序區域＜40或＞40，我們都認為其發生了拷貝數變異。

2.1.2 仿真數據

由于受到目前測序儀器和水平的限制，測序所得堿基序列的reads數會不一致，reads數可能會上下波動，但仍然處于相同水平。這里選reads∈［39，40，41］來模擬實驗數據。

Simulation的過程：

（1）隨機構建一個染色體位點數為2 000的樣本，并對每個位點編號。

（2）任取多個區域如100－149，500－529，900－919，1 600－1 650，對其進行信號加強/減弱處理，模擬成這幾段標記區域發生reads數變化（即拷貝數改變）［12－13］。

為了使實驗數據更逼近實際測得序列，減少誤差，需要對仿真數據進行加噪聲處理。這里主要是利用高斯噪聲處理，并對隨機其他位點噪聲處理。

經過上述步驟，產生了一個包含2 000個位點的樣本。而在統計實驗中一個樣本不能證明任何實驗問題，需要大量的樣本才能減少誤差，得出結論。因此我們重復上述步驟，產生了50個樣本用于實驗。

2.2 置換檢驗

2.2.1 置換檢驗概述

通常顯著性檢驗可以確定一個觀測值是否有效［16］，如假設檢驗中檢測兩組樣本的均值是否有相等（或者檢測哪一均值更大）。本次實驗仿真出一些小樣本結果（這里是50個小樣本），借助于Permutation test置換檢驗來分析小樣本的總體分布。

Permutation test是20世紀30年Fisher提出的基于大量計算，根據對樣本中的數據隨機（或全）排列，統計并推斷的一種方法。算法公布之初，由于它的運算量沒能得到重視與應用。近年來隨著計算機的性能提高，我們可以借助計算機的計算能力來實現置換檢驗來解決問題。它是基于樣本本身的，對樣本的總體分布要求自由，因此應用相對較廣泛，尤其適合用于對總體分布未知的小樣本數據分析，以及一些用常規方法難以分析的假設檢驗問題。置換檢驗的過程一般是：首先對樣本內的數據進行順序置換，然后重新計算檢驗統計量，并構造出經驗分布，最后求出P?value來推斷結果。

2.2.2 算法設計與實現

假設設計一個實驗來驗證仿真實驗中樣本位點數100－149，500－529，900－919，1 600－1 650的區域發生了拷貝數改變。（實驗數據被保存在merge1.txt中。）

用假設檢驗的方法完成驗證，選取樣本位點對應的reads數構造為檢驗統計量。零假設為：樣本堿基對應的reads數沒有發生改變（即堿基序列對應的copy number未發生改變，是正常的）。在這個檢驗中，最終計算出2 000個位點對應的p?value值若＜0.05的區域，則表明小概率事件發生，而原假設是正常的，因此原假設錯誤，此區域（位點）發生了拷貝數變異。

采用置換實驗，計算每個位點對應的p?value（P ［j］，j＝1，2，…，2 000）：

P［j］＝f［j］/1 000

其中，f［j］為每個位點對應的頻數。

觀察每個位點對應的p?value，并繪制見圖2。

圖2 各位點對應P值Fig 2 The p?value of every site

2.2.3 實驗結論

由圖1放大可以直觀看出在100－149，500－529，900－919，1 600－1 650區域的p?value值大小明顯＜0.01，說明在這些區域小概率事件發生，原假設錯誤，而是在這些區域發生了拷貝數變異。這與仿真數據時的變異區域相同，因此本算法可以檢測拷貝數變異。

3 真實數據CNV檢測

3.1 數據來源與處理

為檢測上述方法的適用性，本文從 1 000 genomes project數據庫中獲得真實數據，為了保證數據可用性，下載真實數據要確保控制單一變量reads數變化，其他如read depth、read coverage、測序儀器等要控制一致［17］。這里采用 HG00096. mapped.ILLUMINA.bwa.GBR.low＿coverage.20120522. bam.中chrom20的數據作為數據應用上述檢測方法。同時為了證明在high coverage數據同樣適用，還處理了 HG00096.mapped.ILLUMINA.bwa.high＿coverage.bam數據。

利用Samtools軟件對真實數據進行處理，CBS方法去除噪聲，提取reads數，統計頻數最多的reads值。考慮真實數據噪聲和測量誤差，可確定實驗數據區域正常情況下reads數在［39，43］，并以此為基準檢測該區域內是否發生了拷貝數變異。若區域內位點對應的reads值小于或大于這個區間，認為對應位點發生了拷貝數變異。

3.2 CNV值計算

正常情況下，人類基因拷貝數變異的值為2。研究表明，某一位點拷貝數變異的數目與對應的reads值成如下的關系［19］：其中R0為測序深度、覆蓋度一致時正常情況下區域或位點對應的reads值，R1為待測區域或位點對應的reads值，x即為待測區域或位點的拷貝數的值。因此可以計算任意位點的拷貝數。

2/R0＝x/R1

3.3 數據結果分析

本實驗HG00096.mapped.ILLUMINA.bwa.GBR. low＿coverage.20120522.bam.中chrom20上的seq1：1－1 569位點和seq2：37－1 567位點上的reads數據應用上述算法，并繪制如下圖3。圖3為chrom20的seq1：1－1 569和seq2：37－1 567位點對應的reads數分布情況，圖4和圖5分別為seq1和seq2相應位點reads數分布圖，其中紅線部分表示被測區域內reads數出現最多的數值。大多數位點都在紅線附近上下波動，當位點對應的reads數距離紅線越遠時，我們認為該位點可能發生了拷貝數變異。如圖5中seq1：1－220點附近，圖5中seq2：190－250位點附近等，我們可以很直觀地推測這些區域可能發生了拷貝數變異。還可根據數據確定變異邊界，利用公式計算各位點對應的拷貝數值。

圖3 Chrom20 seq1：1－1 569和seq2：37－1 567上位點對應的reads數分布圖Fig.3 Reads distribution map of the sitechrom20 seq1：1－1 569 and seq2：37－1 567

圖4 Chrom20 seq1：1－1 569上位點對應的reads數分布圖Fig.4 Reads distribution map of the sitechrom20 seq1：1－1 569

圖5 Chrom20 seq2：37－1 567上位點對應的reads數分布圖Fig.5 Reads distribution map of the site chrom20 seq2：37－1 567

3.5 算法的性能與評價

3.5.1 仿真代碼實現上

本算法程序代碼基于R語言相對容易實現，對于涉及的數據預先分配空間，大大降低了時空復雜度。但是在permutation、merge data以及做test時會涉及到雙層for循環，再加之數據樣本自身很大，因此增加了時間復雜度。為了減少時間消耗，提升速度，在編寫代碼時除采用了向量化避免for循環，加入并行運算方法。

3.5.2 算法應用上

算法基于新一代測序技術測序數據，與基于芯片的檢測技術相比，本算法無需參考樣本，數據來源更真實，使得檢測的拷貝數也更真實，大大減少了誤差，同時也最大的降低了檢測費用。

本算法在雙核x86 32 bit的處理器中執行，經檢驗，在內存占有量相差無幾時，時間復雜度降低明顯（本算法樣本計算時間120.2 s，CNV－seq計算時間251.5 s，，FREEC計算時間319.6 s），如圖6所示。同時在檢測邊界也具有相當高的靈敏度，直接從比對后的數據處理，也降低了從raw data到mapped data中產生的各種誤差。

本算法能夠檢測出拷貝數變異，但是對拷貝數變異的類型不能很清晰的界定，這一方面有待改善。它對測序數據的格式等要求比較嚴格，要保證實驗數據序列的read coverage，read depth等一致，還要保證數據是基于同一測序技術測得的。同時，它只對新一代測序的數據有效，隨著第三代測序技術的萌芽，在檢測拷貝數變異時可能會出現瑕疵，但可以借鑒思想，在未來很長時間仍然受用無窮。

圖6 算法時間復雜度對比Fig.6 Algorithm time complexity contrast

4 結束語

CNV作為基因結構變異的一種重要形式，對人類遺傳進化、疾病和藥物研究等具有重要的意義［20］。在眾多檢測方法中，目前急需效率高和準確性高算法。本文提出了一種基于新一代測序數據的CNAs檢測算法，無需額外的參考樣本序列，利用置換檢驗的方法檢驗，降低假陽性率，增強結論的真實性，提高了準確度。實驗表明，這種基于新一代測序的拷貝數變異檢測算法，可快捷方便地找出由新一代測序技術測得的染色體上可能發生拷貝數變異的位點，大大降低了時間復雜度。這對今后拷貝數與疾病的研究具有重要意義。

［1］ FREEMAN J L，PERRY G H，FEUK L，et al.Copy number variation：new insights in genome dicersity［J］. Genome Res，2006，16：949－961.

［2］ SHENDURE J，JIH.Next?generation DNA sequencing ［J］.Nat Biotechnical，2008，26：1135－45.

［3］ SCHUSTER SC.Next?generation DNA sequencing transforms today's biology［J］.Nat Methods，2008，5：16－8.

［4］ IAFRATE A J，FEUK L，RIVERA M N，et al.Detec?tion of large?scale variation in the human genome［J］.Nat Genet，2004，36（9）：949－951.

［5］ XIE C，TAMMIM T.CNV?seq，a new method to detect copy number variation using high?throughput sequencing ［J］.BMC Bioinformatics，2009，10：80.

［6］ BOEVA V，ZINOVYEV A，BLEAKLEY K，et al.Con?trol?freecalling of copy number alterations in deep?sequen?cing data using GC?contentnormalization［J］.Bioinformat?ics，2011，27（2）：268－269.

［7］ REDON R，ISHIKAWA S，FITCHK R，et al.Global variation in copy number in the human genome［J］. Nature，2006，444：444－454.

［8］ COOPER G M，NICKERSON D A，EICHLER E E. Mutational and selective effects on copy?number variants in the human genome［J］.Nature Genetics，2007，39：S22－29.

［9］ CHIANG D Y，GETZ G，JAFFE D B，et al.High?resolution mapping of copy?number alterations with massively parallel sequencing［J］.Nat Methods，2008，6 （1）：99－103.

［10］MILLER C A，HAMPTON O，COARFA C，et al. ReadDepth：a parallel R package for detecting copy number alterations from short sequencing reads［J］.PLoS ONE，2011，6：16327.

［11］YOON S，XUAN Z，MAKAROV V，et al.Sensitive and accurate detection of copy number variants using read depth of coverage［J］.Genome Res，2009，19：1586－1592.

［12］VANCE A.Data analysts captivates by R's power［J］. New York Times，2009，6：22－29.

［13］VENABLESW N，SMITH D M，TEAM R D C.An introduction to R［M］.Network Theory，2006：34－38..

［14］LIH，HANDSAKER B，WYSOKERA，et al：The sequence alignment/map format and SAMtools［J］. Bioinformatics，2009，25（16）：2078－2079.

［15］MEDVEDEV P，FIUME M，DZAMBA M，et al. Detecting copy number variation with mated short reads ［J］.Genome Res，2010，20（11）：1613－1622.

［16］MAGIS，TATTINI L，PIPPUCCI T，et al.Read count approach for dna copy number variants detection［J］. Bioinformatics，2012，28（4）：470－478.

［17］WANG J，WANG W，LI R，et al.The diploid genome sequence of an Asian indicidual［J］.Nature，2008，456：60－65.

［18］TIERNEY L，ROSSINI A J，LI N.Snow：A parallel computing framework for the R system［J］.Int JParallel Program，2009，37（1）：78－90.

［19］ABYZOV A，URBAN A E，SNYDER M，et al.An approach to discover，genotype，and characterize typical and atypical cnvs from family and population genome sequencing［J］.Genome Res，2011，21（6）：974－984.

［20］KORBEL JO，URBAN A E，AFFOURTIT J P，et al. Paired?end mapping reveals extensive structural variants detection in the human genome［J］.Science，2007，318 （5849）：420－426.

An algorithm for detecting copy number alteration from next generation sequencing of human genome

Based on different sequencing technologies，the detection methods of gene copy number variation are available.However，with the development of new generation sequencing technology，a new field for researchingcopy number variations has been opened up.Through the simulation experimentand the replacement test，this paper designs a new copy number variation detection algorithm based on the new generation of sequencing.Unlike other algorithms，our algorithm doesn't need thereference samples，butuses themapped datafrom next generation sequencing platforms and the relationship between reads and gene copy number to detect gene copy number variations in the genome.The experimental results show that the performance in time complexity can be improved bymore than 50％，indicating the important significance for the further study of gene copy number and disease in the future.

Next generation sequencing；Copy number variations；Simulation；Permutation test

TP301.6

A

1672－5565（2015）03－186－06

2015－06－19；

2015－07－14.

黑龍江省教育廳科學技術研究項目（12541565）。

李燕，女，教授，研究方向：數據庫與數據挖掘；E?mail：qliyan＠163.com.