生信转录组

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

2023-02-22  本文已影响0人  生信分析笔记

limma差异表达分析

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

[TOC]


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

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

环境部署与安装

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

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

输入数据准备

差异表达分析过程

准备环节

载入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,row.names = 1,sep = ",") #输入文件TPM原始值,行名是基因,列名是样本
expr_data <- expr_data[which(rowSums(expr_data)!=0),] #删除表达量为0的基因
expr_data = log2(expr_data) #log化处理
expr_data[expr_data == -Inf] = 0 #将log化后的负无穷值替换为0
group<-read.csv("sampleinfo.csv",header = T,row.names = 1,sep = ",") #输入文件,样本信息表,包含分组信息

构建分组矩阵

根据样本的分组信息,构建分组矩阵,最终得到的design矩阵由0和1构成,为斜对角矩阵。

# 构建分组矩阵--design ---------------------------------------------------------

design <- model.matrix(~0+factor(group$group))
colnames(design) <- levels(factor(group$group))
rownames(design) <- colnames(expr_data)

构建比较矩阵

设置样本的比较方式,这里为CK对照比HT处理,该步骤生成的文件为1和-1构成的矩阵。

# #构建比较矩阵——contrast -------------------------------------------------------

contrast.matrix <- makeContrasts(CK-HT,levels = design) #根据实际的样本分组修改,这里对照组CK,处理组HT

线性混合模拟

该步骤是limma包的核心步骤,首先使用lmFit函数进行非线性最小二乘法分析,然后用经验贝叶斯调整t-test中方差部分,得到差异表达结果。

# #线性拟合模型构建 ---------------------------------------------------------------

fit <- lmFit(expr_data,design) #非线性最小二乘法
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)#用经验贝叶斯调整t-test中方差的部分
DEG <- topTable(fit2, coef = 1,n = Inf,sort.by="logFC")
DEG <- na.omit(DEG)

最终生成的DEG文件包含以下几列信息:

> colnames(DEG)
[1] "logFC"     "AveExpr"   "t"         "P.Value"   "adj.P.Val" "B"        
[7] "regulate"  "Genes"    

差异基因标注

本步骤中对P.ValuelogFC进行筛选,并对每个基因进行标注,显示其表达变化情况。

DEG$regulate <- ifelse(DEG$P.Value > 0.05, "unchanged",
        ifelse(DEG$logFC > 1, "up-regulated",
        ifelse(DEG$logFC < -1, "down-regulated", "unchanged")))

结果保存

生成两个文件,保存差异表达结果。

write.table(table(DEG$regulate),file = paste0(job,"_","DEG_result_1_005.txt"),
            sep = "\t",quote = F,row.names = T,col.names = T)
write.table(data.frame(gene_symbol=rownames(DEG),DEG),file = paste0(job,"_","DEG_result.txt"),
            sep = "\t",quote = F,row.names = F,col.names = T)

区分上下调基因

该步骤中DEG$P.Value<0.05&abs(DEG$logFC)>1为参数,筛选了差异倍数和显著性,并将基因按照上调和下调进行区分。

DE_1_0.05 <- DEG[DEG$P.Value<0.05&abs(DEG$logFC)>1,]
upGene_1_0.05 <- DE_1_0.05[DE_1_0.05$regulate == "up-regulated",]
downGene_1_0.05 <- DE_1_0.05[DE_1_0.05$regulate == "down-regulated",]
write.csv(upGene_1_0.05,paste0(job,"_","upGene_1_005.csv"))
write.csv(downGene_1_0.05,paste0(job,"_","downGene_1_005.csv"))

差异基因名称提取

提取差异倍数最大的1000个基因,按照顺序保存到txt文本中,只保留基因名称ID,用于后续验证。

tem1 <- head(rownames(upGene_1_0.05),1000) #以logFC差异倍数从大到小为序,提取前1000个基因名称
tem2 <- as.data.frame(tem1) # 转化数据类型为数据框
write.table(tem2$tem,file = paste0(job,"_","upgene_head100_name.txt"),
          row.names = FALSE,col.names = FALSE)

