GSEA教程

写在前面的话:GSEA对于生物信息学技能的要求不高,以至于不需要太多的生信知识就能学会,这也成为了生物研究生的必备技能,如果你还不会的话,这篇教程将让你100%学会,只要你从头到尾跟着我的步骤走。

我按照网上的教程学习GSEA的时候绕了不少弯路,参考了很多不同的教程,并且不断地Debug,花了一周时间才跑出结果。之后我再做GSEA就如行云流水般顺畅。这篇教程,保证看完包会,文末我甚至给出了可以直接拿来跑的文件,不耽误大家时间,100%能重复出结果!

下面进入正题:

GSEA(Gene Set Enrichment Analysis),基因集富集分析,与GO分析类似,是一种基因差异表达分析。但是在原理上,与GO分析有很大区别。GO分析是针对差异基因,即看差异表达的基因在哪些GO分类或者pathway中富集,从而说明这些差异基因影响了哪些pathway。GSEA分析是针对所有基因的表达(不能用差异基因来跑GSEA,否则GSEA就失去了意义),这些基因的表达矩阵称为dataset。dataset一般是经过处理的RNA测序结果的基因表达矩阵,或者是从数据库中下载得到的表达矩阵。这些表达矩阵必须手动分组,比如分为control组与treat组。然后GSEA会自动分析出control组与treat组之间每个基因的差异并按照差异大小给每个基因算分,然后GSEA会根据这些基因的分值以及他们在某些基因集(geneset)中的富集程度来计算geneset富集分数(EnrichmentScore)并绘图(geneset由GSEA软件自带也可自行到数据库下载(后面由链接)),这样就得到了geneset的富集信息,从而推断出control组与treat组究竟在哪些geneset所代表的通路中存在差异。以上为GSEA简单的原理,算法什么的我不太懂,也不需要去了解,只要会用GSEA进行数据分析就可以了。

一. GSEA的前期准备

1.电脑一台,华为matebook,win10,8G内存。

2.GSEA软件:

http://software.broadinstitute.org/gsea/downloads.jsp

