本發(fā)明涉及生物信息學領域，更具體地涉及一種基于基因組環(huán)境確定遺傳變異功能影響的方法。

背景技術：

隨著以深度測序為代表的高通量遺傳變異檢測技術的快速發(fā)展，目前已可以快速鑒定個體基因組上的遺傳變異。然而，如何準確確定這些遺傳變異對生物分子功能的影響，從而為后續(xù)的個性化醫(yī)療、分子育種等應用提供線索、指導與支持，仍是目前該領域面臨的重大挑戰(zhàn)。

目前在變異注釋領域常用的方法(例如，VEP^[1]、ANNOVAR^[2])通常是以變異為單位，基于參考基因模型(reference gene model)獨立處理每個變異產(chǎn)生的影響。顯然，這種假定每個變異獨立工作產(chǎn)生影響的做法忽略了變異所在的基因組環(huán)境，是不符合生物學實際情況的。大規(guī)模人群基因型數(shù)據(jù)的分析結(jié)果顯示，有大量的變異影響被這種獨立處理每個變異的注釋方法錯誤處理。

另外，在生物學通路水平，目前的注釋方法(例如，DAVID^[3])仍然是以富集分析為主，其利用統(tǒng)計顯著性檢驗方法，根據(jù)用戶提交的大量受變異影響的基因列表中找出在某些通路中顯著富集的結(jié)果。然而，統(tǒng)計檢驗分析旨在找出那些被變異顯著影響的生物學通路，并不能準確直接地指出變異對生物學通路的具體影響。

因此，需要有克服上述缺陷的準確確定遺傳變異對生物分子功能的影響的方法。

技術實現(xiàn)要素：

本發(fā)明針對上述缺陷，提供了一種基于基因組環(huán)境確定遺傳變異影響的方法，所述方法以每個基因為單位注釋該基因上所有變異共同的影響，其包括：1)將所有變異根據(jù)其坐標位置映射到給定基因模型的各個基因上；2)根據(jù)各基因上的所有變異重構(gòu)出各基因的個體化序列；3)對所得個體化序列進行分析以得到變異對該基因的影響。

本發(fā)明還提供了一種基于基因組環(huán)境確定遺傳變異對生物學通路影響的方法，所述方法包括如下步驟：1)將基因/蛋白相互作用通路抽象成有向無環(huán)圖；2)刪除功能缺失基因?qū)膱D節(jié)點以及相應的邊；3)找出因節(jié)點刪除造成的最遠的不連通路徑。

本發(fā)明的方法充分考慮了變異所在的基因組環(huán)境，避免了大量的注釋錯誤，提高了注釋變異影響的準確性。

附圖說明

圖1示出了使用本發(fā)明方法對蛋白編碼基因進行注釋的流程圖。

圖2示出了使用本發(fā)明方法確定變異對基因CHD7影響的結(jié)果。

圖3示出了使用本發(fā)明方法重新分析1000基因組的基因組數(shù)據(jù)的結(jié)果。

圖4示出了使用本發(fā)明方法對轉(zhuǎn)錄因子結(jié)合位點進行注釋的流程圖。

圖5示出了使用本發(fā)明方法確定變異對轉(zhuǎn)錄因子結(jié)合位點TFAP2結(jié)合位點影響的結(jié)果。

圖6示出了利用1000基因組和GTEx項目的基因組數(shù)據(jù)對TFBS注釋沖突進行具體分析的結(jié)果。

圖7示出了使用本發(fā)明方法對microRNA進行注釋的流程圖。

圖8示出了使用本發(fā)明方法所證明的SNP rs56301829與SNP rs2276448變異不會導致microRNA MIMAT0027683失去對基因ZNF716的調(diào)控的示意圖。

圖9示出了利用1000基因組數(shù)據(jù)分析來自TargetScan和miRanda的轉(zhuǎn)錄因子結(jié)合點上的互補突變的結(jié)果。

圖10示出了使用本發(fā)明的方法確定變異對生物學通路的影響的一般性流程圖。

具體實施方式