tem1 <- head(rownames(downGene_1_0.05),1000) #以logFC差异倍数从大到小为序,提取前1000个基因名称
tem2 <- as.data.frame(tem1) # 转化数据类型为数据框
write.table(tem2$tem,file = paste0(job,"_","downgene_head100_name.txt"),
          row.names = FALSE,col.names = FALSE)

自定义筛选阈值

之前的结果均为默认设置,如果你需要修改,仅需更改下面开头两行参数即可,运行后可以得到3个文件,分别是差异基因集、上下调过滤所得基因信息。

foldChange = 2 # 自定义修改筛选参数
padj = 0.05 # 自定义修改筛选参数
All_diffSig <- DEG[(DEG$adj.P.Val < padj & (DEG$logFC > foldChange | DEG$logFC < (-foldChange))),]
#dim(All_diffSig)
write.csv(All_diffSig, paste0(job,"_","all_diffsig_filtered.csv"))  ##输出差异基因数据集

### 自定义筛选上调和下调的基因 ===================================================================

diffup <-  All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig$logFC > foldChange)),]
write.csv(diffup, paste0(job,"_","diffup_filtered.csv"))
diffdown <- All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig$logFC < -foldChange)),]
write.csv(diffdown, paste0(job,"_","diffdown_filtered.csv"))

利用limma分析后,正常情况下应该会生成下面这些数据文件,可以用于验证过程中是否有问题。

image-20230222202113721

结果可视化

火山图绘制方法一

有了上面的DEG差异分析结果,根据gene、logFC、adj.P.Val三列信息即可绘制火山图,第一种方法使用ggvolcano,绘图代码如下:

# ggvolcano绘制火山图 ----------------------------------------------------------

pdf(paste0(job,"_","volcano1.pdf"),width = 10,height = 10)
ggvolcano(data = DEG,x = "logFC",y = "P.Value",output = FALSE,label = "Genes",
          fills = c("#00AFBB", "#999999", "#FC4E07"),
          colors = c("#00AFBB", "#999999", "#FC4E07"),
          x_lab = "log2FC",
          y_lab = "-Log10P.Value",
          legend_position = "UR") #标签位置为up right
dev.off()
image-20230222202352090

火山图绘制方法二

第二种方法使用ggplot2,得到另外一种形式的火山图,绘图代码如下:

# 火山图的绘制 ------------------------------------------------------------------

DEG$Genes <- rownames(DEG)
pdf(paste0(job,"_","volcano2.pdf"),width = 7,height = 7)
ggplot(DEG,aes(x=logFC,y=-log10(P.Value)))+ #x轴logFC,y轴adj.p.value
  geom_point(alpha=0.5,size=2,aes(color=regulate))+ #点的透明度,大小
  ylab("-log10(P.Value)")+ #y轴的说明
  scale_color_manual(values = c("blue", "grey", "red"))+ #点的颜色
  geom_vline(xintercept = c(-1,1),lty=4,col ="black",lwd=0.8)+ #logFC分界线
  geom_hline(yintercept=-log10(0.05),lty=4,col = "black",lwd=0.8)+ #adj.p.val分界线
  theme_bw()  #火山图绘制
dev.off()
image-20230222202415207

热图绘制

根据基因的表达变化信息,绘制热图并展示聚类树,详细代码如下:

# 热图的绘制 -------------------------------------------------------------------

DEG_genes <- DEG[DEG$P.Value<0.05&abs(DEG$logFC)>1,]
DEG_gene_expr <- expr_data[rownames(DEG_genes),]
#DEG_gene_expr[is.infinite(DEG_gene_expr)] = 0
#DEG_gene_expr[DEG_gene_expr == -Inf] = 0
pdf(paste0(job,"_","pheatmap.pdf"))
pheatmap(DEG_gene_expr,
         color = colorRampPalette(c("blue","white","red"))(100), #颜色
         scale = "row", #归一化的方式
         border_color = NA, #线的颜色
         fontsize = 10, #文字大小
         show_rownames = F)
dev.off()
image-20230222202505089
参考资料:https://zhuanlan.zhihu.com/p/437712423

本文由mdnice多平台发布

上一篇下一篇

猜你喜欢

热点阅读