右边选择内存大小,最好选小于电脑内存的最大内存,比如我电脑8G内存,我就选择4GB(for 64-bit Java only)。内存大小决定了GSEA能处理的数据大小以及计算的精细程度。这里有两个坑我曾经栽过——请注意:第一是Java环境是必须的,请到java官网下载64位java8 (https://www.java.com/zh_CN/,官网只会推荐最适合你系统的java,没有其他版本的下载选项,但这不一定是我们所需要的版本!)或者用我文末分享的网盘链接进行下载,GSEA软件不支持java9或更高版本。第二是:如果选择2GB以上内存的java需要64位系统,如果你的电脑是32位系统,请选择1GB,win7-8系统一般会安装64与32两个系统,而win10不再有64与32位的区别,win10用户可放心下载。点launch进行下载,下载完成后可直接运行,如无法运行或没有相应请检查java是否安装正确。

下载完成后的图标为下图所示,如果图标未显示,请检查java是否安装正确。


这样,软件就下载完成了。

3. 实验数据(dataset)的准备。

NCBI的GEO数据库收录的绝大多数已发表文章的测序数据,而且下载起来很方便。本次使用的是人肝脏单细胞测序的数据,GEO编号GSE115469。

由于该测序数据集太大,EXECEL无法打开,需要用R语言进行整理,所以我将样本数缩小为100个,基因缩小为5000个,重新构建了一个表达矩阵,可以直接在EXCEL里进行排版(注意啊,我是方便教大家用excel操作才这样搞的,实际操作千万不能随意缩小数据啊,如果有几万基因几十个样本就用R处理)。这个矩阵我放在网盘里,文件名“human liver dataset”,文末有下载链接。打开之后如下图所示,第一列是基因名,第一行是样本名。这种基因名(ID)的专业术语叫做Gene symbol,我们平时接触到的Actin等基因名都是Gene symbol,此外还有很多其他的基因名,比如Ensemble ID,不同的测序平台对基因的编号也不同,这里不详细讨论Gene ID的转换问题,本篇主要带大家重复一遍GSEA的流程。下图是normalize后的单细胞测序表达矩阵,可用作GSEA的dataset。事实上,如果不是单细胞的数据,普通转录组测序的数据样本数一般只有10个以内,基因数至少2万个,完全可以用EXCEL进行处理。

对dataset的下一步处理是分组。我们这次就研究一下APOA2基因高表达的细胞与低表达的细胞有什么差异吧,首先我们将数据按照APOA2的表达水平排序。搜索到APOA2,在1454行

选择除基因名外的所有列,然后点击工具栏的“排序”

在新的对话框中,点击选项,选择按行排序,确定。

选择主要关键字“更多行”

点击1454行标,点击确定。

然后升序排列,点击确定。

排序后的APOA2

我们在心中把前10列定义为APOA2low,把后10列定义为APOA2high,保存EXCEL。

下一步,将.csv格式的ECXEL转换为.gct格式——可被GSEA软件识别的格式。这里使用既可以绘制热图又可以转换格式的在线软件:Morpheus。

https://software.broadinstitute.org/morpheus/

界面如下

然后上传上文处理好的表达矩阵。

点OK,会出现一个热图,无视这个热图,对我们没有帮助,我们是利用这个软件转换格式而不是画热图。然后点击工具栏File-Save dataset

命名文件,切记选择GCT1.2, 点OK,然后保存文件到本地

到这里我们对数据集(dataset)的处理就全部完成了。接下来进行最后一步准备工作。

4. 数据的Annotation文件的准备。

打开excel,输入

第一行20表示20个样本,2表示分为2组,1为识别字符(必须为1)

第二行表示分组信息,分为APOA2low,APOAhigh两组,注意#号不要丢。

第三行对应于dataset中每一列的分组,前10列输入APOA2low,后10列输入APOA2high。

然后点击保存,文件名结尾加上.cls,可被GSEA识别,储存为制表符分隔的文本文件(.txt)

打开文件检查一下,应该是这样的。

至此,GSEA前期准备就完成了。

下一步运行GSEA软件。

二. GSEA软件的运行

上传文件,点击load data, 使用Method1依次上传.gct与.cls结尾的txt文件。

成功会有以下提示,如果失败则可能是文件格式不对或者数据的排列不对。

然后我们成功上传了两个文件,如下图。

接下来我们设置GSEA output的储存位置,file-preferences。根据自己习惯设置。

2.设置参数。

我们进入run GSEA界面

第一行选择上传的dataset,第二行genesets选择c5.bp(biological process)。genesets详细介绍见MsigDB:http://software.broadinstitute.org/gsea/msigdb.

下一个选项,permutations内存够大CPU够强,就选1000,这里选100。越大越好,越大结果越精确。

phenotype二选一均可。

这里选false,因为我们的基因名已经是gene symbol了,不需要转换。

然后选phenotype,最后一行不选,如果基因名是测序平台的ID的话最后一行要选对应的测序平台并将上面的false改成ture

最后点击最下面的run.

成功的提示:

三. 结果解读

打开储存的文件夹,可以看到很多GSEA结果图

值得注意的是一个名为index的网页文件

打开它我们看以看到有446个基因sets在APOA2low中上调,1008个在APOA2high中上调,

点击enrichment results in html,得到下图的列表,GS为基因集的名字,SIZE代表该基因集下的基因总数,ES代表Enrichment score, NES代表归一化后的Enrichment score, NOM p-val代表pvalue,表征富集结果的可信度,FDR q-val代表qvalue, 是多重假设检验矫正后的p值,注意GSEA采用pvalue < 5%, qvalue < 25% 对结果进行过滤。Rank at max表示ES分数达到最大值时最后一个基因所处的位置。说的通俗一点,比如说在我们此次GSEA的5000个基因里,有3个基因对于对于cell cycle通路的ES分数有正贡献,三个基因分别位于第1,2,3位(本文开头说过,GSEA会根据基因表达差异大小对基因进行排序打分),他们对于ES分数的贡献分别为0.1,0.01,0.001,ES分数在第3位达到最大值0.111, 此时的Rank at max就为3,第三位往后的基因对cell cycle通路的ES分数是负贡献。

GSEA的结果图:

峰值的顶点为最大富集分数,下面的黑色竖条表示基因的位置,在最大富集分数对应横坐标位置左侧的黑线表示leading edge subset,意思是这些gene的排序靠前,对ES分数有主要贡献。总而言之,这张图主要是为了展示GSEA结果,说明我们确实通过GSEA分析发现了low和high两个表型在这条通路存在差异,这张图上的其他参数不必细究。

补充一下:

1. GSEA软件自带的geneset只有人的,如果是鼠的测序数据则要用鼠的geneset,把它下载下来再导入GSEA,很遗憾,我没有找到鼠的geneset。或者将鼠的基因转换为人的进行GSEA分析,这种转换方法我后续会进行介绍。

2. GSEA软件还有很多可供调整的参数,文中没有提到,还有一些其他功能,如leading edge analysis,这些功能的介绍我放在了文末的云盘中。

最后祝大家都能跑出结果!


资料链接:

数据集(dataset)以及可直接用来跑的的.cls与.gct文件,不过我还是建议用dataset从头到尾自己分析一遍:

链接:https://pan.baidu.com/s/1mTlXDXSUe9CRMZKu8P74uw

提取码:oyl2

java64位:

链接:https://pan.baidu.com/s/1qXoT_ocDA1F8zWlxbIJpRg

提取码:i1xz

参考资料与深度学习资料:

链接:https://pan.baidu.com/s/1jNc3Ho8lRoONuKyCDbgJHw

提取码:lvds


最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 161,780评论 4 369
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 68,424评论 1 305
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 111,397评论 0 254
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 44,576评论 0 218
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,997评论 3 295
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,945评论 1 224
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 32,107评论 2 317
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,850评论 0 208
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,625评论 1 250
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,804评论 2 253
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,285评论 1 265
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,613评论 3 261
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,291评论 3 242
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,164评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,963评论 0 201
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 36,096评论 2 285
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,886评论 2 278

推荐阅读更多精彩内容