利用2b-RAD技術(shù)檢測(cè)基因組區(qū)段缺失變異的應(yīng)用潛力評(píng)價(jià)*

程陶然，李語(yǔ)麗，劉平平，楊志輝，張玲玲，胡曉麗，包振民，王 師

(中國(guó)海洋大學(xué)海洋生命學(xué)院，海洋生物遺傳學(xué)與育種教育部重點(diǎn)實(shí)驗(yàn)室，山東 青島 266003)

結(jié)構(gòu)變異是個(gè)體基因組結(jié)構(gòu)上的微觀和亞微觀變異，通常包括缺失、重復(fù)、插入、拷貝數(shù)變異、倒置和易位等，一般結(jié)構(gòu)變異片段的大小介于1 kb到3 Mb之間[1]。從核苷酸的數(shù)量來看，結(jié)構(gòu)變異所含的核苷酸總數(shù)遠(yuǎn)遠(yuǎn)超過SNPs，從而豐富了DNA遺傳變異的多樣性[2]。缺失是結(jié)構(gòu)變異中一種常見的變異形式，是重要的變異來源之一。缺失變異是同正常染色體相比，個(gè)體染色體某一區(qū)段丟失，丟失區(qū)域大小從單個(gè)堿基到整條染色體。導(dǎo)致缺失變異的來源多種多樣，如易位丟失，減數(shù)分裂過程中染色體交叉錯(cuò)誤等，但總的歸結(jié)為是由不等的互換所導(dǎo)致，如圖1所示，由于不等的互換，將序列1上B區(qū)段轉(zhuǎn)移至1’序列上，造成序列1上B/B’區(qū)段的缺失。

缺失變異按照對(duì)個(gè)體的影響程度可分為3種：(1)影響微小。缺失片段處于非基因區(qū)，對(duì)個(gè)體的生長(zhǎng)發(fā)育影響微小。(2)導(dǎo)致疾病，但不影響個(gè)體生存。缺失片段處于基因區(qū)，造成部分基因丟失，對(duì)個(gè)體的生長(zhǎng)發(fā)育造成一定影響，如威廉斯綜合征[3-4]。此外，很多遺傳疾病也與缺失變異有關(guān)，如杜氏肌萎縮癥[5]等。(3)造成胚胎或幼苗的致死性危害。缺失片段處于重要的基因區(qū)，使個(gè)體生存的重要基因丟失，造成個(gè)體無法存活，如SMN編碼基因的缺失導(dǎo)致脊髓性肌肉萎縮癥(Spinal muscular atrophy, SMA)，是嬰兒死亡的最常見的遺傳原因[6]。目前全基因組缺失變異的分析平臺(tái)(如重測(cè)序、SNP芯片等)價(jià)格昂貴，難于實(shí)現(xiàn)對(duì)大量個(gè)體的低成本分析。簡(jiǎn)化基因組測(cè)序(Reduced-representation genome sequencing，RRGS)[7]將高通量測(cè)序技術(shù)[8]和限制性內(nèi)切酶的使用結(jié)合起來，通過對(duì)基因組部分位點(diǎn)進(jìn)行深度測(cè)序，大大降低了測(cè)序成本。其獲取的遺傳變異信息目前主要用于SNPs分型，對(duì)于結(jié)構(gòu)變異的探索還非常少。簡(jiǎn)化基因組測(cè)序的代表技術(shù)有RAD-seq(Restriction-site associated DNA sequencing)[9-10]、2b-RAD(IIb restriction site-associated DNA)[11]、ddRAD-seq(“Double-digest” restriction site-associated DNA sequencing)[12]和GBS(Genotyping-by-sequencing)[13]，其中，由于RAD-seq酶切所產(chǎn)生的限制性片段長(zhǎng)短不一，影響測(cè)序位點(diǎn)的基因組代表性及位點(diǎn)間測(cè)序深度的均勻度[14]，從而增大結(jié)構(gòu)變異檢測(cè)的難度。而2b-RAD技術(shù)一方面繼承了RAD-seq的大部分優(yōu)點(diǎn)，同時(shí)通過改進(jìn)其所使用的酶(BsaXI)來獲得片段長(zhǎng)度均一的酶切標(biāo)簽，保證了測(cè)序位點(diǎn)的基因組代表性和測(cè)序深度的均勻性，從而為結(jié)構(gòu)變異研究提供了可能性。

