转录组不求人系列(七)：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介绍完之后，我们一一详细解决这些问题!