如何快速下载GEO数据并获取其表达矩阵与临床信息 | 附完整代码 + 注释

本文主要是介绍如何快速下载GEO数据并获取其表达矩阵与临床信息 | 附完整代码 + 注释,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

GEO数据库可以说是大家使用频率贼高的数据库啦!那它里面的数据怎么下载大家知道嘛!今天给大家展示一种快速获取它的表达矩阵和临床信息的方法!

话不多说!咱们直接开始!

GEO编号获取

在GEO数据库中,你找到了你需要的数据,接下来怎么办嘞!下载它!处理它!

比如,咱们今天需要的数据是这个:

GEO编号这不就有咯!咱们开始下载!噢对,提前吱一声,芯片数据和非芯片数据进行ID转换的时候方法不太一样,大家要注意噢!我今天会都给大家演示一下!

GEO数据库今天咱就不介绍了哈!我后面再专门出一期超级详细地介绍它!迫切的小伙伴可以在后台或群里催我哈哈哈哈哈哈哈哈!

另外注意:不是所有的GEO数据集都可以使用这种方法下载。不过通常来说,大部分的GEO数据集都可以通过getGEO函数进行下载,但也有一些例外情况,特别是单细胞数据[无奈脸]!多数情况下,应该是数据集可能不是以类似的格式存储,或者有的作者上传格式不对等等,所以就需要找其他办法或者手动下载啦!

建议就是先试试这种方法,要是里面是空的,就去手动或者其他!

数据下载

芯片数据

#################### GEO数据下载及表达矩阵与临床信息获取 ######################## 加载包,没有安装的记得安装一下哟!
library(tidyverse)
library(GEOquery)
library(tinyarray)# 芯片数据# GEO编号,替换成你自己的就好啦!
geo_number = "GSE66360"# 这一步啊,看运气!网络时好时坏,是个玄学!一次不行不要慌!咱多跑几次!总会成功的!
geo_data <- getGEO(geo_number, destdir = './', getGPL = F) # 自己设置想保存的路径
geo_data0 <- geo_data[[1]]

数据长这样!

# 提取表达矩阵
exp <- exprs(geo_data0)
dim(exp)
# [1] 54675    99exp[1:4, 1:4]
#           GSM1620819 GSM1620820 GSM1620821 GSM1620822
# 1007_s_at   5.866378   4.948736   5.148384   5.064419
# 1053_at     8.225790   7.525865   7.764637   7.198461
# 117_at      6.179732   6.628137   5.859151   5.974135
# 121_at      6.179478   6.582889   6.602135   6.545905# 可以看到这里还是探针ID,咱们后给它面转换为基因symbol!# 可以自行判断是否需要log
exp <- log2(exp + 0.01)# 接下来我们进行ID转换# 探针注释获取
gpl_number <- geo_data0@annotation
ids <- AnnoProbe::idmap(gpl_number)
head(ids)
#         probe_id symbol
# 193731   1053_at   RFC2
# 193732    117_at  HSPA6
# 193733    121_at   PAX8
# 193734 1255_g_at GUCA1A
# 193735   1316_at   THRA
# 193736   1320_at PTPN21# ID转换
exp <- as.data.frame(exp)
exp$probe_id <- rownames(exp)
exp <- inner_join(ids, exp, by = "probe_id")
exp <- exp[!duplicated(exp$symbol), ] # 关于重复基因去重的问题,咱们后面专门出一期介绍!这里先用最大道至简的方法!
rownames(exp) <- exp$symbol
exp <- exp[ , -c(1,2)]
head(exp)[1:4, 1:4]
#        GSM1620819 GSM1620820 GSM1620821 GSM1620822
# RFC2     3.041907   2.913773   2.958775   2.849691
# HSPA6    2.629877   2.730778   2.553152   2.581143
# PAX8     2.629818   2.720911   2.725116   2.712795
# GUCA1A   1.486482   1.436010   1.373053   1.444733# 提取临床信息
cli <- pData(geo_data0)
head(cli)[1:4, 1:6]
#                               title geo_accession                status submission_date last_update_date type
# GSM1620819 Muse_1_DISCOVERY_Control    GSM1620819 Public on Feb 28 2015     Feb 27 2015      Feb 28 2015  RNA
# GSM1620820 Muse_2_DISCOVERY_Control    GSM1620820 Public on Feb 28 2015     Feb 27 2015      Feb 28 2015  RNA
# GSM1620821 Muse_3_DISCOVERY_Control    GSM1620821 Public on Feb 28 2015     Feb 27 2015      Feb 28 2015  RNA
# GSM1620822 Muse_4_DISCOVERY_Control    GSM1620822 Public on Feb 28 2015     Feb 27 2015      Feb 28 2015  RNA# 为了后续分析,咱还可以给它进一步调整,使表达矩阵列名和临床信息的行名完全一致
p <- identical(rownames(cli), colnames(exp))
p
# [1] TRUE# 噢啦!# 噢不!还没结束!那万一不是芯片数据怎么办呢!咱们再演示一个!

非芯片数据