圖1 缺失變異產(chǎn)生的主要來源：不等互換Fig.1 The main generation source of deletion variation: Unequal crossing over

本文選擇穩(wěn)定可靠的2b-RAD測(cè)序數(shù)據(jù)，討論2b-RAD數(shù)據(jù)不同測(cè)序深度和不同缺失區(qū)域大小對(duì)檢測(cè)缺失變異的影響及解決方法，并用擬南芥2b-RAD半模擬數(shù)據(jù)進(jìn)行驗(yàn)證。

1 缺失片段鑒定的原理

全基因組測(cè)序中已證實(shí)基因組中某一序列下讀取reads深度可以反映該位點(diǎn)在基因組的情況，如拷貝數(shù)變異與reads深度成正比[15]，同理該理論可利用每個(gè)位點(diǎn)的深度信息推測(cè)基因組中缺失區(qū)域，即reads深度為零的區(qū)域?yàn)槿笔^(qū)域。然而簡(jiǎn)化基因組測(cè)序的酶切標(biāo)簽在基因組間斷性分布，導(dǎo)致數(shù)據(jù)中僅獲得部分位點(diǎn)的深度信息。理論上對(duì)于某個(gè)體基因組的缺失片段，若該片段中含有2b-RAD的酶切標(biāo)簽，在數(shù)據(jù)比對(duì)時(shí)，該標(biāo)簽下reads深度為零。所以從比對(duì)結(jié)果中基因組酶切標(biāo)簽的reads深度出發(fā)，若某標(biāo)簽下reads深度為零，說明測(cè)序個(gè)體中該標(biāo)簽所在的片段可能是缺失的。若相鄰的幾個(gè)標(biāo)簽下reads深度均為零，則認(rèn)為這些標(biāo)簽所在的片段是缺失片段。但實(shí)際情況中，簡(jiǎn)化基因組數(shù)據(jù)的測(cè)序深度往往對(duì)酶切標(biāo)簽深度有較大影響，如某標(biāo)簽下reads深度為零，并不能說明該標(biāo)簽一定是缺失的，而可能由測(cè)序錯(cuò)誤或測(cè)序深度不足導(dǎo)致。

實(shí)際數(shù)據(jù)中可能的缺失情況如圖2a所示，當(dāng)測(cè)序后由于測(cè)序數(shù)據(jù)量不足時(shí)，可能會(huì)出現(xiàn)3種情況：(1)部分非缺失的區(qū)段由于沒有被測(cè)到而形成缺失片段，這種情況在低深度數(shù)據(jù)中是常見而且不可避免的。(2)部分非缺失片段沒有被測(cè)到，但這些片段的上游或下游區(qū)域存在一個(gè)真實(shí)的缺失片段，從而導(dǎo)致這個(gè)缺失片段的延長(zhǎng)。因?yàn)閺慕M成上來說，被延長(zhǎng)的缺失片段的確含有真實(shí)的缺失區(qū)間，因此在統(tǒng)計(jì)中不會(huì)影響缺失片段的檢出率，但會(huì)對(duì)缺失區(qū)間大小的判斷產(chǎn)生一定的影響。(3)因測(cè)序數(shù)據(jù)量不足導(dǎo)致部分非缺失片段沒有reads覆蓋，恰好將兩個(gè)真實(shí)的缺失片段連成一個(gè)大缺失片段，與情況(2)一樣，不會(huì)影響缺失片段的檢出率，但對(duì)正確判斷缺失區(qū)間的大小帶來影響。

(①Genome;②Sequencing data;③Ref tags of genome;④Reads of sequencing data;⑤Real deletion regions;⑥False deletion regions;⑦Data;⑧Reads mapping;⑨Reads depth;⑩Calculate continuous deletion tags;Unrelible decetion regions;Reliable decetion regions.a.基因組中真實(shí)缺失區(qū)域在測(cè)序數(shù)據(jù)中可能會(huì)出現(xiàn)的3種情況；b.檢測(cè)缺失區(qū)間的流程。a.The three cases of the real deletion regions in the genome may appear in the sequencing data; b.Process for detecting deletion regions.)
圖2 缺失片段的鑒定原理和處理流程
Fig.2 The identification principle and process flow of deletion regions

