基础生信分析——deseq2包

大家好，欢迎来到IT知识分享网。

基础生信分析——deseq2包

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

前言

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

注释信息： 表达矩阵可能还包含了基因的注释信息，比如基因的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 的值 ) 

总结

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