GEO数据下载到差异处理笔记

最新推荐文章于 2025-05-10 13:42:11 发布

Tiaoxiaobai

最新推荐文章于 2025-05-10 13:42:11 发布

阅读量601

点赞数 5

CC 4.0 BY-SA版权

文章标签：数据库笔记

本文链接：https://blog.youkuaiyun.com/Tiaoxiaobai/article/details/143890231

1.进入网站：

https://ftp.ncbi.nlm.nih.gov/geo/series/GSE56nnn/

找到自己要下载的GEO号，比如GSE56427

2.解压文件

得到.CEL结尾文件

tar -xvf GSE56416_RAW.tar

3.处理CEL格式文件

根据网上介绍利用affy包进行处理，但是本次下载数据不适应于affy处理（报错如下），因此更换oligo软件包进行处理

（1）原始数据得到表达矩阵：

.libPaths("D:\\software\\R")
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("oligo", lib="D:\\software\\R")

BiocManager::install("pd.mogene.1.1.st.v1", lib="D:\\software\\R")
BiocManager::install("GEOquery", lib="D:\\software\\R")

library(pd.mogene.1.1.st.v1)
library(oligo)
library(GEOquery)
celFiles <- list.celfiles('./', full.names=TRUE) #读取.CEL格式文件
rawData <- read.celfiles(celFiles, pkgname = "pd.mogene.1.1.st.v1", verbose = F)
exset_rma <- rma(rawData) #利用 RMA 算法进行背景校正、归一化及 summarization
expmatx <- exprs(exset_rma)  ##表达矩阵（探针ID）

（2）探针ID转换为基因ID，Ensembl ID

首先根据提供GSE号进入相应网页，然后点击下载注释文件信息

应该有代码下载方式，为了方便省时间，直接从网页下载然后读取文件，接下来在R中处理

.libPaths("D:\\software\\R")
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("oligo", lib="D:\\software\\R")

BiocManager::install("pd.mogene.1.1.st.v1", lib="D:\\software\\R")
BiocManager::install("GEOquery", lib="D:\\software\\R")

library(pd.mogene.1.1.st.v1)
library(oligo)
library(GEOquery)
library(stringr)
celFiles <- list.celfiles('./', full.names=TRUE) 
rawData <- read.celfiles(celFiles, pkgname = "pd.mogene.1.1.st.v1", verbose = F)
exset_rma <- rma(rawData) #利用 RMA 算法进行背景校正、归一化及 summarization. 得到 ExpressionSet 对象
expmatx <- as.data.frame(exprs(exset_rma))
expmatx$ID <- rownames(expmatx)
## 读取所有注释文件
data <- read.delim("metadata/GPL11533_annot.txt", header = T) # 读取下载的注释文件
annotfiles <- data[which(data$gene_assignment != "---"), c(1,10)]
## gene处理
library(stringr)
annotfiles$gene <- str_split(annotfiles$gene_assignment,'//',simplify = T)[,2]
annotfiles$gene_assignment <- NULL
## 获取symbol名对应的表达矩阵
df <- merge(annotfiles, expmatx, by="ID")
df$ID <- NULL
write.csv(file='express_rma.txt', row.names = F, quote = F)

得到如下表达矩阵，可以进行后续的可视化以及差异分析

4.差异基因分析

##limma进行差异分析
library(limma)
library(dplyr)
names(expr) <- c("geneID","GSM1361232","GSM1361233", "GSM1361234", "GSM1361238", "GSM1361239", "GSM1361240")
rownames(expr) <- expr$geneID
expr$geneID <- NULL
meta <- c(rep("WT_MTX", 3), rep("WT.MTX",3)) %>% factor(., levels = c("WT_MTX", "WT.MTX"), ordered = F)

meta <- model.matrix(~factor(meta)+0)  #把group设置成一个model matrix
colnames(meta) <- c("WT_MTX", "WT.MTX")
expr.fit <- lmFit(expr, meta)  ## 数据与meta进行匹配  

# 差异分析
expr.matrix <- makeContrasts(WT.MTX - WT_MTX , levels = meta)
fit <- contrasts.fit(expr.fit, expr.matrix)
fit <- eBayes(fit)
tempOutput <- topTable(fit,n = Inf, adjust = "fdr")
tempOutput$stat <- "non-significant"
tempOutput$stat[which(tempOutput$logFC >0.5849625 & tempOutput$P.Value < 0.05)] <- 'up-regulated'
tempOutput$stat[which(tempOutput$logFC  < -0.5849625 & tempOutput$P.Value < 0.05)] <- 'down-regulated'