从GEO下载数据预处理+DESeq差异分析+KEGG、GO富集分析

最新推荐文章于 2024-09-15 12:13:20 发布

onethree07

最新推荐文章于 2024-09-15 12:13:20 发布

阅读量2w

点赞数 43

分类专栏：数据处理文章标签： r语言

本文链接：https://blog.csdn.net/dudasha/article/details/125799817

版权

数据处理专栏收录该内容

1 篇文章 8 订阅

订阅专栏

因为自己最近需要用GEO的数据来画火山图和富集分析图，就整理了一下操作流程。

一、从GEO中下载数据

我是用代码直接从GEO中下载数据，也可以自己手动去官网上下载，以GDS1906为例

workdir <- "E:5"
setwd(workdir)
getwd()

###预处理分为两布
# 1、找到探针与基因的对应关系
# 2、根据对应关系把探针转换为基因
# 转换时有两点需要注意：
# 可能会出现多个探针对应同一个基因的情况，这种情况一般取各个探针表达的中位数作为基因的表达。当然，视情况也可以取最大值、最小值、均值等。
# 可能会出现一个探针对应对多个基因，例

# BiocManager::install("GEOquery")
library(GEOquery)
eSet <- getGEO("GSE1906", destdir = ".", getGPL = T)   #eSet为一个集合
# geneNames/sampleNames/pData/exprs(这些都是对expression set 对象的操作函数)
exprSet = exprs(eSet[[1]])#表达量矩阵
fdata = fData(eSet[[1]])#需要getGPL = T,平台信息,注释文件包含探针与基因的对应关系
pdata = pData(eSet[[1]]) #表型信息，里面包含了需要的分组信息
samplenames <- sampleNames(eSet[[1]])
dim(exprSet)  
[1] 12488    12
dim(fdata)
[1] 12488    16
dim(pdata)
[1] 12 27
length(samplenames)
[1] 12

#找出ID对应的基因
genes <- fdata[,c("ID","Gene Symbol")]
exprSet <- data.frame(exprSet)
exprSet$ID <- rownames(exprSet)
exprSet <- merge(exprSet,genes,by.x = "ID",by.y = 'ID')
exprSet <- exprSet[,-1]
dim(exprSet)    #最后一列为基因
[1] 12488    13
#处理重复基因
rowMeans = apply(exprSet[,c(1:12)],1,function(x) mean(as.numeric(x), na.rm = T))####计算行平均值
rowMeans_2 = data.frame(rowMeans)
exprSet = exprSet[order(rowMeans, decreasing = T),] #排序
dim(exprSet)
[1] 12488    13
exprSet_2 = exprSet[!duplicated(exprSet[, dim(exprSet)[2]]),] #去掉重复基因
dim(exprSet_2)
[1] 9175   13
exprSet_na = na.omit(exprSet_2)   #去掉缺失值
explan_final = exprSet_na[exprSet_na$Gene Symbol != "",]
dim(explan_final)
[1] 9174   13

#处理一个探针对应多个基因
explan_final <- data.frame(explan_final[-grep("/",explan_final$"Symbol"),]) #去一对多，grep是包含的意思，-就是不包含
dim(explan_final)
[1] 8738   13
rownames(explan_final) <- explan_final$Gene.Symbol
explan_final <- explan_final[,c(1:12)]
# 此时explan_final为所需文件,行为基因，列为样本

二、DESeq差异分析以及火山图

差异分析

BiocManager::install('DESeq2') 
library(DESeq2)
library(pheatmap)  # 用于作热图的包
library(ggplot2)   # 用于作图的包

#DESeq2需要基因表达矩阵和分组信息，且基因表达矩阵是整数而且是没有经过标准化的
# 构建分组信息：即告诉这个软件哪几个样品是对照组，哪几个样品是处理组！
# 我的数据前6组是对照组，后六组是实验组
# counts为上一步预处理的结果，行为基因列为样本
condition <- factor(c(rep("untreated",6), rep("treated ",6)))    
coldata <- data.frame(row.names = colnames(counts), condition)
coldata    #显示coldata值,看看分组信息与真实数据是否一致，很关键！
         condition
