一、GEO数据库数据下载和初步分析
注释:本文章是本人学习 B站 生信技能树-jimmy 课程 整理所得,如有错误,请见谅;
·
1、R包安装
所使用R 的版本为 4.4.1
备注:低版本的R可能安装不上“GSEABase","GSVA","clusterProfiler”
# GEO数据分析,安装相关包
# GEO数据分析所用到的包如下
#BiocManager
#GEOquery
#KEGG.db
#GSEABase","GSVA","clusterProfiler
#limma","impute
#genefu","org.Hs.eg.db","hgu133plus2.db
#设置当前工作目录
#setwd(r"[D:/Bio_info/GEO Data Explore/CEO_EXPLORE_SWW/GEO]")
rm(list = ls())#删除工作空间中所有的对象
#设置R包下载的镜像
#options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) #对应清华源
#options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") #对应中科大源
#查看目前所使用的镜像
#options()$repos
#options()$BioC_mirror
rm(list = ls())
options()$repos
options()$Bioc_mirror # 查看当前R所使用的镜像
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) #对应清华源
#options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") #对应中科大源
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options()$repos
options()$Bioc_mirror # 查看当前R所使用的镜像
#安装所用的包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("KEGG.db",ask = F,update = F)
BiocManager::install(c("GSEABase","GSVA","clusterProfiler" ),ask = F,update = F)
BiocManager::install(c("GEOquery","limma","impute" ),ask = F,update = F)
BiocManager::install(c("genefu","org.Hs.eg.db","hgu133plus2.db" ),ask = F,update = F)
install.packages('WGCNA')
install.packages(c("FactoMineR", "factoextra"))
install.packages(c("ggplot2", "pheatmap","ggpubr"))
#加载相应的包
library(GSEABase)
library(GSVA)
library(clusterProfiler)
library(genefu)
library(ggplot2)
library(ggpubr)
library(hgu133plus2.db)
library(limma)
library(org.Hs.eg.db)
library(pheatmap)
2、数据下载
rm(list = ls())#删除工作空间中所有的对象
options(stringsASFactors =F) # 在调用as.data.frame 时,将StringAsFactors 设置为F,避免将charactor 转换成factor类型
#注意查看下载文件的大小,检查数据
GEO_file1 = 'GSE42872_eSet.Rdata'
#需要注意的是,如果你没有设置这个选项,GEOquery会默认使用基于curl和wget的自动选择下载方法。
#设置这个选项为'auto',会使得GEOquery在可能的情况下优先使用基于curl的下载方法,因为它通常更快,且对防火墙的要求较低。
## 这个包需要注意两个配置,一般来说自动化的配置是足够的。
#Setting.options('download.file.method.GEOquery' = 'auto')
#setting.options('GEOquery.inmemory.gpl' = FALSE)
#加载GEOquery包
library(GEOquery)
#获取数据
if( !file.exists(GEO_file1) ){
#getGEO函数是GEOquery包中的一个重要函数,用于从GEO数据库下载数据。该函数的参数包括:
#GE:决定下载的数据种类,可以是GDS/GSE/GSM/GPL,具体取决于提供的GEO参数。
#filenam:如果已经下载好了文件,可以直接读取。
#destdi:下载文件存放的地址,默认为工作目录。
#GSElimit:用于限制下载的GSE数据集的数量。
#GSEMatri:若为TRUE,则下载Matrix(GSE数据集的表达矩阵)文件;若为FALSE,则下载SOFT文件。默认为TRUE。
#AnnotGP:若为TRUE,则下载GPL注释文件;若为FALSE,则不下载。默认为FALSE。
#getGP:若为TRUE,则下载GPL注释文件;若为FALSE,则不下载。默认为TRUE。
#parseCharacteristic:未在搜索结果中明确说明,可能是与解析数据特性相关的参数。
file1 <- getGEO("GSE42872",
destdir = ".",
AnnotGPL = F, #注释文件
getGPL = F) # 平台文件
save(file1,file = GEO_file1) ## 保存文件到本地
}
3、数据分析
#因为这个GEO数据集只有一个GPL平台,所以下载到的是一个含有一个元素的list
length(file1)
class(file1)
#使用exprs函数:exprs函数可以从矩阵、数组或数据框等多种数据结构中提取基因表达值。
#例如,如果expression_set是一个包含基因表达数据的数据结构,那么通过exprs(expression_set)可以提取出基因表达值,并返回一个向量、矩阵或数组形式的输出
b <- file1[[1]]
#提取的数据如下
Geodata = exprs(file1[[1]])
class(Geodata) #[1] "matrix" "array"
# GEO数据集的每一列代表一个样本。
# 在GEO数据集中,每一列并不是直接对应于某个特定的生物样本属性,而是代表一个具体的样本。
# 这些数据集通常包含大量的基因表达数据,其中每一行代表一个基因,而每一列则代表一个样本。
# 这样的数据结构使得研究人员能够分析特定基因在不同样本(如不同组织、不同条件下的样本)中的表达情况
# ,从而研究基因表达与特定生物学过程或疾病状态之间的关系。
dim(Geodata)
# [1] 33297 6
# 33297个基因,6个样本
4、下载GPL注释文件
# GPL6244 GPL平台注释
if(F){
library(GEOquery)
#Download GPL file, put it in the current directory, and load it:
gpl <- getGEO('GPL6244', destdir=".")
colnames(Table(gpl))
# [1] "ID" "GB_LIST" "SPOT_ID" "seqname" "RANGE_GB" "RANGE_STRAND" "RANGE_START"
# [8] "RANGE_STOP" "total_probes" "gene_assignment" "mrna_assignment" "category"
# 正常应该有gene-Sympol 这一列,但是现在没有,所以需要切割字符串,获得geneid
class(head(Table(gpl)))
#[1] "data.frame"
# 字符串切割获取id
geneid_list <- strsplit(Table(gpl)$gene_assignment,"//")
#class(geneid_list) #list
gene_Symbol <- lapply(geneid_list,\(x) ifelse(length(x) < 2,x[[1]], x[[2]]))
#获取探针名
probe_id <- Table(gpl)$ID
#class(gene_Symbol) list
#class(probe_id) Integer
#合并两个数据,保存在数据框中
probetogene <- data.frame(
`probe_id` <-probe_id,
`gene_Symbol`<- unlist(gene_Symbol)
)
save(probetogene,file='probetogene.Rdata')
}
class((Table(gpl))
#[1] "data.frame"
Table(gpl)
# 字符串切割获取id
geneid_list <- strsplit(Table(gpl)$gene_assignment,"//")
#class(geneid_list) #list
# 获取基因名
gene_Symbol <- lapply(geneid_list,\(x) ifelse(length(x) < 2,x[[1]], x[[2]]))
Table(gpl)$gene_assignment
注释:本文章是本人学习 B站 生信技能树-jimmy 课程 整理所得,如有错误,请见谅;如有侵权,请联系本人,本人会立刻删除,谢谢。
更多推荐
已为社区贡献1条内容
所有评论(0)