查看原文
其他

一万人陪你学习GEO数据库挖掘知识(公益视频听课笔记分享)

生信技能树 生信技能树 2022-06-06

耗费半年的时间精心制作了成套的生物信息学入门视频教程,并且在生信技能树联盟平台发布了这个长达74个小时全套生物信息学入门视频:生信技能树视频课程学习路径,这么好的视频还免费!在B站看了看,大家学的热火朝天, 接下来我们就一个个知识点进行专题介绍,主要是一些优秀学生的笔记分享,希望大家在学习的过程中也能吸收到我传达的学习经验,人生感悟,只要你发给我笔记(邮箱 jmzeng1314@163.com),就有惊喜!专题历史目录:3个学生的linux视频学习笔记
生信人应该这样学R语言系列视频学习心得笔记分享
接下来介绍GEO数据库挖掘:

【生信技能树】GEO数据库挖掘


首先是厦门大学小板蓝根观看完生信技能树GEO数据库挖掘系列视频后对每讲内容的归纳总结。

1-通用文献阅读及其规律

我们的目标是要从读懂文献到复刻文献实验,再到掌握GEO数据挖掘的能力。首先便是要广泛阅读,在读文献时,提炼脉络,读懂文献使用了哪个或哪些GSE数据集,对数据做了哪些处理。了解清楚后,便可下载相应的数据集,得到表达矩阵,作差异分析,注释等一系列下游分析。

2-了解GEO数据库

GEO数据库的基本介绍:GEO数据库起初是为表达芯片数据准备的,在此之前可以了解一下芯片的基础知识。介绍GSE, GSM, GPL, GDS的含义。引用教程里的一段话,详细讲解了这四者之间的关系。

“一篇文章可以有一个或多个GSE数据集,一个GSE里可以有一个或多个GSM样本。多个研究的GSM样本可以根据研究目的整合为一个GDS,每个数据集有着自己对应的芯片平台(GPL)。”

再补充一点,一个GSE里可能有多个平台测出的数据。

3-数据下载的3种方法

  • GEO主页下载原始数据(RAW.tar)

下载下来是CEL格式,需要自己进行预处理。hgu 95系列和133系列的芯片常用affy包中ReadAffy函数进行读取,有些平台用affy包处理不了,例如[HuGene-1_1-st] Affymetrix Human Gene 1.1 ST Array这个平台,就需要oligo包read.celfiles函数进行读取。illumina的芯片用lumi包来处理。之后可以进行rma或mas5进行归一化操作。

  • GEO主页下载series_matrix.txt文件

exprSet = read.table('GSE42872_series_matrix.txt.gz', sep = '\t', fill = T, comment.char = "!", header = T)
  • GEOquery包直接下载表达量

