ToB企服应用市场:ToB评测及商务社交产业平台

标题: 生信小白菜之GEO数据库简介及数据获取方式 [打印本页]
作者: 前进之路    时间: 2024-6-11 02:27
标题: 生信小白菜之GEO数据库简介及数据获取方式
title: “Biotrainee Note7 GEO Entry Level”
author: “yuluyang”
date: “2024-03-15”
   生信技能树数据挖掘课程条记～小洁老师授课
  重要内容


表达矩阵


数据库

不同数据库来源的数据读取方式大概存在差异
不消硬性规定一定是哪个数据库来源
有什么类型的数据可以挖掘


当你有了一个表达矩阵该怎样分析

图表介绍

热图


散点图和箱线图 (单基因分析表达差异)


多基因用差异分析


主成分分析


   https://mp.weixin.qq.com/s/NR8Ou372l7Ule0lO6Hkpog
  GEO背景知识+表达芯片分析思绪


怎样下载表达矩阵

依旧是先装包

  1. options("repos"="https://mirrors.ustc.edu.cn/CRAN/") # 设置镜像
  2. if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
  3. options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") # # 设置镜像
  5. cran_packages <- c('tidyr',
  6.                    'tibble',
  7.                    'dplyr',
  8.                    'stringr',
  9.                    'ggplot2',
  10.                    'ggpubr',
  11.                    'factoextra',
  12.                    'FactoMineR',
  13.                    'devtools',
  14.                    'cowplot',
  15.                    'patchwork',
  16.                    'basetheme',
  17.                    'paletteer',
  18.                    'AnnoProbe',
  19.                    'ggthemes',
  20.                    'VennDiagram',
  21.                    'tinyarray')  # `cran`里面的R包
  22. Biocductor_packages <- c('GEOquery',
  23.                          'hgu133plus2.db',
  24.                          'ggnewscale',
  25.                          "limma",
  26.                          "impute",
  27.                          "GSEABase",
  28.                          "GSVA",
  29.                          "clusterProfiler",
  30.                          "org.Hs.eg.db",
  31.                          "preprocessCore",
  32.                          "enrichplot") # `Biocductor`里面的R包
  34. for (pkg in cran_packages){                            # `quietly = T`的意思是不要返回`warning`
  35.   if (! require(pkg,character.only=T,quietly = T) ) {  # 放进循环的时候`character.only=T`这句代码一定要记得带上
  36.     install.packages(pkg,ask = F,update = F)  # 为什么用`require()`而不用`library()`呢
  37.     require(pkg,character.only=T)   # 遇到安装失败`library()`会报错，从而导致整个循环终止，而`require()`只会`warning`
  38.   }
  39. }
  42. for (pkg in Biocductor_packages){
  43.   if (! require(pkg,character.only=T,quietly = T) ) {
  44.     BiocManager::install(pkg,ask = F,update = F)
  45.     require(pkg,character.only=T)
  46.   }
  47. }
  49. #前面的所有提示和报错都先不要管。主要看这里
  50. for (pkg in c(Biocductor_packages,cran_packages)){
  51.   require(pkg,character.only=T)
  52. }
  53. #没有任何提示就是成功了，如果有warning xx包不存在，用library检查一下。
  55. #library报错，就单独安装。
复制代码
用代码下载数据

  1. # 下载前建议先清空环境
  2. rm(list = ls())
  3. # 打破下载时间的限制,改前60秒，改后10w秒
  4. options(timeout = 100000)
  5. options(scipen = 20) # 不要以科学计数法表示
  7. # 传统下载方式
  8. library(GEOquery)
  9. eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
  11. # 其他下载方式
  12. ## 1. 从网页上下载/发链接让别人帮忙下，放在工作目录里
  13. ## 2. 试试`geoChina`,只能下载2019年前的表达芯片数据
  14. ## library(AnnoProbe)
  15. ## eSet = geoChina("GSE7305")   # 选择性代替`eSet = getGEO("GSE7305", destdir = '.', getGPL = F)`
  17. # 研究一下这个eSet
  18. class(eSet) # 返回结果`list`
  19. length(eSet) # 长度是1，即只有一个元素的列表
  20. ## 当一个列表里面只放了一个元素，那么就不再需要列表这个外壳
  22. eSet = eSet[[1]] # 取出唯一的一个元素，这个代码只能运行一遍哦，不可以重复
  23. class(eSet) # 返回结果`ExpressionSet`
  24. ## 不认识的东东就`?`问下
  25. ## `ExpressionSet`是`Biobase`包里面规定的数据结构
复制代码
只要末了数据在工作目录里即可，无所谓下载方式
‼️一个工作目录里面不要放两个数据集
R包里面不止有函数，还有它定义的数据类型
由R包作者定义的、以某种方式组织起来的、存放特定数据的数据结构
提取GSE内的信息

  1. # 提取表达矩阵exp
  2. exp <- exprs(eSet) # `exprs()`函数提取表达矩阵
  3. dim(exp) # 查看表达矩阵维度
  4. range(exp) # 看数据范围决定是否需要`log`，是否有负值、异常值
  6. exp = log2(exp+1) # 由上一步决定需要`log`才`log`
  7. # 如果`log`出来所有值都在2-4，可能你取了两次或者这个数据本身就不需要取`log`
  8. boxplot(exp,las = 2) # 通过箱线图看是否有异常样本
  9. # 如果存在某一个样本值异常低，可以考虑剔除或者运行标准化代码
  10. exp=limma::normalizeBetweenArrays(exp)
  12. # 提取临床信息
  13. pd <- pData(eSet)
  14. # 阅读理解环节，把分组信息提取出来
  15. ## 只可意会不可言传，仔细读样本名称来确定分组
复制代码
如果一个数据集有一半左右是负值，中位数在0左右，且所有数据集中在-5～5或-10～10，则不正常
如果一个数据集出现这个情况是不适合做差异分析的，必要从原始数据开始处理
如果不处理原始数据，那么这个数据集最多只能用来作热图和生存分析
仅有少量负值为正常
表达矩阵与临床信息对应

  1. # 让exp列名与pd的行名顺序完全一致
  2. # 即：表达矩阵的行名和临床信息的列名（样品名）一一对应
  3. # 为了准确传递分组信息
  4. p = identical(rownames(pd),colnames(exp));p # `identical()`判断是否一致
  5. if(!p) {                                    # 一致的话`(!p)`就是`False`,那么后面的代码就不运行啦
  6.   s = intersect(rownames(pd),colnames(exp)) # `intersect()`取交集 ## 非常重要的一个函数
  7.   exp = exp[,s]                             # 表达矩阵按列取
  8.   pd = pd[s,]                               # 临床信息按行取
  9. }
复制代码
提取芯片平台编号

  2. gpl_number <- eSet@annotation;gpl_number  # 有些是用`@`,有些是用`$`
  3. save(pd,exp,gpl_number,file = "step1output.Rdata") # 知道平台，后续才能知道探针注释
复制代码
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。



欢迎光临 ToB企服应用市场:ToB评测及商务社交产业平台 (https://dis.qidao123.com/) Powered by Discuz! X3.4