deseq与无重复的差异表达分析

下面是完整的 R 脚本，包含 RNA-seq 数据分析和结果可视化的步骤：

# 加载 DESeq 库
library(DESeq)

# 读取计数数据
countdata <- read.table("1516_1_all_gene_count.tsv", row.names = 1, header = TRUE, check.names = TRUE)

# 设置样本名称
#读取到的列名前有X
# 判断列名前是否有 X，并替换掉
names(countdata) <- ifelse(grepl("^X", names(countdata)),
                           gsub("^X", "", names(countdata)),
                           names(countdata))
#names(countdata) <- c("1516.714W.1", "1516.714T.2")

# 创建实验设计表
DesignTable <- 	data.frame(condition=c("W", "T"),
				libType = c("paired-end", "paired-end"),
				row.names=names(countdata))


# 定义条件
conds <- as.factor(DesignTable$condition)

# 创建 CountDataSet 对象
cds <- newCountDataSet(countdata, conds)

# 估算样本大小因子
cds <- estimateSizeFactors(cds)

# 估算离散度
cds <- estimateDispersions(cds, method = "blind", sharingMode = "fit-only", fitType = "local")

# 进行差异表达分析
res <- nbinomTest(cds, "W", "T")

# 调整 p 值
res$padj <- p.adjust(res$pval, method = "BH")

# 统计显著性结果
cat("Number of significant genes (pval < 0.05):", length(which(res$pval < 0.05)), "\n")
cat("Number of significant genes (padj < 0.05):", length(which(res$padj < 0.05)), "\n")

# 保存差异表达结果
write.table(res, file = "1516_1.T2W.DE.res.txt", quote = FALSE, sep = "\t", col.names = TRUE, row.names = FALSE)

# 加载必要的函数脚本
source("~/yiyaoran/rna-seq/my_rna_seq/scripts/deg_plot_funcs.r")

# 读取 DE 结果数据
#res <- read.table("~/yiyaoran/rna-seq/my_rna_seq/1516-ran/4.expression/1516_1/1516_1.T2W.DE.res.txt", 
#                   row.names = 1, header = TRUE, check.names = FALSE)

# 替换无限值为 NA
res$log2FoldChange[is.infinite(res$log2FoldChange)] <- NA
res$padj[is.infinite(res$padj)] <- NA

# 查看缺失值数量
cat("Number of NA in log2FoldChange:", sum(is.na(res$log2FoldChange)), "\n")
cat("Number of NA in padj:", sum(is.na(res$padj)), "\n")

# 移除包含 NA 的行
res <- na.omit(res)

# 统计上调、下调和正常调节基因的数量
nup <- nrow(res[(res$padj < 0.05) & (res$log2FoldChange > 1), ])
cat("Number of up-regulated genes:", nup, "\n")

ndown <- nrow(res[(res$padj < 0.05) & (res$log2FoldChange < -1), ])
cat("Number of down-regulated genes:", ndown, "\n")

nnormal <- nrow(res[(res$padj >= 0.05) | (res$log2FoldChange <= 1) & (res$log2FoldChange >= -1), ])
cat("Number of normal-regulated genes:", nnormal, "\n")

# 定义显著性标记
significant <- rep("Normal", nrow(res))  # 初始化、初始标记为 "Normal"
significant[(res$padj < 0.05) & (res$log2FoldChange > 1)] <- "Up"
significant[(res$padj < 0.05) & (res$log2FoldChange < -1)] <- "Down"
significant[(res$padj >= 0.05) | ((res$log2FoldChange <= 1) & (res$log2FoldChange >= -1))] <- "Normal"

# 调用火山图函数
pdf("ma_vo_plot.pdf")
ma_vo <- plot_MA_Volcano(ma_file = "1516_1.FC_count.png",
                         vo_file = "1516_1.FC_FDR.png",
                         log10exp = log10(res$baseMean),
                         log2FC = res$log2FoldChange,
                         FDR = res$padj,
                         Significant = significant,
                         yline = -log10(0.05),
                         xline = c(1, -1))
pdf("ma_vo_plot.pdf")
write.table(res, file = "6463_3.DE.res.txt", quote = FALSE, sep = "\t", col.names = TRUE, row.names = FALSE)
# 显示结果
print(ma_vo)

脚本说明

1. 加载 DESeq 库:

   • library(DESeq)：加载 DESeq 包。

2. 读取计数数据:

   • read.table()：从文件中读取计数数据，并设置样本名称和实验设计。

3. 创建 CountDataSet 对象并估算参数:

   • newCountDataSet()：创建 DESeq 数据集对象。
   • estimateSizeFactors()：估算样本大小因子。
   • estimateDispersions()：估算离散度。

