RNA_seq(1)植物轉錄組實戰(下)之DESeq2進行差異基因分析

Candlelight_yujia發表於2018-08-22

四、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實戰,對這些工具有了更深的理解,以後自己如果要開發相關軟體,也要多用用,比較工具之間的優劣。

相關文章