如上所述，本發(fā)明提供了一種基于基因組環(huán)境確定遺傳變異功能影響的方法，所述方法以每個基因為單位注釋該基因上所有變異共同的影響，其包括：1)將所有變異根據(jù)其坐標位置映射到給定基因模型的各個基因上；2)根據(jù)各基因上的所有變異重構(gòu)出各基因的個體化序列；3)對所得個體化序列進行分析以得到變異對該基因的影響。

在本發(fā)明的方法中，術語“給定基因模型”是指用戶指定的一套完整的基于參考基因組的基因結(jié)構(gòu)描述，包括基因中可變剪切轉(zhuǎn)錄本外顯子、內(nèi)含子在參考基因組上的坐標。

在本發(fā)明的方法中，術語“個體化序列”是指根據(jù)個人基因型得到的個體基因組中的真實序列。

在一個實施方案中，所述基因為蛋白編碼基因。在一個具體的實施方案中，所述蛋白編碼基因為CHD7。在蛋白編碼基因的情況下，本發(fā)明方法的步驟2)通過根據(jù)各個基因上的所有變異推斷該基因的蛋白編碼區(qū)，并將其翻譯成蛋白序列來進行，并且步驟3)通過將所得蛋白序列與已知的參考蛋白序列比較來進行。

在一個實施方案中，所述基因為轉(zhuǎn)錄因子結(jié)合位點。在一個具體的實施方案中，所述轉(zhuǎn)錄因子結(jié)合位點為轉(zhuǎn)錄因子TFAP2結(jié)合位點。

在一個實施方案中，所述基因為microRNA。具體地，在確定變異對microRNA靶基因影響的情況下，本發(fā)明方法包括確定變異對microRNA生成和microRNA靶位點的影響。

本發(fā)明還提供了一種基于基因組環(huán)境確定遺傳變異對生物學通路影響的方法，所述方法包括如下步驟：1)將基因/蛋白相互作用通路抽象成有向無環(huán)圖；2)刪除功能缺失基因?qū)膱D節(jié)點以及相應的邊；3)找出因節(jié)點刪除造成的最遠的不連通路徑。

以下通過具體實施例來說明本發(fā)明的內(nèi)容。應理解，所述具體實施例僅為說明目的，并不意味著本發(fā)明的內(nèi)容僅限于具體實施例。

實施例1：對蛋白編碼基因的注釋

在本實施例中，使用本發(fā)明的方法對蛋白編碼基因進行注釋。

對于蛋白編碼基因的注釋，可通過以下步驟進行：1)將變異根據(jù)其坐標位置映射到給定基因模型的各個基因上；2)根據(jù)各個基因上的所有變異推斷該基因的蛋白編碼區(qū)；具體地，對于含有影響剪切的變異的轉(zhuǎn)錄本，在給定區(qū)間內(nèi)(默認為+/-100bp)的范圍內(nèi)尋找隱藏的剪切位點，以及由其他變異造成的新的剪切位點。對于無法找到可替代的剪切位點的轉(zhuǎn)錄本則外顯子跳過和內(nèi)含子保留都會被考慮。3)根據(jù)所得的蛋白編碼區(qū)序列將其翻譯成蛋白序列；4)將所得蛋白序列與已知的參考蛋白序列比較以確定變異對該基因的影響。圖1示出了使用本發(fā)明方法對蛋白編碼基因進行注釋的流程圖。

以千人基因組項目中HG02861個體基因組中的變異為例進行具體說明。對于SNP rs549508773，VEP會將它注釋為發(fā)生在基因CHD7上的終止密碼子(gag→tag)獲得突變(STOP-gained variant)。然而，本發(fā)明(COPE，基于基因組環(huán)境的變異注釋工具)會根據(jù)HG02861個體所包含的所有變異，發(fā)現(xiàn)該SNP的旁邊還存在另一個SNP rs567756521，兩個SNP的綜合結(jié)果是導致基因CHD7中的一個氨基酸替換(密碼子gag→ttg)(如圖2)。因而，我們的結(jié)果發(fā)現(xiàn)在考慮變異基因組環(huán)境后，相對于現(xiàn)有的VEP來說，可以更準確地注釋變異影響。

