生信小白菜之GEO数据库简介及数据获取方式

---

title: “Biotrainee Note7 GEO Entry Level”
author: “yuluyang”
date: “2024-03-15”

---

生信技能树数据挖掘课程笔记～小洁老师授课

主要内容

• 介绍公共数据库的挖掘方法
• 基因筛选的方法
• 讲解象限图的意义
• 介绍logFC的意义
• 火山图范围及数据分析的意义
• 讲解数据的主成分分析
• 讲解实验设计的目的
• 介绍数据挖掘的方法
• 数据下载的方法及热图
• 简单介绍数据下载方式
• 关于数据赋值的问题

表达矩阵

• 行是某一个基因在不同样本中表达量；列是某一个样本的不同基因表达量

数据库

不同数据库来源的数据读取方式可能存在差异
不用硬性规定一定是哪个数据库来源

1. GEO
   GENE EXPRESSION OMNIBUS，基因表达数据库
   收录了世界各国研究机构提交的基因表达数据，主要包括肿瘤、非肿瘤、芯片、NGS、差异分析、分子验证等各种公开数据

2. NHANES
   National Health and Nutrition Examination Survey，美国国家健康与营养调查
   一项基于全美各层次人群的横断面调查，旨在收集有关美国家庭人口健康、营养和社会学信息

3. TCGA
   The Cancer Genome Atlas，肿瘤基因组图谱
   涵盖癌症病人各种各样的测序数据，如RNA sequencing、MicroRNA sequencing、DNA sequencing、SNP-based platforms、Array-based DNA methylation sequencing、Reverse-phase array

4. ICGC
   International Cancer Genome Consortium，国际肿瘤基因组联盟
   与TCGA数据库类似之处：都进行多组学研究；不同之处：ICGC数据库的数据来源更加丰富，涵盖了多个国家和地区的数据，包括亚洲人的数据

5. CCLE
   Cancer Cell Line Encyclopedia，癌症细胞系百科全书
   一个癌症细胞系数据库，内含超过1000种细胞系的mRNA、miRNA表达、外显子突变、H3组蛋白尾部修饰与部分代谢物丰度数据

6. SEER(only WIN)
   The Surveillance, Epidemiology, and End Results Program，监测、统计流行病学和最终结果
   收录了美国约30％人口的癌症诊断、治疗和生存数据，在人群特征、病理、免疫组化、放化疗以及随访等多方面均有统计，以临床回顾性数据总结归纳

有什么类型的数据可以挖掘

• 基因表达芯片
• 转录组
• 单细胞
• 甲基化、拷贝数突变
• 临床数据
• …

当你有了一个表达矩阵该如何分析

1. 数据下载
   GEO、TCGA有专属R包用于下载数据
2. 差异分析
   芯片与转录组分析方法不相同哦
3. WGCNA(加权共表达网络)
   一种能够将基因、样本、基因模块和生物性状相互关联起来的分析方法
4. 富集分析
   两种方法：ORA、GSEA
5. PPI网络(蛋白质互作)
6. 预后分析
   不局限于肿瘤，只要是影响生存的疾病都可以

图表介绍

热图

• 看差异基因
• pheatmap包
• 输入数据是数值型matrix/data.frame
• 颜色变化是数值大小
• 聚类树：层次聚类

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

• 表达的内容、意思相似
• 箱线图输入数据要求
  是一个连续型向量和一个有重复值的离散型的量
  最大值、最小值是去掉离群值后计算出来的
• 箱线图常用于比较单个基因在两组或多组间的表达量差异
• 箱线图的5条线，从下至上依次代表数据’min’、‘25%’、‘50%’、‘75%’、'max’几个位置
  'min’和’max’是去掉了离群值后根据算法得出的结果

多基因用差异分析

• 拿到差异分析的结果，最需要关注的几个值包括：logFC，P.Value，FDR
• Pvalue即P值，因为基因和基因并不是相互独立的,通常都要进行校正来降低结果的假阳性
• adj.P.Val、qvalue是校正后的P值，常用的校正方法为FDR校正
• 火山图的横坐标通常是LogFC，纵坐标通常是-log10(qvalue)——即P值越小，这个值就越大
  原始Pvalue取值0.1、0.01、0.001在数轴上差距其实非常小，用log10转换来放大差异
