基础生信分析——deseq2包

适合入门级生信学习，以及湿实验补充生信分析提分

---

---

前言

测序数据的表达矩阵通常包含了基因组数据分析的关键信息。一般来说，表达矩阵中的数据包括：

基因表达水平：
表达矩阵中的每一行代表一个基因，每一列代表一个样本。矩阵中的每个元素通常表示相应基因在相应样本中的表达水平，可以是原始计数值或标准化后的表达值（如FPKM、TPM等）。

样本信息：
表达矩阵中通常还包含了关于每个样本的信息，比如样本ID、类型（如对照组、实验组）、处理条件等。这些信息有助于进一步的数据分析和解释。

注释信息： 表达矩阵可能还包含了基因的注释信息，比如基因的ID、基因名、功能等。这些信息可以帮助用户理解和解释分析结果。

质控信息： 有时候表达矩阵中还会包含质控信息，比如基因的丰度、缺失值处理等。质控信息有助于评估数据质量并进行相应的数据预处理。

这些数据一起构成了表达矩阵，可以用来进行基因表达分析、差异表达分析、聚类分析等。在进行数据处理和分析时，通常需要结合这些数据来解释和理解实验结果。

---

一、什么是差异分析

差异分析就是从基因组学上对不同分组的样本筛选不同表达水平的基因

二、简单步骤

1.环境

install.packages('DESeq2', force = TRUE)
install.packages('dplyr', force = TRUE)
install.packages('ggplot2', force = TRUE)
install.packages('ggpubr', force = TRUE)
install.packages('ggrepel', force = TRUE)
library(DESeq2)#read count数据即可
library(dplyr)
library(ggplot2)
library(ggpubr)
library(ggrepel)
rm(list=ls(all=TRUE))
options(stringsAsFactors = F)

2.整理矩阵数据

mRNA_count <- read.csv("matrix.csv")
mRNA_count <- aggregate( . ~ symbol,data=mRNA_count, max) #重复基因取最大表达量项
rownames(mRNA_count) <- mRNA_count[, 1]#将第一列向量--基因名，设置为矩阵列名
mRNA_count <- mRNA_count[, -1]#删除第一列

expr <- mRNA_count#转换成操作数据，保留原始数据
expr <- ceiling(expr)#取整数
sapply(expr, class)#查看矩阵向量分布
is.na(expr)#确认无缺失值

expr <- expr[rowMeans(expr) > 1, ]  # 保留行均值大于1的行（排除不重要的极值基因）
write.csv(expr,'matrix.csv',row.names = TRUE)#重新写入整理后的矩阵文件

3.创建分组数据

table(substr(colnames(expr),14,16))#查看14到16位编号确认肿瘤样本和正常样本数量
Tumor <- grep('01A',colnames(expr))
Tumor#查看肿瘤样本所处位置
Normal <- grep('11A',colnames(expr))
Normal#查看正常样本所处位置

Tumor_matrix <- expr[Tumor,]#提取肿瘤样本组矩阵
Normal_matrix <- expr[Normal,]#提取正常样本组矩阵
expr <- rbind(Tumor_matrix,Normal_matrix)#合并矩阵

group <- factor(c(rep("Tumor",times=length(Tumor)),rep("Normal",times=length(Normal))))#创建分组（因子变量）
Data <- data.frame(row.names = colnames(expr), group = group)#创建分组数据框
Data#查看确认
write.csv(Data,"group.csv")#写入

4.开始进行差异表达分析

#第一步：构建DEseq2对象(dds)
dds <- DESeqDataSetFromMatrix(countData = expr,
                              colData = Data,
                              design = ~ group)
#第二步：开始差异分析
dds2 <- DESeq(dds)
res <- results(dds2, contrast=c("group", "Tumor", "Normal"))#后者为对照组
res <- res[order(res$pvalue),]#按P值从小到大排序
summary(res)#检查

my_result <- as.data.frame(res)#转成容易查看的数据框
my_result <- na.omit(my_result)#删除倍数为0的值

#第三步：保存差异分析的结果
my_result$Gene_symbol<-rownames(my_result)
my_result <- my_result %>% dplyr::select('Gene_symbol',colnames(my_result)[1:dim(my_result)[2]-1],everything())
rownames(my_result) <- NULL
write.csv(my_result,file="my_result_deseq2.csv")#写入

