本發(fā)明屬于基因組結(jié)構(gòu)變異檢測領(lǐng)域，具體涉及一種基于三代測序的全基因組結(jié)構(gòu)變異分析方法和系統(tǒng)。
背景技術(shù)：
：基因組結(jié)構(gòu)變異通常是指基因組內(nèi)較大片段的插入、缺失、重復(fù)、倒位、易位以及dna拷貝數(shù)變異(cnv)等。較之于短的序列變異(snp、indel等)，基因組結(jié)構(gòu)變異影響了更多的基因組序列(～13％)，因此也在多種疾病中扮演非常重要的角色。目前，基因組結(jié)構(gòu)變異的檢測主要包括，oligonucleotide-basedarray-cgh、snparray、mlpa、qpcr等一代測序技術(shù)，基于二代測序的breakdancer,readdepth,delly,pindel分析技術(shù)，基于三代測序的pbhoney,sniffles分析技術(shù)。由于一代基于存在價(jià)格高、通量低等弊端，已經(jīng)越來越不適應(yīng)目前的檢測需求；第二代測序技術(shù)的發(fā)展，使得snp、indel等遺傳變異得以廣泛檢測。然而，由于二代測序讀長短(100～150bp左右)的特點(diǎn)，reads不能跨過整個(gè)變異的區(qū)域，盡管使用了多種算法，基因組結(jié)構(gòu)變異的檢測依然存在準(zhǔn)確率低，敏感性低的不足；三代測序技術(shù)具有讀長特別長(最高可達(dá)40k以上)，單堿基錯(cuò)誤率高(15％)，錯(cuò)誤隨機(jī)性好(基本不受gc含量影響)等特點(diǎn)，目前基于第三代的基因組結(jié)構(gòu)變異檢測技術(shù)(pbhoney,sniffles等)雖然大大改善了二代技術(shù)敏感性低的問題，但準(zhǔn)確率低的缺點(diǎn)依然存在。技術(shù)實(shí)現(xiàn)要素：為了解決上述問題，本發(fā)明提供了一種基于三代測序的全基因組結(jié)構(gòu)變異分析方法和系統(tǒng)。所述方法和系統(tǒng)通過整合現(xiàn)有的三代基因組結(jié)構(gòu)變異檢測技術(shù)，能有效提高低覆蓋度下基因組結(jié)構(gòu)變異檢測的準(zhǔn)確性和敏感性，在降低檢測成本的同時(shí)保證檢測結(jié)果的可靠性。本發(fā)明的技術(shù)方案為：一種基于三代測序的全基因組結(jié)構(gòu)變異分析方法，其特征在于,包括以下流程：1)序列拆分，將基因組的測序序列拆分成若干個(gè)用于同步分析的子序列；2)序列比對(duì)，將每個(gè)所述子序列分別通過兩種比對(duì)工具與參考基因組比對(duì)，獲得的比對(duì)結(jié)果分別通過合并工具合并得到兩組比對(duì)序列；3)基因組結(jié)構(gòu)變異初步檢測，將所述兩組比對(duì)序列中每組比對(duì)序列僅通過對(duì)應(yīng)的一種結(jié)構(gòu)變異分析工具進(jìn)行檢測，兩組比對(duì)序列經(jīng)分別檢測后得到兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果；4)基因組結(jié)構(gòu)變異初步檢測結(jié)果合并篩選：4.1)分別將兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果轉(zhuǎn)換成統(tǒng)一格式；4.2)合并兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果：4.2.1)遍歷兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中的缺失序列，如果所述兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中缺失序列重疊部分的長度分別占兩缺失序列長度的比例均大于50％，則判定該兩個(gè)缺失序列為同一個(gè)缺失序列；4.2.2)分別計(jì)算4.2.1)中所述兩個(gè)缺失序列的起始位點(diǎn)和終止位點(diǎn)的均值，所述均值為4.2.1)所述判定的缺失序列的起始位點(diǎn)和終止位點(diǎn)；4.2.3)重復(fù)4.2.1)和4.2.2)中步驟，篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有缺失序列的交集；篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有缺失序列的并集；4.2.4)遍歷兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中的插入序列，判斷如果兩插入序列的距離小于1000bp，則判定該兩個(gè)插入序列為同一個(gè)插入序列；4.2.5)分別計(jì)算4.2.4)中所述兩個(gè)插入序列的起始位點(diǎn)和終止位點(diǎn)的均值，所述均值為4.2.4)所述判定的插入序列的起始位點(diǎn)和終止位點(diǎn)；4.2.6)重復(fù)4.2.4)和4.2.5)中步驟，篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有插入序列的交集；篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有插入序列的并集；4.3)數(shù)據(jù)結(jié)果質(zhì)控：根據(jù)交集和并集中的基因組結(jié)構(gòu)變異檢測結(jié)果的比例以及該區(qū)域的覆蓋度，支持?jǐn)?shù)低于20％的基因組結(jié)構(gòu)變異刪除，得到基因組結(jié)構(gòu)變異最終檢測結(jié)果；5)基因組結(jié)構(gòu)變異功能注釋，利用注釋工具注釋基因組結(jié)構(gòu)變異最終檢測結(jié)果。所述步驟2)中所述兩種比對(duì)工具分別為blasr和bwa；所述步驟2)中合并工具為samtools。所述步驟3)中通過blasr比對(duì)得到的比對(duì)序列對(duì)應(yīng)的結(jié)構(gòu)變異分析工具為pbhoney；所述步驟3)中通過bwa比對(duì)得到的比對(duì)序列對(duì)應(yīng)的結(jié)構(gòu)變異分析工具為sniffles。所述步驟4.1)中的統(tǒng)一格式為bed格式。所述步驟5)中的注釋工具為annovar。一種基于三代測序的全基因組結(jié)構(gòu)變異分析系統(tǒng)，其特征在于，所述基于三代測序的全基因組結(jié)構(gòu)變異分析系統(tǒng)包括以下模塊：序列拆分模塊，用于將基因組的測序序列拆分成若干個(gè)用于同步分析的子序列；序列比對(duì)模塊，包括兩個(gè)并列的比對(duì)單元，所述比對(duì)單元用于所述子序列與參考基因組的比對(duì)，獲得兩組比對(duì)序列；基因組結(jié)構(gòu)變異初步檢測模塊，包括兩個(gè)并列的結(jié)構(gòu)變異分析單元，所述兩個(gè)結(jié)構(gòu)變異分析單元用于同步檢測兩組比對(duì)序列中的基因組結(jié)構(gòu)變異，獲得兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果；基因組結(jié)構(gòu)變異初步檢測結(jié)果合并篩選模塊，包括格式轉(zhuǎn)換單元、數(shù)據(jù)分析單元、交集單元、并集單元和數(shù)據(jù)結(jié)果質(zhì)控單元；所述格式轉(zhuǎn)換單元用于將兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果轉(zhuǎn)換成統(tǒng)一格式；所述數(shù)據(jù)分析單元用于分析基因組結(jié)構(gòu)變異初步檢測結(jié)果，具體為遍歷兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中的缺失序列，如果所述兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中缺失序列重疊部分的長度分別占兩缺失序列長度的比例均大于50％，則判定該兩個(gè)缺失序列為同一個(gè)缺失序列；分別計(jì)算所述兩個(gè)缺失序列的起始位點(diǎn)和終止位點(diǎn)的均值，所述均值為所述判定的缺失序列的起始位點(diǎn)和終止位點(diǎn)；篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有缺失序列的交集結(jié)果，將所述交集結(jié)果置于交集單元中；篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有缺失序列的并集結(jié)果，將所述交集結(jié)果置于并集單元中；遍歷兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中的插入序列，判斷如果兩插入序列的距離小于1000bp，則判定該兩個(gè)插入序列為同一個(gè)插入序列；分別計(jì)算所述兩個(gè)插入序列的起始位點(diǎn)和終止位點(diǎn)的均值，所述均值為所述判定的插入序列的起始位點(diǎn)和終止位點(diǎn)；篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有插入序列的交集結(jié)果，并將交集結(jié)果置于交集單元中；篩選出兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中所有插入序列的并集結(jié)果，并將并集結(jié)果置于并集單元中；所述數(shù)據(jù)結(jié)果質(zhì)控單元根據(jù)交集單元和并集單元中基因組結(jié)構(gòu)變異檢測結(jié)果的比例以及該區(qū)域的覆蓋度，支持?jǐn)?shù)低于20％的基因組結(jié)構(gòu)變異刪除，得到基因組結(jié)構(gòu)變異最終檢測結(jié)果；基因組結(jié)構(gòu)變異功能注釋模塊，包括注釋單元，所述注釋單元用于注釋基因組結(jié)構(gòu)變異最終檢測結(jié)果。所述兩個(gè)比對(duì)單元運(yùn)用的分析工具分別為blasr和bwa，分析后的數(shù)據(jù)均用合并工具samtools合并。所述兩個(gè)結(jié)構(gòu)變異分析單元運(yùn)用的工具分別為pbhoney和sniffles；blasr的運(yùn)用與pbhoney相對(duì)應(yīng)；bwa的應(yīng)用與sniffles相對(duì)應(yīng)。所述格式轉(zhuǎn)換單元轉(zhuǎn)換后的統(tǒng)一格式為bed格式。所述基因組結(jié)構(gòu)變異功能注釋模塊中的注釋工具為annovar。本發(fā)明的有益效果為：基因組的一代測序和二代測序耗時(shí)長，三代測序雖然速度得到大幅提升，但是準(zhǔn)確度低，要得到更精準(zhǔn)的數(shù)據(jù)需要很高的覆蓋深度，成本大大提高。本發(fā)明根據(jù)兩種三代測序工具測序后所得結(jié)果進(jìn)行并集或交集來輸出最終的結(jié)構(gòu)變異分析結(jié)果，來滿足準(zhǔn)確度或靈敏度要求，特別是實(shí)現(xiàn)了低覆蓋深度下，基因組結(jié)構(gòu)變異檢測結(jié)果的可靠性，提升檢測速度的同時(shí)降低了檢測成本。附圖說明圖1為本發(fā)明實(shí)施例1和實(shí)施例2的流程圖。圖2為本發(fā)明的所述系統(tǒng)的結(jié)構(gòu)示意圖。圖3為本發(fā)明所述系統(tǒng)中基因組結(jié)構(gòu)變異初步檢測結(jié)果合并篩選模塊的結(jié)構(gòu)示意圖。圖4為圖2不同軟件在實(shí)施例1樣品中缺失序列檢測準(zhǔn)確率/檢出率比較。圖5為不同軟件在實(shí)施例1樣品中插入序列檢出率比較。圖6為不同軟件在實(shí)施例2樣品中缺失序列檢測準(zhǔn)確率/檢出率比較。圖7不同軟件在實(shí)施例2樣品中插入序列檢測準(zhǔn)確率/檢出率比較。具體實(shí)施方式結(jié)合附圖和具體實(shí)施例，對(duì)本發(fā)明作進(jìn)一步描述。結(jié)合附圖1對(duì)本發(fā)明實(shí)施例所述基于三代測序的全基因組結(jié)構(gòu)變異分析方法的工作流程進(jìn)行說明，詳細(xì)流程如下所示：步驟1，獲得原始bam文件數(shù)據(jù)；步驟2，將bam文件中的序列拆分，將基因組的測序序列拆分成若干個(gè)用于同步分析的子序列，即將原始reads數(shù)拆分成多個(gè)fastq文件；每個(gè)fastq文件進(jìn)入步驟3和步驟4；步驟3和步驟4同步進(jìn)行，將fastq文件中的數(shù)據(jù)進(jìn)行基因比對(duì)，步驟3中fastq文件用blasr比對(duì)，比對(duì)結(jié)果文件用samtools合并；步驟4中fastq文件用bwa比對(duì)，比對(duì)結(jié)果文件用samtools合并；步驟3合并后的數(shù)據(jù)進(jìn)入步驟5用pbhoney做基因組結(jié)構(gòu)變異檢測；步驟4合并后的數(shù)據(jù)進(jìn)入步驟6用sniffles做基因組結(jié)構(gòu)變異檢測；步驟5獲得的基因組結(jié)構(gòu)變異初步檢測結(jié)果進(jìn)入步驟7轉(zhuǎn)化成bed格式；步驟6獲得的基因組結(jié)構(gòu)變異初步檢測結(jié)果進(jìn)入步驟8轉(zhuǎn)化成bed格式；步驟9遍歷兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中的缺失序列，如果所述兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中缺失序列重疊部分的長度分別占兩缺失序列長度的比例均大于50％，則判定該兩個(gè)缺失序列為同一個(gè)缺失序列，進(jìn)入步驟10；步驟9判定該兩個(gè)缺失序列不是同一個(gè)缺失序列時(shí)，進(jìn)入步驟12；步驟10計(jì)算被判定為同一個(gè)缺失序列的兩缺失序列的起始位點(diǎn)和終止位點(diǎn)的均值，所述均值為所述判定的缺失序列的起始位點(diǎn)和終止位點(diǎn)；進(jìn)入步驟11；步驟11將步驟10篩選出的所有缺失序列合并作為缺失序列的交集結(jié)果進(jìn)入步驟12；步驟12將步驟11中缺失序列的交集結(jié)果和步驟9中判定為不是同一缺失序列的缺失序列合并，作為所有缺失序列的并集結(jié)果進(jìn)入步驟13；步驟9遍歷兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中的插入序列，如果所述兩組基因組結(jié)構(gòu)變異初步檢測結(jié)果中兩插入序列的距離小于1000bp，則判定該兩個(gè)插入序列為同一個(gè)插入序列，進(jìn)入步驟10；步驟9判定該兩個(gè)缺失序列不是同一個(gè)插入序列時(shí)，進(jìn)入步驟12；步驟10計(jì)算被判定為同一個(gè)插入序列的兩插入序列的起始位點(diǎn)和終止位點(diǎn)的均值，所述均值為所述判定的插入序列的起始位點(diǎn)和終止位點(diǎn)；進(jìn)入步驟11；步驟11將步驟10篩選出的所有插入序列合并作為交集結(jié)果進(jìn)入步驟12；步驟12將步驟11中交集結(jié)果和步驟9中判定為不是同一插入序列的插入序列合并，作為所有插入序列并集結(jié)果，進(jìn)入步驟13步驟13將步驟11和步驟12得到的基因組結(jié)構(gòu)變異中支持?jǐn)?shù)低于20％的基因組結(jié)構(gòu)變異刪除，得到基因組結(jié)構(gòu)變異最終檢測結(jié)果；進(jìn)入步驟14；步驟14利用注釋工具注釋出基因組結(jié)構(gòu)變異最終檢測結(jié)果中基因組結(jié)構(gòu)變異的不同功能類型以及其他相關(guān)信息，得到最終結(jié)果。從圖2可以看出，本發(fā)明實(shí)施例所述基于三代測序的全基因組結(jié)構(gòu)變異分析系統(tǒng)包括序列拆分模塊10，序列對(duì)比模塊20，基因組結(jié)構(gòu)變異初步檢測模塊30，基因組結(jié)構(gòu)變異初步檢測結(jié)果合并篩選模塊40，基因組結(jié)構(gòu)變異功能注釋模塊50。從圖3可以看出，基因組結(jié)構(gòu)變異初步檢測結(jié)果合并篩選模塊包括格式轉(zhuǎn)化單元41，格式轉(zhuǎn)化單元42，數(shù)據(jù)分析單元43，交集單元44，并集單元45和數(shù)據(jù)結(jié)果質(zhì)控單元46。實(shí)施例中各模塊和單元中采用了適用于三代測序超長讀長的多種生物信息分析軟件，具體如下：1、blasr比對(duì)是一個(gè)非常耗費(fèi)計(jì)算資源和時(shí)間的過程，所以本系統(tǒng)首先將原始的測序數(shù)據(jù)按照原始reads數(shù)拆分成多個(gè)fastq文件，在比對(duì)過程中采用多個(gè)任務(wù)并行的模式，大量的節(jié)省了時(shí)間。2、基因組結(jié)構(gòu)變異pbhoney檢測2.1)fastq文件分別用blasr比對(duì)。2.2)比對(duì)結(jié)果文件用samtools合并，用pbhoney做基因組結(jié)構(gòu)變異檢測。3、基因組結(jié)構(gòu)變異sniffles檢測3.1)fastq文件分別用bwa比對(duì)。3.2)比對(duì)結(jié)果文件用samtools合并，用sniffles做基因組結(jié)構(gòu)變異檢測。4、原始基因組結(jié)構(gòu)變異初步檢測結(jié)果合并篩選4.1)分別將pbhoney、sniffles結(jié)果轉(zhuǎn)換成統(tǒng)一的bed格式，方便后續(xù)的合并與篩選。4.2)合并pbhoney、sniffles結(jié)果。4.2.1)遍歷pbhoney、sniffles結(jié)果中的缺失序列，判斷如果兩個(gè)缺失序列重疊部分的長度占兩缺失序列長度的比例大于50％，則判定該兩個(gè)缺失序列為同一個(gè)缺失序列。4.2.2)分別計(jì)算pbhoney、sniffles缺失序列起始位點(diǎn)和終止位點(diǎn)的均值作為合并后結(jié)果的起始位點(diǎn)和終止位點(diǎn)。4.2.3)將pbhoney、sniffles結(jié)果中的intersection部分輸出到intersection結(jié)果中；將intersection和其他結(jié)果輸出到union結(jié)果中。4.2.4)遍歷pbhoney、sniffles結(jié)果中的插入序列，判斷如果兩插入序列的距離小于1000bp，則認(rèn)為該兩個(gè)插入序列為同一個(gè)插入序列，否則則認(rèn)為兩插入序列不同。4.2.5)分別計(jì)算pbhoney、sniffles缺失序列起始位點(diǎn)和終止位點(diǎn)的均值作為合并后結(jié)果的起始位點(diǎn)和終止位點(diǎn)。4.2.6)將pbhoney、sniffles結(jié)果中的intersection部分輸出到intersection結(jié)果中；將intersection和其他結(jié)果輸出到union結(jié)果中。4.3)數(shù)據(jù)結(jié)果質(zhì)控根據(jù)支持基因組結(jié)構(gòu)變異reads的比例以及該區(qū)域的覆蓋度，支持?jǐn)?shù)低于20％的基因組結(jié)構(gòu)變異刪除。5、基因組結(jié)構(gòu)變異功能注釋本系統(tǒng)利用annovar注釋出基因組結(jié)構(gòu)變異的不同功能類型以及其他相關(guān)信息，方便用戶的進(jìn)一步篩選。本系統(tǒng)結(jié)果分為union(并集)和intersection(交集)兩種模式，union模式敏感性方面非常好，而intersection模式則在準(zhǔn)確性方面具有極大的優(yōu)勢。在10x覆蓋度的情況下，本發(fā)明的union模式對(duì)indel的檢出率達(dá)到75％以上，intersection模式的準(zhǔn)確率接近90％，用戶可以根據(jù)自己的需求選擇適合自己的模式。以下通過具體實(shí)施例對(duì)本發(fā)明的結(jié)果與技術(shù)參數(shù)做詳細(xì)說明。實(shí)施例1.樣品：該樣品來自本公司一個(gè)自愿捐獻(xiàn)者,該樣品有很好的一代和二代測序的研究基礎(chǔ)，所以本實(shí)施例將該樣品作為一個(gè)democase來說明本系統(tǒng)的準(zhǔn)確性。數(shù)據(jù)分析與結(jié)果統(tǒng)計(jì)：原始數(shù)據(jù)統(tǒng)計(jì)表1原始數(shù)據(jù)統(tǒng)計(jì)測序base數(shù)34.28gpolymerread數(shù)3.59mpolymerread平均長度9,441polymerread長度n5016,694subread數(shù)12.88msubread平均長度2,624subread平均n503,208比對(duì)結(jié)果統(tǒng)計(jì)通過blasr比對(duì)，最終有12.85mreads被比對(duì)到基因組(版本號(hào)hg19)上。與標(biāo)準(zhǔn)數(shù)據(jù)比較目前已知本實(shí)施例所用樣品中長度大于200bp的缺失序列和插入序列共有2194和68個(gè)。標(biāo)準(zhǔn)結(jié)果中插入序列數(shù)量較少，該情況應(yīng)該是因?yàn)橐淮约岸鷾y序技術(shù)對(duì)插入序列檢測結(jié)果太差造成的。表2實(shí)施例1與其他軟件對(duì)缺失序列檢測結(jié)果比較表3實(shí)施例1與其他軟件對(duì)插入序列檢測結(jié)果比較實(shí)施例2.樣品：該樣品是本公司利用三代測序技術(shù)完成個(gè)一個(gè)全基因組測序樣品。該樣品的測序深度高達(dá)100x,所以該樣品基因組結(jié)構(gòu)變異的檢測結(jié)果具有較高的可信度。本實(shí)施例將多種系統(tǒng)在高深度條件下檢測出的基因組結(jié)構(gòu)變異作為標(biāo)準(zhǔn)集，并隨機(jī)挑選了10x數(shù)據(jù)做為測試數(shù)據(jù)測試本發(fā)明的準(zhǔn)確性。數(shù)據(jù)分析與結(jié)果統(tǒng)計(jì)：本實(shí)施例測試數(shù)據(jù)統(tǒng)計(jì)結(jié)果如下表表4原始數(shù)據(jù)統(tǒng)計(jì)測序base數(shù)34.22gpolymerread數(shù)2.39mpolymerread平均長度14,344polymerread長度n5012,169subread數(shù)3.03msubread平均長度11,294subread平均n509,954比對(duì)結(jié)果統(tǒng)計(jì)通過blasr比對(duì)，最終有3.03mreads被比對(duì)到基因組(版本號(hào)hg19)上。與標(biāo)準(zhǔn)數(shù)據(jù)比較經(jīng)過檢測，該樣品中共發(fā)現(xiàn)缺失序列和插入序列分別為2978和2950個(gè)，根據(jù)比較結(jié)果intersection的準(zhǔn)確率可以高達(dá)90％。表5實(shí)施例2與其他軟件對(duì)缺失序列檢測結(jié)果比較表6實(shí)施例2與其他軟件對(duì)插入序列檢測結(jié)果比較通過兩個(gè)標(biāo)準(zhǔn)樣品的驗(yàn)證，本發(fā)明在測序深度約為10x的情況下，缺失/插入的準(zhǔn)確率和檢出率分別達(dá)到90％和75％以上，將三代基因組結(jié)構(gòu)變異檢測準(zhǔn)確度提高了1倍。根據(jù)實(shí)施例1和實(shí)施例2我們可以得出本發(fā)明union部分敏感性可以達(dá)到75％以上，intersection部分準(zhǔn)確性可以達(dá)到90％。以上所述僅為本發(fā)明的較佳實(shí)施例，凡在本發(fā)明的原則之內(nèi)所做的任何簡單修改、等同變換與改型，仍應(yīng)屬于本發(fā)明的保護(hù)范圍之內(nèi)。當(dāng)前第1頁12