如何以Excel為基礎做RNA-seq 數據可視化分析

沒有任何計算機語言基礎的小夥伴們該如何分析NGS（Next Generation Sequence）的數據呢？看到文章當中炫酷的，多彩的，數據分析結果，是不是一邊流著口水，一邊望洋興嘆呢？實驗室窮，各種商業軟體Strand RNA, Gene Investigator又貴的離譜，網上魚龍混雜的的免費軟體用哪一款好呢？

當然，可以去學習R語言，https://www.bioconductor.org/help/workflows/rnaseqGene/ 但是，等你學成歸來，估計導師都等的不耐煩了。

作為一個分子生物學的學生，沒有半點計算機語言基礎，又被坑在一個經費緊張的實驗室，該怎麼入門數據可視化分析呢？

我當初面對著好容易得到的一批數據，也是一臉無奈。和導師商量，要不買一個分析軟體？導師問：「多少錢啊？」

我聯繫了軟體代理商以後，跟他說：「大概2萬塊一年使用權。「

「兩萬？一年使用權？我再考慮考慮。」

於是，就一直考慮了三四年……

那時，我面對一堆數據無可奈何，急的就差看著數據一個一個的找了，但是想想，自己好歹擁有一顆Ph.D水平的大腦怎麼能用這種費力不討好的愚蠢方法呢？於是，等不及學Phyton、R，我以EXCEL為基礎，做了一些RNA-seq數據的可視化分析。

現在，我將我的經驗做出一些簡單的總結，給那些和曾經和我一樣迷茫的分子生物學的小白分享一下我的一些經驗。會計算機或者生物信息學大佬們，看到有不足的地方，希望能批評指正。

通過NGS高通量測序，可以看到全基因組的基因差異化表達，這裡的基因差異表達，說的就是mRNA的差異性表達，而這個全基因組的基因差異化表達的實現，就是用以二代測序（NGS）為基礎的mRNA sequence，這也是目前分子生物學最常用的一種方法。mRNA，全稱message RNA就是通常所說的信使RNA，作為DNA與蛋白質之間的媒介，其表達，受到環境、時間等等條件的影響。很多表型，都是基因差異表達造成的，但是如何將這些數據變成論文中易讀的FigureABC呢？這就需要一下方法將所有數據可視化（Visualizating all the data），比如，下面幾張圖，代表了集中典型的RNA-seq data可視化的結果。

• A: Volcano Plot (火山圖)
• B: Venn Diagram (韋恩圖)
• C: Intensity Scatter Plot (密度分散圖)
• D: Heatmap (熱圖)
• E: Gene Ontology (GO) analysis
• F: GO analysis

如果你認真看完本篇文章，你也能夠做出從A到F這些圖。

結合之前文章提到的配色方案(關於配色問題，參考以下鏈接），你會做出很漂亮的A-F，然後開心地放在PPT和論文中！

漆黑的師兄：為什麼 PNAS 的圖畫得往往沒有《Nature》、《Science》精緻？

RNA-seq的一般流程是這樣的

mRNA-seq 測序的最原始結果是核苷酸序列，文件大小在4-10GB之間，分析起來相當的繁瑣，分析也不是一般的小本本能跑起來的，好在一般來說，我們從測序公司拿到的數據，已經是經過比對和Mapping之後的結果，公司給這個樣子的表格：

做好分析，首先要理解每一列每一行代表的意思，那麼每一列代表的是什麼意思呢？

• 1) GeneID: 分析出來的基因在基因資料庫中的名稱。
• 2) 從第二列到第七列,
  分別是WT-raw-counts(-1,-2,-3)和lsd1-counts(-1,-2,-3),可以簡單的理解為改Gene在這次測序中讀到的總次數。
• 3) 從第八列到第十三列，分別是WT-CPM(-1,-2,-3)和lsd1-CPM(-1,-2,-3)，可以簡單的理解成一個是對照組WT另一個是實驗組（Var2）,CPM( counts per million reads)每一百萬次讀數該基因記到的次數。
• 4) LogFC （FC: Fold changes）,一般是以2為底數的對數值，代表這個基因在實驗組的表達是對照組的多少倍。LogFC≥1（upregulated上調）或≤-1（downregulated下調），即表達量大於等於2或者小於等於2，認為是比較顯著的，是很重要的一個參數。
• 5) LogCPM (CPM, counts per million reads)一般好像用不到這個值，不知道為什麼。
• 6) P_value(基因表達差異的可信度)一般取小於0.05。
• 7) FDR（False discovering rate）也是指示基因表達差異性是否可靠的值，一般取小於0.05。
• 8) Gene Type，這是根據資料庫對基因的注釋（Gene Annotation），一個基因可以從三個方面進行闡釋，1）Biological Process; 2) Molecular Function; 3) Cellular Component，即這個基因參與到的生物進程（信號通路），表達的蛋白有什麼樣的分子功能，表達的蛋白亞細胞定位在哪裡，三個方面。

