title: “Biotrainee Note7 GEO Entry Level”
author: “yuluyang”
date: “2024-03-15”
生信技能树数据挖掘课程条记～小洁老师授课
重要内容
- 介绍公共数据库的挖掘方法
- 基因筛选的方法
- 讲解象限图的意义
- 介绍logFC的意义
- 火山图范围及数据分析的意义
- 讲解数据的主成分分析
- 讲解实验计划的目标
- 介绍数据挖掘的方法
- 数据下载的方法及热图
- 简朴介绍数据下载方式
- 关于数据赋值的问题
表达矩阵
- 行是某一个基因在不同样本中表达量；列是某一个样本的不同基因表达量
数据库
不同数据库来源的数据读取方式大概存在差异
不消硬性规定一定是哪个数据库来源
- GEO
GENE EXPRESSION OMNIBUS，基因表达数据库
收录了天下各国研究机构提交的基因表达数据，重要包括肿瘤、非肿瘤、芯片、NGS、差异分析、分子验证等各种公开数据
- NHANES
National Health and Nutrition Examination Survey，美国国家康健与营养观察
一项基于全美各层次人群的横断面观察，旨在收集有关美国家庭人口康健、营养和社会学信息
- TCGA
The Cancer Genome Atlas，肿瘤基因组图谱
涵盖癌症病人各种各样的测序数据，如RNA sequencing、MicroRNA sequencing、DNA sequencing、SNP-based platforms、Array-based DNA methylation sequencing、Reverse-phase array
- ICGC
International Cancer Genome Consortium，国际肿瘤基因组同盟
与TCGA数据库类似之处：都举行多组学研究；不同之处：ICGC数据库的数据来源更加丰富，涵盖了多个国家和地区的数据，包括亚洲人的数据
- CCLE
Cancer Cell Line Encyclopedia，癌症细胞系百科全书
一个癌症细胞系数据库，内含超过1000种细胞系的mRNA、miRNA表达、外显子突变、H3组蛋白尾部修饰与部门代谢物丰度数据
- SEER(only WIN)
The Surveillance, Epidemiology, and End Results Program，监测、统计流行病学和最终结果
收录了美国约30％人口的癌症诊断、治疗和生存数据，在人群特性、病理、免疫组化、放化疗以及随访等多方面均有统计，以临床回顾性数据总结归纳
有什么类型的数据可以挖掘
- 基因表达芯片
- 转录组
- 单细胞
- 甲基化、拷贝数突变
- 临床数据
- …
当你有了一个表达矩阵该怎样分析
- 数据下载
GEO、TCGA有专属R包用于下载数据
- 差异分析
芯片与转录组分析方法不雷同哦
- WGCNA(加权共表达网络)
一种能够将基因、样本、基因模块和生物性状相互关联起来的分析方法
- 富集分析
两种方法：ORA、GSEA
- PPI网络(蛋白质互作)
- 预后分析
不范围于肿瘤，只要是影响生存的疾病都可以
图表介绍
热图
- 看差异基因
- 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报错，就单独安装。
- # 下载前建议先清空环境
- 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`包里面规定的数据结构
‼️一个工作目录里面不要放两个数据集
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)
- # 阅读理解环节，把分组信息提取出来
- ## 只可意会不可言传，仔细读样本名称来确定分组
如果一个数据集出现这个情况是不适合做差异分析的，必要从原始数据开始处理
如果不处理原始数据，那么这个数据集最多只能用来作热图和生存分析
仅有少量负值为正常
表达矩阵与临床信息对应
- # 让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,] # 临床信息按行取
- }
- gpl_number <- eSet@annotation;gpl_number # 有些是用`@`,有些是用`$`
- save(pd,exp,gpl_number,file = "step1output.Rdata") # 知道平台，后续才能知道探针注释
|
