一種測定待測基因組區(qū)域表達(dá)水平的方法及系統(tǒng)的制作方法
【專利摘要】本發(fā)明提供了一種檢測基因組區(qū)域表達(dá)水平(RPKM)的方法和系統(tǒng)，采用本發(fā)明，一方面，可以檢測出整個(gè)基因的表達(dá)水平及其所有外顯子各自的表達(dá)水平；另一個(gè)方面可以檢測出同一個(gè)基因不同的同源異構(gòu)體的表達(dá)水平及其所有外顯子各自的表達(dá)水平；最后還可以檢測出基因組任意指定區(qū)間的表達(dá)水平。
【專利說明】一種測定待測基因組區(qū)域表達(dá)水平的方法及系統(tǒng)
【技術(shù)領(lǐng)域】
[0001]本發(fā)明涉及生物技術(shù)和生物信息學(xué)領(lǐng)域，具體涉及一種測定基因組區(qū)域表達(dá)水平的方法及系統(tǒng)。
【背景技術(shù)】
[0002]生命遺傳信息的表達(dá)調(diào)控既是生物學(xué)研究的重點(diǎn)領(lǐng)域，也是揭示生物學(xué)各種生命現(xiàn)象的重要手段，尤其是隨著21世紀(jì)大量物種基因組序列的測定以及大量測序技術(shù)推陳出新，使得基因表達(dá)定量方面的研究突飛猛進(jìn)。測序技術(shù)也從傳統(tǒng)Sanger測序技術(shù)，迅速發(fā)展為多種第二代高通量測序技術(shù)，如羅氏454、IlluminaHiSeq和AB公司的SOLiD，以及第三代的單分子實(shí)時(shí)DNA測序技術(shù)。其中，Sanger測序技術(shù)和羅氏454測序技術(shù)的測序讀長在700-1000bp, Illumina測序技術(shù)的測序讀長平均IOObp左右,而單分子實(shí)時(shí)DNA測序技術(shù)的讀長達(dá)到了 2500-3000bp。
[0003]第二代測序技術(shù)也被稱為新一代測序技術(shù)(NGS, Next Generation Sequencing),目前主要是Illumina公司出的HiSeq為主，它通過從物種中提取出的RNA轉(zhuǎn)錄本中隨機(jī)進(jìn)行的短片段測序(通常平均讀長50bp、75bp、100bp)獲得所測樣本的整體表達(dá)譜。轉(zhuǎn)錄本是通過以連續(xù)性基因組為模板進(jìn)行轉(zhuǎn)錄，然后剪切去除內(nèi)含子，拼接剩余的外顯子而形成的。測序過程中，如果一個(gè)轉(zhuǎn)錄本的豐度高，則測序后定位基因組區(qū)域的測序讀段也就多，可以通過對定位到基因上的外顯子區(qū)的測序讀段數(shù)來估計(jì)基因表達(dá)水平。測序讀段數(shù)除了與基因真實(shí)表達(dá)水平成正比，還與基因長度成正比，同時(shí)也與測序深度即測序?qū)嶒?yàn)中得到的總讀段數(shù)正相關(guān)。為了保持對不同基因和不同實(shí)驗(yàn)間估計(jì)的基因表達(dá)值的可比性，Mortazavi等人提出了 RPKM(Reads Per Kilo-base per Million reads)的概念，并成為 RNA-seq 應(yīng)用早期估計(jì)基因表達(dá)水平和外顯子表達(dá)水平的主要方法。RPKM是每百萬讀段中來自于某基因每千堿基長度的讀段數(shù)，考慮了測序深度對讀段計(jì)數(shù)的影響。
[0004]新一代測序技術(shù)的廣泛普及，使得RNA測序(RNA-seq)已成為基因表達(dá)和轉(zhuǎn)錄組分析的重要手段。在NGS測序技術(shù)出現(xiàn)之前，不同基因表達(dá)水平測量的主要手段是基因芯片，利用在基因芯片上高密度集成特點(diǎn)的寡核苷酸，可以對不同組織或者不同發(fā)育階段的特定基因表達(dá)差異和模式進(jìn)行分析。但是與基因芯片數(shù)據(jù)相比，RNA-seq得到的是全基因組轉(zhuǎn)錄水平的數(shù)字化信號，具有高靈敏度、高分辨率、無飽和區(qū)等優(yōu)勢。
[0005]隨著新一代測序技術(shù)的不斷進(jìn)步，產(chǎn)生的RNA-seq數(shù)據(jù)通量高、周期短和成本低，越來越多的人選擇轉(zhuǎn)錄組測序作為科學(xué)研究的首選。RPKM在評估基因表達(dá)水平上的作用越來越顯著，人們通過基因包含的外顯子信息，和轉(zhuǎn)錄組測序數(shù)據(jù)在基因組上的定位信息，來計(jì)算出 RPKM值。FPKM(fragments per kilobase of exon per million fragments mapped)也可以用來表示基因表達(dá)水平。FPKM與RPKM計(jì)算方法基本一致。不同點(diǎn)就是FPKM計(jì)算的是片段(fragments),而RPKM計(jì)算的是測序讀段(reads)。目前cuff links軟件包中的cufflinks模塊和cuffdiff模塊及eXpress軟件可以計(jì)算相關(guān)基因表達(dá)水平,具體計(jì)算過程為，首先統(tǒng)計(jì)出映射定位到基因組上的所有測序讀段數(shù)目，然后統(tǒng)計(jì)出定位到各個(gè)基因外顯子區(qū)間上的所有測序讀段的數(shù)目，再計(jì)算出基因包含的外顯子的長度，最后計(jì)算出基因的FPKM值。
[0006]但是，上述軟件存在以下問題:
[0007](I)目前大部分計(jì)算RPKM的程序,僅支持TopHat、Bowtie、bwa等少數(shù)常用的序列比對定位程序，不能支持所有的Illumina/Solexa測序平臺的讀段定位程序；
[0008](2)在選擇注釋文件的時(shí)候，通常僅支持已知的基因注釋文件，不能支持多種文件格式；
[0009](3)在計(jì)算基因表達(dá)水平的時(shí)候，通常計(jì)算的是片段的表達(dá)水平值，而不是整個(gè)基因的表達(dá)水平值；
[0010](4)在計(jì)算表達(dá)水平的時(shí)候，沒有計(jì)算出單個(gè)外顯子的表達(dá)水平；
[0011](5)在計(jì)算表達(dá)水平的時(shí)候，不能夠計(jì)算出基因組任意指定區(qū)間的表達(dá)水平；
[0012](6)在計(jì)算表達(dá)水平的時(shí)候，通常僅支持計(jì)算一個(gè)轉(zhuǎn)錄組測序結(jié)果，不能夠同時(shí)支持多個(gè)轉(zhuǎn)錄測序結(jié)果的基因表達(dá)水平的計(jì)算。
[0013]因此，本領(lǐng)域期 待一種能夠檢測基因表達(dá)水平和基因組任意指定區(qū)間表達(dá)水平的方法。

