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

(中國海洋大學海洋生命學院，海洋生物遺傳學與育種教育部重點實驗室，山東 青島 266003)

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

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

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

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

1 缺失片段鑒定的原理

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

實際數(shù)據(jù)中可能的缺失情況如圖2a所示，當測序后由于測序數(shù)據(jù)量不足時，可能會出現(xiàn)3種情況：(1)部分非缺失的區(qū)段由于沒有被測到而形成缺失片段，這種情況在低深度數(shù)據(jù)中是常見而且不可避免的。(2)部分非缺失片段沒有被測到，但這些片段的上游或下游區(qū)域存在一個真實的缺失片段，從而導致這個缺失片段的延長。因為從組成上來說，被延長的缺失片段的確含有真實的缺失區(qū)間，因此在統(tǒng)計中不會影響缺失片段的檢出率，但會對缺失區(qū)間大小的判斷產(chǎn)生一定的影響。(3)因測序數(shù)據(jù)量不足導致部分非缺失片段沒有reads覆蓋，恰好將兩個真實的缺失片段連成一個大缺失片段，與情況(2)一樣，不會影響缺失片段的檢出率，但對正確判斷缺失區(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.基因組中真實缺失區(qū)域在測序數(shù)據(jù)中可能會出現(xiàn)的3種情況；b.檢測缺失區(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ū)間檢測流程如圖2b所示。為了降低缺失片段檢測的錯誤率，對比對結(jié)果中所有無reads覆蓋的酶切標簽進行篩選，當標簽中無reads覆蓋的連續(xù)標簽數(shù)(continuous deletion tags, cdt)大于某一個閾值cdtset時，認為該處存在缺失變異，反之認為該處是由于外部因素所導致的不可靠的缺失區(qū)域，在數(shù)據(jù)處理過程中cdtset設置為3。為了明確缺失區(qū)域的大小，我們?nèi)藶槎x了缺失區(qū)域的起始和終止位置，其中起始位置為該缺失區(qū)域前最近一個酶切標簽的末尾，終止位置是該缺失片段后第一個酶切標簽的起始，如圖3所示。判斷不同深度測序數(shù)據(jù)對缺失區(qū)域檢測的影響時，我們用錯誤率作為評估參數(shù)，其計算公式如下：錯誤率=(檢測到的缺失片段數(shù)-實際缺失片段數(shù))/檢測到的缺失片段數(shù)%。

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

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

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

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

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

3 效果評價

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

3.1 測序深度的影響

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

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

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

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

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

3.2 缺失片段大小的影響

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

圖4 100 kb大小的缺失文庫在抽取不同測序深度下處理前、后檢測到的缺失片段大小的分布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×測序數(shù)據(jù)中檢測情況Table 2 The detection summary of deletion regions with different sizes in 20× sequencing data

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

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

4 討論

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

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

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