基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法
2023-10-28 20:02:02 2
專利名稱：：基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法
技術領域：
：本發明涉及一種基因組長度多態性標記的引物設計方法。具體為基於酶切建庫雙末端(Pair-end)測序的長度多態性標記的弓I物設計開發方法；在缺少參考序列的情況下，尋找到個體間的Indel標記，並能夠在兩端設計引物。屬於生物信息學
技術領域：
。這對於缺少參考序列的非模式生物的研究具有重要的意義。
背景技術：
：InDel(insertion-deletion)插入缺失標記，指的是兩種親本中在基因組上的差異，相對另一個親本而言，其中一個親本的基因組中有一定數量的核苷酸插入或缺失。Indel位點信息的獲得可以有許多重要的應用，如構建遺傳圖譜，基因分型，分子標記輔助育種，疾病檢測等。如今，第二代DNA測序技術是一種高通量低成本的測序技術，基本原理是邊合成邊測序。以solexa測序方法為例，先用物理方法將DNA鏈隨機打斷，然後在片段兩端加上特定接頭，接頭上有擴增引物序列。測序時，DNA聚合酶合成待測片段的互補鏈，通過檢測新合成鹼基所攜帶的螢光信號讀取鹼基序列，從而獲得待測片段的序列。第二代測序技術已經廣泛應用於生物科學的許多領域，特別是研究一個物種不同個體之間的多態性。傳統CallIndel標記的方法是將測序個體得到的短reads通過比對軟體比對回參考序列，從而得到測序個體的Indel信息。常見的流程有使用BWA軟體將reads比對回參考序列，使用SAMtools軟體處理比對結果尋找Indel位點^2。大體過程如圖I所示。目前，有參考序列的物種都可以很方便的進行Indel標記的查找，並在兩端設計引物進行實驗驗證。但是對於那些非模式生物而言，基本上是不存在參考序列的。而在沒有參考序列的情況下，傳統尋找Indel標記的方法存在著技術上的瓶頸。I.LiH.andDurbinR.(2009)FastandaccurateshortreadalignmentwithBurrows-WheelerTransform.Bioinformatics,25:1754-60.[PMID:19451168]2.LiH.*，HandsakerB.WysokerA.,FennellT.,RuanJ.,HomerN.,MarthG.,AbecasisG.,DurbinR.and1000GenomeProjectDataProcessingSubgroup(2009)TheSequencealignment/map(SAM)formatandSAMtools.Bioinformatics,25,2078-9.[PMID:19505943]RAD-PE測序技術採用了新的建庫方式(酶切建庫)，其測序具體過程如圖2所示，用限制性內切酶切斷DNA特定的位點，再用物理方法將酶切之後的DNA分子隨機打斷，通過瓊脂糖膠DNA分離技術挑選特定長度的DNA分子，然後在挑選出來的DNA末端添加特定的擴增接頭與測序接頭，從而構建上機文庫進行高通量測序。其中RAD測序方法為本領域公知的方法，例如可參考以下文獻(I)MichaelRMiller,TressaSAtwood,BFrankEames,etal,RADmarkermicroarraysenablerapidmappingofzebrafishmutations,GenomeBiology,2007,8(6):R105.1-R105.10;(2)MichaelR.Miller,JosephP.Dunham,AngelAmores,etal,Rapidandcost-effectivepolymorphismidentificationandgenotypingusingrestrictionsiteassociatedDNA(RAD)markers,GenomeResearch,2007,17,240-248；(3)NathanA.Bairdl,PaulD.Etter,TressaS.Atwood,etal,RapidSNPDiscoveryandGeneticMappingUsingSequencedRADMarkers,PLoSONE,2008,3(10),e3376,doi:10.1371/journal,pone.0003376.散列表(Hashtable,或哈希表),是根據關鍵碼值(Keyvalue)而直接進行訪問的數據結構。也就是說，它通過把關鍵碼值映射到表中一個位置來訪問記錄，以加快查找的速度。這個映射函數叫做散列函數，存放記錄的數組叫做散列表。使用哈希表對數據進行索引基本是隨著數據量的上升線性增長，而且由「ATCGN」構成的字符串，鍵值出現衝突的可能性非常低。這樣在處理海量測序數據的時候有著很好的性能。
發明內容本發明的目的是提供一種基於酶切建庫雙末端(Pair-end)測序的長度多態性標記的引物設計開發方法；它是一種通過處理基於酶切建庫pair-end測序(RAD-PE測序技術)得到的測序數據，在兩個個體之間尋找長度多態性位點，並能夠在兩端側翼序列設計引物進行實驗驗證的技術方案。本發明的目的通過以下技術方案來實現基於酶切建庫雙末端(Pair-end)測序的長度多態性標記的引物設計開發方法，其步驟如下I)在獲得RAD高通量測序技術的測序結果後，對RAD雙末端測序序列進行過濾以去除不合格的測序序列。其中，RAD高通量測序技術可以為IlluminaGA測序技術，也可以為現有的其他高通量測序技術。所述的不合格的測序序列為測序質量低於預定的低質量閾值的鹼基個數超過整條序列鹼基個數的50%的序列。2)根據測序個體基因組酶切一端的測序序列，利用序列的全同性生成每個個體堆的信息。例如，將每個個體過濾後的酶切一端的測序序列信息作為哈希的鍵，哈希的值指向一個鍊表，用於存放另一端的序列信息，並計算測序深度信息。可用任何一種程式語言實現該過程。3)過濾掉酶切一端序列測序深度為I的結果(成對過濾)。4)兩個個體內分別將酶切一端的測序序列數據進行不容許空隙的兩兩比對，對堆進行聚類以確定個體內在酶切一端序列上的雜合SNP信息。所述的不容許空隙的兩兩比對是指比對的時候不容許開空位。其中，只有一個堆的聚類結果表明在酶切一端測序片段上不存在雜合位點，只有兩個堆的聚類結果表明在酶切一端測序片段上存在雜合位點，一般情況下這個雜合位點不會處於重複區域。對於那些堆的個數超過兩個的聚類結果，通常由於酶切一端測序序列處於基因組的重複區域所造成的，因此這些聚類結果將會被過濾掉。所使用的比對軟體可以是任何一款序列比對軟體，如blast、blat等。計算堆的個數為一和二的聚類結果的總深度，並進一步過濾掉低深度和高深度的聚類結果。低深度的閾值通常為平均測序深度的四分之一，高深度的閾值通常為平均測序深度的兩倍。5)在兩個個體內部，對每個堆的另一端數據進行局部組裝，採用的組裝可以是任何一款組裝軟體，如基於重疊關係的組裝軟體phrap,如基於DeBruijngraph算法的組裝軟體SOAPDenovo。利用重疊關係進行組裝的時候通常不容許空位的存在，其他參數可採用軟體默認設置。利用DeBruijngraph算法進行組裝的時候,通常kmer的大小要在21到31之間，其他參數可採用軟體默認設置。6)利用兩個個體酶切一端的測序序列信息將兩個個體堆的信息相互進行兩兩對齊，即在個體A和個體B中，個體A的某個堆能夠和個體B的某個堆對齊，若且唯若兩個個體堆中的酶切一端的測序序列完全相同。對能夠對齊的堆，兩個個體之間的另一端的組裝結果序列相互進行比對，來尋找Indel位點信息。比對的時候，兩個個體之間容許的空位數為一、並且公共區域內不存在錯配。使用的比對軟體可以是任何一款序列比對軟體，如blast、blat等。通過以上的步驟，將會得到兩個個體之間，高可信度的Indel位點信息，還得到了在Indel位點周圍的側翼序列信息。7)最後可在Indel位點周圍的側翼序列上設計引物，以便於後續的大規模實驗應用。如構建遺傳圖譜，個體基因分型等。本發明的技術方案採用了生物信息學分析方法，處理RAD(restriction-siteassociatedDNA,whichareshortfragmentsofDNAadjacenttoeachinstanceofaparticularrestrictionenzymerecognitionsite)雙末端測序的測序數據，從而尋找RAD測序片段上的Indel位點信息，以突破非模式生物缺少參考序列的瓶頸，簡化了基因組的複雜度，同時也減少了測序成本。常規的過程針對的物種通常為二倍體生物，但並不僅局限於二倍體生物。圖I為傳統CallIndel標記方法的過程原理圖；圖2為RAD-PE測序技術的測序具體過程原理圖；圖中，(A)限制性內切酶酶切基因組DNA，並加上Pl接頭，每一個Pl接頭含有不同的序列標籤；(B)帶有不同Pl接頭的樣品混合在一起，打斷；(C)加上接頭P2；(D)擴增富集RADtags；圖3為RAD雙末端測序的例圖；圖4為生成堆的聚類過程圖；圖5為堆中酶切一端序列信息示意圖；圖6為堆中另一端序列信息示意圖；圖7為堆中測序數據過濾和整合流程示意圖；圖8為局部組裝和Indel位點查找流程示意圖；圖9為引物設計方式示意圖10為RAD序列測序深度分布圖。具體實施例方式下面結合附圖與具體實施例進一步闡述本發明的技術特點。如圖2所示，獲得兩個個體基因組的RAD雙末端(pair-end)測序序列。圖3為RAD雙末端測序的例圖。在圖3中顯示了用限制性內切酶Ecorl，識別DNA分子上「G~AATTC」的迴文序列，並在G與A之間將DNA分子切斷，將酶切後的DNA分子用物理方法打斷成短的序列片段，並在含有酶切末端序列的DNA片段兩端加上接頭，PCR富集並進行雙末端(pair-end)測序,測序讀長一般為IOOnt,也可以為50nt。通常，回收300bp到500bp的DNA片段進行測序，這樣PE測序的左端由於酶切的原因是對齊的。另一端測序片段由於物理打斷時DNA分子長度的隨機性，片段之間就會存在著overlap關係，可以進行局部組裝，獲得較長(20(T300bp)的DNA片段，並在組裝結果上查找個體間的長度多態性位點，即Indel標記。基於酶切建庫雙末端(Pair-end)測序的長度多態性標記的引物設計開發方法，其步驟如下I)在獲得RAD高通量測序技術的測序結果後，對RAD雙末端測序序列進行過濾以去除不合格的測序序列。其中，RAD高通量測序技術可以為IlluminaGA測序技術，也可以為現有的其他高通量測序技術。所述的不合格的測序序列為測序質量低於預定的低質量閾值的鹼基個數超過整條序列鹼基個數的50%的序列。低質量閾值由具體測序技術及測序環境而定，例如設定為單鹼基測序質量低於20;測序序列中測序結果不確定的鹼基(如IlluminaGA測序結果中的N)個數超過整條測序序列鹼基個數的10%則認為是不合格序列；除樣本接頭序列外，與其它實驗引入的外源序列比對，如各種接頭序列。若序列中存在外源序列則認為是不合格序列；在酶切一端測序序列中，若起始的幾個鹼基不是酶切末端序列則過濾掉(如限制性內切酶Ecorl，測序序列開頭若不是「AATTC」則過濾掉整個測序序列)。2)根據測序個體基因組酶切一端的測序序列，利用序列的全同性生成每個個體堆的信息。例如，將每個個體過濾後的酶切一端的測序序列信息作為哈希的鍵，哈希的值指向一個鍊表，用於存放另一端的序列信息，並計算測序深度信息。可用任何一種程式語言實現該過程。具體過程如圖4所示。堆(Stack)中酶切一端的序列信息可以以圖5的方式保存，在圖5中,第一列表不的是酶切一端的序列信息；第二列表不的是酶切一端序列被測序的次數，即測序深度信息；第三列是該堆的ID，用於唯一確定堆。堆(Stack)中另一端的序列信息可以以圖6的方式保存，類FASTA格式。在圖6中，大於號(>)之後的是堆的ID，用於唯一確定堆。大於號(>)中的內容為堆中另一端的所有測序序列信息。3)過濾掉酶切一端序列測序深度為I的結果(成對過濾)。即過濾掉圖5中第二列為I的堆，同時過濾圖6中對應的數據。深度為I的結果通常是由測序錯誤導致的，過濾掉深度為I的測序序列信息，進一步減少由測序錯誤所帶來的分析困難。4)兩個個體內分別將酶切一端的測序序列數據進行不容許空隙的兩兩比對，對堆進行聚類以確定個體內在酶切一端序列上的雜合SNP信息。所述的不容許空隙的兩兩比對是指比對的時候不容許開空位。即不考慮開空位比對上的情況，例如以下兩條序列的比對結果就不滿足不容許空隙的兩兩比對條件序列IAATTCATCGACo序列2AACATCGTC。比對的時候容許的錯配數隨測序的長度來定，例如在測序長度小於50nt的情況下，容許的錯配數為1，長度在IOOnt的情況下，容許的錯配數為2。具體地，兩條序列之間只有一個鹼基不相同，則這兩條序列歸為一類。如果A序列與B序列之間只有一個鹼基不相同，而B與C之間只有另外一個鹼基不相同，則三條序列歸為一類，以此類推，通過所有測序序列之間的比對，可以將所有滿足比對條件的測序序列進行聚類。挑選出聚類結果中只有一個堆和兩個堆的聚類結果。其中只有一個堆的聚類結果表明在酶切一端測序片段上不存在雜合位點，只有兩個堆的聚類結果表明在酶切一端測序片段上存在雜合位點，一般情況下這個雜合位點不會處於重複區域。對於那些堆的個數超過兩個的聚類結果，通常由於酶切一端測序序列處於基因組的重複區域所造成的，因此這些聚類結果將會被過濾掉。使用的比對軟體可以是任何一款序列比對軟體，如blast、blat等。計算堆的個數為一和二的聚類結果的總深度，並進一步過濾掉低深度和高深度的聚類結果。低深度的閾值通常為平均測序深度的四分之一，高深度的閾值通常為平均測序深度的兩倍。按以上步驟對個體內的數據進行過濾和整合，重新生成如圖5和圖6的結果。酶切一端測序序列存在雜合位點的另一端信息將會被整合到一起，而酶切一端的測序序列信息將會用一致性序列來表示(如某個位點為A/T雜合，則用字母W表示)，並累加深度，重命名ID以利於後續的分析。這樣，處理之後的數據將更利於另一端數據的組裝。過程如圖8所示。5)在兩個個體內部，對每個堆的另一端數據進行局部組裝，採用的組裝可以是任何一款組裝軟體，如基於重疊關係的組裝軟體phrap,如基於DeBruijngraph算法的組裝軟體SOAPDenovo。利用重疊關係進行組裝的時候通常不容許空位的存在，其他參數可採用軟體默認設置。利用DeBruijngraph算法進行組裝的時候,通常kmer的大小要在21到31之間，其他參數可採用軟體默認設置。參數的意義和用法查看相關使用說明。6)利用兩個個體酶切一端的測序序列信息將兩個個體堆的信息相互進行兩兩對齊，即在個體A和個體B中，個體A的某個堆能夠和個體B的某個堆對齊，若且唯若兩個個體堆中的酶切一端的測序序列完全相同。對能夠對齊的堆，兩個個體之間的另一端的組裝結果序列相互進行比對，來尋找Indel位點信息。比對的時候，兩個個體之間容許的空位數為一、並且公共區域內不存在錯配。使用的比對軟體可以是任何一款序列比對軟體，如blast、blat等。過程如圖8所示。通過以上的步驟，將會得到兩個個體之間，高可信度的Indel位點信息，還得到了在Indel位點周圍的側翼序列信息。7)最後可在Indel位點周圍的側翼序列上設計引物，以便於後續的大規模實驗應用。如構建遺傳圖譜，個體基因分型等。如圖9所示。實施例數據瓠瓜F2群體兩個親本的RAD-PE測序數據。(說明父本和母本雜交生成的後代叫F1,F1自交生成的後代叫F2)材料來源浙江省農業科學院。實施例具體操作流程將兩個親本RAD-PE測序得到的測序數據，根據測序質量值，N的含量，以及是否含有酶切末端序列進行過濾，去除不合格的測序序列，得到的有效數據統計如表I所示。權利要求1.基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於其步驟如下1)在獲得RAD高通量測序技術的測序結果後，對RAD雙末端測序序列進行過濾以去除不合格的測序序列；2)根據測序個體基因組酶切一端的測序序列，利用序列的全同性生成每個個體堆的信息；將每個個體過濾後的酶切一端的測序序列信息作為哈希的鍵，哈希的值指向一個鍊表，用於存放另一端的序列信息，並計算測序深度信息；3)過濾掉酶切一端序列測序深度為I的結果；4)兩個個體內分別將酶切一端的測序序列數據進行不容許空隙的兩兩比對，對堆進行聚類以確定個體內在酶切一端序列上的雜合SNP信息；5)在兩個個體內部，對每個堆的另一端數據進行局部組裝；6)利用兩個個體酶切一端的測序序列信息將兩個個體堆的信息相互進行兩兩對齊，即在個體A和個體B中，個體A的某個堆能夠和個體B的某個堆對齊，若且唯若兩個個體堆中的酶切一端的測序序列完全相同；對能夠對齊的堆，兩個個體之間的另一端的組裝結果序列相互進行比對，來尋找Indel位點信息；進而得到兩個個體之間高可信度的Indel位點信息，還得到在Indel位點周圍的側翼序列信息；7)最後在Indel位點周圍的側翼序列上設計引物，應用於後續的大規模實驗。2.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟I)中，RAD高通量測序技術為IlluminaGA測序技術,或現有的其他高通量測序技術。3.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟I)中，所述的不合格的測序序列為測序質量低於預定的低質量閾值的鹼基個數超過整條序列鹼基個數的50%的序列。4.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟3)中，過濾掉酶切一端序列測序深度為I結果的方式為成對過濾。5.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟4)中，只有一個堆的聚類結果表明在酶切一端測序片段上不存在雜合位點；只有兩個堆的聚類結果表明在酶切一端測序片段上存在雜合位點；堆的個數超過兩個的聚類結果被過濾掉。6.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟4)中，計算堆的個數為一和二的聚類結果的總深度，並進一步過濾掉低深度和高深度的聚類結果；其中，低深度的閾值為平均測序深度的四分之一，高深度的閾值為平均測序深度的兩倍。7.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟5)中，組裝選自基於重疊關係的組裝軟體phrap或基於DeBruijngraph算法的組裝軟體SOAPDenovo;利用重疊關係進行組裝的時候不容許空位的存在；利用DeBruijngraph算法進行組裝的時候，kmer的大小要在21到31之間。8.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟7)中，比對的時候，兩個個體之間容許的空位數為一、並且公共區域內不存在錯配。9.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於所使用的比對軟體選自blast或blat序列比對軟體。10.根據權利要求I所述的基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，其特徵在於步驟7)中，應用於構建遺傳圖譜或個體基因分型。全文摘要本發明公開了基於酶切建庫雙末端測序的長度多態性標記的引物設計開發方法，它採用了生物信息學分析方法，處理RAD雙末端測序的測序數據，從而尋找RAD測序片段上的Indel位點信息，以突破非模式生物缺少參考序列的瓶頸，簡化了基因組的複雜度，同時也減少了測序成本。文檔編號G06F19/18GK102831331SQ20121023023公開日2012年12月19日申請日期2012年7月4日優先權日2012年7月4日發明者鄭澤群,任一,陶曄,胡秋萍,黃華生申請人:上海美吉生物醫藥科技有限公司