library(GEOquery)
gset = getGEO('GSE42872'#返回一个list
exprSet = exprs(gset[[1]])

getGEO这个函数,输入不同的参数,下载的东西不一样。输入参数是GDS号时,下载soft文件,参数是GPL号时,下载芯片设计信息,参数是GSE号时,下载series_matrix.txt.gz文件,返回的是ExpressionSet对象,需要掌握geneNames,sampleNames,pData,exprs等对ExpressionSet对象操作的函数

我自己在linux和windows底下用GEOquery包下载表达谱时,会出现得到的表达谱少了第一行(第一个探针对应的表达量奇迹般消失)的情况,咨询了jimmy老师后,发现是readr这个包的问题,详细解决办法见http://www.bio-info-trainee.com/3719.html。

4-ID转换技巧大全

从GEO上下载的表达谱的行名是probe_id探针名,但是不同的平台,探针名不同,我们也无法直观地知道某个样本在某个探针上的表达量是那个基因的表达量,于是就需要将探针名转换为大家公认的NCBI的entrez ID,或者HUGO组织规定的gene symbol以便于后续分析。于是,我们要根据不同的GPL找到该芯片平台有对应的bioconductor注释包来找到探针与基因的对应关系,再进行转换。这里会遇到,“一个探针对应着多个基因”或者“一个基因对应多个探针”或者“探针没有对应基因”的情况,这就需要过滤整合表达矩阵,处理方法不尽相同。

5-了解你的表达矩阵

表达矩阵描述的就是各个基因在各个样本上的表达量。这讲主要是表达矩阵的可视化。无论是芯片表达数据或是转录组高通量测序数据,下载完表达谱需要根据生物学背景验证一下表达谱是不是正确的。只有确定了所得表达谱是正确的,之后的差异分析等一系列分析手段才是有意义的。这里提到的方法是看看管家基因的表达量是不是在表达谱中处于高表达。也可以用boxplot看看每个样本的表达量分布图,看看是否有批次效应等等。这里就需要去了解一些画图函数的使用方法。

6-差异分析

差异分析的目的就是找到在不同分组之间表达差别较大的基因。 得到表达矩阵及分组信息后,便可开始用limma包做差异分析。 这个包需要三个数据:表达矩阵,分组矩阵,差异比较矩阵。再经过lmFit,eBayes,topTable三个步骤便可得到差异分析结果。关键在于用手头上的表达矩阵和分组信息去生成分组矩阵和差异比较矩阵,得到的差异分析结果,可以用火山图及热图进行可视化。差异分析的包不止limma这一个,不变的做法就是看说明书,知道每个步骤的输入输出,生成相应的输入即可。

7-8-KEGG-GO等数据库的注释及GSEA分析

介绍基因的注释及富集分析。差异分析通过自定义的阈值挑选了有统计学显著的基因列表后我们其实是需要对它们进行注释才能了解其功能,最常见的就是GO/KEGG数据库注释,当然也可以使用Reactome和Msigdb数据库来进行注释。而最常见的注释方法就是超几何分布检验。当然还有其他的注释方法。超几何分布检验,运用到通路的富集概念就是“总共有多少基因(这个地方值得注意,主流认为只考虑那些在KEGG等数据库注释的背景基因),你的通路有多少基因,你的通路被抽中了多少基因(在差异基因里面属于你的通路的基因)。”目的就是知道,哪些通路中的哪些基因的表达因为药物或者某些操作的作用发生了较大的变化,导致通路有较大改变。

本讲涉及到clusterProfiler这个包,通过阅读说明书掌握包里面的函数包括函数所需的参数的能力是举一反三的关键。例如,KEGG输入的基因是EntrezID,在此之前需要进行转换。这些需要注意的点,都可以通过阅读说明得到。

GSEA是另一个常用的富集分析方法,目的是看看基因全局表达量的变化是否有某些特定的基因集合的倾向性。使用这些包或函数的方法都是看说明,把自己的数据做成软件要求的那种格式。

9-收尾的几点建议

这里包含了视频中的建议和个人学习之后的一些体会。

  • 学会创建Rproject来组织项目,通过打开Rproject来打开Rstudio。

  • 学会用if(F){},if(T){}来对代码进行折叠,清理视图。

  • 多看文献,了解进阶分析方法。

  • 分析时有意识地按步骤整理好代码,以便后续重复使用,提高效率。(这一点深有体会)

  • 学好R语言是GEO数据挖掘以及一切分析的基础,只有在懂得R语言的基础上,看懂代码,才能在有idea的时候,迅速实现。

  • 对于初学者,需要知道自己想实现某种功能时,应该上哪去找代码,对代码的归纳整理十分重要,绝大部分问题的解决办法和代码基本上都能在jimmy老师的github和生信技能树论坛找到。

10-批量生存分析

谈生存分析,就是说的KM方法估计生存函数,并且画出生存曲线,然后还可以根据感兴趣的基因表达量的高低进行分组,检验一下它们的生存曲线是否有显著的差异!这里介绍了一个做生存分析的R包survival,具体用法详见jimmy老师github上的代码:https://github.com/jmzeng1314/GEO

然后是南医大吴同学的:

P1.绪论--GEO数据挖掘是生物信息学一个重要部分,
将其作为入门可以更快的了解生信操作和背后逻辑。这一讲为系列课程的绪论,
整体介绍GEO数据挖掘的具体步骤,总结数据分析的主要流程和关键问题。
最后是一个完整过程示范,相信大家遇到的问题都可以在后续的教程中找到答案。

P2.从文献中学习是掌握数据挖掘思路最直接的办法,这一讲教程主要是如何快速的阅读文献
用最快的方式找到文献中的关键信息,从而完整get到文献的流程图。通过大量的文献阅读,
就可以知道哪些数据挖掘思路比较高级,更易发高分文章。

P3.如何找到我们需要的数据,这一讲会和大家详细介绍GEO数据库以及如何准确查找我们需要
的数据集。GEO数据有特殊的命名规则,平台号、数据集号、样本号、探针号,了解这些名称的含义
是有效利用GEO数据库的关键。

P4.GEO数据下载的三种方式,包括1.直接从官网上下载原始数据rawdata 2.下载作者预处理过的表达矩阵 
3.直接从R中下载表达矩阵。需要注意的是,第三种方法由于存在网络墙,国内下载时会出现不稳定的现象,
可以尝试从不同的镜像进行下载。

P5.进行ID转换的办法很多,这一讲给大家介绍非常实用的一种-----采用官方包进行ID转换。
每个芯片平台在Bioconductor上都有相应的注释包,一般以db后缀。ID转换中,对于探针和基因多对一
的情况,一般采用取平均值或者最大值的办法。

P6.表达矩阵精进篇,这一讲主要介绍如何对表达矩阵进行质量评估,包括聚类分析,PCA分析,或管家基因进行
分组间的比较。一个简单的箱线图一般就可以判断样品间的表达基线水平。

P7.获取差异基因一般是进行GEO数据挖掘的第一步,这一讲介绍最常用的limma包进行差异分析,对于
已经采用对数化预处理的表达矩阵,根据分组进行t检验计算组间差,从而获得分组间的差异情况,
这就是limma包的原理。获取的差值为LogFC,是差异倍数FC以2为底的对数。

P8.富集分析包括对差异基因的GO注释,KEGG注释和GSEA注释,这一讲具体介绍这3种注释的使用办法。
GO的input文件一般是genesymbol,而KEGG的input文件一般是geneENTREZID。输入文件的标准化,
再进行下一步分析时,可以避免报错。GSEA分析略有不同,不需要差异基因,只需要输入表达矩阵即可。

P9.这一讲和大家分享一些项目管理的小Tips,好的操作习惯事倍功半!

P10.一个小惊喜-----生存分析,可以给文章加分的数据处理办法。这一讲包括进行生存分析
常用的代码和一些做图时的调整办法,好的数据集是进行生存分析的必要条件。

接着分享 linchyu 【生信技能树】GEO数据挖掘导学 --

Part1:生物信息学资源、学习方法介绍,GEO数据挖掘简要浏览

生物信息学资源的获取方式有很多,比如在生信技能树博客中,就基本囊括了你想要的所有知识点。所以,如果你有任何问题,可先自行尝试在生信技能树(http://www.biotrainee.com/)中独自解决。另外,教程中所涉及的数据及代码均可在github网站(https://github.com/jmzeng1314)免费获取。当然,作为新手的话,还是建议自己按照教程手敲代码,以加深理解。

Part2:阅读文献的技巧

通过阅读文献,快速提取我们所需要的信息,包括但不限于数据的GEO号,文章的分析流程,log2FoldChange阈值的设定、P值的设定,文章所使用的R包。
生信菜鸟团网页中生信人的20题介绍,建议有时间的同学可以尝试着做一做。

Part3:熟悉GEO数据库

登录GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)官网主页,我们可以简单了解该平台所收录的数据数量及测序平台数量,也可以了解一下GEO数据库保存数据的方式;当然,如果有兴趣的同学,可以进一步探索各板块的内容,代表的意义;也带领大家探索了常用的芯片平台。当然,这里更多的是点到为止,更多更细的知识点是需要自己用心去探索的哟。

Part4:GEO数据下载

数据的三种下载方式:①直接下载raw数据(对新手来说后续处理太复杂,不推荐);②下载matrix.txt文件(容易受网速的限制,次推荐);③使用getGEO下载(简单易学,墙裂推荐)。getGEO下载方法:

library(GEOquery) 
eSet <- getGEO("GSE42872",destdir = ".",AnnotGPL = F,getGPL = F
save(eSet,file = "GSE42872_eSet.Rdata")

如果上述方式都不行(可能性很小的啦),只能这样:旧版的affymetrix芯片:ReadAffy;新版的affymetrix芯片:read.celfiles;Illumine相关芯片:lumiR.batch。

Part5:ID转换

由探针ID转换为pubmed所存储的通用ID:方法也很简单,我们只需要下载包含有探针ID与基因symbol关系的注释文件,将他们对应起来即可。
首先,将下载的数据读入R,提取下载文件的表达矩阵:

ob <- eSet[[1]] 
exprSet <- exprs(ob) 
samples <- sampleNames(ob) 
pdata <- pData(ob) 
group_list <- as.character(pdata[,2]) 

其次,下载GPL注释信息:①与上面类似,使用getGEO下载:

    library(GEOquery)
    gpl <- getGEO("GPL16570",destdir = "."
    probe2gene <- Table(gpl)[,c(1,11)]
    save(probe2gene,file = "probe2gene.Rdata")

②如果可以找到该平台所对应的R包(http://www.bio-info-trainee.com/1399.html),则可以这样下载:

    library(BiocInstaller)
    BiocInstaller::biocLite("hugene10sttranscriptcluster.db"
library(hugene10sttranscriptcluster.db) 

最后,提取注释信息中探针ID与基因symbol的对应关系。

Part6:表达矩阵的介绍

首先,需要检验表达矩阵在经过我们的多次处理之后,基因symbol是否与其表达量存在明显不合理的地方,即是否出现基因symbol与基因表达量之间的不对应:
① 查看内参表达量:exprSet[“GAPDH”,]、boxplot(exprSet[“GAPDH”,])
② hcluster图:

colnames(new_exprSet) <- paste(group_list,1:ncol(new_exprSet)) 
nodepar <- list(lab.cex=0.6,pch=c(NA,19),
                cex=0.7,col="blue")
hc <- hclust(dist(t(new_exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc),nodePar=nodepar,horiz=T)

③ PCA图:

library(ggfortify)
library(ggplot2)
df=as.data.frame(t(new_exprSet))
df$group=group_list
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data
=df,colour = 'group')

设置分组信息,方便我们后续使用:group_list

group_list <- factor(c(rep("control",3), rep("experiment",3)))  

Part7:差异分析

使用limma包进行差异基因的分析:首先需要准备表达矩阵和分组信息,使得我们的数据满足limma包的要求,进而计算p值:

suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(new_exprSet)
contrast.matrix <- makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
##step1
fit <- lmFit(new_exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2) ## default no trend !!!
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)

画热图:

library(pheatmap)
choose_gene=names(tail(sort(apply(new_exprSet,1,mad)),50))
choose_matrix=new_exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

Part8:火山图、散点图等的制作,高阶美化

在有一定R语言背景下,使用R绘图比较简单,就是一串代码的问题,如果一串不行,那就再来一串。当然,各种绘图代码不用我们自己写,直接谷歌搜索火山图,就会出现相对应的代码(如http://www.sthda.com/english/),然后你在其基础上修改你自己的参数即可。

Part9:GO-KEGG功能注释

当你找到差异基因后,可以继续对其进行功能注释,以便于了解该基因是在哪条通路中发挥作用,方便从生物学角度进行解释。
超几何分布检验(GSEA):上述的功能注释使用的是p值或矫正p值,但由于现在出现一些声音开始质疑p值的临床意义;且由于生物体不同组织的差异,机械的设定阈值筛选差异基因可能存在一定的假阴/阳性(比如,在脑的某些区域,神经递质的变化可能并没有达到设定阈值,但却能表现出临床差异)。因此,GSEA的存在则有一定的必要性(想要详细了解GSEA原理,请转至https://www.jianshu.com/p/b409a5576ce1)。

Part10:来自jimmy的进阶建议

当你自己跑完整个流程,建议分类将代码保存,方便后续寻找和稍加修改后使用,便于简化流程;另外,就如在part1所提到的一样,如果你需要更多的R代码,可以在github网站(https://github.com/jmzeng1314)进行资源获取。

Part11:最后的彩蛋

走过路过,不拿白不拿的生存分析代码免费送:GitHub(https://github.com/jmzeng1314)。

最后是杰根妹妹 的,原文链接:https://www.jianshu.com/p/13d4289d5bc9

1.项目总览及Github介绍

介绍整个项目
简介Github,注册,下载其中代码

  • R语言用文件夹+project方式组织,定位所有数据和代码

2.通用文献阅读及规律

  • identification of the interaction network of hub genes for melanoma treated with vemurafenib based on microarray data
    文献导读,注意文章中找到差异性的方法
    例子中是p值<0.01,|logFC|≥2

  • 差异基因要通过阈值控制,~200多差异比较正常
    -至少看20篇相关文章,提炼脉络,选择GSE- 表达矩阵-差异分析-5大数据库的注释-PPI等网络

了解GEO数据库(生新技能树公众号,解读GEO)

GSE号-修改URL即可到数据库
refseq_id, GEOquery
芯片基础知识(生信技能树论坛)
HG-U133_Plus_2(经典芯片)

3. 数据下载的3种方式

  • ①下载rawdata(不推荐)

  • ②下载表达矩阵(matrix)

  • ③ R语言直接读取GSE号 (GEOquery)
    getGEO("GSE42549", GSEMatrix = TRUE, AnnotGPL = FALSE, getGPL= FALSE)

  • 不同芯片用不同的R包

4.ID转换技巧大全

downGSE

  • geneID,探针和基因不是一一对应的,且基因本身就是多种多样(entrez ID和symbol是最重要的)

  • ID转换,library(hgu95av2.db )
    不同平台对应不同R包,可谷歌

  • ID转换实操,

5. 了解你的表达矩阵

实操,跑代码,了解PCA,hclust图等

6. 差异分析

limma对芯片数据做差异分析
需要

  • 表达矩阵

  • 分组矩阵

  • 差异比较矩阵
    实现步骤

  • lmFit

  • eBayes

  • topTable



■   ■   ■




生信基础知识大全系列:生信基础知识100讲   

史上最强的生信自学环境准备课来啦!! 7次改版,11节课程,14K的讲稿,30个夜晚打磨,100页PPT的课程。   

如果需要组装自己的服务器;代办生物信息学服务器

如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等),点我?

如果需要线下辅导及培训,看招学徒 

如果需要个人电脑:个人计算机推荐

如果需要置办生物信息学书籍,看:生信人必备书单

如果需要实习岗位:实习职位发布

如果需要售后:点我

如果需要入门资料大全:点我

欢迎大家留言选出最优秀笔记

 1. 厦门大学小板蓝根

 2.南医大吴同学

 3.linchyu

 4. 杰根妹妹 

也欢迎提交自己的学习笔记,我会记住你的

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存