徐安琪 蔡 磊 高敬陽(yáng)

(北京化工大學(xué) 信息科學(xué)與技術(shù)學(xué)院， 北京 100029)

引 言

拷貝數(shù)變異(copy number variation, CNV)是人類基因組中一種重要的結(jié)構(gòu)變異(structural variation，SV)，它的位點(diǎn)突變率遠(yuǎn)高于單核苷酸變異(single nucleotide polymorphism，SNP)[1-2], 是人類疾病的重要致病因素之一。多項(xiàng)研究表明，拷貝數(shù)變異與多種人類腫瘤的發(fā)生密切相關(guān)，例如乳腺癌、肝癌、腎癌、食管癌等[3-5]。

現(xiàn)有的二代/三代測(cè)序技術(shù)使用的測(cè)序數(shù)據(jù)是多種細(xì)胞的混合樣本，對(duì)這樣的數(shù)據(jù)進(jìn)行分析得到的結(jié)果只是一群細(xì)胞中信號(hào)的平均值，單個(gè)細(xì)胞的特異性被忽略，導(dǎo)致細(xì)胞的異質(zhì)性不能得到體現(xiàn)[6]，而異質(zhì)性恰恰是癌癥發(fā)生、發(fā)展和復(fù)發(fā)的基礎(chǔ)[7]。作為在單個(gè)細(xì)胞水平上對(duì)基因組、轉(zhuǎn)錄組、表觀組進(jìn)行高通量測(cè)序以揭示單個(gè)細(xì)胞的基因結(jié)構(gòu)和基因表達(dá)狀態(tài)的新技術(shù),單細(xì)胞測(cè)序(single cell sequencing，SCS)是研究細(xì)胞異質(zhì)性的有力工具[8]，對(duì)于探索癌癥細(xì)胞異質(zhì)性有著重要的意義。但由于單細(xì)胞測(cè)序是在單個(gè)細(xì)胞水平上進(jìn)行的測(cè)序，所以需要對(duì)單個(gè)細(xì)胞微量全基因組DNA進(jìn)行擴(kuò)增。一般使用全基因組擴(kuò)增(whole genome amplification，WGA)技術(shù)進(jìn)行擴(kuò)增，這種擴(kuò)增技術(shù)可以大幅度增加DNA總量，但同時(shí)也會(huì)使得到的數(shù)據(jù)具有較低的基因覆蓋度和較高的擴(kuò)增偏好性，從而導(dǎo)致定量分析結(jié)果不準(zhǔn)確[9]。

對(duì)于這種具有低覆蓋度和高擴(kuò)增偏好性的單細(xì)胞數(shù)據(jù)，大多數(shù)用于高通量測(cè)序(next generation sequencing，NGS)的CNV檢測(cè)方法并不能取得良好的檢測(cè)效果，只有少數(shù)可適用于單細(xì)胞數(shù)據(jù)CNV檢測(cè)，這些方法從原理上基本可以分為3類，即滑動(dòng)窗口、構(gòu)造函數(shù)和隱馬爾可夫模型[10]。這3種方法各有其優(yōu)勢(shì)和局限性[11-12]：滑動(dòng)窗口能合并相同狀態(tài)區(qū)間以提高斷點(diǎn)檢測(cè)的準(zhǔn)確性，但忽略了單細(xì)胞數(shù)據(jù)覆蓋度低和分布不均衡的數(shù)據(jù)特點(diǎn)，導(dǎo)致產(chǎn)生大量假陽(yáng)性變異；構(gòu)造函數(shù)雖然能很好地適應(yīng)單細(xì)胞的數(shù)據(jù)特征，但對(duì)變異區(qū)間和斷點(diǎn)的檢測(cè)不考慮區(qū)域性，無(wú)法很好地合并公共斷點(diǎn)； 隱馬爾可夫模型雖然擅長(zhǎng)處理異常數(shù)據(jù)點(diǎn)，但易出現(xiàn)過(guò)擬合或擬合程度不高的情況。