my_result <- read.csv("my_result_deseq2.csv",row.names = 1)#读入

5.筛选DEGs（差异基因）

#定义差异基因，参数可根据需求调整
my_result$regulate <- ifelse(my_result$padj > 0.05, "unchanged",
                             ifelse(my_result$log2FoldChange > 2, "up-regulated",
                                    ifelse(my_result$log2FoldChange < -2, "down-regulated", "unchanged")))
table(my_result$regulate)#查看DEGs数量
#可以把上调基因和下调基因取出放在一块
DEG_deseq2 <-subset(my_result, padj < 0.05 & abs(log2FoldChange) > 2) 
upgene <- DEG_deseq2[DEG_deseq2$regulate=='up-regulated',]
downgene <- DEG_deseq2[DEG_deseq2$regulate=='down-regulated',]

write.csv(DEG_deseq2,file= "DEG_deseq2.csv")#写入
DEG_deseq2 <- read.csv("DEG_deseq2.csv",row.names = 1)#读入

7.火山图无标签

plot(DEG_deseq2$log2FoldChange,-log2(DEG_deseq2$padj))#基础函数查看火山图
#ggplot2美化
p <- ggplot(data=DEG_deseq2, aes(x=log2FoldChange, y=-log10(padj),color=regulate)) + 
  geom_point(shape = 16, size=2) + 
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2 fold change") + #X轴标题
  ylab("-log10 p-value") + #Y轴标题
  theme(plot.title = element_text(size=15,hjust = 2.5)) + 
  theme_classic()+
  scale_colour_manual(values = c('#86D47D','#DFE0DF','#C34B99'))+#颜色自定义
  geom_vline(xintercept = c(-2,2),lty=4,col ="gray",lwd=0.8)+ #FC分界线
  geom_hline(yintercept=-log10(0.05),lty=2,col = "gray",lwd=0.6)+#P值分界线
  annotate("text",label="886 genes are up-regulated",x = 5.5,y = 30)+ #可以去掉
  annotate("text",label="623 genes are down-regulated",x = 5.7,y = 29)+ #可以去掉
  labs(title='TCGA-ESCA')+#标题
  annotate("text",x=upgene$log2FoldChange[1:3],y=(-log10(upgene$padj[1:3])),label=upgene$Gene_symbol[1:3], size=5.0)+
  annotate("text",x=downgene$log2FoldChange[1:3],y=(-log10(downgene$padj[1:3])),label=downgene$Gene_symbol[1:3], size=5.0)
plot(p)#出图

7.火山图带标签

DEG_deseq2$log10padj <- -log10(DEG_deseq2$padj)#生成新的一列v
non_zero_values <- DEG_deseq2$v[DEG_deseq2$v != Inf]# 提取padj列中所有非零值
DEG_deseq2$v[DEG_deseq2$v == 0] <- sample(non_zero_values, sum(DEG_deseq2$v == 0), replace = TRUE)# 替换padj列中的0值

#挑选要展示的基因（2种方式）
my_select <- subset(DEG_deseq2, padj < 0.0001 & abs(log2FoldChange) > 2)#高要求
my_select <- my_select$Gene_symbol#提取向量
#DIY
my_select <- c("ACHE","CDS1","PLA2G12A","LPCAT4","PEMT",'PLD1','PTDSS2','ARSA','GALC','SPHK1','UGT8',"CERS1",'PLA2G10')
#绘图
library(ggrepel)
ggscatter(DEG_deseq2,
          x = "log2FoldChange",
          y = "log10padj",
          ylab = "-log10(adjust p-value)",
          size = 2,
          color = "regulate",
          palette = c('#86D47D', '#C34B99')) +
  geom_text_repel(data = subset(DEG_deseq2, Gene_symbol %in% my_select),
                  aes(label = Gene_symbol),
                  color = "black",
                  box.padding = 0.5,
                  point.padding = 1,
                  segment.color = "black",
                  show.legend = FALSE, max.overlaps = 3000  # 增加 max.overlaps 的值
  )

---

总结

差异分析是生信分析的起点，思维可以不用局限在不同分组的样本，也可以是单细胞测序中的不同细胞簇