本發明涉及生物，具體涉及一種應用于二代測序數據的str分型方法。

背景技術：

1、str(short?tandem?repeats)又名微衛星(microsatellites)序列，是由一系列2-6個核苷酸的dna重復單元構成的長度可達100個核苷酸的序列。研究表明，與只有兩個等位基因型的snv(single?nucleotide?variations)不同，str在人類基因組中具有高度多態性的特點，因此被廣泛應用于醫學遺傳學、群體遺傳學和法醫科學領域。在醫學領域中，str被證明在超過65種孟德爾遺傳病中起著關鍵作用，并且參與基因調控以及復雜性狀的形成。在法醫科學中，str則被廣泛應用于個體身份識別、親子鑒定以及姓氏推斷等工作中。

2、str等位基因分型主要基于毛細管電泳(capillary?electrophoresis,ce)和熒光標記多重擴增的str分型方法，該類方法操作簡單、經濟實惠、檢測快速，因此，在數十年來都依舊被作為法醫dna實驗室的“金標準”。然而，基于毛細管電泳的str分型技術受光譜分辨率的限制，導致單次只能分型數十個str位點，從而無法應對大規模樣本對str分型高吞吐量的需求。

3、隨著高通量測序技術的出現和發展，科研人員可以通過多重擴增子技術(multiplex?amplicon?sequencing,mas)同時對成百上千的str基因座進行擴增，很大程度上可以消除傳統ce技術對str基因座分型數目的限制。此外，基于高通量測序技術的str分型技術具有更高的靈敏度和準確性，研究表明62pg的基因組dna便可以成功分型出97.7％的等位基因。然而，由于測序深度、側翼序列重復性、pcr影峰等多重因素影響，從低覆蓋度的測序數據中獲得準確的str分型結果具有較大難度。目前已經存在許多基于高通量測序數據的str分型方法，如strait?razor、hipstr、strling、expansionhunter、strinngs等，但是這些方法在進行str分型時存在兼容性差、性能低下、準確率低等問題。strait?razor和hitstr無法很好處理側翼序列具有較高重復性的str基因座；strling和expansionhunter專為長片段str位點設計，無法很好工作在單端測序數據中；strinngs通過正則表達式匹配str重復單元，導致其分型速度緩慢且分型準確率低。

4、如中國專利cn113362892a公開了一種短串聯重復序列重復數的檢測和分型方法，該方法對側翼序列的保守性依賴較高，如果側翼序列存在變異，將會對比對產生很大的影響，造成str的邊界識別問題，進一步影響str的分型。中國專利cn117037906a公開了一種基于二代測序的短串聯重復序列的分型方法，該方法極大依賴參考基因組，如果檢測樣本與參考基因組差異較大，reads將無法難以比對到參考基因組上，即使比對上，其邊界也無法準備匹配。

5、基于上述問題，研發一種應用于二代測序數據的str分型方法，提高其分型效率和準確率，為海量組學數據的準確分型提供基礎支撐，是本領域研究人員的研究重點。

技術實現思路

1、本發明針對上述問題，提供了一種應用于二代測序數據的str分型方法，能夠從測序數據中分型獲得指定str基因座的等位基因及其似然概率。

2、為實現上述目的，本發明采用的技術方案如下：

3、一方面，本發明提供了一種應用于二代測序數據的str分型方法，包含以下步驟：

4、s1：定制str基因座配置文件；

5、s2：構建str基因座側翼序列k-mer索引：為配置文件中每個str基因座的側翼序列構建k-mer索引庫；

6、s3：針對步驟s2所述str基因座配置文件中的每個str基因座，分別執行步驟s4-s7；

7、s4：提取str基因座等位基因序列：通過結合cigar值和k-mer索引兩種方法從序列比對文件的每條read中提取str等位基因序列，生成str基因座的等位基因集合；