# 非芯片数据# 比如这个GEO数据集
geo_number = "GSE16561"# 这一步啊,看运气!网络时好时坏,是个玄学!一次不行不要慌!咱多跑几次!总会成功的!
geo_data <- getGEO(geo_number, destdir = './', getGPL = F)
geo_data0 <- geo_data[[1]]# 提取表达矩阵
exp <- exprs(geo_data0)
dim(exp)
# [1] 54675    99exp[1:4, 1:4]
#                GSM416528  GSM416529  GSM416530   GSM416531
# ILMN_1343291  0.22550488  0.7348099  0.7431908  0.63261986
# ILMN_1651209 -0.03326082  0.5184302  0.5990343  0.25144434
# ILMN_1651228 -0.20895290 -0.9036322 -0.3511963 -0.06179523
# ILMN_1651229  0.18835354  0.3693209  0.3938241  0.14767408# # 可以自行判断是否需要log
# exp <- log2(exp + 1)# 大家可以看到这玩意儿不是探针,那咱们怎么办呢!

接下来我们进行非芯片数据的ID转换!

在你的GEO数据集的主页,可以看见下面这个红框框的地方!这里代表测序平台(咱们后面会在超详细介绍GEO数据库的时候给大家进行介绍这部分内容),直接点进去就好!

长这样!别急!往下滑滑!

哎对!就是这里!点红框框的地方!下载下来!

来!咱们继续看看怎么搞!

# 文件下载地址:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6883
ids <- data.table::fread("./GPL6883-11606.txt")# ID转换
exp <- as.data.frame(exp)
exp$ID <- rownames(exp)
ids <- ids[ , c(1,12)] # 注意:这里选择的是与表达矩阵行名相匹配的列和基因symbol列
exp <- inner_join(ids, exp, by = "ID")
exp <- exp[!duplicated(exp$Symbol), ]
exp <- as.data.frame(exp)
rownames(exp) <- exp$Symbol
exp <- exp[ , -c(1,2)]
head(exp)[1:4 ,1:4]
#              GSM416528  GSM416529   GSM416530  GSM416531
# JMJD1A    -0.401840200 -0.2110777 -0.42790460 -0.4261866
# NCOA3      0.296299930 -0.3155098 -0.01228809  0.5331011
# LOC389834 -0.001213074  0.5300436  0.55930424  0.2009177
# SPIRE2     0.127047540  0.5759945  0.45181800  0.4050984# 噢啦!# 老步骤!# 提取临床信息
cli <- pData(geo_data0)
head(cli)[1:4, 1:6]
#                    title geo_accession                status submission_date last_update_date type
# GSM416528 3100083_Stroke     GSM416528 Public on Jun 01 2010     Jun 11 2009      Jun 11 2009  RNA
# GSM416529 3100191_Stroke     GSM416529 Public on Jun 01 2010     Jun 11 2009      Jun 11 2009  RNA
# GSM416530 3100068_Stroke     GSM416530 Public on Jun 01 2010     Jun 11 2009      Jun 11 2009  RNA
# GSM416531 3100060_Stroke     GSM416531 Public on Jun 01 2010     Jun 11 2009      Jun 11 2009  RNA# 为了后续分析,咱还可以给它进一步调整,使表达矩阵列名和临床信息的行名完全一致
p <- identical(rownames(cli), colnames(exp))
p
# [1] TRUE

这回是真的噢啦!

一般使用GEO数据库很可能是多个数据集合并使用,那这个时候,由于数据来源不一样,肯定会有批次效应存在,那大家一定要记得进行批次处理噢!有兴趣的小伙伴们可以查看:数据挖掘 | 批次效应的鉴定与处理 | 附完整代码 + 注释 | 看完不会来揍我

在获取表达矩阵和临床信息后,我们还可以进行生存分析:看完不会来揍我 | 生存分析详解 | 从基础概念到生存曲线绘制 | 代码注释 + 结果解读

还有差异分析:看完还不会来揍/找我 | 差异分析三巨头 —— DESeq2、edgeR 和 limma 包 | 附完整代码 + 注释

还有富集分析:看完还不会来揍我 | GSEA富集分析详解(一)—— 代码实操 | MSigDB数据库介绍 | 附完整代码 + 注释

哎呀,还有好多好多!大家需要什么,都可以随时搜索!像这样:

暂时还没有的,也可以催更我!!!

文末碎碎念

那今天的分享就到这里啦!我们下期再见哟!

最后顺便给自己推荐一下嘿嘿嘿!

如果我的分享对你有用的话,欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴!这可是我的第一原动力!

蟹蟹你们的喜欢和支持!!!

如果小伙伴们有需求的话,可以加入我们的交流群:一定要知道 | 永久免费的生信交流群终于来啦!!在这里,你可以稍有克制地畅所欲言!

超级建议大家在入群前或入群后可以看一下这个:干货满满 | 给生信小白的入门小建议 | 掏心掏肺版!绝对干货满满!让你不虚此看!

如果有需要个性化定制分析服务的小伙伴,可以看看这里:你要的个性化生信分析服务今天正式开启啦!定制你的专属解决方案!全程1v1答疑!!绝对包你满意!