GSM34219 untreated
GSM34220 untreated
GSM34221 untreated
GSM34222 untreated
GSM34223 untreated
GSM34224 untreated
GSM34225  treated 
GSM34226  treated 
GSM34228  treated 
GSM34229  treated 
GSM34230  treated 
GSM34231  treated 
# 构建dds矩阵
dds <- DESeqDataSetFromMatrix(countData = round(counts), colData = coldata, design= ~condition)
dds
class: DESeqDataSet 
dim: 8738 12 
metadata(1): version
assays(1): counts
rownames(8738): Rpl41 Actb ... Vmn2r60 Tex21
rowData names(0):
colnames(12): GSM34219 GSM34220 ... GSM34230
  GSM34231
colData names(1): condition
# 低质量数据过滤,非必要
keep <- rowSums(counts(dds) >= 10) >= 3  #过滤低表达基因，至少有3个样品都满足10个以上的reads数  
dds <- dds[keep, ] 

# 利用DESeq（）函数标准化dds矩阵
dds1 <- DESeq(dds,fitType = 'mean', minReplicatesForReplace = 7, parallel = FALSE)    # 将数据标准化，必要步骤！！！
resultsNames(dds1)    # 查看结果的名称。
dds1$condition        #默认后者的处理组比前面的对照组。
#将结果用result()函数来获取
res <- results(dds1)  
summary(res)          #看一下结果的概要信息，p值默认小于0.1。
out of 8505 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 43, 0.51%
LFC < 0 (down)     : 37, 0.44%
outliers [1]       : 49, 0.58%
low counts [2]     : 494, 5.8%
(mean count < 14)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

