本發(fā)明涉及基因組學及生物信息學技術領域，具體涉及一種基因表達的定量方法及裝置。

背景技術：
轉(zhuǎn)錄組測序技術（RNA-seq，RNAsequencing）是把小RNA（RibonucleicAcid，核糖核酸）、mRNA和非編碼RNA等或者其中一些用高通量測序技術把它們的序列測出來。目前RNA-seq測序平臺有多種，包括Hiseq、RocheFLX、IlluminaSolexa、ABIsolid等。不同測序平臺的測序原理有所不同，但測序步驟基本包括文庫制備，聚合酶鏈式反應（PCR，PolymeraseChainReaction）擴增等。通過RNA-seq，科研工作者能夠獲得生物中基因表達的情況，研究不同個體、不同時期、不同形態(tài)的組織的基因表達水平的差異。中國專利申請（申請?zhí)枺?01110283718.2，名稱：一種分析基因表達定量的方法）基于Illumina平臺公開一種分析基因表達定量的方法，可以克服數(shù)字基因表達譜（DGE，DigitalGeneExpression）技術對CATG位點和參考基因完整性依賴性強的缺點。但是，該方法測序分析需時較長，勞動效率有待提高。

技術實現(xiàn)要素：
本發(fā)明提供一種基因表達的定量方法及裝置，可以快速地完成基因表達的定量。依據(jù)本發(fā)明的一方面提供一種基因表達的定量方法,包括：獲取含有核酸序列信息的讀段序列；將讀段序列與所有參考基因進行比對，獲取比對上的讀段序列；對比對上的讀段序列進行過濾，舍去軟剪切比例超過第一預設值，序列長度小于第二預設值，以及比對得分小于第三預設值的讀段序列，軟剪切比例是指沒有比對上的堿基數(shù)目占該讀段序列總堿基數(shù)目的比例；比對得分是按照每個讀段序列與參考基因的匹配程度以及讀段序列的長度而確定的數(shù)值；對于已過濾的讀段序列，使用每百萬讀段序列中來自目標基因每千堿基長度的讀段序列數(shù)目RPKM對所述目標基因表達進行定量，定義為RPKM=（比對到目標基因?qū)膮⒖蓟虻淖x段序列的數(shù)目）*109/（比對到所有參考基因的讀段序列的數(shù)目*目標基因的長度）。優(yōu)選地，比對到目標基因?qū)膮⒖蓟虻淖x段序列的數(shù)目是指只能比對到目標基因?qū)膮⒖蓟蛏?，而且能夠比對到所述參考基因的至少一個轉(zhuǎn)錄本的讀段序列的數(shù)目；目標基因的長度是指目標基因的所有轉(zhuǎn)錄本中最長的轉(zhuǎn)錄本的長度。依據(jù)本發(fā)明的另一方面提供一種基因表達的定量裝置，包括：數(shù)據(jù)輸入單元，用于輸入數(shù)據(jù)；數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)；存儲單元，用于存儲數(shù)據(jù)，其中包括可執(zhí)行的程序；處理器，與數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元及存儲單元數(shù)據(jù)連接，用于執(zhí)行存儲單元中存儲的可執(zhí)行的程序，該程序的執(zhí)行包括完成上述基因表達的定量方法。本發(fā)明的有益效果是：通過將讀段序列與參考基因進行比對，而不是現(xiàn)有的與參考基因組進行比對，可以簡化比對過程，提高比對效率。特別地，比對到目標基因?qū)膮⒖蓟虻淖x段序列的數(shù)目是指只能比對到目標基因?qū)膮⒖蓟蛏?，而且能夠比對到所述參考基因的至少一個轉(zhuǎn)錄本的讀段序列的數(shù)目，則不會認為這部分讀段序列是重復比對而需要被過濾，從而提高RPKM和QPCR的相關性，即提高基因表達定量的準確性。附圖說明圖1為現(xiàn)有技術中RNA-seq的流程圖；圖2為本發(fā)明實施例一的流程圖（A）；圖3為本發(fā)明實施例一的流程圖（B）；圖4為本發(fā)明實施例一的讀段序列選擇示意圖；圖5是本發(fā)明實施例一的HBRR標準品和QPCR標準的相關性結果圖；圖6是本發(fā)明實施例一的HBRR標準品的重復性結果圖。具體實施方式下面通過具體實施方式結合附圖對本發(fā)明作進一步詳細說明?，F(xiàn)有的高通量測序平臺有多種，包括Roche454，IonPGM和IonProton等。本發(fā)明中的實施例以IonProton測序平臺作說明，其他測序平臺亦同樣適用本發(fā)明所提供的方法，測序平臺并不構成本發(fā)明的限制。RNA樣本的文庫構建一般包括將RNA反轉(zhuǎn)錄為DNA來進行文庫構建，RNA的提取、構建文庫等均可利用現(xiàn)有技術進行，測序文庫構建步驟一般包括打斷、末端修復、加proton接頭、擴增等，請參考圖1，測序步驟及參數(shù)可以根據(jù)不同測序平臺的建議操作說明、測試樣本種類進行調(diào)整，不構成對本發(fā)明的限制。實施例中未注明具體條件的，按照常規(guī)條件或制造商建議的條件進行；所用試劑或儀器未注明生產(chǎn)廠商的，均為可以通過市面購買獲得的常規(guī)產(chǎn)品。實施例一：本實施例采用RNA樣本構建文庫。RNA樣本使用人組織混合液RNA的微陣列質(zhì)量控制標準品（UHRR-MAQC，UniversalHumanReferenceRNA-MicroArrayQualityControl）和人腦混合液RNA微陣列質(zhì)量控制標準品（HBRR-MAQC，HumanBrainReferenceRNA-MicroArrayQualityControl），其中UHRR-MAQC標準品采購自安捷倫公司（AgilentTechnologies,Inc.），HBRR-MAQC購自Ambion公司。在其他具體實施方式中，亦可以使用其他種類的RNA標準品，或是采購自其他公司所生產(chǎn)的RNA標準品，對本發(fā)明不構成限制。本實施例構建文庫的過程如下：取總RNA樣品，用DEPC（diethylpyrocarbonate，焦碳酸二乙酯）水稀釋，混勻，65℃變性，使用dT（DynalbeadsOligo）25磁珠將總RNA中的信使RNA（mRNA）調(diào)取出來并純化；將所得mRNA與打斷試劑混合得到打斷的mRNA，再與試劑I混合進行一鏈合成反應；將一鏈合成反應后的體系與試劑II混合，進行二鏈合成反應，反應完成后，用AmpureXP磁珠純化二鏈產(chǎn)物；所得二鏈產(chǎn)物與試劑III混合進行末端修復，并用AmpureXP磁珠純化末端修復產(chǎn)物；所得末端修復產(chǎn)物與試劑IV混合進行加接頭，并用AmpureXP磁珠純化加接頭產(chǎn)物；采用PCR儀擴增，并用AmpureXP磁珠純化PCR產(chǎn)物，獲得測序文庫。構建轉(zhuǎn)錄本文庫或其它RNA文庫亦可利用現(xiàn)有方法，文庫構建并不構成本發(fā)明的限制。試劑I：0.5μl的100mM二硫蘇糖(DTT，DL-Dithiothreitol)、0.5μl的10mM脫氧核糖核苷三磷酸(dNTPMix，deoxy-ribonucleosidetriphosphate)、0.5μl的RNases抑制劑（RNaseInhibitor)。試劑II：10μlGEXSecondStrandBuffer、2μl10mMdNTPMix，0.2μl逆轉(zhuǎn)錄酶RNaseH、2.5μlDNA聚合酶I（DNAPolI）。試劑III：5μl10X末端修復緩沖液（EndRepairBuffer）、0.4μl25mMdNTPMix、1.2μlT4DNA聚合酶（T4DNAPolymerase）、0.2μlKlenowDNA聚合酶（KlenowDNAPolymerase）、1.2μlT4多聚核苷酸激酶（T4PNK）。試劑IV：2μlT4DNA連接酶（T4DNALigase）、2μlprotonAdapterOligoMix(12um)、25μl2XRapidT4DNALigaseBuffer。利用Agilent2100質(zhì)檢構建得的文庫，上機測序，獲得測序序列，即獲得讀段序列（reads）。請參考圖2至圖6，本實施例提供一種基因定量表達方法，可以快速地完成定量表達。其中在先步驟如文庫制備、PCR擴增等采用前述步驟與參數(shù)。本實施例具體包括：S100：獲取含有核酸序列信息的讀段序列readsS101：對讀段序列進行修剪（trimming）Trimming可以減少堿基序列在拼接之后產(chǎn)生的錯誤。在其他具體實施方式中，亦可以不對讀段序列進行修剪，直接進行后續(xù)步驟；或者使用校正（correct），或修剪與校正結合的方式，以進一步提高測序分析的準確率。Trimming針對讀段序列的開頭和末尾的的3到4bp，這幾個bp通常帶有測序接頭。包括低質(zhì)量reads，接頭（adapter），基因組3’端位置相同的reads。在高通量測序中，每測一個堿基會給出一個相應的質(zhì)量值（Q-Value)，可以參考公開號為CN102653784A，名稱為《用于多重核酸測序的標簽及其使用方法》的中國專利申請。質(zhì)量值可以反映測序質(zhì)量的好壞，數(shù)值越高表示測序質(zhì)量越好。因此，低質(zhì)量reads是指質(zhì)量值低于y1的堿基的數(shù)目超過該reads總堿基數(shù)目的y2%，y1的取值范圍為15＜y1≤20，y2的取值范圍為15＜y2≤25，本實施例取y1為17，y2為20。本領域人員知道，譬如Q20是指質(zhì)量值大于20的堿基在所有堿基中所占的比例，取值范圍為[0,1]，Q20數(shù)值越接近1，質(zhì)量值大于20的堿基在所有堿基中所占的比例越大。因此，低質(zhì)量reads可以描述為Q(y1)小于（100-y2）%的reads，或其他等同描述方式。譬如本實施例的低質(zhì)量reads，亦可以描述為Q17小于80%的reads，其中80來源于100-y2=100-20。譬如對于Hiseq測序平臺，y1優(yōu)選設置為20，y2優(yōu)選設置為20，則低質(zhì)量reads可以描述為Q20小于80%的reads。y1和y2的取值之間沒有必然的數(shù)值聯(lián)系，可以相同或是不同的數(shù)值。在其他具體實施方式中，y1及y2的取值可以根據(jù)樣品、測試平臺等有所調(diào)整，y1、y2越高，被篩選的reads越多，即留下的reads越少；y1、y2越低，則被篩選的reads越少，處理效率越慢。S102：將讀段序列與參考基因進行比對，獲取比對上的讀段序列基因組作圖(genomemapping)是應用界標或遺傳學標記對基因組進行精細的劃分，進而標示出堿基序列或基因排列。本實施例中利用reads與參考基因進行比對，而不是現(xiàn)有的reads和參考基因組比對，從而提高比對準確性及比對效率。對于真核生物，基因是由基因組中的外顯子拼接而成，而測序平臺測出來的是拼接之后的序列，直接和參考基因進行比對可以較為直接、準確。另外，在輸出比對結果時，本實施例是輸出所有的比對匹配結果，即如果有兩條以上的讀段序列都與參考基因比對匹配，則這兩條以上的讀段序列都會輸出，而不是只輸出唯一匹配的reads。一個基因包括多個轉(zhuǎn)錄本，很多轉(zhuǎn)錄本是來自外顯子的不同組合方式，所以有些轉(zhuǎn)錄本會有許多同源序列，所以有許多序列會比對到多個轉(zhuǎn)錄本上，因此保留所有這些堿基序列，用來判斷這些序列是否來自同一個基因。在本實施例中，應用tmap比對工具。tmap是一款適用proton測序平臺的商業(yè)比對軟件，由LifeTech.公司開發(fā)。比對的過程主要通過比對得分進行，利用設置基礎比對分值，比如本實施例設置基礎分為0，reads上的一個堿基位置匹配上參考基因加一分，一個位置錯配減一分，該位置缺失計0分等，由此對該read的比對情況進行打分，一般地，一條reads越長，與參考基因匹配程度越高，則其得分越高。在其他具體實施方式中，計分的規(guī)則可以根據(jù)實現(xiàn)的程序進行調(diào)整，譬如基礎分為100，每匹配上一個參考基因加0.1分，具體的計分規(guī)則不構成本發(fā)明的限制。在其他實施方式中，亦可以根據(jù)測序平臺的不同使用合適的商用比對軟件，比如Bowtie、SOAP2、BWA-SW等，或者是自編程序，只要該程序可以達到reads與參考基因進行比對并輸出所有的比對匹配結果的目的即可，因此具體的設置參數(shù)及比對工具并不構成本發(fā)明的限制。S103：對比對上的讀段序列進行過濾對步驟S102得出的比對讀段序列過濾，去掉含軟剪切比例超過第一預設值x1的reads，序列長度小于第二預設值x2的reads，以及比對得分小于第三預設值x3的reads。軟剪切是指沒有比對匹配的reads段，例如一條100bp的reads，共有90bp與參考序列比對匹配，但剩下的10bp沒有比對匹配，則這10bp稱為軟剪切，該reads的軟剪切比例為10%。在本實施例中，第一預設值x1為自然數(shù)，取值范圍是[10%，30%]，優(yōu)選為20%；x1越大，被過濾的reads數(shù)目越多，可能導致后面檢測到的基因數(shù)目偏少，x1如果過小，則可能導致部分錯誤的reads沒有被過濾掉。第二預設值x2為正整數(shù)，取值范圍是[15，25]，優(yōu)選為20，對于過短的序列，如10bp的reads，由于長度較短，可能會比對到參考基因的多個區(qū)域。第三預設值x3為正整數(shù)，取值范圍為[20，50]，x3過低則說明比對匹配的程度過差，易引入錯誤，x3過高則會導致reads過多的被去掉。值得注意的是，x3的取值范圍必然根據(jù)步驟S102的比對得分規(guī)則而調(diào)整，對于本實施例的proton測序平臺以及比對得分規(guī)則而適用于為[20，50]的取值范圍。在其他具體實施方式中，x1、x2、x3的具體數(shù)值可以根據(jù)測試平臺、測試樣品進行調(diào)整。x1、x2、x3之間沒有必然的數(shù)值聯(lián)系，可以相同或是不同的數(shù)值。S104：對基因表達進行定量本實施例用RPKM來定量，RPKM（readsperkilobaseofexonmodelpermillionmappedreads）是目前通用的定量歸一化的方法，定義為：RPKM=（比對到目標基因?qū)膮⒖蓟虻淖x段序列的數(shù)目）*109/（比對到所有參考基因的讀段序列的數(shù)目*目標基因的長度）。選取唯一比對到參考基因上的read作為比對到目標基因?qū)膮⒖蓟虻膔ead。對于比對到多個參考基因的read，無法區(qū)分來自哪個參考基因，因此將比對到多個參考基因的read都去掉。對于一條read比對到一個參考基因的多個同源轉(zhuǎn)錄本，或者一個參考基因的多個位置的情況，則認為只比對到該參考基因一次。當一條read比對到多個轉(zhuǎn)錄本時判斷所有比對上的轉(zhuǎn)錄本是否來自同一個基因，即是否所有比對上的轉(zhuǎn)錄本同源，如果判斷結果為是，即所有比對上的轉(zhuǎn)錄本都是來自同一個基因，則這條read并不是重復比對（multiplemap）而不需要去除；如果判斷為否，則該條read是multiplemap而需要去除，不能作為唯一比對到參考基因上的read。在本實施例中，步驟S102的顯示結果可以包括reads比對上哪些轉(zhuǎn)錄本，可以有multiplemap的顯示提示，因此可以利用基因和轉(zhuǎn)錄本對應的數(shù)據(jù)庫，來對multiplemap的reads進行過濾。然后，統(tǒng)計比對到該參考基因的總reads數(shù)目，一個基因可以存在多個轉(zhuǎn)錄本或者多個位置，但是這些read都來自同一個參考基因，不會干擾基因的定量，選取該基因的最長轉(zhuǎn)錄本代表該基因的長度?；虻拈L度越長，在同等表達水平下產(chǎn)生的read會比長度短的基因要多。因此在計算RPKM的時候除以基因的長度，能夠盡量避免基因長度對定量的影響。請參考圖4，以基因A（GeneA）為例進行說明。圖4是基因A三個轉(zhuǎn)錄本（transcript）的覆蓋（coverage）情況，分別是transcript1、transcript2、transcript3。在計算RPKM的時候，覆蓋到基因A的read數(shù)目為3，包括read1，read2，read3，其中基因的長度我們用最長的轉(zhuǎn)錄本3（transcript3）的長度來當做該基因的總長度。對于本實施例中的RPKM計算公式，由于前述步驟中比對、過濾的設置，以及本步驟中對參數(shù)的限制選擇，使得基因表達的定量變得快速、簡單。本實施例提供的基因表達的定量準確性用QPCR的相關性作評價。這里以皮爾森相關性系數(shù)（pearsoncorrelation）作說明。皮爾森相關系數(shù)是用來反映兩個變量線性相關程度的統(tǒng)計量，皮爾森相關性系數(shù)越高，QPCR的相關性越強，基因表達的定量越準確。相關系數(shù)用r表示，其中n為樣本量，分別為兩個變量的觀測值和均值。r描述的是兩個變量間線性相關強弱的程度，絕對值越大表明相關性越強，具體公式為在其他具體實施方式中，亦可以與其他相關性系數(shù)聯(lián)合評價，如斯皮爾曼相關性系數(shù)(spearmanrelativity)等。圖5是HBRR標準品和QPCR標準的相關性結果圖，其中橫坐標是HBRR標準品的proton測序結果計算出來的RPKM值的以10為底的對數(shù)值，縱坐標是QPCR值的以10為底的對數(shù)值，一個黑點代表一個基因。該標準品的QPCR基因為1000個，即genenum為1000。經(jīng)計算，pearsoncorrelation可達到0.917，spearmanrelativity亦可達到0.868。圖6是HBRR標準品的重復性結果圖，分別使用了兩個HBRR標準品，分別命名為proton_A和proton_B以作說明上的區(qū)分，實質(zhì)并無區(qū)別。橫坐標是proton_A用proton測序得到的RPKM值的以10為底的對數(shù)值，縱坐標是重復proton_B用proton測序得到的RPKM值的以10為底的對數(shù)值。基因數(shù)目genenum為17463表示，在proton_A和proton_B中均能夠檢測到的基因個數(shù)為17463。圖5中的genenum數(shù)目與圖6的genenum數(shù)目不同是因為圖5中的genenum中的QPCR結果是標準品RNA提供方Agilent公司提供的經(jīng)過驗證的1000個，而圖6中的genenum是proton_A和proton_B都能測出來的基因，但是其中很大部分基因仍未有經(jīng)過驗證的QPCR結果?？梢钥闯觯瑘D6的pearsoncorrelation可達到0.997，用spearmanrelativity亦能達到0.985，說明對于不同的樣品的定量結果具有很好的重復性。對于UHRR的標準品，QPCR的相關性亦達到0.86以上，詳細的結果請見表1。以8個樣本為例，其中UHRR為4個，HBRR為4個，其中樣本的名稱不具有實質(zhì)意義，只是作為不同樣本的區(qū)分之用。表1不同樣本的基因定量表達評價然后，可以根據(jù)國際標準化的基因功能分類體系GeneOntology全面描述基因的屬性，其中包括基因的分子功能molecularfunction、所處的細胞位置cellularcomponent、參與的生物過程biologicalprocess。亦可以通過比較不同樣本間的數(shù)據(jù)從而篩選出差異表達的基因，后續(xù)分析中的差異基因表達模式聚類分析，GeneOntology功能顯著性富集分析，Pathway顯著性富集分析，蛋白互作網(wǎng)絡分析均是基于差異表達基因。本領域技術人員可以理解，上述實施方式中各種方法的全部或部分步驟可以通過程序來指令相關硬件完成，該程序可以存儲于一計算機可讀存儲介質(zhì)中，存儲介質(zhì)可以包括：只讀存儲器、隨機存儲器、磁盤或光盤等。依據(jù)本發(fā)明的另一方面還提供一種基因定量表達的裝置，包括：數(shù)據(jù)輸入單元，用于輸入數(shù)據(jù)；數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)；存儲單元，用于存儲數(shù)據(jù)，其中包括可執(zhí)行的程序；處理器，與上述數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元及存儲單元數(shù)據(jù)連接，用于執(zhí)行存儲單元中存儲的可執(zhí)行的程序，該程序的執(zhí)行包括完成上述實施方式中各種方法的全部或部分步驟。以上內(nèi)容是結合具體的實施方式對本發(fā)明所作的進一步詳細說明，不能認定本發(fā)明的具體實施只局限于這些說明。對于本發(fā)明所屬技術領域的普通技術人員來說，在不脫離本發(fā)明構思的前提下，還可以做出若干簡單推演或替換。