【原】轉(zhuǎn)錄組不求人系列(七)：DESeq2分析轉(zhuǎn)錄組測(cè)序數(shù)據(jù)

TS的美夢(mèng) 2021-12-24

展開全文

DESeq2是一個(gè)比較常用的轉(zhuǎn)錄組分析R包，包的使用非常簡(jiǎn)單，與之前的limma包不一樣，DESeq2需要的數(shù)據(jù)是Row counts矩陣，這點(diǎn)非常重要。

首先我們讀入下載的數(shù)據(jù)，使用GSE169758高通量轉(zhuǎn)錄組數(shù)據(jù)：

setwd("F:/生物信息學(xué)")A <- read.csv("GSE169758_markdup.featco.2.counts.csv",header = T,row.names = 1)# 預(yù)處理，過(guò)濾低豐度的數(shù)據(jù)A <- A[apply(A, 1, sum) > 0, ]#過(guò)濾前#[1] 63677 12#過(guò)濾后#[1] 35072 12

counts數(shù)據(jù)其實(shí)有很多基因都是在樣品中都沒有檢測(cè)到，所以counts為0，在進(jìn)行分析之前將低豐度的數(shù)據(jù)刪除。

加載R包，并構(gòu)建分組：

library(DESeq2)metadata <- data.frame(condition = factor(rep(c('Mcc', 'Pan'), each = 6), levels = c('Mcc', 'Pan')))

分析過(guò)程很簡(jiǎn)單，只有幾句代碼，因?yàn)楹芏喙δ苷系揭粋€(gè)函數(shù)中了。首先是構(gòu)建DESeqDataSet對(duì)象，對(duì)輸入矩陣數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化等。之后DESeq()函數(shù)進(jìn)行因子大小估計(jì)、離散度估計(jì)、負(fù)二項(xiàng)模型擬合、Wald統(tǒng)計(jì)等多步在內(nèi)的過(guò)程，將得到的結(jié)果將返回至DESeqDataSet對(duì)象。獲得差異列表。

最終的所有差異分析結(jié)果都存儲(chǔ)在res中。

dds <- DESeqDataSetFromMatrix(countData = A, colData = metadata, design= ~condition)dds1 <- DESeq(dds, fitType = 'mean', minReplicatesForReplace = 7, parallel = FALSE)res <- results(dds1, contrast = c('condition', 'Pan', 'Mcc'))head(res)#log2 fold change (MLE): condition Pan vs Mcc #Wald test p-value: condition Pan vs Mcc #DataFrame with 6 rows and 6 columns# baseMean log2FoldChange lfcSE stat# <numeric> <numeric> <numeric> <numeric>#ENSG00000000003 1.62508e+02 9.388941 0.814107 11.532803#ENSG00000000005 6.10349e-02 0.401448 2.866978 0.140025#ENSG00000000419 1.74896e+03 -0.155045 0.126955 -1.221255#ENSG00000000457 2.95673e+02 0.196623 0.143674 1.368538#ENSG00000000460 4.97165e+02 -0.131487 0.190582 -0.689924#ENSG00000000938 6.50361e-01 -0.521221 1.041933 -0.500245# pvalue padj# <numeric> <numeric>#ENSG00000000003 9.01594e-31 4.34219e-29#ENSG00000000005 8.88640e-01 NA#ENSG00000000419 2.21990e-01 3.25837e-01#ENSG00000000457 1.71144e-01 2.64359e-01#ENSG00000000460 4.90242e-01 6.03378e-01#ENSG00000000938 6.16903e-01 7.15152e-01

保存結(jié)果，差異基因可以在excel表中手動(dòng)篩選。

df <- as.data.frame(res)write.csv(df,file = "DEGs.csv")

也可以通過(guò)代碼篩選：

df <- df[order(df$padj, df$log2FoldChange, decreasing = c(FALSE, TRUE)), ]#對(duì)數(shù)據(jù)進(jìn)行排序#按實(shí)際情況以不同的閾值去定義差異條件df[which(df$log2FoldChange >= 5 & df$padj < 0.01),'sig'] <- 'up'df[which(df$log2FoldChange <= -5 & df$padj < 0.01),'sig'] <- 'down'df[which(abs(df$log2FoldChange) <= 5 | df$padj >= 0.01),'sig'] <- 'none'#可以選擇性的導(dǎo)出差異結(jié)果df_sig <- subset(df, sig %in% c('up', 'down'))

到這里我想大家應(yīng)該有兩個(gè)疑問(wèn)。第一counts矩陣可以直接把值當(dāng)作表達(dá)量嗎？答案是不可以，需要進(jìn)行轉(zhuǎn)化計(jì)算FPKM才行。第二個(gè)問(wèn)題是數(shù)據(jù)行名是基因ID，不知道這具體是什么基因，沒有g(shù)ene symbol看起來(lái)直觀，這就需要對(duì)基因進(jìn)行注釋。

這兩個(gè)問(wèn)題暫且耐心留在這里，下節(jié)另外一個(gè)差異分析R包edgeR介紹完之后，我們一一詳細(xì)解決這些問(wèn)題!