Bulk RNA Seq Analysis

本文档使用DESeq2包对组织的转录组测序进行下游分析，从raw counts数据到差异表达基因（Differentially Expressed Genes, DEGs），绘制火山图，再进一步做通路富集分析。

数据读取和清洗

library(dplyr)
library(stringr)
library(ggplot2)
library(DESeq2)

读取counts矩阵，注意此处不能是reads，也不能是基因表达矩阵。 清洗数据，使每行为一个基因，每列为一个样本，并设置行名为基因名，列名为样本名。

f <- file.path('_')
countData <- readxl::read_xlsx(f)
countData

设置分组信息。是一个只有一列的数据框，行名为样本名。此处演示前8个样本为实验组，后10个样本为对照组。

coldata <- data.frame(
  row.names = colnames(countData),
  condition = c(rep('treatment', 8), rep('control', 10))
)

DESeq2流程

# 构建 DESeq2 数据对象
dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = coldata,
                              design = ~ condition)

# 过滤低表达基因，此步骤可有可无，对结果几乎无影响。
dds <- dds[rowSums(counts(dds)) > 10, ]

# 差异表达分析
dds <- DESeq(dds)

# 获取结果。
# 注意此处contrast参数的顺序，表示分组信息储存在condition中，计算treatment相对于 control的变化，也就是 > 0 表示在实验组上调。
res <- results(dds, contrast = c("condition", "treatment", "control"))

head(res)

明确上调和下调的判断标准，一般log2FoldChange设置为1.5或1.2或1。

# calculate statistical difference ----
res_sig <- res |>
  tibble::rownames_to_column(var = 'gene') |> 
  dplyr::mutate(signif = case_when(
    log2FoldChange > 1 & pvalue < 0.05 ~ 'up',
    log2FoldChange < -1 & pvalue < 0.05 ~ 'down',
    .default = 'ns'
  ))
count(res_sig, signif)

# 保存数据
readr::write_csv(res_sig, file = '_.csv')