李 燕，李垚垚（哈爾濱醫(yī)科大學大慶校區(qū)，黑龍江大慶163319）

doi：10.3969/j.issn.1672－5565.2015.03.07

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

李 燕?，李垚垚
（哈爾濱醫(yī)科大學大慶校區(qū)，黑龍江大慶163319）

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

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

LIYan?，LIYaoyao

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

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

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

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

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

1 拷貝數(shù)變異概述

1.1 拷貝數(shù)變異含義

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

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

1.2 目前檢測方法

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

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

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

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

2 算法介紹

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

2.1 仿真實驗

2.1.1 檢驗標準

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

2.1.2 仿真數(shù)據(jù)

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

Simulation的過程：

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

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

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

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

2.2 置換檢驗

2.2.1 置換檢驗概述

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

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

2.2.2 算法設計與實現(xiàn)

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

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

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

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

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

觀察每個位點對應的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區(qū)域的p?value值大小明顯＜0.01，說明在這些區(qū)域小概率事件發(fā)生，原假設錯誤，而是在這些區(qū)域發(fā)生了拷貝數(shù)變異。這與仿真數(shù)據(jù)時的變異區(qū)域相同，因此本算法可以檢測拷貝數(shù)變異。

3 真實數(shù)據(jù)CNV檢測

3.1 數(shù)據(jù)來源與處理

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

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

3.2 CNV值計算

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

2/R0＝x/R1

3.3 數(shù)據(jù)結果分析

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

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

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

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

3.5 算法的性能與評價

3.5.1 仿真代碼實現(xiàn)上

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

3.5.2 算法應用上

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

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

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

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

4 結束語

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

［1］ FREEMAN J L，PERRY G H，F(xiàn)EUK 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，F(xiàn)EUK 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，F(xiàn)ITCHK 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，F(xiàn)IUME 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）。

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