R語言分析菌相資料 By李孟展

以下使用的Ubuntu版本為18.10

R語言介紹
一種自由軟體程式語言與操作環境,主要用於統計分析、繪圖、資料探勘,由來自紐西蘭奧克蘭大學的羅斯·伊哈卡和羅伯特·傑特曼開發。
R內建多種統計學及數字分析功能,並透過安裝由用戶撰寫的套件(packages),大幅增強特殊功能包括繪圖、統計技術,以及編程介面和資料輸出/輸入功能等等。



在linux上安裝R、Rstudio

安裝R:
如果只是安裝R的話非常簡單,只要在terminal輸入以下指令就好:
sudo apt-get install r-base

在terminal內輸入指令:R,即可進入R語言環境,如圖:         這時的R環境為交談式介面


安裝Rstudio:
RStudio是R語言的一種整合開發環境,是免費的軟體,個人較習慣在Rstudio上編寫
直接進入Rstudio網站下載符合自己作業系統的安裝檔安裝,網址如下:
https://www.rstudio.com/products/rstudio/download/

安裝完後會在應用程式中看到Rstudio #請忽略肥宅的桌布#


打開後介面如圖,左上角可直接進行R script的編寫,左下為交談式介面,可以顯示文字的執行結果,右上會顯示現有的已儲存變數,右下會顯示執行的圖像化結果


使用Rstudio開啟OTU table進行繪圖

在這次報告中,我將使用表格形式的NGS菌相資料(OTU table)
1.利用ggplot2套件繪製barchart
2.計算出菌之間的correlation,再繪製成heatmap
OTU table表格:


讀取excel檔
在報告中我使用的讀檔指令為read.csv("檔案路徑", sep="分割符號"),因為我習慣將table存為csv格式,比較好處理
有個指令叫read.table,可以讀取大多數ASCII資料,語法格式一樣;另有套件為"xlsx",可以讀取新版本的excel檔,指令同樣為read.xlsx,只是需要額外安裝套件

使用read.csv讀入的檔案,會形成list形式的變數,如圖
#is.list()為判斷是否為list的函式#


使用ggplot2套件進行繪圖
ggplot2是一個功能相當強大的繪圖套件,可繪製各種統計圖,結果較R內建的函式更為美觀且功能更為多樣。
套件安裝的指令為 install.packages("想安裝的套件"),安裝完成後再用 library(套件) 引入套件。

雖說ggplot2的功能強大,但同時也必須給予正確的檔案格式才能畫出想要的圖
ggplot2的基本繪圖指令為:
ggplot(變數名稱, aes(x軸, y軸, fill=填入顏色的項目)) + 其他功能性函數

在畫圖之前需要注意幾件事:
1.變數形式是否正確
2.變數內容的排版是否符合函數的運作

這邊使用的變數形式必須為list或是data.frame,一般來說data.frame是大多數都支援的格式,可使用函式as.data.frame(變數)進行更改
設定ggplot2的XY軸時一定要知道ggplot2是如何進行繪圖的,不然會發生系統錯誤或是出來的圖是悲劇
ggplot2的運作模式是將給定的函數中的XY作為data中的column name,將該行的data內容作為X跟Y來繪制,舉例來說:
如下圖,假如我在函數X軸Y軸的地方填入aes(x, y),程式會將X下的10、100、10作為3個sample,而非數值,而在稍微研究了一下後會發現,用這種排版的data會無法將每一個sample(x,y,z)內的value全部導入到圖上,因為Y軸的column一次只能給一行


輸入:ggplot(a, aes(a[,1],y, fill = y)) + geom_bar(stat = "identity"),這句程式代表說我將a,b,c那行作為sample,將y那行的值作為Y值帶入,會發現這種格式的data一次只能畫出一行的data
更麻煩的是,一般來說NGS data給定的table,橫向為sample,縱向才是菌的數量,因此硬要畫的話還必須要將list轉換,十分麻煩
這邊解釋一下前面沒提到的函數:
geom_bar()代表要繪制的圖為長條圖
stat = "identity"代表將Y軸的內容視為已經統計好的數字,如果沒有加的話系統會將其中的數字視為一個元素去計算數量而非一個值



