一、GEO数据库数据下载和初步分析
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"")
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[]
https://i-blog.csdnimg.cn/direct/477ffc75df454a0ca563e2cb15482d27.png
#提取的数据如下
Geodata = exprs(file1[]) https://i-blog.csdnimg.cn/direct/61d36fc649f5430e82eaa44531dfe449.png
class(Geodata) # "matrix" "array"
# GEO数据集的每一列代表一个样本。
# 在GEO数据集中，每一列并不是直接对应于某个特定的生物样本属性，而是代表一个具体的样本。
# 这些数据集通常包含大量的基因表达数据，其中每一行代表一个基因，而每一列则代表一个样本。
# 这样的数据结构使得研究人员能够分析特定基因在不同样本（如不同组织、不同条件下的样本）中的表达情况
# ，从而研究基因表达与特定生物学过程或疾病状态之间的关系。
dim(Geodata)
# 33297 6
# 33297个基因，6个样本 4、下载GPL注释文件
# GPL6244GPL平台注释
if(F){
library(GEOquery)
#Download GPL file, put it in the current directory, and load it:
gpl <- getGEO('GPL6244', destdir=".")
colnames(Table(gpl))
# "ID" "GB_LIST" "SPOT_ID" "seqname" "RANGE_GB" "RANGE_STRAND" "RANGE_START"
# "RANGE_STOP" "total_probes" "gene_assignment" "mrna_assignment" "category"
# 正常应该有gene-Sympol 这一列，但是现在没有，所以需要切割字符串，获得geneid
class(head(Table(gpl)))
# "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[], x[]))
#获取探针名
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))
# "data.frame"
Table(gpl) https://i-blog.csdnimg.cn/direct/77e17676be3a4202aa98e69e27eb0481.png
# 字符串切割获取id
geneid_list <- strsplit(Table(gpl)$gene_assignment,"//")
#class(geneid_list) #list
# 获取基因名
gene_Symbol <- lapply(geneid_list,\(x) ifelse(length(x) < 2,x[], x[]))
Table(gpl)$gene_assignment https://i-blog.csdnimg.cn/direct/87cdac9a0b234488bd2b784ecfd12f21.png
注释：本文章是本人学习 B站 生信技能树-jimmy 课程 整理所得，如有错误，请见谅；如有侵权，请接洽本人，本人会立刻删除，谢谢。
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。
页:
[1]