目前已有的這些方法各有利弊,且由于單細(xì)胞數(shù)據(jù)的特異性從而普遍存在檢測(cè)精度不高、變異類型識(shí)別不準(zhǔn)確、假陽(yáng)性變異數(shù)量偏高等問(wèn)題[13]。為了提高檢測(cè)的精度，減少單細(xì)胞數(shù)據(jù)CNV檢測(cè)中的假陽(yáng)性變異結(jié)果，本文提出了一種針對(duì)單細(xì)胞測(cè)序數(shù)據(jù)的拷貝數(shù)變異檢測(cè)方法(FL- CNV)，通過(guò)對(duì)不同檢測(cè)方法的優(yōu)勢(shì)進(jìn)行互補(bǔ)，以期解決現(xiàn)有檢測(cè)中遇到的精度問(wèn)題和假陽(yáng)性問(wèn)題，并在后續(xù)實(shí)驗(yàn)中驗(yàn)證了本文方法用于單細(xì)胞CNV檢測(cè)時(shí)精度的提高。

1 拷貝數(shù)變異檢測(cè)方法

1.1 工作流程

FL- CNV工作流程如圖1所示。首先該方法綜合了冷泉港實(shí)驗(yàn)室不定長(zhǎng)窗口的數(shù)據(jù)預(yù)處理思想[10]，以避免由于刪除低映射質(zhì)量序列而出現(xiàn)的假陽(yáng)性變異；其次利用構(gòu)造函數(shù)方法中負(fù)二項(xiàng)建模的降噪優(yōu)勢(shì)進(jìn)行數(shù)據(jù)處理，以修正在GC矯正中不能有效消除的數(shù)據(jù)偏好性；最后通過(guò)計(jì)算單細(xì)胞測(cè)序數(shù)據(jù)降噪后的當(dāng)前窗口及周邊一定數(shù)量窗口的讀段深度(read depth，RD)數(shù)據(jù)[14]，來(lái)確定數(shù)據(jù)異常范圍并得到斷點(diǎn)的位置及變異類型。確定斷點(diǎn)和變異類型的步驟執(zhí)行兩次，第一次粗略劃定可能存在變異的區(qū)域，第二次進(jìn)一步細(xì)化最終的變異區(qū)間。

圖1 FL- CNV工作流程Fig.1 The workflow for FL- CNV

檢測(cè)流程如下：1)數(shù)據(jù)預(yù)處理，使用模擬數(shù)據(jù)將參考基因分割成變長(zhǎng)窗口，通過(guò)Bowtie2將單細(xì)胞reads比對(duì)到參考基因，統(tǒng)計(jì)窗口真實(shí)數(shù)據(jù)RD，對(duì)窗口RD進(jìn)行GC矯正；2)RD降噪，對(duì)步驟1)輸出的RD數(shù)據(jù)進(jìn)行矯正，消除擴(kuò)增帶來(lái)的噪聲；3)變異檢測(cè)，綜合當(dāng)前窗口與前后相鄰若干窗口的RD值，通過(guò)兩次窗口掃描確定變異。

1.2 數(shù)據(jù)預(yù)處理

數(shù)據(jù)預(yù)處理的步驟如下。