綜合考慮以上因素，研究中制定的缺失區(qū)間檢測(cè)流程如圖2b所示。為了降低缺失片段檢測(cè)的錯(cuò)誤率，對(duì)比對(duì)結(jié)果中所有無reads覆蓋的酶切標(biāo)簽進(jìn)行篩選，當(dāng)標(biāo)簽中無reads覆蓋的連續(xù)標(biāo)簽數(shù)(continuous deletion tags, cdt)大于某一個(gè)閾值cdtset時(shí)，認(rèn)為該處存在缺失變異，反之認(rèn)為該處是由于外部因素所導(dǎo)致的不可靠的缺失區(qū)域，在數(shù)據(jù)處理過程中cdtset設(shè)置為3。為了明確缺失區(qū)域的大小，我們?nèi)藶槎x了缺失區(qū)域的起始和終止位置，其中起始位置為該缺失區(qū)域前最近一個(gè)酶切標(biāo)簽的末尾，終止位置是該缺失片段后第一個(gè)酶切標(biāo)簽的起始，如圖3所示。判斷不同深度測(cè)序數(shù)據(jù)對(duì)缺失區(qū)域檢測(cè)的影響時(shí)，我們用錯(cuò)誤率作為評(píng)估參數(shù)，其計(jì)算公式如下：錯(cuò)誤率=(檢測(cè)到的缺失片段數(shù)-實(shí)際缺失片段數(shù))/檢測(cè)到的缺失片段數(shù)%。

(①Genome;②Restriction digestion tags;③Artifically specificed deletion region size;④Detected deletion region)
圖3 缺失片段大小的計(jì)算
Fig.3 Calculation of deletion regions’size

2 擬南芥半模擬數(shù)據(jù)準(zhǔn)備

首先利用RADTyping軟件[16]中的Extract_cut_site.pl腳本在擬南芥基因組序列(https://www.ncbi.nlm.nih.gov/genome/?term=Arabidopsis%20thaliana)中提取出39 678個(gè)BsaXI酶的基因組酶切標(biāo)簽(ref標(biāo)簽)用作后續(xù)序列比對(duì)的參考序列。

本文采用的擬南芥2b-RAD數(shù)據(jù)來源于Wang等[17]的5標(biāo)簽串聯(lián)的2b-RAD實(shí)驗(yàn)，且該擬南芥的Multi-isoRAD文庫(kù)(SRP068382)含有同一個(gè)擬南芥?zhèn)€體的5個(gè)串聯(lián)文庫(kù)，依次命名為AL1-AL5，每個(gè)文庫(kù)的測(cè)序深度均在200×左右。隨機(jī)挑選AL4文庫(kù)作為測(cè)試數(shù)據(jù)，為了避免文庫(kù)自身由于測(cè)序不足產(chǎn)生的零reads覆蓋標(biāo)簽對(duì)檢測(cè)的影響，人為將AL4文庫(kù)中59個(gè)零覆蓋標(biāo)簽補(bǔ)齊，標(biāo)簽深度按照文庫(kù)平均深度的泊松分布隨機(jī)分配。

對(duì)于已不存在零覆蓋標(biāo)簽的AL4文庫(kù)，分別做以下處理：(1)分別構(gòu)建5個(gè)不同大小(5、10、50、100和500 kb)的缺失片段文庫(kù)。對(duì)于每種規(guī)格的缺失片段，分別在擬南芥基因組上隨機(jī)選取50個(gè)無重疊的區(qū)域(每條染色體選10個(gè))，并要求這些區(qū)域內(nèi)至少含有3個(gè)或3個(gè)以上的ref標(biāo)簽。利用SOAP2軟件[18]的比對(duì)結(jié)果，剔除AL4文庫(kù)中比對(duì)到人為選擇的片段含有的ref標(biāo)簽下的reads，分別生成含有真實(shí)缺失片段的5個(gè)缺失文庫(kù)。2)在1)所得數(shù)據(jù)中按一定梯度對(duì)測(cè)序數(shù)據(jù)進(jìn)行抽樣。將上述5個(gè)缺失文庫(kù)分別抽取5×-100×生成梯度測(cè)序深度數(shù)據(jù)，并用SOAP2軟件與基因組ref標(biāo)簽比對(duì)，通過比對(duì)結(jié)果獲取ref標(biāo)簽的reads覆蓋深度信息，統(tǒng)計(jì)連續(xù)無reads覆蓋的標(biāo)簽區(qū)域。

3 效果評(píng)價(jià)