當初我拿到這個含有上千甚至上萬個基因的Excel表格是，一陣莫名激動，感覺自己馬上就要發文章成為大佬受萬人膜拜，前途亮的能瞎眼。然而，一股莫名的憂傷湧現，那麼多基因，我該怎麼辦呢？總不能把基因都貼滿論文吧？好像沒有這樣搞的啊！

後來才知道，第一步，先什麼都不管，直接把所有基因進行可視化處理，就是上文提到的圖A，火山圖； 圖D, 熱圖； 根據情況用C。

Volcano Plot

繪製起來比較簡單，用到的值就是上文EXCEL表格中的P_value和LogFC。如果要繪製成像圖A，灰色的是不顯著的基因，紅色的是顯著上調的基因,就需要在EXCEL里的根據logFC 和P_value,把基因們分為三類（這點，請熟練運用EXCEL中的篩選，和排序功能）

• 1）P_value大於0.05，不顯著的，
• 2）P_value小於0.05， LogFC大於等於1，顯著上調的，
• 3）P-_value小於0.05， logFC小於等於-1，顯著下調的。

**選中以上三組數據中的任意一組，比如1）P_value大於0.05，不顯著的，

LogFC和P_value兩列，X軸設置為LogFC, Y軸設置為P_value

點擊 插入圖標>選擇XY 散點圖>

是不是得到一個很奇怪的圖，沒關係，

• a) 可以點擊Y軸，設置Y軸格式，選中 「對數刻度」和「逆序刻度值」
• b) 點擊X軸，設置X軸格式，選擇「標籤」>「高」，
• c) 將縱坐標交叉，自定義設置成最小值
• d) 選中圖中的點，設置一下數據系列格式，找到「標記」，設置一下標記的大小和顏色（文中設置成灰色）
• e) 右擊散點圖，選擇「添加數據」分別把其他兩組數據添加進去，Volcano Plot就做好了。

Heatmap

繪製起來要稍微麻煩一些，我暫時還沒找到可以用Excel繪製Heatmap比較簡單的方法或者程序包。為了簡便起見，我選擇了一個免費軟體MeV (Multiple experiment viewer) 下載地址：http://www.mybiosoftware.com/mev-4-6-2-multiple-experiment-viewer.html

這個軟體功能很強大，也可以用作Microarray data analysis，有包括HCL，K-mean Clustering 在內的好幾種聚類分析方案，重要的是一鍵可用。好用雖然好用，但是前期的數據準備還是要在EXCEL裡面先做好，這次用的是CPM，即上文提到的WT-CPM(-1,-2,-3)和lsd1-CPM(-1,-2,-3)。

• a) 首先，算一下均值， 「=mean(wt_cpm-1, wt_cpm-2,wt-cpm-3)」和「=mean(lsd1_cpm-1, lsd1_cpm-2, lsd1-cpm-3)」。如果，你的數據很多，比如在圖D中，還有一個lsd1NahG,也可以算好均值，但是要注意的是，將每一個基因一一對應好，不要有錯漏，有些基因可能在一個樣本里沒有檢測到，而在另一個樣本里有。那麼基因樣本的預設值，可設置0.00000001，意即表達量極低，以至於不能檢測到。也可以設置能NA，意為錯誤或預設值。兩種處理方式，在文章中都是可以接受的，但為了美觀，我通常會選擇前者。後者，用MeV作圖時會出現一塊灰色的區域。

• b) 接著，把整理好的數據，複製到txt文檔中，因為MeV只支持txt格式的文本文檔數據。
• c) 打開軟體,File>load data ,注意軟體左下方的紅字：Click the upper-leftmost expression
  value. Click the load button to finish.

• d) 依次點擊Adjust data 選擇
• e)

• f)

• g)