為了進一步說明本發(fā)明方法的優(yōu)勢，我們使用本發(fā)明方法重新分析了千人基因組的基因組數(shù)據(jù)(下載于ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/)。我們選擇了其中被VEP注釋為功能缺失型的9728個終止密碼子獲得突變，2092個移碼突變和5559個剪切突變重新用本發(fā)明方法進行分析，結(jié)果如圖3所示。從圖3中可以看出，有大約23.2％影響剪切的變異，6.5％移碼突變以及2.1％終止密碼子獲得的變異是被以前算法錯誤注釋的(圖3a)。富集分析發(fā)現(xiàn)，這些被錯誤注釋為功能缺失的基因與特定的生物功能相關(圖3b)。另外，本發(fā)明還發(fā)現(xiàn)了38對由兩個突變造成新的終止密碼子獲得的變異(被VEP忽略了的功能缺失變異)(圖3c)。由此進一步說明了本發(fā)明方法的優(yōu)勢。

實施例2：對轉(zhuǎn)錄因子結(jié)合位點的注釋

在本實施例中，使用本發(fā)明的方法對轉(zhuǎn)錄因子結(jié)合位點進行注釋。

對轉(zhuǎn)錄因子結(jié)合位點的注釋分為轉(zhuǎn)錄因子結(jié)合位點的丟失(TFBS loss)和轉(zhuǎn)錄因子結(jié)合位點獲得(TFBS gain)。對于轉(zhuǎn)錄因子結(jié)合位點丟失的預測，通過對比變異前后序列對應的位點權(quán)重矩陣得分的大小判斷是否存在轉(zhuǎn)錄因子結(jié)合位點丟失。另一方面，對于轉(zhuǎn)錄因子結(jié)合位點獲得的預測，首先找出所有位于啟動子內(nèi)的變異，進而重構(gòu)啟動子序列。然后，根據(jù)位點權(quán)重矩陣預測轉(zhuǎn)錄因子結(jié)合位點獲得。圖4示出了使用本發(fā)明方法對轉(zhuǎn)錄因子結(jié)合位點進行注釋的流程圖。

以SNP rs9274541和rs9274542為例進行具體說明，這兩個SNP都是位于轉(zhuǎn)錄因子TFAP2結(jié)合位點上的兩個SNP。以往的注釋方法，如SNP2TFBS[⁴]，會分別處理這兩個變異，得到結(jié)果是rs9274541能降低TFAP2的結(jié)合強度(ΔPMW分數(shù)<0)，而rs9274542能增加它的結(jié)合強度(ΔPMW分數(shù)>0)(如圖5a)，兩者沖突。而本發(fā)明的方法綜合考慮兩者的影響，根據(jù)變異重構(gòu)出個體化的轉(zhuǎn)錄因子結(jié)合位點，然后根據(jù)位點權(quán)重矩陣打分，比較參考基因組得分和個體化序列得分，最后判斷為降低了TFAP2的結(jié)合強度(ΔPMW分數(shù)<0)(如圖5b)。

接著，本發(fā)明人利用1000基因組和GTEx項目的基因組數(shù)據(jù)具體分析了上述注釋出現(xiàn)沖突的情況，結(jié)果如圖6所述，從圖中可以看出大約有28％的含有多個變異的TFBSs會出現(xiàn)這類注釋沖突，進一步說明了本發(fā)明方法的優(yōu)勢。

實施例3：對microRNA的注釋

在本實施例中，使用本發(fā)明的方法對microRNA進行注釋。

