四、DESeq2差異基因分析

獲得reads-counts之後，我們就可以開展差異基因分析了。我們以subread中的featureCounts工具得到的counts_id.txt為例，來進行後續的差異基因分析。
目前常見的差異基因分析工具有DESeq2、limma包等等，此處以DESeq2為工具來進行差異基因的篩選。

First step:獲取表達矩陣和分組資訊

進行差異基因分析之前，首先要獲取表達矩陣和分組資訊。我們的表達矩陣是剛才用featureCounts定量得到的counts_id.txt ，經過格式處理之後，是這樣（部分擷取）：

第一列是基因ID，之後的列都是樣本ID
每一行代表不同的基因在不同樣本中的表達量.

我們選day0和day1做比較，
為了方便，分組矩陣的製作我在R裡面完成，輸入如下程式碼：

us_count<-read.csv("C:\\Users\\admin\\Desktop\\RNA_seq_Ara_counts_day1_day0.csv",head=T,row.names=1) #輸入表達矩陣資料路徑
us_count<-round(us_count,digits=0) #將輸入資料取整

#準備
us_count<-as.matrix(us_count) #將資料轉換為矩陣格式
condition<-factor(c("day0","day0","day0"
 
,"day0","day1","day1","day1","day1")) ## 設定分組資訊，建立環境(8個樣本，2組處理)
coldata<-data.frame(row.names=colnames(us_count),condition)  

coldata #展示coldata內容
condition  #展示環境

上面的程式碼執行之後，我們的分組資訊就是這樣噠：

現在，就可以正式使用DESeq2做差異基因分析了，總共其實只有三步：

• 構建dds矩陣
• 對dds矩陣進行標準化
• 提取結果並繪製火山圖

程式碼如下：

library(DESeq2) #使用library函式載入DEseq2包
 

##構建dds矩陣 
dds<-DESeqDataSetFromMatrix(us_count,coldata,design=~condition)
head(dds) #檢視構建好的矩陣

##進行差異分析
dds<-DESeq(dds) #對原始的dds進行標準化
resultsNames(dds)  #檢視結果名稱
res<-results(dds) #用results函式提取結果，並賦值給res變數
summary(res) #檢視結果
plotMA(res,ylim=c(-2,2)) 
mcols(res,use.names=TRUE)

plot(res$log2FoldChange,res$pvalue) #繪製火山圖

#提取差異基因
res <- res[order(res$padj),]
resdata <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
deseq_res<-data.frame(resdata)
up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) #提取上調差異表達基因
down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -1)) #提取下調差異表達基因

write.csv(up_diff_result,"C:\\Users\\admin\\Desktop\\上調_day0_VS_day1_diff_results.csv") #輸出上調基因
write.csv(down_diff_result,"C:\\Users\\admin\\Desktop\\下調_day0_VS_day1_diff_results.csv") #輸出下調基因

至此，差異基因就成功提取了，看看火山圖（火山圖繪製是這句程式碼起效果：plot(res$log2FoldChange,res$pvalue) #其實就是對log2foldchange和pvalue作圖）：

五、 總結

自己之前處理線蟲資料主要是用的RSEM（加bowtie2）來做RNA-seq的序列比對和生物學定量，但是沒有接觸過salmon和subread。這次使用這兩個工具完成了植物的RNA-seq實戰，對這些工具有了更深的理解，以後自己如果要開發相關軟體，也要多用用，比較工具之間的優劣。