转录组分析差异表达

使用DEseq：

今天我们主要讲一下如何使用R包，选择 DESeq2，这是一个常用的包来进行差异表达分析。我们之所以选择 DESeq2，一个重要的优点是它能够处理有无样品重复的情况。如果没有样品重复，其他软件比如 edgeR，它是不支持的。不过，edgeR 也有一些解决方法，可以通过手动设置样品之间的标准差（SD值）来进行处理。但是 DESeq2 就非常灵活，它可以支持有无样品重复的情况。

把序列比对中的两个表格合并就得到了input.txt

不过在处理数据时，我们还需要注意，如果某个基因在两个样品中的表达量都为零，我们应该去掉这样的数据。根据我们的经验，一般来说，当每个样品中的基因表达量小于5时，我们会去除这些基因，因为它们在两个样本中都没有表达。这是因为我们使用的是测试数据，数据集只包含了几万条基因，因此会得到一些表达量为零的结果。如果测序的深度达到一定要求，表达量为零的基因会非常少。

因此，我们会去掉在两个样品中都为零的基因，之后再进行分析。关于如何去除这部分数据，这里就不再详细讲解了，大家可以使用任何方法来去除表达量都很小的基因。我们将使用以前的数据进行处理，教大家如何构建一个矩阵，并展示如何用输入文件进行操作。

首先是安装Deseq：

install.packages("BiocManager")

BiocManager::install("DESeq2")

（下面代码只是用来作解释，是老代码，不要使用）

在做差异表达分析时，首先我们需要查看每个组内的标准差（SD）。计算完组内的标准差之后，我们再来比较组间的差异是否在组内允许的标准差范围之内。如果组间的差异在组内允许的标准差之内，那么p值会较大，意味着差异性较小，显著性较低。如果组间的差异超出了组内的标准差范围，p值会较小，意味着差异性较大，显著性也会提高。因此，计算标准差（SD）本质上就是在衡量组内变异与组间差异的大小关系。

上面代码只是用来做解释：

真正运行的时候用这个：

# 设置工作目录

setwd("C:/Users/admin/Desktop/生物/复试资料/生信复试内容/转录组/差异表达/DESeq")

# 加载 DESeq2 包

library("DESeq2")

# 读取输入数据

data = read.table("input.txt", row.names=1 ,header=T, check.names=F)

# 创建设计因子（根据实验设计选择对照组与实验组）

design = factor( c( rep("con", 2), rep("treat", 2)) )

# 使用 DESeq2 创建数据集

dds = DESeqDataSetFromMatrix(countData = data, colData = data.frame(condition = design), design = ~ condition)

# 归一化数据并进行差异分析

dds = DESeq(dds)

# 提取标准化的计数数据

norm_counts = counts(dds, normalized = TRUE)

# 绘制箱型图 - 原始数据

boxplot(data, col = "blue", xaxt = "n", outline = F)

# 绘制箱型图 - 标准化数据

boxplot(norm_counts, col = "red", xaxt = "n", outline = F)

# 估计离散度（DESeq2会自动完成此步骤，不需要手动调用）

# 你不需要再次运行 estimateDispersions，DESeq() 已经处理了

# 执行差异分析

res = results(dds, contrast = c("condition", "con", "treat"))

# 检查结果是否成功生成

if (is.null(res)) {

  stop("差异分析结果为空，请检查DESeq()是否正确执行。")

}

# 绘制 p 值分布直方图

hist(res$padj, breaks = 100, col = "skyblue", border = "slateblue", main = "P-value distribution")

# 输出差异分析结果

write.table(res, file = "DESeqOut.xls", sep = "\t", quote = F, row.names = F)

# 筛选显著差异基因

resSig = res[!is.na(res$padj),]

resSig = resSig[(resSig$padj < 0.05 & (resSig$log2FoldChange > 1 | resSig$log2FoldChange < -1)),]

resSig = resSig[order(resSig$padj),]

# 输出显著差异基因

write.table(resSig, file = "DESeqSig.xls", sep = "\t", quote = F, row.names = F)

# 热图

hmExp = log10(norm_counts[rownames(resSig),] + 1e-5)  # 加上一个小常数避免log(0)

library('gplots')

hmMat = as.matrix(hmExp)

# 保存热图为 PDF

pdf(file = "heatmap.pdf")

par(oma = c(3, 3, 3, 5))  # 调整外边距

heatmap.2(hmMat, col = 'greenred', trace = "none", cexCol = 1)  # 绘制热图

dev.off()

# 火山图

pdf(file = "volcano.pdf")

allDiff = res[!is.na(res$padj),]

xMax = max(-log10(allDiff$padj)) + 1

yMax = 10

plot(-log10(allDiff$padj), allDiff$log2FoldChange, xlab = "-log10(padj)", ylab = "log2FoldChange", main = "Volcano", xlim = c(0, xMax), ylim = c(-yMax, yMax), yaxs = "i", pch = 20, cex = 0.4)

diffSub = subset(allDiff, allDiff$padj < 0.05 & allDiff$log2FoldChange > 1)

points(-log10(diffSub$padj), diffSub$log2FoldChange, pch = 20, col = "red", cex = 0.4)

diffSub = subset(allDiff, allDiff$padj < 0.05 & allDiff$log2FoldChange < -1)

points(-log10(diffSub$padj), diffSub$log2FoldChange, pch = 20, col = "green", cex = 0.4)

abline(h = 0, lty = 2, lwd = 3)  # 添加水平虚线

dev.off()

结果分析：

原始箱线图：

标准箱线图：

可以看出校正确实有一定的效果的。

差异表达结果：

热图：

火山图：