PRO-seq數據分析

背景知識

大多數RNA-seq都是研究不同條件下細胞內mRNA變化。除了基因的編碼區（CDS）可以轉錄成mRNA，基因組上的其他區域也能不同程度地轉錄（例如poly A,下游區域以及Enhancer），Enhancer可以產生短的且不穩定的RNA來調控轉錄，而這種調控的錯誤會引發多種疾病，因此，理解這種調控機理十分重要，然而傳統RNA-seq技術在檢測這種不穩定的RNA方面效率很低。
而PRO-seq技術就是對傳統RNA-seq技術在這方面的改進，它可以富集并且測出剛剛被RNA聚合酶轉錄出來的新生RNA，并且精度達到一個堿基對。

相關文獻：Nature protocol Base-pair-resolution genome-wide mapping of active RNA polymerases using precision nuclear run-on (PRO-seq)
https://www.ncbi.nlm.nih.gov/pubmed/?term=27442863
http://www.docin.com/p-1691424856.html

數據來源

文章標題：Nascent RNA sequencing reveals a dynamic global transcriptional response at genes and enhancers to the natural medicinal compound celastrol
數據來源：2017年5月23日冷泉港實驗室更新的PRO-seq表達譜
實驗設計：
K562細胞系在加入雷公藤紅素（中藥的一種）后，于0min,10min,20min,40min,60min,160min共六個時間點進行測量，每次2個重復，共計12個數據。
數據下載網址：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96869

創新亮點

傳統的RNA-seq研究都是在測量mRNA的量，而通過PRO-seq，可以測出新合成的RNA，并且檢測到幾分鐘后轉錄水平的變化，這可以更精確地分辨出調控的信號傳導通路。

數據預處理

由于此數據原始數據sra太大，沒有表達矩陣，只提供了測序序列reads在染色體上分布的位置文件(bw文件)，所以需要進行數據預處理，統計每個基因上reads的數量，作為表達量，此處調用了R語言的rtracklayer包讀入bw文件，接下來使用GenomicRanges包統計每個基因上的reads數。

差異表達基因篩選

由于樣本量較少，故考慮專為小樣本設計的T檢驗方法。又因為是時間序列，不能輕易劃分成兩個大組，考慮到時間是一個連續性因素，設計了如下的篩選方法：
1、對0min和10min的兩組四個樣本進行t檢驗；接下來依次進行10min和20min；20min和40min；40min和60min；60min和160min的t檢驗，相鄰兩個時間點進行t檢驗共計5次。
2、在上述5次t檢驗中，如果有4次發生顯著性差異，且p<0.1，則說明細胞在加入雷公藤紅素后，該基因表達有著顯著性變化。
經過這樣的篩選后，共有19個差異表達的基因，詳見diffgene.txt，第一列為EntrezID。

表達量變化圖

對于上述的19個差異表達基因，隨機選取幾個繪制表達量變化圖，縱坐標為表達量，橫坐標為不同時間。
可以發現兩個基因表達量都在降低，并且在40-60之間有個轉錄反應的峰，這與文獻摘要的This transcriptional response occurred in two major waves, one within 10 minutes, and a second 40-60 minutes after treatment.相對應。

圖表 1BTBD2基因

圖表 2PEAR1基因

表達譜繪制

首先，對于差異表達基因繪制表達譜，先是只對基因聚類，可以看出來從左到右，顏色由紅到綠；這表明隨著時間增長，大部分基因的表達量都是由高變低，這與文獻摘要中提到的“雷公藤紅素會抑制大部分的基因轉錄”相吻合。

圖表 3表達譜單向聚類
接下來，對表達譜雙向聚類，可以發現同一時間測得的兩個重復試驗都能聚類到一起去。

圖表 4表達譜雙向聚類

分析與討論

1、對于時間序列數據的處理，這個相鄰兩組t檢驗的模型顯然還是太過于簡單，
2、GEO下載下來的PRO-seq數據是有作為內參對照的spikein數據，可以利用這些內參對照數據對數據進行歸一化，將預處理做得更精細。

總結

以上是生活随笔為你收集整理的PRO-seq数据分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。