4. 差异表达分析:

   • nbinomTest()：进行负二项检验以找出差异表达的基因。
   • p.adjust()：调整 p 值。

5. 统计显著性结果:

   • 计算显著性基因的数量，并将结果保存到文件。

6. 读取和处理 DE 结果数据:

   • 读取 DE 结果数据，替换无限值，移除缺失值，统计基因调节状态。

7. 绘制火山图:

   • 调用自定义的 plot_MA_Volcano 函数生成 MA 图和火山图。

使用说明

• 请确保文件路径、文件名和函数定义与实际情况相符。
• 根据需要调整输出文件名和参数设置。

nbinomTes函数

nbinomTest 是 DESeq 包中的一个函数，主要用于检测两个条件（如 RNA-Seq 实验中的两种处理条件）之间基因表达的差异。具体来说，它可以用于比较两种条件下基因表达的均值差异。

具体解释：

• 功能描述: 这个函数用于测试两种条件（如 A 和 B）下基因表达均值的差异，用于差异基因表达分析。
• 使用方法: nbinomTest(cds, condA, condB, pvals_only = FALSE, eps=NULL)
  • cds: 包含原始计数数据和归一化因子的对象（通常是 RNA-Seq 数据）。
  • condA 和 condB: cds 中的两个条件，例如条件 A 和条件 B。
  • pvals_only: 如果设置为 TRUE，只返回未调整的 p 值，而不是完整的数据框。

函数返回值：

这个函数返回一个数据框，其中包含以下信息：

• id: 基因或其他观察对象的 ID（通常是行名）。
• baseMean: 两个条件下计数的均值，经过归一化后的均值。
• baseMeanA: 条件 A 下计数的归一化均值。
• baseMeanB: 条件 B 下计数的归一化均值。
• foldChange: 条件 B 相对于条件 A 的折叠变化（baseMeanB/baseMeanA）。
• log2FoldChange: 折叠变化的对数值（log2 值）。
• pval: 原假设为两个条件的均值相等（meanA == meanB）的 p 值。
• padj: 使用 Benjamini-Hochberg 方法校正的 p 值。

示例代码:

cds <- makeExampleCountDataSet()
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
head(nbinomTest(cds, "A", "B"))

这个示例展示了如何使用 nbinomTest 函数来进行差异表达分析。在分析过程中，首先生成了一个示例数据集，接着对数据进行了归一化和方差估计，最后比较了 A 和 B 两个条件下的基因表达差异。

estimateDispersions 函数:

estimateDispersions 是 DESeq 包中的一个函数，用于估计 RNA-seq 数据中基因的离散度（dispersion）。离散度描述了基因表达量的变异性，它对于后续的差异表达分析非常重要。

函数的核心步骤和参数：

1. 目的：

   • 该函数用于估计每个基因的离散度，以便更好地拟合模型并进行差异表达分析。

2. 主要参数：

   • object: 一个包含计数数据的 CountDataSet 对象，它已经通过 estimateSizeFactors 计算出样本的大小因子。

   • method: 定义如何计算每个基因的经验离散度，有以下几种选择：

     • pooled: 使用所有有生物重复的样本来估计一个单一的经验离散度，并将其分配给所有样本。
     • pooled-CR: 使用 Cox-Reid 调整的轮廓似然进行估计，适合复杂的实验设计。
     • per-condition: 对每个条件单独估计离散度，适合有生物重复的情况。
     • blind: 忽略样本标签，将所有样本视为同一个条件的重复样本来估计离散度，适合没有生物重复的情况，但可能会导致统计功效降低。

   • sharingMode: 决定如何使用基因的经验值和拟合值来确定最终的离散度估计值：

     • fit-only: 仅使用拟合值，适合没有太多生物重复的情况。
     • maximum: 使用经验值和拟合值中的最大值，适合有 3-4 个重复的情况。
     • gene-est-only: 仅使用经验值，适合有足够生物重复的情况。

   • fitType: 指定如何拟合离散度与均值的关系：

     • parametric: 使用参数模型拟合，通常是更为稳健的选择。
     • local: 使用局部回归进行拟合。

不同方法的区别：

• method="blind"：忽略生物学条件的区别，适用于没有生物重复的情况，但可能导致功效损失。
• method="pooled"：基于所有有重复的条件进行估计，适合有生物重复的实验。
• fitType="parametric"：参数模型通常更为稳健，适用于大多数情况。
• fitType="local"：使用局部回归，适合某些特殊的表达模式。

示例：

cds <- makeExampleCountDataSet()  # 创建示例数据集
cds <- estimateSizeFactors(cds)   # 估计大小因子
cds <- estimateDispersions(cds)   # 估计离散度

这个过程用于调整和估计每个基因的变异性，以便后续的差异表达分析变得更加可靠。