本文利用擬南芥半模擬數(shù)據(jù)分別討論了不同測(cè)序深度和不同大小缺失片段對(duì)缺失變異區(qū)域檢測(cè)的準(zhǔn)確性的影響。

3.1 測(cè)序深度的影響

以片段大小為100 kb的缺失文庫(kù)的梯度測(cè)序數(shù)據(jù)為例，研究不同測(cè)序深度對(duì)檢測(cè)缺失變異準(zhǔn)確性的影響。檢測(cè)到的連續(xù)無reads比對(duì)的酶切標(biāo)簽區(qū)域(即缺失片段)及與根據(jù)已知缺失片段計(jì)算的檢測(cè)錯(cuò)誤率的統(tǒng)計(jì)結(jié)果如表1所示。

表1 擬南芥100 kb缺失文庫(kù)的梯度測(cè)序數(shù)據(jù)缺失片段檢測(cè)的統(tǒng)計(jì)結(jié)果Table 1 The detection statistics of deletion regions in the gradient sequencing data extracted from A.thaliana 100 kb deletion library

從表1中可以看出，未處理時(shí)，測(cè)序深度越淺，數(shù)據(jù)量越少，大部分標(biāo)簽未被測(cè)到，導(dǎo)致了較高的缺失片段檢測(cè)的錯(cuò)誤率，如5×、10×?xí)r，其錯(cuò)誤率高達(dá)90%以上。當(dāng)測(cè)序深度增大時(shí)，檢測(cè)到的缺失片段數(shù)目逐漸趨近于真實(shí)缺失片段數(shù)，錯(cuò)誤率有所降低，但當(dāng)測(cè)序深度達(dá)到100×?xí)r錯(cuò)誤率仍高達(dá)45%以上。這說明直接將連續(xù)無reads比對(duì)的ref標(biāo)簽所在區(qū)域視為缺失區(qū)域，會(huì)帶來很高的錯(cuò)誤率。

為了保證檢測(cè)缺失片段的準(zhǔn)確性，對(duì)連續(xù)無reads比對(duì)的ref標(biāo)簽區(qū)域進(jìn)行過濾篩選，要求連續(xù)無reads比對(duì)的標(biāo)簽數(shù)cdt大于某一個(gè)閾值cdtset時(shí)，認(rèn)為該處存在的缺失變異是可靠的。當(dāng)cdtset=3，即cdt≥3時(shí)，檢測(cè)到的缺失片段數(shù)接近真實(shí)缺失片段數(shù)，錯(cuò)誤率也隨之降低。例如5×測(cè)序數(shù)據(jù)未過濾時(shí)檢測(cè)到2 813個(gè)缺失片段，錯(cuò)誤率為98.84%，cdt≥3過濾處理后，檢測(cè)到的片段數(shù)目為70，錯(cuò)誤率降低至28.57%。表1梯度測(cè)序深度下檢測(cè)缺失片段的結(jié)果顯示，當(dāng)測(cè)序深度達(dá)到20×以上時(shí)，cdt≥3條件下的錯(cuò)誤率均降到0，說明cdt≥3條件下20×的測(cè)序深度即可正確檢測(cè)測(cè)序數(shù)據(jù)中缺失變異的區(qū)間。

理論上真實(shí)缺失片段大小約為100 kb，所以梯度測(cè)序深度數(shù)據(jù)統(tǒng)計(jì)的缺失片段大小范圍應(yīng)均在100 kb左右。圖4-a中顯示，當(dāng)未作任何處理時(shí)，造成大量無reads比對(duì)的區(qū)域被誤判成缺失片段，導(dǎo)致箱型圖波動(dòng)較大，尤其是當(dāng)測(cè)序深度過低時(shí)，檢測(cè)到的缺失片段大小的中位值被拉低，遠(yuǎn)遠(yuǎn)低于100 kb。當(dāng)要求連續(xù)缺失標(biāo)簽數(shù)cdt≥3時(shí)，如圖4b中所示，除了5×和10×數(shù)據(jù)的結(jié)果偏離100 kb之外，20×及以上的數(shù)據(jù)的缺失片段大小基本在100 kb左右，表明檢測(cè)出來的缺失片段基本和真實(shí)的缺失片段相符合，從而進(jìn)一步證實(shí)了方法的準(zhǔn)確性。

3.2 缺失片段大小的影響

