转录组丨limma差异表达分析，绘制火山图和热图

limma差异表达分析

本篇笔记的内容是在R语言中利用limma包进行差异表达分析，主要针对转录组测序得到的基因表达数据进行下游分析，并将分析结果可视化，绘制火山图和热图

基因表达差异分析是我们做转录组最关键根本的一步，不管哪种差异分析，其本质都是广义线性模型，limma也是广义线性模型的一种，其对每个gene的表达量拟合一个线性方程。

limma包是2015年发表在Nucleic Acids Resarch一个做差异分析的工具，目前引用次数高达七千多次，最流行的差异分析软件之一就是limma。

环境部署与安装

• 安装limma包

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("limma")

• 安装ggVolcano包

# install.packages("devtools")
devtools::install_github("BioSenior/ggVolcano")

• 安装ggplot2包

#直接安装tidyverse，一劳永逸（推荐,数据分析大礼包）
install.packages("tidyverse")
#直接安装ggplot2
install.packages("ggplot2")
#从Github上安装最新的版本，先安装devtools(如果没安装的话)
devtools::install_github("tidyverse/ggplot2")

输入数据准备

• 样本信息表：第一列为样本名称ID，第二列为样本的分组（CK、HT），sampleinfo.csv格式如下

  这一步需要注意：每个分组下至少有两个生物学重复，即至少4行数据，如果没有重复的话在后续贝叶斯检验会出错，原因是该模型利用统计假设检验，多个重复能够评估变异性，而如果仅有一组数据，则无法检验。

  sample	group
  A01	CK
  A02	HT
  …	…

• 表达矩阵：每行是一个基因，每列是一个样本，表达数据为TPM，data.csv格式如下

  	A01	A02	…
  gene1	xx	xx	…
  gene2	xx	xx	…
  …	…	…	…

差异表达分析过程

准备环节

载入R包，设置参数，其中job变量用于项目输出文件前缀标识，可以自定义修改。

setwd("D:/LABdata")
options(stringsAsFactors = F)
rm(list=ls()) #清空变量
job <- "test" #设定项目名称
library(limma)
library(ggplot2) #用于绘制火山图
library(ggVolcano)

数据导入

导入样本信息和表达量数据，然后进行删除表达量之和为0基因、log2化、替换异常值等步骤，得到原始数据矩阵。

# 输入表达矩阵和分组文件 -------------------------------------------------------------

expr_data<-read.csv("data.csv",header = T