• logFC即log2(foldchange)差异倍数，Foldchange处理组表达量平均值/对照组表达量平均值 *用原始的foldchange描述上调方便[1,+∞)，下调不方便(0, 1)，绘制图片中上调占的空间太多而下调空间太少，所以一般会做log2`转换*
  ‼️处理和对照一定不能弄反
• 运算一下就是：log2(x/y) = log2(x) - log2(y)，logFC>0即处理组表达量上升，<0即表达量下降
  log2(x/y) = 5是上调了2^5即32倍，log2(x/y) = -4是下调2^4即1/16
• log2(x/y)的正常范围应该是0～20之间 （2^20=1048576）
  芯片差异分析的起点是一个取过log的表达矩阵，转录组、单细胞没有这个要求
  如果芯片数据VALUE值大于24，则需要自己取过log后再分析
• 通常上调/下调基因的筛选条件：|logFC|>2，P<0.01（依此绘制虚线）
  当差异基因数量过多的时候，可适当缩窄筛选标准以获得相对适合数量的差异基因再进行下游的富集分析，比如要求’qvalue<0.001’，或’|logFC|>5’等
  当差异基因数量过少的时候，可适当放宽筛选标准，比如要求’qvalue<0.05’，或’|logFC|>1.5’等
  当Foldchange没有那么明显的时候，横轴可以直接选择展示Foldchange，不取log2

主成分分析

• PCA样本聚类图
• PCA中心思想是降维，把多指标转化为少数几个综合指标（覆盖大部分信息），即主成分dim1、dim2…（或称为PC1、PC2）
  不用刻意找出每个主成分具体是什么
• 根据这些主成分对样本进行聚类，图上的点代表样本，点与点之间的距离代表样本与样本之间的差异
• 在坐标轴上距离越远，说明样本差异越大；理想情况是同一分组聚成一簇（组内重复好），而组间中心点存在差距（组间差异大）
• PCA的优点
  简化运算
  去掉数据噪音
  利用散点图实现多维数据可视化
  发现隐性相关变量

https://mp.weixin.qq.com/s/NR8Ou372l7Ule0lO6Hkpog

GEO背景知识+表达芯片分析思路

• 官网首页有GEO2R，是它自带的分析代码
• 官网首页’Browse Content’模块
  • 一个’Series’是一个独立的实验设计
  • 'Platforms’是测序平台
  • 'Samples’是涉及到的样本
• GSM是用户提供到GEO的样本单位
• GPL用户测定表达量使用的芯片/平台
  部分数据没有提供’gene symbol’，需要根据芯片编号找到二者的对应关系
• GSE是完成的研究单位，包括研究描述、数据描述、总结分析等
  挖掘数据就找GSE
• 芯片测序数据GSE名称Expression profiling by array
• 单细胞/转录组测序数据GSE名称Expression profiling by high throughput sequencing
• 芯片的表达矩阵下载位置
  • 页面下方Download family栏目下Series Matrix File点进去看见GSE000000 series matrix.txt.gz
• 转录组/单细胞表达矩阵下载位置
  • 页面下方Supplementary file栏目下GSE000000 RAW.tar
• 下载前先判断样本合不合格
  • 首先看样本量以及文件大小(至少>500k)
  • 再点进一个样本看看VALUE值，如果有一半是负值可能有问题
• 下载哪些内容
  • 表达矩阵（列名是探针id，行名是样本名称GSMxxx）
  • 临床信息（分组信息）
    后续分析都是基于分组而非单个的样本
  • GPL编号（探针对应的基因注释，有些有GENE-SYMBOL就更简单）

如何下载表达矩阵

依旧是先装包

options("repos"="https://mirrors.ustc.edu.cn/CRAN/") # 设置镜像
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") # # 设置镜像

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'paletteer',
                   'AnnoProbe',
                   'ggthemes',
                   'VennDiagram',
                   'tinyarray')  # `cran`里面的R包
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot") # `Biocductor`里面的R包

for (pkg in cran_packages){                            # `quietly = T`的意思是不要返回`warning`
  if (! require(pkg,character.only=T,quietly = T) ) {  # 放进循环的时候`character.only=T`这句代码一定要记得带上 
    install.packages(pkg,ask = F,update = F)  # 为什么用`require()`而不用`library()`呢
    require(pkg,character.only=T)   # 遇到安装失败`library()`会报错，从而导致整个循环终止，而`require()`只会`warning`
  }
}


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T,quietly = T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}
#没有任何提示就是成功了，如果有warning xx包不存在，用library检查一下。

#library报错，就单独安装。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

用代码下载数据

# 下载前建议先清空环境
rm(list = ls())
# 打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) 
options(scipen = 20) # 不要以科学计数法表示

# 传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)

# 其他下载方式
## 1. 从网页上下载/发链接让别人帮忙下，放在工作目录里
## 2. 试试`geoChina`,只能下载2019年前的表达芯片数据
## library(AnnoProbe) 
## eSet = geoChina("GSE7305")   # 选择性代替`eSet = getGEO("GSE7305", destdir = '.', getGPL = F)`

# 研究一下这个eSet
class(eSet) # 返回结果`list`
length(eSet) # 长度是1，即只有一个元素的列表
## 当一个列表里面只放了一个元素，那么就不再需要列表这个外壳 

eSet = eSet[[1]] # 取出唯一的一个元素，这个代码只能运行一遍哦，不可以重复
class(eSet) # 返回结果`ExpressionSet` 
## 不认识的东东就`?`问下
## `ExpressionSet`是`Biobase`包里面规定的数据结构
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

只要最后数据在工作目录里即可，无所谓下载方式
‼️一个工作目录里面不要放两个数据集
R包里面不止有函数，还有它定义的数据类型
由R包作者定义的、以某种方式组织起来的、存放特定数据的数据结构

提取GSE内的信息

# 提取表达矩阵exp
exp <- exprs(eSet) # `exprs()`函数提取表达矩阵
dim(exp) # 查看表达矩阵维度
range(exp) # 看数据范围决定是否需要`log`，是否有负值、异常值

exp = log2(exp+1) # 由上一步决定需要`log`才`log` 
# 如果`log`出来所有值都在2-4，可能你取了两次或者这个数据本身就不需要取`log`
boxplot(exp,las = 2) # 通过箱线图看是否有异常样本
# 如果存在某一个样本值异常低，可以考虑剔除或者运行标准化代码
exp=limma::normalizeBetweenArrays(exp)

# 提取临床信息
pd <- pData(eSet)
# 阅读理解环节，把分组信息提取出来
## 只可意会不可言传，仔细读样本名称来确定分组
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

如果一个数据集有一半左右是负值，中位数在0左右，且所有数据集中在-5～5或-10～10，则不正常
如果一个数据集出现这个情况是不适合做差异分析的，需要从原始数据开始处理
如果不处理原始数据，那么这个数据集最多只能用来作热图和生存分析
仅有少量负值为正常

表达矩阵与临床信息对应

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

提取芯片平台编号

gpl_number <- eSet@annotation;gpl_number  # 有些是用`@`,有些是用`$`
save(pd,exp,gpl_number,file = "step1output.Rdata") # 知道平台，后续才能知道探针注释
1
2
3