基因表達的定量方法及裝置與流程
2023-06-12 17:51:56
本發明涉及基因組學及生物信息學技術領域，具體涉及一種基因表達的定量方法及裝置。
背景技術：
轉錄組測序技術（RNA-seq，RNAsequencing）是把小RNA（RibonucleicAcid，核糖核酸）、mRNA和非編碼RNA等或者其中一些用高通量測序技術把它們的序列測出來。目前RNA-seq測序平臺有多種，包括Hiseq、RocheFLX、IlluminaSolexa、ABIsolid等。不同測序平臺的測序原理有所不同，但測序步驟基本包括文庫製備，聚合酶鏈式反應（PCR，PolymeraseChainReaction）擴增等。通過RNA-seq，科研工作者能夠獲得生物中基因表達的情況，研究不同個體、不同時期、不同形態的組織的基因表達水平的差異。中國專利申請（申請號：201110283718.2，名稱：一種分析基因表達定量的方法）基於Illumina平臺公開一種分析基因表達定量的方法，可以克服數字基因表達譜（DGE，DigitalGeneExpression）技術對CATG位點和參考基因完整性依賴性強的缺點。但是，該方法測序分析需時較長，勞動效率有待提高。
技術實現要素：
本發明提供一種基因表達的定量方法及裝置，可以快速地完成基因表達的定量。依據本發明的一方面提供一種基因表達的定量方法,包括：獲取含有核酸序列信息的讀段序列；將讀段序列與所有參考基因進行比對，獲取比對上的讀段序列；對比對上的讀段序列進行過濾，捨去軟剪切比例超過第一預設值，序列長度小於第二預設值，以及比對得分小於第三預設值的讀段序列，軟剪切比例是指沒有比對上的鹼基數目佔該讀段序列總鹼基數目的比例；比對得分是按照每個讀段序列與參考基因的匹配程度以及讀段序列的長度而確定的數值；對於已過濾的讀段序列，使用每百萬讀段序列中來自目標基因每千鹼基長度的讀段序列數目RPKM對所述目標基因表達進行定量，定義為RPKM=（比對到目標基因對應的參考基因的讀段序列的數目）*109/（比對到所有參考基因的讀段序列的數目*目標基因的長度）。優選地，比對到目標基因對應的參考基因的讀段序列的數目是指只能比對到目標基因對應的參考基因上，而且能夠比對到所述參考基因的至少一個轉錄本的讀段序列的數目；目標基因的長度是指目標基因的所有轉錄本中最長的轉錄本的長度。依據本發明的另一方面提供一種基因表達的定量裝置，包括：數據輸入單元，用於輸入數據；數據輸出單元，用於輸出數據；存儲單元，用於存儲數據，其中包括可執行的程序；處理器，與數據輸入單元、數據輸出單元及存儲單元數據連接，用於執行存儲單元中存儲的可執行的程序，該程序的執行包括完成上述基因表達的定量方法。本發明的有益效果是：通過將讀段序列與參考基因進行比對，而不是現有的與參考基因組進行比對，可以簡化比對過程，提高比對效率。特別地，比對到目標基因對應的參考基因的讀段序列的數目是指只能比對到目標基因對應的參考基因上，而且能夠比對到所述參考基因的至少一個轉錄本的讀段序列的數目，則不會認為這部分讀段序列是重複比對而需要被過濾，從而提高RPKM和QPCR的相關性，即提高基因表達定量的準確性。附圖說明圖1為現有技術中RNA-seq的流程圖；圖2為本發明實施例一的流程圖（A）；圖3為本發明實施例一的流程圖（B）；圖4為本發明實施例一的讀段序列選擇示意圖；圖5是本發明實施例一的HBRR標準品和QPCR標準的相關性結果圖；圖6是本發明實施例一的HBRR標準品的重複性結果圖。具體實施方式下面通過具體實施方式結合附圖對本發明作進一步詳細說明。現有的高通量測序平臺有多種，包括Roche454，IonPGM和IonProton等。本發明中的實施例以IonProton測序平臺作說明，其他測序平臺亦同樣適用本發明所提供的方法，測序平臺並不構成本發明的限制。RNA樣本的文庫構建一般包括將RNA反轉錄為DNA來進行文庫構建，RNA的提取、構建文庫等均可利用現有技術進行，測序文庫構建步驟一般包括打斷、末端修復、加proton接頭、擴增等，請參考圖1，測序步驟及參數可以根據不同測序平臺的建議操作說明、測試樣本種類進行調整，不構成對本發明的限制。實施例中未註明具體條件的，按照常規條件或製造商建議的條件進行；所用試劑或儀器未註明生產廠商的，均為可以通過市面購買獲得的常規產品。實施例一：本實施例採用RNA樣本構建文庫。RNA樣本使用人組織混合液RNA的微陣列質量控制標準品（UHRR-MAQC，UniversalHumanReferenceRNA-MicroArrayQualityControl）和人腦混合液RNA微陣列質量控制標準品（HBRR-MAQC，HumanBrainReferenceRNA-MicroArrayQualityControl），其中UHRR-MAQC標準品採購自安捷倫公司（AgilentTechnologies,Inc.），HBRR-MAQC購自Ambion公司。在其他具體實施方式中，亦可以使用其他種類的RNA標準品，或是採購自其他公司所生產的RNA標準品，對本發明不構成限制。本實施例構建文庫的過程如下：取總RNA樣品，用DEPC（diethylpyrocarbonate，焦碳酸二乙酯）水稀釋，混勻，65℃變性，使用dT（DynalbeadsOligo）25磁珠將總RNA中的信使RNA（mRNA）調取出來並純化；將所得mRNA與打斷試劑混合得到打斷的mRNA，再與試劑I混合進行一鏈合成反應；將一鏈合成反應後的體系與試劑II混合，進行二鏈合成反應，反應完成後，用AmpureXP磁珠純化二鏈產物；所得二鏈產物與試劑III混合進行末端修復，並用AmpureXP磁珠純化末端修復產物；所得末端修復產物與試劑IV混合進行加接頭，並用AmpureXP磁珠純化加接頭產物；採用PCR儀擴增，並用AmpureXP磁珠純化PCR產物，獲得測序文庫。構建轉錄本文庫或其它RNA文庫亦可利用現有方法，文庫構建並不構成本發明的限制。試劑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逆轉錄酶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質檢構建得的文庫，上機測序，獲得測序序列，即獲得讀段序列（reads）。請參考圖2至圖6，本實施例提供一種基因定量表達方法，可以快速地完成定量表達。其中在先步驟如文庫製備、PCR擴增等採用前述步驟與參數。本實施例具體包括：S100：獲取含有核酸序列信息的讀段序列readsS101：對讀段序列進行修剪（trimming）Trimming可以減少鹼基序列在拼接之後產生的錯誤。在其他具體實施方式中，亦可以不對讀段序列進行修剪，直接進行後續步驟；或者使用校正（correct），或修剪與校正結合的方式，以進一步提高測序分析的準確率。Trimming針對讀段序列的開頭和末尾的的3到4bp，這幾個bp通常帶有測序接頭。包括低質量reads，接頭（adapter），基因組3』端位置相同的reads。在高通量測序中，每測一個鹼基會給出一個相應的質量值（Q-Value)，可以參考公開號為CN102653784A，名稱為《用於多重核酸測序的標籤及其使用方法》的中國專利申請。質量值可以反映測序質量的好壞，數值越高表示測序質量越好。因此，低質量reads是指質量值低於y1的鹼基的數目超過該reads總鹼基數目的y2%，y1的取值範圍為15＜y1≤20，y2的取值範圍為15＜y2≤25，本實施例取y1為17，y2為20。本領域人員知道，譬如Q20是指質量值大於20的鹼基在所有鹼基中所佔的比例，取值範圍為[0,1]，Q20數值越接近1，質量值大於20的鹼基在所有鹼基中所佔的比例越大。因此，低質量reads可以描述為Q(y1)小於（100-y2）%的reads，或其他等同描述方式。譬如本實施例的低質量reads，亦可以描述為Q17小於80%的reads，其中80來源於100-y2=100-20。譬如對於Hiseq測序平臺，y1優選設置為20，y2優選設置為20，則低質量reads可以描述為Q20小於80%的reads。y1和y2的取值之間沒有必然的數值聯繫，可以相同或是不同的數值。在其他具體實施方式中，y1及y2的取值可以根據樣品、測試平臺等有所調整，y1、y2越高，被篩選的reads越多，即留下的reads越少；y1、y2越低，則被篩選的reads越少，處理效率越慢。S102：將讀段序列與參考基因進行比對，獲取比對上的讀段序列基因組作圖(genomemapping)是應用界標或遺傳學標記對基因組進行精細的劃分，進而標示出鹼基序列或基因排列。本實施例中利用reads與參考基因進行比對，而不是現有的reads和參考基因組比對，從而提高比對準確性及比對效率。對於真核生物，基因是由基因組中的外顯子拼接而成，而測序平臺測出來的是拼接之後的序列，直接和參考基因進行比對可以較為直接、準確。另外，在輸出比對結果時，本實施例是輸出所有的比對匹配結果，即如果有兩條以上的讀段序列都與參考基因比對匹配，則這兩條以上的讀段序列都會輸出，而不是只輸出唯一匹配的reads。一個基因包括多個轉錄本，很多轉錄本是來自外顯子的不同組合方式，所以有些轉錄本會有許多同源序列，所以有許多序列會比對到多個轉錄本上，因此保留所有這些鹼基序列，用來判斷這些序列是否來自同一個基因。在本實施例中，應用tmap比對工具。tmap是一款適用proton測序平臺的商業比對軟體，由LifeTech.公司開發。比對的過程主要通過比對得分進行，利用設置基礎比對分值，比如本實施例設置基礎分為0，reads上的一個鹼基位置匹配上參考基因加一分，一個位置錯配減一分，該位置缺失計0分等，由此對該read的比對情況進行打分，一般地，一條reads越長，與參考基因匹配程度越高，則其得分越高。在其他具體實施方式中，計分的規則可以根據實現的程序進行調整，譬如基礎分為100，每匹配上一個參考基因加0.1分，具體的計分規則不構成本發明的限制。在其他實施方式中，亦可以根據測序平臺的不同使用合適的商用比對軟體，比如Bowtie、SOAP2、BWA-SW等，或者是自編程序，只要該程序可以達到reads與參考基因進行比對並輸出所有的比對匹配結果的目的即可，因此具體的設置參數及比對工具並不構成本發明的限制。S103：對比對上的讀段序列進行過濾對步驟S102得出的比對讀段序列過濾，去掉含軟剪切比例超過第一預設值x1的reads，序列長度小於第二預設值x2的reads，以及比對得分小於第三預設值x3的reads。軟剪切是指沒有比對匹配的reads段，例如一條100bp的reads，共有90bp與參考序列比對匹配，但剩下的10bp沒有比對匹配，則這10bp稱為軟剪切，該reads的軟剪切比例為10%。在本實施例中，第一預設值x1為自然數，取值範圍是[10%，30%]，優選為20%；x1越大，被過濾的reads數目越多，可能導致後面檢測到的基因數目偏少，x1如果過小，則可能導致部分錯誤的reads沒有被過濾掉。第二預設值x2為正整數，取值範圍是[15，25]，優選為20，對於過短的序列，如10bp的reads，由於長度較短，可能會比對到參考基因的多個區域。第三預設值x3為正整數，取值範圍為[20，50]，x3過低則說明比對匹配的程度過差，易引入錯誤，x3過高則會導致reads過多的被去掉。值得注意的是，x3的取值範圍必然根據步驟S102的比對得分規則而調整，對於本實施例的proton測序平臺以及比對得分規則而適用於為[20，50]的取值範圍。在其他具體實施方式中，x1、x2、x3的具體數值可以根據測試平臺、測試樣品進行調整。x1、x2、x3之間沒有必然的數值聯繫，可以相同或是不同的數值。S104：對基因表達進行定量本實施例用RPKM來定量，RPKM（readsperkilobaseofexonmodelpermillionmappedreads）是目前通用的定量歸一化的方法，定義為：RPKM=（比對到目標基因對應的參考基因的讀段序列的數目）*109/（比對到所有參考基因的讀段序列的數目*目標基因的長度）。選取唯一比對到參考基因上的read作為比對到目標基因對應的參考基因的read。對於比對到多個參考基因的read，無法區分來自哪個參考基因，因此將比對到多個參考基因的read都去掉。對於一條read比對到一個參考基因的多個同源轉錄本，或者一個參考基因的多個位置的情況，則認為只比對到該參考基因一次。當一條read比對到多個轉錄本時判斷所有比對上的轉錄本是否來自同一個基因，即是否所有比對上的轉錄本同源，如果判斷結果為是，即所有比對上的轉錄本都是來自同一個基因，則這條read並不是重複比對（multiplemap）而不需要去除；如果判斷為否，則該條read是multiplemap而需要去除，不能作為唯一比對到參考基因上的read。在本實施例中，步驟S102的顯示結果可以包括reads比對上哪些轉錄本，可以有multiplemap的顯示提示，因此可以利用基因和轉錄本對應的資料庫，來對multiplemap的reads進行過濾。然後，統計比對到該參考基因的總reads數目，一個基因可以存在多個轉錄本或者多個位置，但是這些read都來自同一個參考基因，不會干擾基因的定量，選取該基因的最長轉錄本代表該基因的長度。基因的長度越長，在同等表達水平下產生的read會比長度短的基因要多。因此在計算RPKM的時候除以基因的長度，能夠儘量避免基因長度對定量的影響。請參考圖4，以基因A（GeneA）為例進行說明。圖4是基因A三個轉錄本（transcript）的覆蓋（coverage）情況，分別是transcript1、transcript2、transcript3。在計算RPKM的時候，覆蓋到基因A的read數目為3，包括read1，read2，read3，其中基因的長度我們用最長的轉錄本3（transcript3）的長度來當做該基因的總長度。對於本實施例中的RPKM計算公式，由於前述步驟中比對、過濾的設置，以及本步驟中對參數的限制選擇，使得基因表達的定量變得快速、簡單。本實施例提供的基因表達的定量準確性用QPCR的相關性作評價。這裡以皮爾森相關性係數（pearsoncorrelation）作說明。皮爾森相關係數是用來反映兩個變量線性相關程度的統計量，皮爾森相關性係數越高，QPCR的相關性越強，基因表達的定量越準確。相關係數用r表示，其中n為樣本量，分別為兩個變量的觀測值和均值。r描述的是兩個變量間線性相關強弱的程度，絕對值越大表明相關性越強，具體公式為在其他具體實施方式中，亦可以與其他相關性係數聯合評價，如斯皮爾曼相關性係數(spearmanrelativity)等。圖5是HBRR標準品和QPCR標準的相關性結果圖，其中橫坐標是HBRR標準品的proton測序結果計算出來的RPKM值的以10為底的對數值，縱坐標是QPCR值的以10為底的對數值，一個黑點代表一個基因。該標準品的QPCR基因為1000個，即genenum為1000。經計算，pearsoncorrelation可達到0.917，spearmanrelativity亦可達到0.868。圖6是HBRR標準品的重複性結果圖，分別使用了兩個HBRR標準品，分別命名為proton_A和proton_B以作說明上的區分，實質並無區別。橫坐標是proton_A用proton測序得到的RPKM值的以10為底的對數值，縱坐標是重複proton_B用proton測序得到的RPKM值的以10為底的對數值。基因數目genenum為17463表示，在proton_A和proton_B中均能夠檢測到的基因個數為17463。圖5中的genenum數目與圖6的genenum數目不同是因為圖5中的genenum中的QPCR結果是標準品RNA提供方Agilent公司提供的經過驗證的1000個，而圖6中的genenum是proton_A和proton_B都能測出來的基因，但是其中很大部分基因仍未有經過驗證的QPCR結果。可以看出，圖6的pearsoncorrelation可達到0.997，用spearmanrelativity亦能達到0.985，說明對於不同的樣品的定量結果具有很好的重複性。對於UHRR的標準品，QPCR的相關性亦達到0.86以上，詳細的結果請見表1。以8個樣本為例，其中UHRR為4個，HBRR為4個，其中樣本的名稱不具有實質意義，只是作為不同樣本的區分之用。表1不同樣本的基因定量表達評價然後，可以根據國際標準化的基因功能分類體系GeneOntology全面描述基因的屬性，其中包括基因的分子功能molecularfunction、所處的細胞位置cellularcomponent、參與的生物過程biologicalprocess。亦可以通過比較不同樣本間的數據從而篩選出差異表達的基因，後續分析中的差異基因表達模式聚類分析，GeneOntology功能顯著性富集分析，Pathway顯著性富集分析，蛋白互作網絡分析均是基於差異表達基因。本領域技術人員可以理解，上述實施方式中各種方法的全部或部分步驟可以通過程序來指令相關硬體完成，該程序可以存儲於一計算機可讀存儲介質中，存儲介質可以包括：只讀存儲器、隨機存儲器、磁碟或光碟等。依據本發明的另一方面還提供一種基因定量表達的裝置，包括：數據輸入單元，用於輸入數據；數據輸出單元，用於輸出數據；存儲單元，用於存儲數據，其中包括可執行的程序；處理器，與上述數據輸入單元、數據輸出單元及存儲單元數據連接，用於執行存儲單元中存儲的可執行的程序，該程序的執行包括完成上述實施方式中各種方法的全部或部分步驟。以上內容是結合具體的實施方式對本發明所作的進一步詳細說明，不能認定本發明的具體實施只局限於這些說明。對於本發明所屬技術領域的普通技術人員來說，在不脫離本發明構思的前提下，還可以做出若干簡單推演或替換。