零基础入门转录组数据分析——limma差异分析

零基础入门转录组数据分析——limma差异分析

---

---

1. 转录组分析基础知识

1.1 什么是转录组？
转录组（transcriptome） 是指在某一生理条件下，细胞内所有转录产物的集合，这些产物主要包括信使RNA（mRNA）、核糖体RNA（rRNA）、转运RNA（tRNA）以及非编码RNA（ncRNA）。其中，编码RNA可以通过转录和翻译过程产生蛋白质，从而直接调控生命活动；而非编码RNA具有调节基因表达的作用，它们可以对染色体结构、RNA加工修饰及稳定性、转录和翻译、甚至蛋白质的转运及稳定性产生影响。
（简单来说：转录组是指一个活细胞在特定时间和环境下，所能转录出来的所有RNA的总和。因此，转录组分析就是为了了解基因的表达调控机制和功能。）

1.2 差异分析是什么？为什么要做差异分析？
转录组差异分析是一种生物信息学分析方法，旨在比较不同样本或条件下转录组数据的差异，以揭示基因表达的变化。
（简单来说：就是从上w个基因中找出不同样本中表达具有差异的基因，从而了解这些差异基因如何影响生物体的生理功能和疾病发生发展。）

1.3 在R中对于转录组数据都有哪些差异分析的方法？
包括但不局限于DESeq2、limma和edgeR等。
DESeq2： 基于负二项分布模型进行差异分析，考虑了测序深度和基因长度对count数据的影响，因此在处理高度变异的基因时表现较好。DESeq2还提供了多种标准化方法，可以帮助减少批次效应和其他技术噪声。
limma： 最初是为微阵列数据设计的，后来也扩展到RNA-Seq数据分析。它使用线性模型进行差异分析，并且可以轻松处理多个因素和协变量。limma的优势在于其速度和灵活性，可以处理大量的基因和样本。
edgeR： 也是基于负二项分布模型进行差异分析，但与DESeq2不同的是，它使用了经验贝叶斯方法来估计分散参数，从而提高了差异检验的准确性和稳定性。edgeR在处理大型转录组数据集时表现出色，尤其是在识别低表达水平的差异基因方面。

1.4 三种方法对于输入数据的要求？
DESeq2： 需要原始的count矩阵（如通过htseq-count工具获得），并且只支持重复样品。
limma： 可以接受原始的count矩阵，但需要用户自行进行标准化（通常是log转换），并且也支持重复样品。
edgeR： 要求输入原始的count矩阵，既支持单个样品也支持重复样品。

1.5 在做转录组分析的时候该选用哪种方法？（关键）
GEO芯片数据——用limma差异分析，因为GEO芯片数据底层已经对数据做了标准化处理。
（数据参考：零基础入门转录组分析——数据处理（GEO数据库——芯片数据））

GEO高通量数据——用DESeq2差异分析，因为可以获取到count数据。
（数据参考：零基础入门转录组分析——数据处理（GEO数据库——高通量测序数据））

TCGA_count数据——用DESeq2差异分析，因为可以获取到count数据。
（数据参考：零基础入门转录组分析——数据处理（TCGA数据库））

自测序数据——用DESeq2差异分析，因为可以获取到count数据。
（数据参考：零基础入门转录组分析——数据处理（自测序数据））

注：目前了解到的很少用到edgeR，一般来说有重复的样本通常limma和DESeq2就能解决（edgeR这种方法不做过多介绍）。

---

---

2. limma差异分析（Rstudio）

本项目以GSE103552数据集(芯片数据)展示limma差异分析过程
实验分组：疾病组（21例），对照组（16例）
R版本：4.2.2
R包：tidyverse，lance，limma

数据处理过程参考之前的教程：零基础入门转录组分析——数据处理（GEO数据库——芯片数据）

设置工作空间：

rm(list = ls()) # 删除工作空间中所有的对象
setwd('/XX/XX/XX') # 设置工作路径
if(!dir.exists('./01_DEGs')){
  dir.create('./01_DEGs')
} # 判断该工作路径下是否存在名为00_rawdata的文件夹，如果不存在则创建，如果存在则pass
setwd('./01_DEGs/') # 设置路径到刚才新建的00_rawdata下

加载包：

library(tidyverse)
library(lance)
library(limma)

导入处理好的表达矩阵：

# 表达矩阵
train_data <- read.csv('../../00_rawdata/GSE103552/dat.GSE103552.csv', row.names = 1) %>% lc.tableToNum()

train_data 如下图所示行名为基因symbol，列名为样本ID，表中数字是经过标准化的基因表达量

判断标准化数据的两个简要特征：

• 数据类型是小数，不是整数
• 表达量数值在20以内（比较小）

为啥说是在20以内，因为在前面处理过程中提到不做标准化处理的数据一般是呈现偏态分布，而用log2处理后数据就会基本处于正态分布，打个比方说有个基因的比对reads是1048576，换算成log2(1048576) = 20，所以说标准化后的数值一般都会小，毕竟是2的对数。