8、s5：構建str基因座的pcr影峰模型，使用二代測序樣本數據集和最大似然估計法估計pcr影峰模型的參數；

9、s6：估計str基因座的候選等位基因頻率：使用二代測序樣本數據集估計候選等位基因頻率；

10、s7：推斷序列比對文件中str基因座最有可能的等位基因及其后驗概率：利用s4提取的等位基因集合計算所有候選等位基因的后驗概率，以此推斷最有可能的等位基因。

11、優選地，s1中，所述str基因座配置文件用于定義str位點的主要信息和在參考基因組上的精確位置。

12、優選地，s1中，所述str基因座配置文件包含如下字段：str基因座名稱、染色體、起始位置、結束位置、重復單元長度、排除堿基數目和是否為單倍型。

13、優選地，所述重復單元長度用于計算等位基因大小。

14、優選地，所述排除堿基數目可以排除不需要參與計算等位基因的序列。

15、優選地，所述是否為單倍型是用于分型性染色體中由兩個片段組成的str基因座。

16、優選地，s2中，所述側翼序列為str核心序列5’和3’端兩側相連的堿基序列。

17、優選地，s2中，所述k-mer索引為局部唯一k-mer子序列，可以提高str基因座側翼序列匹配準確率。

18、優選地，s2中，所述構建str基因座側翼序列k-mer索引，采用動態k-mer索引構建策略，即根據str基因座側翼序列的重復性選擇最佳長度k作為該基因座k-mer子序列長度，然后將k-mer子序列長度k、k-mer在側翼序列中的相對位置p存儲到哈希表中，為后續等位基因的提取做準備。

19、優選地，s4中，所述比對文件為比對到參考基因組后的bam文件。

20、優選地，s4中，所述提取str等位基因序列的具體操作為：

21、s4-1、以序列比對文件為輸入，首先利用每條read的cigar值和相應str位點在參考基因組上的位置計算str側翼序列位置；

22、s4-2、隨后利用配置文件中的側翼序列對計算得到的read側翼序列進行驗證，以此直接計算得到str等位基因；

23、s4-3、當cigar值中插入或缺失出現在側翼序列與核心序列交界處且堿基數目不超過2倍str重復單元長度時，則利用cigar值校正方法將5’和3’端側翼序列中與重復單元相同且比對錯誤的插入或缺失移回str核心區域，并最終使用校正后的cigar值計提取str等位基因序列；

24、s4-4、當樣本在str基因組上的等位基因長度超過參考基因組在該位點的等位基因長度的40％時，序列比對軟件往往會引入完全錯誤的cigar值，為了彌補cigar方法的缺陷，啟用k-mer索引方法進行等位基因的提取：針對每條read序列，從s2構建的k-mer索引庫中提取對應str位點信息，然后匹配側翼序列并計算得到str等位基因。

25、優選地，s5中，所述pcr影峰模型包含參數u和d，u代表與真實等位基因相比多一個重復單元的概率，而d則表示與真實等位基因相比少一個重復單元的概率。

26、優選地，s5中，所述pcr影峰模型的分布如下所示：

27、對于給定的str基因座，x表示觀察到的等位基因與真實等位基因之間相差的重復單元的個數，共計分為三種情況：

28、與真實等位基因相同(x＝0)；

29、與真實等位基因相比多一個重復單元(x＝1)；

30、以及與真實等位基因相比少一個重復單元(x＝-1)；

31、因此，可以寫出pcr影峰模型的分布：

32、

33、優選地，s5中，所述使用二代測序樣本數據集和最大似然估計法估計pcr影峰模型的參數的具體操作為：

34、s5-1、對于一個給定的個體i在當前str基因座上具有ni條完全橫跨該基因座的reads，則令pcr影峰集合為其中表示個體i第j條橫跨當前基因座的read中觀察到的等位基因與真實等位基因的之差；

35、然后，對于二代測序樣本數據集中的n個個體的似然函數如下所示：

36、