HCL clustering可以用軟體的默認設置，具體每一條什麼意思，大家可以查詢一下，解釋起來很麻煩，而且也不是本文的重點。如果樣品不是很多可以不用選Sample Tree. 具體的軟體使用，可以參考軟體使用手冊，一定要讀一遍手冊！！！然後調整一下字體什麼的，就做成圖D了。

韋恩圖

的意思一般是，找到兩個或者三個樣本之間相同的差異表達基因，比如與wt比較，樣本lsd1中有很多LogFC>1的基因，lsd1NahG也有很多LogFC>1的基因，那麼這些基因中有沒有重疊的（Overlapped Genes）,這時就可以用韋恩圖了。這時，用的就是GeneID,可以

a) 同時選中兩列GeneID,用EXCEL中的「條件格式」>「突出顯示單元格規則」「重複值」把重複值標記成紅色選項,手動作圖。

b) 也可以用現成的工具 http://bioinfogp.cnb.csic.es/tools/venny/

這樣 圖B就做好了。

GO analysis

現在，我們已經把這些基因分類了，也Visualize了，其實這些炫酷的圖基本沒啥卵用的。也就是為了好看，告訴讀者，我做了RNA-seq結果是這樣的，給你看給你看，多好看啊！然而，讀者肯定會一頭霧水！這是個啥，啥，啥！！！！！為了解釋這是個啥啥啥！我們就要對這些基因進行功能性分析，也就是Gene Ontology（基因本體論）從，三個方面

• 1) Biological Process （BP);
• 2) Molecular Functionv(MF);
• 3) Cellular Component (CC)，

對這些基因進行注釋，然後找到一些線索，為接下來的研究、設計實驗做些準備。一些人開發了Excel分析包，我沒用過，因為用於GO analysis的分析軟體特別多。一款比較常用的就是 gprofiler: http://biit.cs.ut.ee/gprofiler/

這裡用到的是GeneID

把得到的數據，以excel的方式存起來，打開文檔把corrected p-value 用-Log10作一下轉化得到-log10(p-value)，簡單的按從小到大排一下序,然後選中term name，插入條形圖，就可以做出圖E。按顏色，把BP，MF，CC標出來。

圖F表述的意思其實和圖E差不多，根據自己的需要可以做出選擇，作圖F用到的也是GeneID,不同的是，它是用一個軟體Cytoscape:http://www.cytoscape.org中的App BiNGO做的，操作起來很簡單，有興趣的可以去下載軟體自己去試一試，這裡就不詳細說了。

最後，還差一個圖C,圖C是幹嘛的呢？圖C學名Ratio -intensity plot

假如，有一些基因在mutant lsd1 裡面是上調的（圖C中全部為紅色），那麼這些上調的基因，在敲出一個關鍵基因以後，表達量會不會受到影響呢？這時候就要用到圖C：密度散點圖（Intensity scatterplot）,這裡用到的值是「mean(wt_raw counts-1, wt_ raw counts -2,wt_raw counts -3)」和「mean(lsd1_raw counts -1, lsd1_ raw counts -2, lsd1_raw counts -3)」

首先，將兩組數據的均值用log2做一下轉換。然後選擇轉換後的數據，插入XY scatter，編輯一下X、Y軸就可以了。

還有一個類似的，叫做R-I (MA)plot,看起來像是把上述密度散點圖旋轉45°，也是一種比較初級的Normalization的方法，一般的用於全基因組基因可視化繪製，用的也是mean(wt_raw counts-1, wt_ raw counts -2,wt_raw counts -3)」和「mean(lsd1_ raw counts -1, lsd1_ raw counts -2, lsd1_raw counts -3)」，也做轉換 不過不同的是：

I=log2(wt_mean*lsd1_mean)：作為X軸

R=log2(wt_mean/lsd1_mean)：作為Y軸

大概是這個樣子的。

好了，如果你認真看完這篇貌似技術貼的水貼，那麼就可以解決你的燃眉之急，解決一大半類似的數據分析問題了。

我現在也在學習R，以及Python。因為，要真正掌握生物信息學的精髓，這些是遠遠不夠的。雖然我只是一隻分子植物學喵！如果你能熟練運用R語言，那麼做出以上圖，將會更加容易。

最後，歡迎大家提出各種意見，歡迎各種姿勢騷擾，交流學習經驗。

歡迎大家轉載，但需註明出處哦～

推薦閱讀：