(1)窗口劃分 為了避免單細(xì)胞數(shù)據(jù)覆蓋率低且分布不均勻的特性對(duì)變異檢測(cè)的影響，采用可變長(zhǎng)的窗口來(lái)劃分參考基因。首先以參考基因hg19為模版使用art_illumina軟件[15]生成數(shù)據(jù)預(yù)處理階段所需的模擬數(shù)據(jù)，再將模擬數(shù)據(jù)比對(duì)回參考基因(此處使用比對(duì)工具Bowtie2，參數(shù)設(shè)置與后續(xù)的真實(shí)數(shù)據(jù)比對(duì)相同)，根據(jù)比對(duì)回參考基因的模擬讀數(shù)的分布，將染色體劃分為具有相同數(shù)量的模擬讀數(shù)的窗口，即窗口劃分規(guī)則是每個(gè)窗口內(nèi)比對(duì)的模擬讀數(shù)相同。已有研究表明窗口劃分的最適數(shù)量為50 000左右[11]，故設(shè)定比對(duì)到每個(gè)窗口的模擬讀數(shù)數(shù)量為3 621，將全基因劃分為50 042個(gè)窗口。因?yàn)閱渭?xì)胞數(shù)據(jù)分布不均衡，所以劃分出的窗口長(zhǎng)短不一，格式如表1所示，其中，Chr為窗口所在染色體號(hào)，Start為窗口的起始位置，End為窗口終止位置，Size為窗口大小，Reads為落在該窗口的reads數(shù)。

表1 窗口文件Table 1 Window file

(2)比對(duì)得窗口讀段數(shù)量(read count，RC) 使用Bowtie2將單細(xì)胞的真實(shí)數(shù)據(jù)比對(duì)到參考基因，將窗口文件轉(zhuǎn)化為bed文件，將bed文件和比對(duì)后得到的Bam文件作為輸入，通過(guò)bedtools計(jì)算得到每個(gè)窗口的RC及GC含量。

(3)GC矯正 單細(xì)胞數(shù)據(jù)不均衡性的主要來(lái)源是整個(gè)基因組中GC含量的變化，為了對(duì)此進(jìn)行調(diào)整，使用deeptools根據(jù)GC含量對(duì)窗口內(nèi)的RC進(jìn)行矯正。

(4)歸一化 將GC矯正后得到的RC通過(guò)歸一化得到RD。

1.3 負(fù)二項(xiàng)分布消除RD噪聲

僅使用GC矯正不足以消除WGA帶來(lái)的偏差，還需要進(jìn)一步的矯正。常用的矯正方式有兩種：第一種是設(shè)置對(duì)照組，將對(duì)照組樣本測(cè)序數(shù)據(jù)作為參考，對(duì)每個(gè)窗口進(jìn)行矯正；第二種是通過(guò)分析噪聲特點(diǎn)有針對(duì)性地消除噪聲。已有研究表明，來(lái)自單細(xì)胞測(cè)序數(shù)據(jù)的RD信號(hào)的固有噪聲可以通過(guò)負(fù)二項(xiàng)(NB)分布準(zhǔn)確地表征[11]。現(xiàn)有的利用負(fù)二項(xiàng)分布建模的方式雖然在確定斷點(diǎn)和變異類型時(shí)效率低下，但在消除RD噪聲方面性能良好，對(duì)比之下選擇使用負(fù)二項(xiàng)分布對(duì)窗口RD進(jìn)行矯正。

1.4 變異檢測(cè)

構(gòu)造函數(shù)具有良好的數(shù)據(jù)矯正能力，但其對(duì)斷點(diǎn)的判斷缺乏全局性考量，在變異區(qū)域的融合和共享斷點(diǎn)的合并上存在明顯缺陷。為達(dá)到更好的邊界檢測(cè)效果，本文摒棄通過(guò)目標(biāo)函數(shù)求導(dǎo)確定斷點(diǎn)的方法，通過(guò)兩次掃描全局窗口的線性計(jì)算，高效準(zhǔn)確地劃分可能存在變異的區(qū)域，具體做法如圖2所示，步驟如下。

圖2 變異檢測(cè)模塊Fig.2 Variation detection module

(1)第一次窗口掃描 從頭依次掃描每個(gè)窗口，計(jì)算該窗口及其前10個(gè)窗口和后10個(gè)窗口矯正后RD值的三均值M3，若M3落在正常閾值外，則認(rèn)為當(dāng)前窗口內(nèi)可能存在發(fā)生拷貝數(shù)變異的區(qū)域。這里所述的正常閾值范圍為(N-N·m·n，N+N·m·n)，其中N為樣本倍性，n為所有窗口數(shù)據(jù)的標(biāo)準(zhǔn)差，m取1、2或3。在依次掃描全部窗口后，將變異類型相同的窗口標(biāo)記為一個(gè)變異區(qū)域，初步得到可能存在變異的區(qū)域。