导入处理好的分组信息表：

# 分组信息
group <- read.csv('../../00_rawdata/GSE103552/group.GSE103552.csv')
table(group$group)
group$group <- factor(group$group, levels = c("disease","control")) 

group 如下图所示，第一列为样本ID，第二列为样本分组

先构建分组矩阵design

# 分组矩阵design
design <- cbind(control = ifelse(group$group == "control", 1, 0), 
                disease = ifelse(group$group == "disease", 1, 0))

design 如下图所示，行名中1，2，3…与之前的group分组信息表中行名数字是对应的

注：control组被编码为0，disease组被编码为1。这意味着该对比方式将以control组为参考（基准）

之后构建对比矩阵contrast.matrix，makecontrasts函数是R语言中的线性假设检验函数。它提供了一种方便的方式，可以生成所需的对比矩阵，这里就可以用到上一步构建的分组矩阵design。

contrast.matrix <- makeContrasts(contrasts = 'disease-control', levels = design)

contrast.matrix 如下图所示：control为-1，disease为1。

注：-1和1的意思是control是用来被比的，disease是用来比的

构建拟合模型并且进行差异分析：（每一行代码后面都有注释，解释代码是干什么的）

fit <- lmFit(train_data, design) # 构建拟合模型，用来描述train_data中的数据是如何与design中的变量相关联的
fit2 <- contrasts.fit(fit, contrast.matrix) # 拟合模型的对比。对比方式是之前构建的对比矩阵：contrast.matrix（即：比较拟合模型中不同分组条件下基因表达是否存在差异）
fit2 <- eBayes(fit2 , trend = TRUE) # 是一种贝叶斯统计方法，为了让模型的参数估计更加准确，使它们更接近真实值

将差异分析的结果转换成数据框：

tempOutput <- topTable(fit2, coef = 1, n = nrow(fit2), adjust="fdr")  # 将差异结果转换成数据框
tempOutput <- na.omit(tempOutput) ## 去掉数据中有NA的行或列

tempOutput 如下图所示：

• 行名——基因名
• logFC——代表基因表达的对数2倍变化（log2 fold change），它表示在两组样本之间基因表达水平的变化量。一个正值表示基因在某一组样本中表达上调，而负值表示表达下调。
• AveExpr——每个基因的平均表达量。
• t值——t值是从t检验中得到的统计量，用于衡量基因表达变化的显著性。一个大的正值或负值表示基因表达的变化是否显著。
• P.Value——P值是从假设检验中得到的，用于评估观察到的效应（如基因表达的变化）是否可能是由于随机误差引起的，通常认为P < 0.05则结果是显著的。
• adj.P.Val——是经过某种调整后的P值，用于控制假阳性率（简单说就是矫正后的P值）。
• B值——是从贝叶斯方法中得到的后验对数比率，它是对基因表达变化的另一种度量，通常用于对基因进行排序或可视化。
  

给差异分析结果添加上下调标签

DEG <- tempOutput
DEG$change <- as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) >= 0.5, 
                               ifelse(DEG$logFC > 0.5, 'Up', 'Down'), 'Not'))
DEG_write <- cbind(symbol = rownames(DEG), DEG)

处理好的 DEG_write 如下图所示，除了新添加的symbol和change，其余与tempOutput一样，change列中: （筛选条件：|logFC| > 0.5 & p.value < 0.05）

• Up 表示上调基因
• Down 表示下调基因
• Not 表示未通过筛选的基因

筛选出表达具有差异的基因

sig_diff <- subset(DEG, DEG$P.Value < 0.05 & abs(DEG$logFC) >= 0.5)
sig_diff_write <- cbind(symbol = rownames(sig_diff), sig_diff)

sig_diff_write 与DEG_write的差别就是少了那些没有差异的基因。

最后，保存差异分析的结果：

write.csv(DEG_write, file = paste0('DEG_all.csv'))
write.csv(sig_diff_write, file = paste0('DEG_sig.csv'))

---

---

3. 结语

以上就是limma差异分析的所有过程，如果有什么需要补充或不懂的地方，大家可以私聊我或者在下方评论。

如果觉得本教程对你有所帮助，点赞关注不迷路！！！

---

注意：有的小伙伴lance包下不下来，可以考虑直接删除这个包，只需要把管道符%>%及其后面的lc.tableToNum()删除即可

%>% lc.tableToNum()

这个包之前是用来将数据框里的内容转成数值型的，但是看很多人反馈说是下载不了，那就直接不要了（但是要注意导入数据后看看是否是数值型）

---

---

与教程配套的原始数据+代码+处理好的数据见配套资源

注：配套资源只要改个路径就能运行，本人已检测过可以跑通，请放心食用，食用过程遇到问题，可先自行百度，实在解决不了可以私信

---

• 目录部分跳转链接：零基础入门生信数据分析——导读