為了研究不同缺失片段大小對(duì)檢測(cè)缺失變異準(zhǔn)確性的影響，將5、10、50、100和500 kb的缺失文庫(kù)的梯度測(cè)序數(shù)據(jù)分別用soap軟件進(jìn)行比對(duì)，獲取擬南芥基因組ref標(biāo)簽下reads覆蓋深度信息，并統(tǒng)計(jì)連續(xù)無reads比對(duì)的標(biāo)簽區(qū)間。表2是真實(shí)缺失片段5、10、50、100和500 kb在測(cè)序深度約為20×左右時(shí)，檢測(cè)到缺失片段的統(tǒng)計(jì)情況。未處理時(shí)缺失片段檢測(cè)的錯(cuò)誤率隨著真實(shí)缺失片段大小的增大略有減小，但都在85%以上，而cdt≥3時(shí)，檢測(cè)的缺失片段數(shù)均為50，錯(cuò)誤率降至0。而圖5顯示，20×的測(cè)序數(shù)據(jù)cdt≥3條件下檢測(cè)到的缺失片段大小的分布與標(biāo)準(zhǔn)缺失片段大小基本吻合。

圖4 100 kb大小的缺失文庫(kù)在抽取不同測(cè)序深度下處理前、后檢測(cè)到的缺失片段大小的分布Fig.4 The size distribution of deletion regions detected before treatment and after treatment at different sequencing depth extract by 100 kb deletion data

表2 不同缺失片段大小下20×測(cè)序數(shù)據(jù)中檢測(cè)情況Table 2 The detection summary of deletion regions with different sizes in 20× sequencing data

Note:①Unfitter;②Detected deletion regons;③Error rate

圖5 20×測(cè)序數(shù)據(jù)cdt≥3時(shí)檢測(cè)出的缺失片段大小的變化范圍Fig.5 Size distribution of detected deletion regions at 20×sequencing depth with cdt≥3

4 討論

缺失變異是近年來熱門的研究之一，但傳統(tǒng)的基于SNP芯片和重測(cè)序方法的分析平臺(tái)成本過高，不適于大量個(gè)體的缺失變異分析。簡(jiǎn)化基因組測(cè)序?qū)⒏咄繙y(cè)序技術(shù)和限制性內(nèi)切酶的使用結(jié)合起來，通過對(duì)基因組部分位點(diǎn)進(jìn)行深度測(cè)序，大大降低了測(cè)序成本，而以往簡(jiǎn)化基因組測(cè)序數(shù)據(jù)主要應(yīng)用于SNP分型，其中的基因組遺傳結(jié)構(gòu)變異等信息沒有被充分發(fā)掘。

本研究選擇了測(cè)序深度均一的2b-RAD數(shù)據(jù)，并克服了簡(jiǎn)化基因組測(cè)序數(shù)據(jù)中酶切標(biāo)簽在基因組分布上的間斷性給缺失變異研究帶來的影響。以擬南芥2b-RAD測(cè)序數(shù)據(jù)為研究對(duì)象，人為生成包含已知缺失區(qū)域的半模擬數(shù)據(jù)，通過對(duì)不同測(cè)序深度和不同缺失片段大小的缺失文庫(kù)分析發(fā)現(xiàn)，當(dāng)測(cè)序深度達(dá)到20×左右，設(shè)定缺失區(qū)域所含連續(xù)標(biāo)簽數(shù)cdt≥3時(shí)，即可準(zhǔn)確地鑒定測(cè)序數(shù)據(jù)中的缺失區(qū)域，且對(duì)于不同大小的缺失片段的檢測(cè)均有效。雖然該方法能大大降低由于抽樣造成的假陽(yáng)性，但是由于我們方法設(shè)置的問題，對(duì)于小片段(<10 kb)缺失變異的檢測(cè)仍有一定難度，因此小片段缺失的檢測(cè)仍需要借助傳統(tǒng)方法如重測(cè)序等。

本文通過對(duì)半模擬數(shù)據(jù)的分析證實(shí)了2b-RAD在缺失變異研究中的可行性，為2b-RAD簡(jiǎn)化基因組數(shù)據(jù)檢測(cè)測(cè)序個(gè)體中的缺失區(qū)域提供了理論指導(dǎo)。在不增加測(cè)序成本的條件下，既節(jié)約了分析缺失變異的成本，又實(shí)現(xiàn)了數(shù)據(jù)的充分利用，為2b-RAD技術(shù)的應(yīng)用開辟了新方向。