【發(fā)明內(nèi)容】

[0014]本發(fā)明的目的是提供一種檢測基因組區(qū)域表達(dá)水平(RPKM)的方法和系統(tǒng)。
[0015]本發(fā)明的第一方面提供了一種測定待測基因組區(qū)域表達(dá)水平的方法，包括以下步驟:
[0016](I)對待測樣本進(jìn)行測序，獲得包含待測基因組區(qū)域轉(zhuǎn)錄本的轉(zhuǎn)錄組測序數(shù)據(jù)；
[0017](2)將獲得的轉(zhuǎn)錄組測序數(shù)據(jù)與同一物種的基因組序列進(jìn)行比對；
[0018](3)對定位到基因組的轉(zhuǎn)錄組測序讀段進(jìn)行篩選，所述篩選包括去除測序質(zhì)量(99.9%的轉(zhuǎn)錄組測序讀段；
[0019](4)將篩選后的轉(zhuǎn)錄組測序讀段，按照其定位到基因組上的起始位置進(jìn)行排序，并對排序結(jié)果建立索引；
[0020](5)根據(jù)待測基因組區(qū)域的位置信息，構(gòu)建出計(jì)算RPKM的基因注釋文件；
[0021](6)計(jì)算能夠映射到基因組上的所有測序讀段的總數(shù)M ;
[0022](7)根據(jù)上述步驟(5)構(gòu)建的基因注釋文件計(jì)算出定位至待測DNA區(qū)間上所有測序讀段的總數(shù)R ；
[0023](8)根據(jù)上述步驟(5)構(gòu)建的基因注釋文件，計(jì)算出待測DNA區(qū)間所有被測序讀段定位的序列長度L ;和
[0024](9)根據(jù)上述步驟(6)-(8)的計(jì)算結(jié)果，將步驟(7)得到的R除以步驟(6)得到的M與步驟(8)得到的L乘以109，得待測基因組區(qū)域的RPKM值，即為待測基因組區(qū)域的表達(dá)水平，計(jì)算公式如下，

D
[0025]RPKli=^IjX IO9




O
[0026]在另一優(yōu)選例中，所述待測基因組區(qū)域包含N個(gè)同源異構(gòu)體，且N≥2。如N可以為 2、3、4、5、6、7、8、9、10 或大于 10。
[0027]在另一優(yōu)選例中，所述方法還包括結(jié)果驗(yàn)證步驟:提取待測樣品的總RNA，經(jīng)過反轉(zhuǎn)錄得到其cDNA，以cDNA作為模板進(jìn)行PCR檢測，驗(yàn)證待測基因組區(qū)域的表達(dá)水平。
[0028]在另一優(yōu)選例中，所述步驟(5)中所述注釋文件整合有已有的基因注釋信息、新預(yù)測的基因注釋信息和/或基因組任意指定區(qū)間的注釋信息。
[0029]在另一優(yōu)選例中，所述待測基因組區(qū)域表達(dá)水平，可以為單個(gè)基因的表達(dá)水平、同一個(gè)基因不同的同源異構(gòu)體的表達(dá)水平、所有外顯子的表達(dá)水平、單個(gè)外顯子的表達(dá)水平以及基因組任意指定區(qū)間的表達(dá)水平。
[0030]在另一優(yōu)選例中，當(dāng)所述待測基因組區(qū)域中包含兩個(gè)以上的同源異構(gòu)體基因序列時(shí)，在測定過程中還包括步驟:將各同源異構(gòu)體的所有外顯子進(jìn)行整合，對于重復(fù)的序列區(qū)間，僅保留單一序列，從而將同一待測基因組區(qū)域中的不同同源異構(gòu)體的外顯子整合成單一序列，將該單一序列的長度作為計(jì)算該基因組區(qū)域表達(dá)水平時(shí)的序列長度L。
[0031]在另一優(yōu)選例中，所述步驟(I)中，所述轉(zhuǎn)錄組序列數(shù)據(jù)由羅氏454測序技術(shù)、Illumina測序技術(shù)、AB公司的SOLiD技術(shù)、或者第三代的單分子實(shí)時(shí)DNA測序技術(shù)獲得。
[0032]在另一優(yōu)選例中，所述步驟(2)中，序列比對程序?yàn)閠ophat2，以程序默認(rèn)參數(shù)進(jìn)行比對。
[0033]在另一優(yōu)選例中，所述步驟(2)中，比對結(jié)果存儲為SAM (Sequence Alignment/Map)格式或其二進(jìn)制版本BAM格式的定位文件。
[0034]在另一優(yōu)選例中，所述步驟⑷中，所述排序方法為:
[0035]a.按照每條測序讀段定位到基因組的起始位置進(jìn)行排序；
[0036]b.如果測序讀段在基因組位置中的起始位置相同，按照其定位到基因組的先后順序進(jìn)行排序，并且保留所有的測序讀段；
[0037]最后對排序結(jié)果建立索引。
[0038]在另一優(yōu)選例中，所述步驟(5)中，基因注釋文件存儲格式為refFlat或者bed格式。
[0039]在另一優(yōu)選例中，所述步驟(7)中，計(jì)算定位至待測DNA區(qū)間上所有測序讀段的總數(shù)R時(shí)，如果一個(gè)轉(zhuǎn)錄組測序讀段定位到兩個(gè)外顯子上，每個(gè)外顯子都會(huì)對這個(gè)測序讀段進(jìn)行統(tǒng)計(jì)，以保證RPKM計(jì)算的準(zhǔn)確性。
[0040]在另一優(yōu)選例中，所述基因組區(qū)域選自如下的組:癌基因基因組區(qū)域、遺傳疾病基因組區(qū)域和/或長非編碼基因區(qū)域或其他任意指定基因組區(qū)域。
[0041]本發(fā)明的第二方面提供了一種檢測基因組區(qū)域表達(dá)水平的系統(tǒng)，所述系統(tǒng)包括:
[0042](I)比對單元，用于轉(zhuǎn)錄組測序讀段與基因組序列進(jìn)行比對；
[0043](2)篩選單元，用于對定位到基因組的轉(zhuǎn)錄組測序讀段進(jìn)行篩選；
[0044](3)排序單元，用于對轉(zhuǎn)錄組測序讀段，按照其定位到基因組上的起始位置進(jìn)行排序;
[0045](4)基因注釋文件構(gòu)建單元，用于構(gòu)建和整合基因注釋文件；和，
[0046](5)計(jì)算單元，包括:
[0047]a.第一模塊，用于計(jì)算能夠映射到基因組上的所有測序讀段的總數(shù)M ;
[0048]b.第二模塊，用于計(jì)算定位至待測DNA區(qū)間上所有測序讀段的總數(shù)R ；[0049]c.第三模塊，用于計(jì)算待測基因組表達(dá)區(qū)域的序列長度之和L ;和，
[0050]d.第四模塊，用于計(jì)算待測基因組區(qū)域的RPKM值，計(jì)算公式為，
[0051]
【權(quán)利要求】
1.一種測定待測基因組區(qū)域表達(dá)水平的方法，其特征在于，包括以下步驟: (1)對待測樣本進(jìn)行測序，獲得包含待測基因組區(qū)域轉(zhuǎn)錄本的轉(zhuǎn)錄組測序數(shù)據(jù)； (2)將獲得的轉(zhuǎn)錄組測序數(shù)據(jù)與同一物種的基因組序列進(jìn)行比對； (3)對定位到基因組的轉(zhuǎn)錄組測序讀段進(jìn)行篩選，所述篩選包括去除測序質(zhì)量(99.9%的轉(zhuǎn)錄組測序讀段； (4)將篩選后的轉(zhuǎn)錄組測序讀段，按照其定位到基因組上的起始位置進(jìn)行排序，并對排序結(jié)果建立索引； (5)根據(jù)待測基因組區(qū)域的位置信息，構(gòu)建出計(jì)算RPKM的基因注釋文件； (6)計(jì)算能夠映射到基因組上的所有測序讀段的總數(shù)M; (7)根據(jù)上述步驟(5)構(gòu)建的基因注釋文件計(jì)算出定位至待測DNA區(qū)間上所有測序讀段的總數(shù)R ； (8)根據(jù)上述步驟(5)構(gòu)建的基因注釋文件，計(jì)算出待測DNA區(qū)間所有被測序讀段定位的序列長度L ;和 (9)根據(jù)上述步驟(6)-⑶的計(jì)算結(jié)果，將步驟(7)得到的R除以步驟(6)得到的M與步驟(8)得到的L乘以109，得待測基因組區(qū)域的RPKM值，即為待測基因組區(qū)域的表達(dá)水平，計(jì)算公式如下，
2.如權(quán)利要求1所述的方法，其特征在于，所述方法還包括結(jié)果驗(yàn)證步驟，優(yōu)選地，所述結(jié)果驗(yàn)證步驟包括:提取待測樣品的總RNA，經(jīng)過反轉(zhuǎn)錄得到其cDNA，以cDNA作為模板進(jìn)行PCR檢測，驗(yàn)證待測基因組區(qū)域的表達(dá)水平。
3.如權(quán)利要求1所述的方法，其特征在于，所述待測基因組區(qū)域包含N個(gè)同源異構(gòu)體，且N≥2。
4.如權(quán)利要求3所述的方法，其特征在于，在測定過程中還包括步驟:將各同源異構(gòu)體的所有外顯子進(jìn)行整合，對于重復(fù)的序列區(qū)間，僅保留單一序列，從而將同一待測基因組區(qū)域中的不同同源異構(gòu)體的外顯子整合成單一序列，將該單一序列的長度作為計(jì)算該基因組區(qū)域表達(dá)水平時(shí)的序列長度L。
5.如權(quán)利要求1所述的方法，其特征在于，所述步驟(1)中，所述轉(zhuǎn)錄組序列數(shù)據(jù)由羅氏454測序技術(shù)、Illumina測序技術(shù)、AB公司的SOLiD技術(shù)、或者第三代的單分子實(shí)時(shí)DNA測序技術(shù)獲得。
6.如權(quán)利要求1所述的方法，其特征在于，所述步驟(4)中，所述排序方法為: a.按照每條測序讀段定位到基因組的起始位置進(jìn)行排序； b.如果測序讀段在基因組位置中的起始位置相同，按照其定位到基因組的先后順序進(jìn)行排序，并且保留所 有的測序讀段； 最后對排序結(jié)果建立索引。
7.如權(quán)利要求1所述的方法，其特征在于，所述基因組區(qū)域選自如下的組:癌基因基因組區(qū)域、遺傳疾病基因組區(qū)域和/或長非編碼基因區(qū)域或任意指定的基因組區(qū)域。
8.一種測定待測基因組區(qū)域表達(dá)水平的系統(tǒng)，其特征在于，所述系統(tǒng)包括:(1)比對單元，用于轉(zhuǎn)錄組測序讀段與基因組序列進(jìn)行比對； (2)篩選單元，用于對定位到基因組的轉(zhuǎn)錄組測序讀段進(jìn)行篩選； (3)排序單元，用于對轉(zhuǎn)錄組測序讀段，按照其定位到基因組上的起始位置進(jìn)行排序； (4)基因注釋文件構(gòu)建單元，用于構(gòu)建和整合基因注釋文件；和， (5)計(jì)算單元,包括: a.第一模塊，用于計(jì)算能夠映射到基因組上的所有測序讀段的總數(shù)M; b.第二模塊，用于計(jì)算定位至待測DNA區(qū)間上所有測序讀段的總數(shù)R； c.第三模塊，用于計(jì)算待測基因組表達(dá)區(qū)域的序列長度之和L;和， d.第四模塊，用于計(jì)算待測基因組區(qū)域的RPKM值，計(jì)算公式為，
9.如權(quán)利要求8所述的系統(tǒng)，其特征在于， 所述篩選單元中,所述篩選包括去除測序質(zhì)量< 99.9%的轉(zhuǎn)錄組測序讀段；和/或， 所述排序單元的排序方法為: a.按照每條測序讀段定位到基因組的起始位置進(jìn)行排序； b.如果測序讀段在基因組位置中的起始位置相同，按照其定位到基因組的先后順序進(jìn)行排序，并且保留所有的測序讀段； 最后對排序結(jié)果建立索引。
10.如權(quán)利要求8所述的系統(tǒng)，其特征在于，所述基因組區(qū)域選自如下的組:癌基因基因組區(qū)域、遺傳疾病基因組區(qū)域和/或長非編碼基因區(qū)域或任意指定的基因組區(qū)域。
【文檔編號】G06F19/22GK103984879SQ201410096063
【公開日】2014年8月13日 申請日期:2014年3月14日 優(yōu)先權(quán)日:2014年3月14日
【發(fā)明者】楊力, 朱閃閃, 薛尉 申請人:中國科學(xué)院上海生命科學(xué)研究院