對miRNA的注釋，同樣根據(jù)給定的變異重構(gòu)出個體化序列綜合判斷變異對miRNA生成和miRNA靶位點的影響。對于microRNA生成的注釋，旨在預測位于pre-microRNA上的基因組變異對pre-microRNA二級結(jié)構(gòu)最小自由能的影響。這里用于計算pre-microRNA最小自由能的工具是RNAfold。首先，找出所有位于同一個pre-microRNA上的變異。然后，重構(gòu)出真實的pre-microRNA序列。最后，計算基因組變異發(fā)生前后，該pre-microRNA的最小自由能的變化，并以此作為對microRNA生成的影響。microRNA靶標結(jié)合注釋是指預測變異對microRNA與3’UTR的相互作用的影響。為了避免預測假陽性過高，使用了兩種常用的工具是TargetScan^[5]和miRanda^[6]。具體地，將靶標結(jié)合缺失定義為在參考序列下TargetScan和miRanda都可以預測到的靶標基因；而在變異發(fā)生后，TargetScan和miRanda都不能預測到的靶標基因。靶標結(jié)合獲得則是指在參考序列下TargetScan和miRanda都不能預測到的靶標基因；而在變異發(fā)生后，TargetScan和miRanda都能預測到的靶標基因。具體地，首先同樣根據(jù)變異重構(gòu)出所有3’UTR和所有microRNA序列。然后，利用miRanda和TargetScan預測新的結(jié)合位點。最后，根據(jù)上面的定義找出所有的靶標丟失和靶標獲得位點。圖7示出了使用本發(fā)明方法對microRNA進行注釋的流程圖。

以SNP rs2276448和rs56301829為例進行具體說明，其分別是在microRNA MIMAT0027683和基因ZNF716上的兩個變異。microRNA對靶位點的識別主要依靠microRNA種子區(qū)(seed region)內(nèi)的幾個堿基與靶位點堿基的互補配對。SNP rs2276448是一個發(fā)生在microRNA種子區(qū)的突變(可以影響microRNA靶位點)，在不考慮其他因素的情況下，基因ZNF716會被錯誤地認為不再受到該microRNA的調(diào)控。本發(fā)明的方法在預測microRNA上的變異對其靶位點的影響時會考慮靶位點上的序列環(huán)境，發(fā)現(xiàn)rs56301829恰好是一個發(fā)生在基因ZNF716上且與SNP rs2276448互補的突變。因此，這對互補的變異不會導致microRNA MIMAT0027683失去對基因ZNF716的調(diào)控(如圖8)。利用1000基因組數(shù)據(jù)分析了來自TargetScan^[5]和miRanda^[6]的轉(zhuǎn)錄因子結(jié)合點上的互補突變，結(jié)果如圖9所示，從圖中可以看出超過半數(shù)的microRNA的靶位點會存在這種互補的變異，進一步說明了本發(fā)明方法的優(yōu)勢。

實施例4：確定變異對生物學通路的影響

在本實施例中，使用本發(fā)明的方法確定變異對生物學通路的影響，具體步驟如下：1)將基因/蛋白相互作用通路抽象成有向無環(huán)圖；2)刪除功能缺失基因?qū)膱D節(jié)點以及相應的邊；3)找出因節(jié)點刪除造成的最遠的不連通路徑。圖10示出了使用本發(fā)明的方法確定變異對生物學通路的影響的一般性流程圖。

以基因GNAS發(fā)生功能缺失為例，本發(fā)明首先找到含有此基因的生物學通路，然后轉(zhuǎn)化為有向無環(huán)圖。最后我們根據(jù)刪除基因GNAS對應的節(jié)點后，找出因該節(jié)點刪除造成的最遠不連通路徑(如對胰腺分泌(Pancreatic secretion)通路造成的最遠不連通路徑是從SCTR(Secretin Receptor，分泌素受體)到ADCY1(Adenylate Cyclase 1，腺苷酸環(huán)化酶1))。本發(fā)明會直接指出變異影響的基因在發(fā)生功能缺失后對生物學通路造成的直接影響，相比于傳統(tǒng)的富集分析來說，這樣的注釋更直觀準確，有助于遺傳診斷分析。

參考文獻：

[1]Deriving the consequences of genomic variants with the Ensembl API and SNP Effect Predictor(McLaren et al.,2010).

[2]ANNOVAR:functional annotation of genetic variants from high-throughput sequencing data(Wang et al.,2010).

[3]DAVID Bioinformatics Resources:expanded annotation database and novel algorithms to better extract biology from large gene lists(Huang,D.W.et al.,2007).

[4]SNP2TFBS-a database of regulatory SNPs affecting predicted transcription factor binding site affinity(Kumar,S.et al.,2017).

[5]Human MicroRNA targets(John,B.et al.,2004).

[6]Conserved seed pairing,often flanked by adenosines,indicates that thousands of human genes are microRNA targets(Lewis,B.P.et al.,2004).