這里使用的三均值是概率分布中的概念，計(jì)算公式如下。

M3=Q1/4+Q2/2+Q3/4

(1)

式中Q1、Q3為數(shù)據(jù)的兩個(gè)四分位點(diǎn)，Q2為中位數(shù)。三均值是一個(gè)有25%崩潰點(diǎn)的具有統(tǒng)計(jì)學(xué)抗性的L-估計(jì)，既對(duì)當(dāng)前窗口進(jìn)行了估計(jì)，又考慮到全局特性并排除了總體中的異常數(shù)值。

(2)第二次窗口掃描 以第一次劃分的某個(gè)變異區(qū)域?yàn)槔?，首先從變異區(qū)域的第一個(gè)窗口開(kāi)始，計(jì)算該窗口及其前后5個(gè)窗口RD值的平均值S2，定義第一個(gè)S2落在正常閾值外的窗口為第一斷點(diǎn)，依次向后掃描計(jì)算，當(dāng)出現(xiàn)連續(xù)兩個(gè)窗口變異類型與之前窗口的不同時(shí)，定義出現(xiàn)變異類型轉(zhuǎn)換的窗口的前一個(gè)窗口為第二斷點(diǎn)。這樣對(duì)某個(gè)初次劃分的可能存在變異的區(qū)域進(jìn)行二次掃描后，將其粗略估計(jì)的變異區(qū)域進(jìn)一步細(xì)化，得到最后的斷點(diǎn)和變異區(qū)間，如圖2(b)所示。

2 實(shí)驗(yàn)及結(jié)果分析

為了驗(yàn)證本文方法的優(yōu)越性，通過(guò)仿真數(shù)據(jù)實(shí)驗(yàn)和真實(shí)數(shù)據(jù)實(shí)驗(yàn)來(lái)評(píng)估方法的檢測(cè)效果。首先比較了FL- CNV與nbCNV[11]、control-FREEC[16]、CNVnator[17]在SCSsim[18]生成的模擬數(shù)據(jù)上的性能，然后在100個(gè)導(dǎo)管癌的真實(shí)數(shù)據(jù)上以倍體信息及aCGH的拷貝數(shù)劃分作為方法性能的真實(shí)性驗(yàn)證[19]。

2.1 模擬數(shù)據(jù)實(shí)驗(yàn)

2.1.1數(shù)據(jù)生成

在模擬實(shí)驗(yàn)方面，本文選取了千人基因組的人類hg19參考基因的21號(hào)染色體為模版，人為地隨機(jī)引入拷貝數(shù)變異事件。由于單細(xì)胞噪聲問(wèn)題突出，所以有效地模擬單細(xì)胞的數(shù)據(jù)特征對(duì)于單細(xì)胞測(cè)序工具的基準(zhǔn)測(cè)試來(lái)說(shuō)是非常必要的。不同于數(shù)據(jù)預(yù)處理階段，在模擬實(shí)驗(yàn)階段本文不再使用傳統(tǒng)仿真工具art_illumina，而是選擇專用于生成單細(xì)胞測(cè)序數(shù)據(jù)的仿真工具SCSsim來(lái)生成模擬實(shí)驗(yàn)所需的數(shù)據(jù)。

模擬數(shù)據(jù)的生成分為3步：首先通過(guò)SCSsim的simuvars功能定義拷貝數(shù)的變異位點(diǎn)及變異類型，每次模擬實(shí)驗(yàn)都隨機(jī)引入15個(gè)互不重疊的拷貝數(shù)變異事件生成fasta文件；之后基于由Illumina儀器生成的正常細(xì)胞的實(shí)際測(cè)序數(shù)據(jù)推斷測(cè)序概況，構(gòu)建profile配置文件，此處使用導(dǎo)管癌T10的一個(gè)單細(xì)胞數(shù)據(jù)，即索引號(hào)為SRR052047的單細(xì)胞測(cè)序樣本；最后基于fasta文件及profile文件得到fastq格式的序列片段文件。

