【工具】GSEA分析RNA-seq數據

相信很多人一提到分析RNA-seq數據就頭大。

今天來教大家用一個R package來處理RNA-seq數據，並進行GSEA分析。

有人可能會問，幹嘛要這麼麻煩用R來處理RNA-seq？直接取實驗組和對照組的平均值做一個fold change不就可以去運行GSEA了嗎？

然而問題並非取一個fold change這麼簡單。

• 在RNA-seq中，很多值都極高，或者極低。如果我的對照組為0，實驗組為100，那麼它的fold change將得不到。
• 而且GSEA的pre-rank只能運行一組值，如果單純的去用fold change, 將不能考慮到p值。舉個例子，如果兩個基因的fold change 都是9，然而其中一個的p值是0.5 另一個是0.05，那麼肯定後者的比重要更大，然而這都無法在fold change上體現。

因此，在運行GSEA前，我們需要處理RNA-seq數據，目標是得到一個值，它既能反應基因的上調和下調，同時還能融入p值。

目前市面上處理RNA-seq的這種R package有很多，這裡我們選擇的package叫FCROS，感興趣的同學可以在網上搜與它有關的論文。

首先，打開RNA-seq的excel，把格式設置成如下：

讓第一列是基因名稱，第一行是樣本名稱。

接下來另存為，格式選擇Tab Delimited Text.txt

打開R studio：

• install.packages("fcros") #安裝fcros#
• library(fcros)#運行fcros#
• getwd() #找到所在目錄#
• setwd(「/Users/XXX/Documents/Computational
  Bio/XXX」)#設置要去的目標目錄下#
• fdata=fcrosRead(「XXXX.txt」)#導入數據#

• cont = c("NR_Pt1", "NR_Pt10", "NR_Pt12"); #設置對照組，注意名字和數據中的樣品名字一一對應#
• test= c("R_Pt13", "R_Pt15", "R_Pt19"); #設置實驗組#
• log2.opt = 0;
• trim.opt = 0.25;
• af = fcros(fdata, cont, test, log2.opt, trim.opt); #運行fcros#
• fcrosWrite(af, file = "test.txt", values = TRUE, thr = 1) #導出fcros#

導出完之後，用excel打開文件。

這裡得到ri值。ri值大於0.5，為轉錄上調；小於0.5的，為轉錄下調。為了保證數據對稱，我們設置將所有小於0.5的值都減1

操作方法為 另起一列，設置為 」 =IF(X<0.5,X-1,X) 」， 如下圖

我們把這個值設為ri2

並將格式調整成如下圖。

另存為，格式為Tab Delimited Text.txt

之後找到文件，把後綴由.txt改成.rnk

這樣就能上GSEA的pre-rank上面運行了。

打開GSEA，接著選擇Load Data。 找到剛才的XXX.rnk文件。

Load成功之後，找到Pre-rank. 如下圖

Gene Set Database可以選擇Hallmark，也可以選擇其它的。

之後點擊Run.

Hallmark的Gene Set和其它相比比較少，所以運行的塊。

看到綠色的Success之後，點擊它。

就會出來完整的報告了。這裡注意到有兩個Enrichment in phenotype

上面的那個是上調的通路

下面的是下調的通路。

點擊Snapshot可以看見更詳細的圖。

這樣一個RNA-seq的數據就用GSEA分析完成了。是不是很簡單呢？

後面一章我們將介紹用IPA來分析RNA-seq數據。

推薦閱讀：