GEO学习笔记-P1-P2

学习材料：【生信技能树】公共数据库挖掘实例（基于R语言） bilibili版本以及后续更新课程中的github材料为基础。

本章节是以：【生信技能树】公共数据库挖掘实例（基于R语言）为基础，进行的代码复现与注释，章节标注以之为准。

【生信技能树】公共数据库挖掘实例（基于R语言）_哔哩哔哩_bilibili今天小编给大家带来的是由我们jimmy大大亲自录制的公共数据库挖掘实例~纸上得来终觉浅，绝知此事要躬行，一起跟着大大来实践吧~https://www.bilibili.com/video/BV1Lt411S786?p=2其他资料：
【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1is411H7Hq?p=1

数据分析github源代码https://github.com/jmzeng1314/GEOhttps://github.com/jmzeng1314/GEO主要是survival文件夹中的代码。

一、绪论

现代处理数据思路：通过筛选进行寻找上调/下调的基因集。

【ps 以功能相关的基因集合为单位，而不是以单个基因通路为单位。因此此处的基因集就是与该类功能相关的所有基因的集合】

【ps 例如，化疗处理前后获得上下调的基因组合，即是一个基因集，因此可以作为一个功能/应激的改变，即应对什么的一个基因功能集合】

eg： Msigdb ，常用的GO/KEGG都是该基因集，并且包含免疫等等基因集。

【根据p值或者变化程度进行筛选，利用t检验等来筛选出变化最大的基因集-差异分析？基因识别？所以合适数量的基因集是很有意义的】

通过p或阈值的差异获取表达矩阵的差异基因，然后构建差异富集分析，同时利用GO/KEGG等进行功能注释。

在获取差异基因时，除了使用p直接比较，也可以使用其他统计学方法，例如生存分析等分析出有差异的部分。

二、GEO库

2.1、GEO处理的基本思路

        1、选择GSE:基于GSE号

        2、表达矩阵：下载表达矩阵

        3、差异分析得到基因集 - limma包就可以用 - 多个差异分析的方法可以使用

5大数据库的注释-对  差异基因进行注释 - misgdb/GO/KEGG

构建PPI网络或其他

GEO主要包含4种类型的数据：平台（Platform,GPL）、样本（Sample,GSM）、系列（Series，GSE）、数据集(Datasets,GDS)，以及表达谱。

四种编号GPL，GDS，GSE和GSM可以获得完整的平台，数据集，系列以及样本的信息.

其中：

原始数据：GPL(Platform),GSM(Sample),GSE(Series)

数据集：GDS(DataSets), 表达谱(Profiles).

GPL:平台，主要是探针列表。

GSM:

每个GSM数据都是会有相应的平台信息-经哪个技术获取该信息。

2.2、GEO数据接口:

直接修改GEO数据后面的数字编号即可直接定位到该文章中的GSM数据库。可以搜索GSM在文中获取响应信息。

2.3、GEO数据下载与安装（R语言）

        2.3.1        R包：GEOquery安装

# 安装GEOquery包
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("GEOquery")

library(GEOquery)

这里遇到一个问题，安装了部分的包后，仍无法library(GEOquery)包。

方法1：update部分包 some

方法2：翻墙试一下（trying）。

解决方案：error可能会提示调用的包版本有问题，需要重新在package栏中重新删除旧的version，并安装新的version。部分包可以通过install.package("*")或者手动操作直接完成，也会有部分新的version需要biocManager::install("*")语句重新安装包。（具体见上）

        2.3.2        数据下载

【方法1】利用R包GEOquery中的语句getGEO完成

if(!file.exists(f)){
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE42872_eSet.Rdata')  ## 载入数据

AnnotGPL和getGPL是两个注释文件，如果不想要其对应的是soft文件。

此时下载的就是表达矩阵，如果无法下载就是被墙了，需要想办法翻才行。

需注意此时在工作目录下会生成相应的gz压缩文件，此时可以load相应文件即可。

如果此时下载好了也可以使用read.table()读取这个文件。

read有很多模式，可以读进不同的参数进而获取不同的文件格式，例如使用不同的分隔符、不同的read后缀等获取不同的文件，需要去看help或者相关书籍重新学习！需要转换的话，可以用expr等语句。

ps：当然也存在直接读这个gz文件进来就是matrix的格式，需要注意。

此处应该常规检查文件格式。str()看内部结构，class()看文件类型。

---

【方法2】网页直接下载压缩文件

利用该GEO数据网页中的supplementary file中download下custom右键进入，或者download family中的series matrix files中下载原始TXT文件。

或者这里有下载的matrix也可以直接下下来，省事。

代码后可以直接格式化

通过class及length可以获得该gset的性质为list

那么gset[[1]]之后可以查看该list中的对象。这个对象的性质是ExpressionSet。

换而言之，这个gset的list中只有一个objective就是这个gset[[1]]，其性质为ExpressionSet，因此此时它还不是一个表达矩阵，需要时可通过exprs()转化成为一个表达矩阵。

此表达矩阵需检查一下是否取对数，一般信号获取的是1k-10k的信号强度，需通过log来降级。

# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]]  #提取gset的list中的第一个元素对象
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
             #exprs是一个将元素转化为表达矩阵的函数，此时得到的exprset是一个matrix。
dim(dat)     #看一下dat这个矩阵的维度
# ？expressionset 可以查看相关部分的性质及介绍

质控表达矩阵是否符合下游分析。

samples=sampleNames(ob) # sample name就是看有多少GSM样本
pdata=pData(ob) 
# pData是看样本怎么分组，相应的生物学分类，下载整理好后的包 -GEOmeta？
# pdata此时是个data.frame
group_list=as.character(pdata[,2])
dim(exprset) #查看exprset矩阵的维度，有33279行（基因表达量），6列（样本数）
exprset[1:5,1:5]

一个是检测samplename里是否一致，其次，pData转换格式以及获取分组信息，相当于直接将ob转化成为了data.frame格式。后续可以dim检测维度，或者将你认为有意义的列提出来，转化为字符或者怎样。

---

【方法3】直接处理raw data

从GEO数据库下载矩阵数据-可以直接进行下游分析 | 生信菜鸟团 http://www.bio-info-trainee.com/941.html公众号可以搜索：从GEO数据库下载表达矩阵 一文搞定。

旧版：affy包

readaffy读取原始文件

新版：oligo包

不同公司测序的读取的包也要不一样，具体可以见公众号中的总结

从GEO数据库下载得到表达矩阵 一文就够https://mp.weixin.qq.com/s/HoRUzx0UJxgkgQDxouj8rw