2.1.2評(píng)價(jià)指標(biāo)

模擬實(shí)驗(yàn)通過(guò)精度Pre、敏感度Sen和F1系數(shù)來(lái)評(píng)估每個(gè)方法的性能，各參數(shù)的計(jì)算公式如下。

Pre=TP/(TP+FP)

(2)

Sen=TP/(TP+FN)

(3)

F1=(2Pre·Sen)/(Pre+Sen)

(4)

式中，TP表示真陽(yáng)性，F(xiàn)P表示假陽(yáng)性，F(xiàn)N表示假陰性。

2.1.3實(shí)驗(yàn)結(jié)果

本文進(jìn)行了100次模擬實(shí)驗(yàn)，共獲得了1 500個(gè)拷貝數(shù)基準(zhǔn)變異，使用FL- CNV、control-FREEC、nbCNV和CNVnator這4種方法對(duì)模擬數(shù)據(jù)進(jìn)行變異檢測(cè)，結(jié)果如表2所示。

表2 模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果Table 2 Experimental results with simulated data

由表2中數(shù)據(jù)可知，F(xiàn)L- CNV在模擬實(shí)驗(yàn)中表現(xiàn)良好，檢測(cè)精度明顯高于其他3種方法，與此同時(shí)敏感度也是4種方法里最高的?？梢?jiàn)本文方法在保證不損失敏感度的同時(shí)，顯著提高了檢測(cè)的精度，達(dá)到了較好的檢測(cè)效果。

2.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

由于單細(xì)胞缺少CNV基準(zhǔn)數(shù)據(jù)，所以現(xiàn)有的單細(xì)胞CNV檢測(cè)方法在進(jìn)行性能評(píng)估時(shí)，對(duì)方法檢測(cè)效果的驗(yàn)證大都使用倍體信息作為類標(biāo)簽對(duì)細(xì)胞個(gè)體數(shù)據(jù)進(jìn)行聚類分析，但這樣的驗(yàn)證方法并不能直觀地體現(xiàn)變異檢測(cè)效果，因?yàn)楸扼w信息僅代表了樣本序列拷貝數(shù)的平均概況，并不能量化單個(gè)細(xì)胞的變異檢測(cè)水平。為了更直觀地比較變異檢測(cè)效果，在真實(shí)數(shù)據(jù)實(shí)驗(yàn)部分，我們不但使用真實(shí)數(shù)據(jù)的倍體信息來(lái)衡量方法的好壞，還將Navin等[10]提出的比較基因組雜交的拷貝數(shù)劃分作為金標(biāo)準(zhǔn)來(lái)直觀比較方法性能。

2.2.1數(shù)據(jù)概況

真實(shí)數(shù)據(jù)為來(lái)自美國(guó)國(guó)家生物技術(shù)信息中心(National Center of Biotechnology Information, NCBI)的100個(gè)導(dǎo)管癌(T10)單細(xì)胞，編號(hào)為SRP002535(https:∥www.ncbi.nlm.nih.gov/sra/?term=SRP002535)。該數(shù)據(jù)集有47個(gè)二倍體D，24個(gè)亞二倍體H和29個(gè)異倍體，其中異倍體又可分為異倍體AA(3N)和異倍體AB(3.3N)。用于檢測(cè)效果驗(yàn)證的比較基因組雜交(array-based comparative genomic hybridization，aCGH)基準(zhǔn)數(shù)據(jù)來(lái)自NCBI，編號(hào)為GSE16607。

2.2.2實(shí)驗(yàn)結(jié)果

