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介绍完之后，我们一一详细解决这些问题!
