【原】转录组不求人系列(七)：DESeq2分析转录组测序数据

DESeq2是一个比较常用的转录组分析R包，包的使用非常简单，与之前的limma包不一样，DESeq2需要的数据是Row counts矩阵，这点非常重要。

首先我们读入下载的数据，使用GSE169758高通量转录组数据：

setwd("F:/生物信息学")A <- read.csv("GSE169758_markdup.featco.2.counts.csv",header = T,row.names = 1)# 预处理，过滤低丰度的数据A <- A[apply(A, 1, sum) > 0, ]#过滤前#[1] 63677 12#过滤后#[1] 35072 12

counts数据其实有很多基因都是在样品中都没有检测到，所以counts为0，在进行分析之前将低丰度的数据删除。

加载R包，并构建分组：

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

分析过程很简单，只有几句代码，因为很多功能整合到一个函数中了。首先是构建DESeqDataSet对象，对输入矩阵数据进行标准化等。之后DESeq()函数进行因子大小估计、离散度估计、负二项模型拟合、Wald统计等多步在内的过程，将得到的结果将返回至DESeqDataSet对象。获得差异列表。

最终的所有差异分析结果都存储在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

保存结果，差异基因可以在excel表中手动筛选。

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

也可以通过代码筛选：

df <- df[order(df$padj, df$log2FoldChange, decreasing = c(FALSE, TRUE)), ]#对数据进行排序#按实际情况以不同的阈值去定义差异条件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'#可以选择性的导出差异结果df_sig <- subset(df, sig %in% c('up', 'down'))

到这里我想大家应该有两个疑问。第一counts矩阵可以直接把值当作表达量吗？答案是不可以，需要进行转化计算FPKM才行。第二个问题是数据行名是基因ID，不知道这具体是什么基因，没有gene symbol看起来直观，这就需要对基因进行注释。

这两个问题暂且耐心留在这里，下节另外一个差异分析R包edgeR介绍完之后，我们一一详细解决这些问题!