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

doi：10.3969/j.issn.1672－5565.2015.03.07

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

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

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

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

LIYan?，LIYaoyao

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

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

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

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

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

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

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

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

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

1.2 目前檢測(cè)方法

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

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

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

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

2 算法介紹

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

2.1 仿真實(shí)驗(yàn)

2.1.1 檢驗(yàn)標(biāo)準(zhǔn)

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

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

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

Simulation的過程：

（1）隨機(jī)構(gòu)建一個(gè)染色體位點(diǎn)數(shù)為2 000的樣本，并對(duì)每個(gè)位點(diǎn)編號(hào)。

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

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

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

2.2 置換檢驗(yàn)

2.2.1 置換檢驗(yàn)概述

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

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

2.2.2 算法設(shè)計(jì)與實(shí)現(xiàn)

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

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

采用置換實(shí)驗(yàn)，計(jì)算每個(gè)位點(diǎn)對(duì)應(yīng)的p?value（P ［j］，j＝1，2，…，2 000）：

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

其中，f［j］為每個(gè)位點(diǎn)對(duì)應(yīng)的頻數(shù)。

觀察每個(gè)位點(diǎn)對(duì)應(yīng)的p?value，并繪制見圖2。

圖2 各位點(diǎn)對(duì)應(yīng)P值Fig 2 The p?value of every site

2.2.3 實(shí)驗(yàn)結(jié)論

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

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

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

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

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

3.2 CNV值計(jì)算

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

2/R0＝x/R1

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

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

圖3 Chrom20 seq1：1－1 569和seq2：37－1 567上位點(diǎn)對(duì)應(yīng)的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上位點(diǎn)對(duì)應(yīng)的reads數(shù)分布圖Fig.4 Reads distribution map of the sitechrom20 seq1：1－1 569

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

3.5 算法的性能與評(píng)價(jià)

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

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

3.5.2 算法應(yīng)用上

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

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

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

圖6 算法時(shí)間復(fù)雜度對(duì)比Fig.6 Algorithm time complexity contrast

4 結(jié)束語

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

黑龍江省教育廳科學(xué)技術(shù)研究項(xiàng)目（12541565）。

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