為了直觀地比較不同方法的檢測(cè)效果，將4種方法在這100個(gè)單細(xì)胞數(shù)據(jù)上的全基因拷貝數(shù)檢測(cè)結(jié)果整理為4個(gè)矩陣，生成如圖3所示的熱圖。其中每行代表一個(gè)細(xì)胞個(gè)體，每列代表一個(gè)檢測(cè)窗口，游標(biāo)為拷貝數(shù)與熱圖中每個(gè)窗口顏色的映射關(guān)系。

由圖3可以看出，F(xiàn)L- CNV和nbCNV能將100個(gè)單細(xì)胞分為明顯的3類，而control-FREEC與CNVnator對(duì)100個(gè)導(dǎo)管癌單細(xì)胞的聚類效果則不太理想，遠(yuǎn)不如FL- CNV和nbCNV的分類效果清晰明了。

圖3 不同方法的熱圖Fig.3 Heatmap displays of different methods

為進(jìn)一步直觀比較其他檢測(cè)方法與本文方法的優(yōu)劣，選擇在熱圖中表現(xiàn)良好的nbCNV與本文方法進(jìn)行皮爾遜相關(guān)系數(shù)(Pearson correlation coefficient)的比較，結(jié)果如圖4所示。由圖4可以看到，綜合不同亞群的相關(guān)性情況，F(xiàn)L- CNV整體來(lái)說(shuō)具有比nbCNV更高的皮爾遜相關(guān)系數(shù)，即FL- CNV檢測(cè)的拷貝數(shù)變異情況更符合真實(shí)數(shù)據(jù)的拷貝數(shù)變異情況，也即FL- CNV檢測(cè)出的拷貝數(shù)變異具有更高的可信度。

圖4 皮爾遜相關(guān)系數(shù)比較Fig.4 Pearson correlation coefficient comparison

綜合以上模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果可知，本文方法在用于拷貝數(shù)變異檢測(cè)時(shí)具有更高的精度和敏感度；在真實(shí)數(shù)據(jù)實(shí)驗(yàn)部分，本文不僅通過(guò)熱圖定性地比較了不同方法對(duì)單細(xì)胞數(shù)據(jù)的分類能力，并且通過(guò)皮爾遜相關(guān)系數(shù)定量地比較了所提方法和nbCNV對(duì)單細(xì)胞拷貝數(shù)變異的檢測(cè)能力，實(shí)驗(yàn)結(jié)果更加可靠。

3 結(jié)論

本文提出了一種基于單細(xì)胞測(cè)序數(shù)據(jù)的拷貝數(shù)變異檢測(cè)方法FL- CNV，根據(jù)單細(xì)胞測(cè)序數(shù)據(jù)的特點(diǎn)，動(dòng)態(tài)劃分變長(zhǎng)窗口并矯正窗口數(shù)據(jù)，對(duì)矯正后的窗口進(jìn)行兩次數(shù)據(jù)估算，這兩次數(shù)據(jù)估算既顧及到全局?jǐn)?shù)據(jù)概況又考慮了局部數(shù)據(jù)變化，有效地避免了異常數(shù)據(jù)點(diǎn)對(duì)變異檢測(cè)帶來(lái)的影響，減少假陽(yáng)性樣本的產(chǎn)生，提高了變異檢測(cè)方法的精度和敏感度，解決了現(xiàn)有檢測(cè)方法精度不高的問(wèn)題。

模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果驗(yàn)證了本文方法的優(yōu)越性，并且避免了單細(xì)胞實(shí)驗(yàn)驗(yàn)證因缺少基準(zhǔn)數(shù)據(jù)而采用倍體標(biāo)簽來(lái)衡量方法好壞的不準(zhǔn)確的結(jié)果檢驗(yàn)方式。以比較基因組雜交的拷貝數(shù)劃分作為金標(biāo)準(zhǔn)，通過(guò)熱圖和皮爾遜相關(guān)系數(shù)的比較進(jìn)一步證明了本文所提FL- CNV具有優(yōu)于其他檢測(cè)方法的精度、敏感度和可信度。