数据读入

使用airway中的数据，SummarizedExperiment格式。

data("airway")
airway

#定义condition
airway$dex <- relevel(airway$dex, "untrt")

#创建DEseq对象
dds <- DESeqDataSet(airway, design = ~ cell + dex)

class: RangedSummarizedExperiment 
dim: 64102 8 
metadata(1): ''
assays(1): counts
rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
rowData names(0):
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample

如果没有SummarizedExperiment格式也可以使用DESeqDataSetFromMatrix导入。

对于运行DESeq2模型，您可以使用R的公式表示法来表达任何固定效果的实验设计。请注意，DESeq2使用与baseR的lm函数相同的公式表示法。如果研究目的是确定各组的治疗效果不同的基因，则可以包括相互作用项并使用设计进行测试如 ~ group + treatment + group:treatment。

countdata <- assay(airway)
coldata <- colData(airway)
dds <- DESeqDataSetFromMatrix(countData = countdata,
                                  colData = coldata,
                                  design = ~ cell + dex)

标准化

探索性多维数据分析的许多常用统计方法，例如聚类和主成分分析（PCA），最适合通常在不同平均值范围内具有相同方差范围的数据。当方差的预期量为约跨越不同的平均值是相同的，该数据被认为是同方差。但是，对于RNA-seq计数，预期方差随平均值而增加。例如，如果直接在计数或标准化计数的矩阵上执行PCA（例如，校正测序深度的差异），则结果图通常主要取决于具有最高序列的基因计数是因为它们在样本之间显示出最大的绝对差异。避免这种情况的一种简单且经常使用的策略是取标准化计数值的对数加上1的伪计数。但是，根据伪计数的选择，现在计数最低的基因将对结果图产生很大的干扰，因为取小计数的对数实际上会夸大其方差。

#vst标准化，blind参数指定实验设计不直接用于转换
vsd <- vst(dds, blind = FALSE)
head(assay(vsd), 3)

#rlog标准化
rld <- rlog(dds, blind = FALSE)
head(assay(rld), 3)

差异分析

#预处理
dds <- dds[ rowSums(counts(dds))>1, ] ##将所有样本基因表达量之和小于1的基因过滤掉

#差异分析
dds <- DESeq(dds)
#根据contrast参数提取结果，alpha 指定p值，lfcThreshold指定logfoldchange
res <- results(dds, contrast=c("dex","trt","untrt"),alpha = 0.05,lfcThreshold = 1)
summary(res)

out of 33469 with nonzero total read count
adjusted p-value < 0.05
LFC > 1.00 (up)    : 142, 0.42%
LFC < -1.00 (down) : 80, 0.24%
outliers [1]       : 0, 0%
low counts [2]     : 13588, 41%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

基因注释

library("AnnotationDbi")
library("org.Hs.eg.db")

res$symbol <- mapIds(org.Hs.eg.db,
                     keys=row.names(res),
                     column="SYMBOL",
                     keytype="ENSEMBL",
                     multiVals="first")
res$entrez <- mapIds(org.Hs.eg.db,
                     keys=row.names(res),
                     column="ENTREZID",
                     keytype="ENSEMBL",
                     multiVals="first")

与基因表达矩阵合并并导出结果

res_deseq <- res[order(res$padj),]
#提取差异基因
diffSig <- subset(res_deseq, padj<0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diffSig <- row.names(diffSig)

#合并结果并导出
res <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
diffmRNAexp<- res %>% 
  remove_rownames() %>% 
  column_to_rownames(var = 'Row.names')

diffmRNAexp<- diffmRNAexp[diffSig,]

结果可视化

#jitter
topGene <- rownames(res)[which.min(res$padj)]
plotCounts(dds, gene = topGene, intgroup=c("dex"))

#MAplot
res <- lfcShrink(dds, coef="dex_untrt_vs_trt")
plotMA(res, ylim = c(-5, 5))

#MAplot中添加label
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

多因素差异分析

DESeq2可用于分析时程实验，例如，找到与一组基线样本相比随时间以特定于条件的方式发生反应的那些基因。在这里，我们展示了使用fission数据包进行的基本时程分析，该数据包包含裂变酵母的RNA-seq时程的基因计数。酵母暴露于氧化应激，一半的样品含有atf21基因的缺失。我们使用一个设计公式来模拟时间0处的应变差异，随时间变化的差异以及随时间变化的任何应变特定的差异（交互作用项strain:minute）。

library("fission")
data("fission")
ddsTC <- DESeqDataSet(fission, ~ strain + minute + strain:minute)

下面的代码块执行似然比测试，其中我们去除了随时间变化的应变特定差异。此测试中p值小的基因是那些在时间0后一个或多个时间点显示出菌株特异性作用的基因。因此请注意，这不会给两个菌株中以相同方式随时间上下移动的基因提供小的p值。

#test指定似然比检验，reduced指定要与之进行比较的简化公式，即去掉相关术语的完整公式。
ddsTC <- DESeq(ddsTC, test="LRT", reduced = ~ strain + minute)
resTC <- results(ddsTC)
resTC$symbol <- mcols(ddsTC)$symbol
head(resTC[order(resTC$padj),], 4)

这只是可以应用于时间序列数据的测试之一。另一个选择是将计数建模为时间的平滑函数，并包括条件与平滑函数的交互项。可以使用R中的样条基函数建立这样的模型，另一种更现代的方法是使用高斯过程（Tonner et al.2017）。

我们可以使用ggplot2随时间变化的组计数，对于调整后的p值最小的基因，测试条件相关的时间曲线，并说明时间0的差异（下图）。

fiss <- plotCounts(ddsTC, which.min(resTC$padj), 
                   intgroup = c("minute","strain"), returnData = TRUE)
fiss$minute <- as.numeric(as.character(fiss$minute))
ggplot(fiss,
  aes(x = minute, y = count, color = strain, group = strain)) + 
  geom_point() + stat_summary(fun.y=mean, geom="line") +
  scale_y_log10()

DEseq2-差异分析