差异分析

发表于 2022-05-18 更新于 2023-04-06 分类于生物信息学，富集分析阅读次数： Waline：

一、什么是差异分析

差异分析简单来说就是我们通过一组实验后，实验组中所提取的总mRNA 较对照组提取的总mRNA有显著性提高或者是降低的部分mRNA，寻找这部分有差异的mRNA就被称为差异分析

在常见的差异分析中，有几个比较常用的库下面来介绍一下

1.GEO数据库

GEO数据库可以理解为一个实验原始数据的储存库。

如图所示，这是一个研究某基因对酮酸脱氢酶活性是否有促进性的实验。

GSE号：是这组实验的编号，在GEO每个不同的实验编号唯一，用于检索。

GPL号：实验芯片平台，在GEO数据库的实验信息，往往是通过基因芯片高通量测序的方式测量转录组。简单来讲就是基因芯片中包含了这个实验测量了哪些基因的表达量

GSM号：是每一种样品处理方式的编号，一组实验中（一个GSE号）往往有多种处理（GSM号），例如，在研究pH压迫某微生物生长时，当你设置了5个梯度，则在一个GSE号中至少会有5个GSM号，即五种处理方式。如果设置了3个重复，则会有15个GSM号，GSM在GEO数据库中也是唯一的，可以通过GSM来唯一检索。

表达矩阵：表达矩阵中包含了每种处理下，样品被提取出的总mRNA的表达情况。简单来讲就是测量的“值域”，而平台信息就是测量的“定义域”

二、如何进行差异分析

1.下载原始文件

下载表达矩阵文件

下载实验芯片平台信息

在处理样本量比较小的情况下，可以手动记录一下，每个GSM号对应的处理方式

2.将原始数据导入你的R中并且进行格式化处理

分析过程中依赖的包：limma，stringr

library(limma)
library(stringr)


# 首先设置你的工作目录
setwd("E:/Project/R/差异分析") 
  
# 导入表达矩阵
expmatrix_data <- read.table("GSE115269_series_matrix.txt", header = TRUE, skip = 64)

# 导入平台文件 
pla_data <- read.table("GPL18233-21178.txt", sep = "\t", header = T)

在表达矩阵导入中，因为我用的表达矩阵从第65行开始才是矩阵数据，所以设置了skip。

# 格式化平台文件（将转录组编号列提出，只保留ENSMUST编号）
mRNA_assigment <- pla_data[, "mrna_assignment"]   

pattern <- "ENSMUST\\d*"

for (i in 1:length(mRNA_assigment)){ 
  # 用正则表达式循环替换单元格内的内容
  mRNA_assigment[i] <- str_match(string = mRNA_assigment[i], pattern = pattern)
}

pla_data[, "mrna_assignment"] <- mRNA_assigment

# 让平台文件以ID列重排序
pla_data <- pla_data[order(pla_data[, 'ID']), ]

格式化的目的是从上图到下图，

即让多个转录组编号转成一个，并且让平台文件按基因编号排序，这样处理是为了后续富集分析的方便。可能会出现转录组缺失的情况，因为多个编号转成了一个。

# 将样品处理方式导入R --手动 这里写样品的处理方式，不同实验有所不同
samp_treat <- c("DJ-1+/+_01", "DJ-1+/+_02", "DJ-1+/+_03", "DJ-1-/-_01", "DJ-1-/-_02", "DJ-1-/-_03")

# 以csv文件保存格式化后的平台注释文件
write.table(pla_data, "pla_data.csv", row.names = T, col.names = T, sep = ",")

因为R使用不熟练的原因，后续富集部分操作可能需要在excel中进行，故导出格式化后的平台注释文件。

# 格式化表达矩阵
rownames(expmatrix_data) <- expmatrix_data[, "ID_REF"]

expmatrix_data <- expmatrix_data[, -1]

colnames(expmatrix_data) <- samp_treat

expmatrix_data <- log2(expmatrix_data)

格式化目的是上图到下图

3.用limma包进行差异分析


# 根据样本处理方式设计design矩阵
group_list <- c(rep("Control", 3), rep("Treatment", 3)) 

design <- model.matrix(~0 + factor(group_list))

colnames(design) = levels(factor(group_list))

rownames(design) = colnames(expmatrix_data)

设置design矩阵的目的是为了对样本分组，但model.matrix函数的构造原理不是很清楚，这个需要进一步学习。总之就是构造出下图一样的分组。

这个design矩阵是含有level信息的。

# 模型拟合及分析
fit <- lmFit(expmatrix_data, design)

# 前面是分子(实验组)，后面是分母(对照组)
contrast.matrix <- makeContrasts("Treatment-Control", levels = design) 
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2) 
expDiff_data <- topTable(fit2, coef = 1, n = Inf) %>% na.omit()

# 保存差异分析结果文件
write.table(expDiff_data, "expDiff.csv", row.names = TRUE, col.names = TRUE, sep = ",")

分析后得到的expDiff_data如图所示