# res格式转化：用data.frame转化为表格形式
res1 <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)
# 依次按照pvalue值log2FoldChange值进行排序
res1 <- res1[order(res1$pvalue, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]
# 获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。
res1_up<- res1[which(res1$log2FoldChange >= 1 & res1$pvalue < 0.05),]      # 表达量显著上升的基因
res1_down<- res1[which(res1$log2FoldChange <= -1 & res1$pvalue < 0.05),]    # 表达量显著下降的基因
res1_total <- rbind(res1_up,res1_down)
# 保存
write.table(res1,file = "res1.txt",row.names = TRUE,col.names = TRUE,sep = '\t')
write.table(res1_up,file = "res1_up.txt",row.names = TRUE,col.names = TRUE,sep = '\t')
write.table(res1_down,file = "res1_down.txt",row.names = TRUE,col.names = TRUE,sep = '\t')
# 或者
# res1[which(res1$log2FoldChange >= 1 & res1$pvalue < 0.05),'sig'] <- 'up'
# res1[which(res1$log2FoldChange <= -1 & res1$pvalue < 0.05),'sig'] <- 'down'
# res1[which(abs(res1$log2FoldChange) <= 1 | res1$pvalue >= 0.05),'sig'] <- 'none'
# res1_up <- subset(res1, sig == 'up')
# res1_down <- subset(res1, sig == 'down')

差异分析的结果

火山图

#加载包
library(ggplot2)
library(ggrepel)
#数据
Dat <- res1  #res1为DESeq2差异分析结果最终结果
Dat$Gene <- rownames(res1)
#确定是上调还是下调，用于给图中点上色）
Dat$threshold = factor(ifelse(Dat$padj < 0.05 & abs(Dat$log2FoldChange) >= 1, ifelse(Dat$log2FoldChange>= 1 ,'Up','Down'),'NoSignifi'),levels=c('Up','Down','NoSignifi'))


ggplot(Dat,aes(x=log2FoldChange,y=-log10(padj),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("#DC143C","#00008B","#808080"))+#确定点的颜色
  geom_text_repel(
    data = Dat[Dat$padj<0.05&abs(Dat$log2FoldChange)>1,],
    aes(label = Gene),
    size = 3,
    segment.color = "black", show.legend = FALSE )+#添加关注的点的基因名
  theme_bw()+#修改图片背景
  theme(
    legend.title = element_blank()#不显示图例标题
  )+
  ylab('-log10 (p-adj)')+#修改y轴名称
  xlab('log2 (FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-1,1),lty=3,col="black",lwd=0.5) +#添加横线|FoldChange|>2
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=0.5)#添加竖线padj<0.05

画出的火山图

三、KEGG富集分析

##富集分析
BiocManager::install("clusterProfiler")
BiocManager::install("topGO")
BiocManager::install("Rgraphviz")
BiocManager::install("pathview")
BiocManager::install("org.Mm.eg.db")
#BiocManager::install("org.Hs.eg.db")

library(clusterProfiler)
library(topGO)
library(Rgraphviz)
library(pathview)
library(org.Mm.eg.db)

# 找到基因对应的ENTREZID，也可以自己手动对应，在第一步从GEO中下载的平台信息里就有对应的ENTREZID
# control_DESeq2_up为上一步的res1_up 即上调基因
# control_DESeq2_down为上一步的res1_down 即下调基因
df1 <- bitr(unique(rownames((control_DESeq2_up))), fromType = "SYMBOL",
           toType = c( "ENTREZID"),
           OrgDb = org.Mm.eg.db)
df2 <- bitr(unique(rownames((control_DESeq2_down))), fromType = "SYMBOL",
            toType = c( "ENTREZID"),
            OrgDb = org.Mm.eg.db)
# 增加列
control_DESeq2_up$symbol <- rownames(control_DESeq2_up)
control_DESeq2_down$symbol <- rownames(control_DESeq2_down)
deg_up =merge(control_DESeq2_up,df1,by.y='SYMBOL',by.x='symbol')
deg_down =merge(control_DESeq2_down,df2,by.y='SYMBOL',by.x='symbol')
#把数据DEG,df通过，DEG的'symbol'列，df的'SYMBOL'列连接在一起，转化ID
# head(DEG)
# save(DEG,file = 'anno_DEG.Rdata')

gene_up= deg_up[,'ENTREZID'] #选出上调基因ID
gene_down=deg_down[,'ENTREZID'] #选出下调基因ID
gene_diff=c(gene_up,gene_down)#得出上下调基因ID

##上调基因的富集分析
library(DO.db)
# R.utils::setOption("clusterProfiler.download.method",'auto') 
kk.up <- enrichKEGG(gene = gene_up,
                    organism = "mmu",   #看自己数据是什么类型我的是小鼠
                    keyType = "kegg",
                    pvalueCutoff = 0.9,
                    qvalueCutoff = 0.9)
head(kk.up)
dotplot(kk.up);
ggsave('kk.up.dotplot.pdf', height = 10, width = 8)
barplot(kk.up)
ggsave('kk.up.barplot.pdf', height = 10, width = 8)
kegg_up_dt <- as.data.frame(kk.up)
##下调基因同上
kk.down <- enrichKEGG(gene = gene_down,
                      organism = "mmu",
                      keyType = "kegg",
                      pvalueCutoff = 0.9,
                      qvalueCutoff = 0.9)
head(kk.down)
dotplot(kk.down );
ggsave('kk.down.dotplot.pdf',height = 10,width = 8)
barplot(kk.down)
ggsave('kk.down.barplot.pdf', height = 10, width = 8)
kegg_down_dt <- as.data.frame(kk.down)
# 差异基因的富集分析
kk.diff <- enrichKEGG(gene = gene_diff,
                      organism = "mmu",
                      keyType = "kegg",
                      pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
dotplot(kk.diff );
ggsave('kk.diff.dotplot.pdf',height = 10,width = 8)

上调基因的两张图

四、GO富集分析

library(clusterProfiler)
library(topGO)
library(Rgraphviz)
library(pathview)
library(org.Mm.eg.db)
##GO分析
#Cellular component解释的是基因存在在哪里
#Biological process是在说明该基因参与了哪些生物学过程
#Molecular function在讲该基因在分子层面的功能是什么
#可以分别来富集分析也可以一起富集分析， ont参数来选择

#Biological process， BP 生物学过程
erich.go.BP_up<-enrichGO(gene = gene_up,   #gene_up为上调基因的信息，在这里是已经找了对应的ENTREZID
                      OrgDb = org.Mm.eg.db,
                      keyType = "ENTREZID",
                      ont = "BP",
                      pvalueCutoff = 0.01,
                      qvalueCutoff = 0.05,
                      readable = T)
#颜色越红，生物过程越显著
barplot(erich.go.BP_up,title = "BP")
ggsave('erich.go.BP_barplot_up.pdf', height = 10, width = 8)
#点越大、颜色越红，相应生物过程的基因所占比越显著重要
dotplot(erich.go.BP_up,title = "BP")
ggsave('erich.go.BP_dotplot_up.pdf', height = 10, width = 8)
#将生物过程以树的形式进行展示
plotGOgraph(erich.go.BP_up)



#Cellular component，CC 细胞成分
erich.go.CC_up<-enrichGO(gene = gene_up,
                      OrgDb = org.Mm.eg.db,
                      keyType = "ENTREZID",
                      ont = "CC",
                      pvalueCutoff = 0.01,
                      qvalueCutoff = 0.05,
                      readable = T)

#颜色越红，细胞构成越显著
barplot(erich.go.CC_up,title = "CC")
ggsave('erich.go.CC_barplot_up.pdf', height = 10, width = 8)
#点越大、颜色越红，相应细胞构成的基因所占比越显著重要
dotplot(erich.go.CC_up,title = "CC")
ggsave('erich.go.CC_dotplot_up.pdf', height = 10, width = 8)
#将细胞构成以树的形式进行展示
##树状图很大，所以我们用代码把它存成pdf
pdf(file="enrich.go.CC_up.tree.pdf",width = 10,height = 15)
plotGOgraph(erich.go.CC_up)
dev.off()


#Molecular function，MF 分子功能
erich.go.MF_up<-enrichGO(gene = gene_up,
                      OrgDb = org.Mm.eg.db,
                      keyType = "ENTREZID",
                      ont = "MF",
                      pvalueCutoff = 0.01,
                      qvalueCutoff = 0.05,
                      readable = T)

#颜色越红，分子功能越显著
barplot(erich.go.MF_up,title = "MF")
ggsave('erich.go.MF_barplot_up.pdf', height = 10, width = 8)
#将细胞构成以树的形式进行展示
#点越大、颜色越红，相应分子功能的基因所占比越显著重要

dotplot(erich.go.MF_up,title = "MF")
ggsave('erich.go.MF_dotplot_up.pdf', height = 10, width = 8)

#将分子功能以树的形式进行展示
pdf(file="enrich.go.MF_up.tree.pdf",width = 10,height = 15)
plotGOgraph(erich.go.MF_up)
dev.off()


###########所有
erich.go_up<-enrichGO(gene = gene_up,
                         OrgDb = org.Mm.eg.db,
                         keyType = "ENTREZID",
                         ont = "ALL",
                         pvalueCutoff = 0.05,
                         qvalueCutoff = 0.05,
                         readable = T)

#颜色越红，分子功能越显著
barplot(erich.go_up)
ggsave('erich.go_barplot_up.pdf', height = 10, width = 8)
#将细胞构成以树的形式进行展示
#点越大、颜色越红，相应分子功能的基因所占比越显著重要

dotplot(erich.go_up)
ggsave('erich.go_dotplot_up.pdf', height = 10, width = 8)

#将分子功能以树的形式进行展示

# pdf(file="enrich.go_up.tree.pdf",width = 10,height = 15)
# plotGOgraph(erich.go_up)
# dev.off()

图示