入群链接后续可能会不定期更新,主要是因为群满换码或是其他原因,如果小伙伴点开它之后发现,咦,怎么失效啦!不要慌!咱们辛苦一下动动小手去主页的要咨询那里,点击进交流群即可入群!

这篇关于如何快速下载GEO数据并获取其表达矩阵与临床信息 | 附完整代码 + 注释的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/851042

相关文章

乐鑫 Matter 技术体验日|快速落地 Matter 产品,引领智能家居生态新发展

随着 Matter 协议的推广和普及,智能家居行业正迎来新的发展机遇,众多厂商纷纷投身于 Matter 产品的研发与验证。然而,开发者普遍面临技术门槛高、认证流程繁琐、生产管理复杂等诸多挑战。  乐鑫信息科技 (688018.SH) 凭借深厚的研发实力与行业洞察力,推出了全面的 Matter 解决方案,包含基于乐鑫 SoC 的 Matter 硬件平台、基于开源 ESP-Matter SDK 的一

uniapp接入微信小程序原生代码配置方案(优化版)

uniapp项目需要把微信小程序原生语法的功能代码嵌套过来,无需把原生代码转换为uniapp,可以配置拷贝的方式集成过来 1、拷贝代码包到src目录 2、vue.config.js中配置原生代码包直接拷贝到编译目录中 3、pages.json中配置分包目录,原生入口组件的路径 4、manifest.json中配置分包,使用原生组件 5、需要把原生代码包里的页面修改成组件的方

公共筛选组件(二次封装antd)支持代码提示

如果项目是基于antd组件库为基础搭建,可使用此公共筛选组件 使用到的库 npm i antdnpm i lodash-esnpm i @types/lodash-es -D /components/CommonSearch index.tsx import React from 'react';import { Button, Card, Form } from 'antd'

Lipowerline5.0 雷达电力应用软件下载使用

1.配网数据处理分析 针对配网线路点云数据,优化了分类算法,支持杆塔、导线、交跨线、建筑物、地面点和其他线路的自动分类;一键生成危险点报告和交跨报告;还能生成点云数据采集航线和自主巡检航线。 获取软件安装包联系邮箱:2895356150@qq.com,资源源于网络,本介绍用于学习使用,如有侵权请您联系删除! 2.新增快速版,简洁易上手 支持快速版和专业版切换使用,快速版界面简洁,保留主

17.用300行代码手写初体验Spring V1.0版本

1.1.课程目标 1、了解看源码最有效的方式,先猜测后验证,不要一开始就去调试代码。 2、浓缩就是精华,用 300行最简洁的代码 提炼Spring的基本设计思想。 3、掌握Spring框架的基本脉络。 1.2.内容定位 1、 具有1年以上的SpringMVC使用经验。 2、 希望深入了解Spring源码的人群,对 Spring有一个整体的宏观感受。 3、 全程手写实现SpringM

【服务器运维】MySQL数据存储至数据盘

查看磁盘及分区 [root@MySQL tmp]# fdisk -lDisk /dev/sda: 21.5 GB, 21474836480 bytes255 heads, 63 sectors/track, 2610 cylindersUnits = cylinders of 16065 * 512 = 8225280 bytesSector size (logical/physical)

通过高德api查询所有店铺地址信息

通过高德api查询所有店铺地址电话信息 需求:通过高德api查询所有店铺地址信息需求分析具体实现1、申请高德appkey2、下载types city 字典值3、具体代码调用 需求:通过高德api查询所有店铺地址信息 需求分析 查询现有高德api发现现有接口关键字搜索API服务地址: https://developer.amap.com/api/webservice/gui

代码随想录算法训练营:12/60

非科班学习算法day12 | LeetCode150:逆波兰表达式 ,Leetcode239: 滑动窗口最大值  目录 介绍 一、基础概念补充: 1.c++字符串转为数字 1. std::stoi, std::stol, std::stoll, std::stoul, std::stoull(最常用) 2. std::stringstream 3. std::atoi, std

Android我的二维码扫描功能发展史(完整)

最近在研究下二维码扫描功能,跟据从网上查阅的资料到自己勉强已实现扫描功能来一一介绍我的二维码扫描功能实现的发展历程: 首页通过网络搜索发现做android二维码扫描功能看去都是基于google的ZXing项目开发。 2、搜索怎么使用ZXing实现自己的二维码扫描:从网上下载ZXing-2.2.zip以及core-2.2-source.jar文件,分别解压两个文件。然后把.jar解压出来的整个c

Android多线程下载见解

通过for循环开启N个线程,这是多线程,但每次循环都new一个线程肯定很耗内存的。那可以改用线程池来。 就以我个人对多线程下载的理解是开启一个线程后: 1.通过HttpUrlConnection对象获取要下载文件的总长度 2.通过RandomAccessFile流对象在本地创建一个跟远程文件长度一样大小的空文件。 3.通过文件总长度/线程个数=得到每个线程大概要下载的量(线程块大小)。