【RNAseq】差异分析

---

前言

对于二代测序的count值（也就是没有标准化后的数据）通常有三个包可以进行差异分析：

• DESeq2
• edgeR
• limma

下面是对整理好的表达矩阵进行下游分析，不是从上游分析开始

---

一.环境设置

代码如下（示例）：

Sys.setenv(language = "en") # 英文环境
options(stringsAsFactors = F) # 全局设置，默认不转化为因子
rm(list = ls()) # 清空环境变量

二.加载R包

代码如下（示例）：

library(GEOquery)
library(GenomicRanges)
library(GenomeInfoDb)
library(DESeq2)
library(dplyr)
library(tibble)
library(edgeR)
library(limma)

三、分析

1、DESeq2

DEseq2针对有生物学重复的样本

代码如下（示例）：

# 读取数据
# 表达矩阵
cds <- read.csv("data/airway_scaledcounts.csv",row.names = "ensgene") %>% 
  as.matrix() 
# 分组信息
coldata <- read.csv("data/airway_metadata.csv",row.names = "id") 
coldata$dex <- factor(coldata$dex,levels = c("control","treated")) #两个分组实验组比对照组 （因子水平要设置contro在前，否则默认按字母排序）
all(colnames(cds) %in% rownames(coldata))
if(!all(colnames(cds) == rownames(coldata))){  
  coldata <- coldata[colnames(cds),]  # 将表达矩阵和分组信息匹配
}

# 创建DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = cds,
                              colData = coldata,
                              design = ~ dex) # 如果是多个分组 design = ~ 0 + dex 

# 过滤
## 减少对象的内存，提高运行速度
keep <- rowSums(counts(dds)) > 1 #过滤行和1以下的
dds <- dds[keep,]

# 质控
## 原始表达矩阵的
## 样本间距离的热图
## 样本的PCA图

# 差异分析
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds,name = "dex_treated_vs_control")
degDeseq2 <- res %>% as.data.frame() %>% 
  na.omit()

2.edgeR

针对没有重复的

# 读取数据
cds <- read.csv("data/airway_scaledcounts.csv",row.names = 1) %>% 
  as.matrix()
View(cds)
coldata <- read.csv("data/airway_metadata.csv",row.names = 1)

# 匹配和分组
all(colnames(cds) == rownames(coldata))
group <- coldata$dex
group <- factor(group,levels = c("control","treated"))

# 创建dge对象
dgeData <- DGEList(counts = cds,
                   group = group)

# 过滤
keep <- filterByExpr(dgeData)
dgeData <- dgeData[keep,,keep.lib.sizes = FALSE]

# 标准化处理
dgeData <- calcNormFactors(dgeData)
design <- model.matrix(~ group)

dgeData <- estimateDisp(dgeData,design)

# 分析
# quasi-likelihood F-tests
# bulk RNA-seq data
fit <- glmQLFit(dgeData,design)
qlf <- glmQLFTest(fit,coef = 2)
degEdgeR <- topTags(qlf,n = Inf) %>% 
  as.data.frame()

3.limma-voom

还有一种针对芯片数据的方式

# 读取数据
cds <- read.csv("data/airway_scaledcounts.csv",row.names = 1) # 表达矩阵
coldata <- read.csv("data/airway_metadata.csv",row.names = 1) # 分组信息
if(!all(rownames(coldata) == colnames(cds))){ # 表达矩阵是否和分组信息一致
  coldata <- coldata[,colnames(cds)]
}

# 分组和设计矩阵
group <- factor(coldata$dex,levels = c("control","treated")) # 因子水平：控制组在前
design <- model.matrix(~ group) # 设计矩阵
colnames(design) <- c("control","controlVStreated")
rownames(design) <- rownames(coldata)

# 创建DGEList
dge <- DGEList(counts = cds) 

# 过滤
keep <- filterByExpr(dge,design) # 去除0低表达的行
dge <- dge[keep,,keep.lib.sizes = F]

# 标准化
dge <- calcNormFactors(dge)

# 差异分析
v <- voom(dge,design,plot = TRUE)
fit <- lmFit(v,design)
fit <- eBayes(fit)
degLimmaVoom <- topTable(fit,coef = 2,number = Inf)

---

总结

以上就是今天学习内容，主要学会了RNA-seq原始数据的三种分析方式，只有两个分组的方式，具体多个分组还是要参考官方R包。之前零零散散的学习方式，对于差异分析有了初步的认识，可能也是可以看懂R包说明文档的一个重要原因。
当然，英文的说明文档看不懂就得硬看
当然，英文的说明文档看不懂就得硬看
当然，英文的说明文档看不懂就得硬看

参考

• DESeq2
• EdgeR
• limma