37、s5-2、基于公式(1)和(2)，可以計算得到關于u和d的似然函數：

38、

39、其中：

40、t表示n個樣本完全橫跨當前str基因座的reads數目之和；

41、m1表示reads中的等位基因與真實等位基因相比多一個重復單元的數目；

42、m2則表示reads中的等位基因與真實等位基因相比少一個重復單元的數目；

43、s5-3、通過公式(3)和(4)，即可估計獲得參數u和d的極大似然估計值：

44、

45、優選地，s5中，對于用戶自定義的pcr影峰模型參數的估計，本發明需要通過嚴格的標準篩選輸入樣本，用以估計獲得可靠的pcr影峰模型參數，其基本規則包括：

46、(1)完全橫跨str基因座的reads數目必須超過設定的最小限制(默認為10)；

47、(2)對于reads支持數最多的等位基因，其reads支持數必須超過該str基因座上總reads數目的70％。

48、優選地，s6中，所述估計候選等位基因頻率，具體為：通過用戶輸入的二代測序樣本數據集，統計當前str基因座等位基因頻率分布。并能在樣本量增加時，重新統計當前基因座的候選等位基因頻率。

49、優選地，s6中，所述估計候選等位基因頻率的標準為：

50、(1)用于估計等位基因頻率分布的測序樣本數不少于20個；

51、(2)候選等位基因必須在至少2個以上樣本中出現。

52、優選地，s6中，對于測序樣本不足以估計自定義str基因座等位基因頻率的數據集，將為每一個出現的潛在str等位基因賦予相同頻率值，以此消除先驗概率對后驗概率估計的影響。

53、優選地，所述s7，s7中，推斷序列比對文件中當前str基因座最有可能的等位基因及其后驗概率，具體操作為：

54、s7-1、如果當前str基因座在配置文件被定義為單倍型。假設其在當前str基因座上有n條完全橫跨該基因座的reads，則這些reads中觀察到的等位基因構成的等位基因集合可以表示為a＝{a1,a2,…,an}，然后，對于每一個可能的候選等位基因a(a∈a)的后驗概率可以通過如下公式進行計算：

55、

56、其中，fa表示候選等位基因a在當前str基因座中的頻率；

57、ai表示該樣本在當前基因座上的第i個觀察到的等位基因；

58、a+1和a-1分別表示與給定候選等位基因a相比增加一個重復單元和減少一個重復單元的影峰等位基因；

59、others則表示觀察到的等位基因ai比給定的候選等位基因a多若干個或少若干個重復單元；

60、s7-2、如果當前str基因座在配置文件被定義為非單倍型。為了估計兩個候選等位基因組合的后驗概率，假定str基因座的兩個等位基因分別為j(j∈a)和k(k∈a)；于是就可以估計出所有兩兩組合的等位基因的聯合后驗概率：

61、p(allele＝(j,k)|a,u,d)∝#

62、

63、其中，fj,fk表示候選等位基因j和等位基因k在當前str基因座中的頻率；

64、ai表示該樣本在當前基因座上的第i個觀察到的等位基因；

65、a+1和a-1分別表示與給定候選等位基因a相比增加一個重復單元和減少一個重復單元的影峰等位基因；

66、s7-3、通過計算，得到當前str基因座的所有候選等位基因或等位基因組合的后驗概率，選取后驗概率最大的等位基因或等位基因組合，作為序列比對文件在當前str基因座上最有可能的等位基因或等位基因組合。

67、相對于現有技術，本發明具有以下有益效果：

68、本發明提出了一種可應用于二代測序數據的str分型技術方案，能夠從測序數據中分型獲得指定str基因座的等位基因及其似然概率，其檢出率和分型準確率更高；除了可應用于全基因組測序(whole?genome?sequencing,wgs)測序數據外，同樣也可以應用于高覆蓋度靶向測序數據或擴增子測序數據；通過算法優化極大提高了str分型速度。