使用reshape2套件進行排版變換
為了能完美畫出我想要的圖,我使用了reshape2套件下的melt()函數來將data重新排版
melt()函數可以將原本的data轉換為ggplot2的繪圖形式,如圖:


在經由melt()變換後,可以發現函數幫我們把data整理成variable跟value兩大行,接下來就是把這兩行填入ggplot函數裡即可
可以發現執行結果就是我想要的形式,其中a$id = a[,1]aa[,1] = factor(aa[,1])只是為了填入顏色而加上的另一行
若想要在barchart顯示的不是絕對數據而是百分比的話只需要在geom_bar()內加入position='fill'即可


進行實際數據的繪圖
在掌握畫圖所需的步驟後就可以完整的用實際數據畫出barchart了,可以在右下方的視窗進行圖片匯出,也可利用函數進行儲存
在script中我加入幾個函數來改圖中的一些顯示:
scale_fill_manual(values=色表,name="圖例名稱")可以更改圖例的顏色及名稱
labs(x="X座標名稱", y = "Y座標名稱")可以更改座標名稱


使用ggpubr套件計算spearman correlation再繪制heatmap

引入ggpubr進行correlation計算
在這邊我使用的函數為corr(),它可以將column name作為目標,將column內的資料作為data,計算每個column name之間的correlation
因為通常NGS table的column names為sample,而correlation要計算的是菌之間的關係,因此,需要先將data轉置
語法為:t(變數)
轉置後需要注意兩件事:
1.若colnames不是OTU菌相,則利用colnames()函數進行更改
2.輸入的data不能有任何數字外的內容,因此若有其他的非數字內容記得進行擷取
更改完成後即可帶入函數中計算,語法如下:
corr(變數, method = "spearman")
可以用round進行四捨五入:
round(corr(變數,method = "spearman"), 取到小數點第幾位)
跑出來的結果會是一個colnames對colnames的matrix,如圖:


接下來畫圖的步驟就如同畫barchart一樣,用melt()函數轉換後用ggplot2畫出heatmap
這邊要使用的繪圖函數為geom_tile(),可以將給定的XY繪製成heatmap
除了基本的的繪圖函數外,同樣可以加入一些美化的功能函數:
scale_fill_gradient2(low = "負相關的顏色", high = "正相關的顏色", mid = "不相關時的顏色", name="圖例名稱")
labs(x="X座標名稱", y = "Y座標名稱")
使用theme()函數改變XY座標標示的角度等等顯示

結果如下圖:


報告後作業

以上資料處理的目的:
1.繪製barplot可以幫助我們快速觀察資料特性,可方便查看高達數千個菌種的菌相資料。
2.計算correlation可以幫助model的建置,以實例來說,菌相分析可以告訴我們土壤初始的菌相,但若要計算環境改變後的菌相,利用correliation建立model,可以利用土壤內少數菌相推知整個菌相,不必花費昂貴的NGS費用。
3.使用spearman correlation的目的是因為其通用性較高,菌相資料常為隨機分布的數據,使用spearman比較合適(假設菌相為兩兩影響的關係)。
4.當correlaion為0的時候代表兩者為獨立,互不影響。

使用matlab匯入excel檔,並轉置table
使用xlsread指令匯入excel檔,檔案會自動匯入成矩陣,再將變數後面加上一撇(')即可轉置。


使用python進行矩陣轉置
python有個package叫做numpy,可以進行多維陣列的運算,其中transpose、T、swapaxes指令為轉置的指令。


R語言能達成的事情十分多樣,同樣可通過迴圈、條件判斷等方法來處理大量資料,同時,可藉由長期以來各個使用者編輯的套件便利的達成分析統計的繁複步驟

Referance:
R軟體:應用統計方法(二版) 作者: 陳景祥
計算correlation並繪制heatmap的教學
關於一些實際菌相correlation的數據