精华内容
下载资源
问答
  • GSEA_Win_4.0.2-installer.exe

    2019-10-21 13:55:23
    GSEA富集分析,1、准备三个文件第一行:#1.2,表示版本号,自己准备文件时照抄就行; 第二行:两个数分别表示gene NAME的数量和样本数量(矩阵列数-2); 矩阵:第一列是NAME;第二列Description,没有的话可以全用...
  • 盖西 用于基因组富集分析的R包
  • GSEA

    千次阅读 2020-04-26 17:07:46
    GSEA分析数据准备KEGG与GO分析GSEA绘图分析 数据准备 我们需要的是带有基因id的log2foldchgange的向量,我们从de_result里面提取出来就行了 de_result$entrezid <-mapIds(org.Mm.eg.db, #.db是这个芯片数据对应的...

    数据准备

    我们需要的是带有基因id的log2foldchgange的向量,我们从de_result里面提取出来就行了

    de_result$entrezid <-mapIds(org.Mm.eg.db, #.db是这个芯片数据对应的注释包
                    keys=de_result$id,
                     column="ENTREZID", #clolumns参数是你要转换的ID类型是什么,只能选择一个。
                     keytype="ENSEMBL" )
    #先将de_result里面的ENSEMBL转换为entrezid
    
    geneList <- de_result$log2FoldChange
    #提取所有基因的log2FoldChange
    names(geneList) <- de_result$entrezid
    #用对应基因的entrezid对geneList命名
    geneList <- sort(geneList, decreasing = T)
    #将其按照降序排列
    

    GSEA分析

    1. KEGG
      这是kegg的富集分析
    require(enrichplot)
    kk <- gseKEGG(geneList, organism = "mmu")
    

    然后出现了如下报错

    preparing geneSet collections...
    GSEA analysis...
    no term enriched under specific pvalueCutoff...
    Warning message:
    In fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize,  :
      There are duplicate gene names, fgsea may produce unexpected results
    

    它说设置pvalueCutoff并且这里有重复的基因名,运行的话会产生意想不到的结果
    于是我就手动设置pvalueCutoff

    kk <- gseKEGG(geneList, organism = "mmu",pvalueCutoff = 0.05)
    GSEA analysis...
    leading edge analysis...
    done...
    Warning message:
    In fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize, :
    There are duplicate gene names, fgsea may produce unexpected results
    

    还是报错说,有重复的基因名,在我的努力下,发现这里的问题和我是一样的
    https://github.com/YuLab-SMU/clusterProfiler/issues/189,于是我按照llrs的fangfa,剔除重复的基因名

    geneList <- geneList[!duplicated(names(ggeneList))]
    

    然后再次运行

    kk <- gseKEGG(geneList, organism = "mmu",pvalueCutoff = 0.05)
    preparing geneSet collections...
    GSEA analysis...
    no term enriched under specific pvalueCutoff...
    

    我的阈值设置得太低,导致没有结果
    我们调整阈值,再来

    kk <- gseKEGG(geneList, organism = "mmu",pvalueCutoff = 1)
    preparing geneSet collections...
    GSEA analysis...
    leading edge analysis...
    done...
    

    现在就成功了!!!!
    然后我们开始画图

    gseKEGG_result<-as.data.frame(kk@result)
    #我们先把结果提取出来,方便查看
    gseaplot2(kk, "mmu04060")
    #然后进行绘图
    

    在这里插入图片描述

    1. GO
      GO富集分析又三种,MF,BP,CC
      虽然可以用命令将他们各自求出来,但是一般情况下,我们是将他们全部求出来
    kk <- gseGO(geneList,ont = "BP",OrgDb = org.Mm.eg.db)	
    one of "BP", "MF", and "CC" subontologies, or "ALL" for all three
    #将ont参数调整为all,就可以了
    #我们可以将结果提取出来。以供查看
    gseGO_result<-as.data.frame(kk@result)
    View(gseGO_result)
    

    这里我们用的是GO分析的方法,所以我们在绘图的时候是使用的GO,作为指定,而不是基因id了

    gseaplot2(kk, "GO:0001816")
    

    在这里插入图片描述
    参考:
    https://www.jianshu.com/p/a1927f75043b
    https://yulab-smu.github.io/clusterProfiler-book.

    展开全文
  • GSEA分析流程

    千次阅读 2021-05-16 09:25:46
    GSEA 是 2003 年提出来的一种对表达谱芯片进行分析的方法,并被编制成软件。它的主要目的就是确定预先定义的基 因集(具有相同或相似的功能,或位于同一染色体相邻位点的一群基因)在表达谱芯片结果中是否有显著性。 ...

    1. 简介
    GSEA 是 2003 年提出来的一种对表达谱芯片进行分析的方法,并被编制成软件。它的主要目的就是确定预先定义的基
    因集(具有相同或相似的功能,或位于同一染色体相邻位点的一群基因)在表达谱芯片结果中是否有显著性。
    GSEA 分析过程分为 5 步:
    1. 基因知识库的获得;
    2. 根据基因表达谱数据对所有基因进行排序;
    3. 计算富集得分(enrichment score,ES);
    4. 估计显著性水平;
    5. 进行多重假设检验。
    GSEA 能够鉴定疾病发生过程中潜在的及起决定作用的遗传改变或信号转导通路,从而把表达谱芯片数据和生物学意
    义偶联起来。一些研究者通过 GSEA 鉴定了糖尿病患者肌肉组织中氧化磷酸化信号通路的下调以及腺癌患者 K-RAS
    的下调等
    2. 下载安装(GSEA v3.0)
    下载地址 http://software.broadinstitute.org/gsea/downloads.jsp (注册后可免费下载)
    下载后安装打开(依托的 java 运行环境按照提示下载安装)
    初始界面如下:

    3. 数据准备
    3.1 数据要求:
    1. 可使用 Excel 或者某种文本编辑器创建或者编辑 GSEA 文件
    2. 文件中以制表符(table)隔开
    3. 输入文件名称不要包含连字符号(hypens, -)
    支持的文件格式类型
    1. Expression Data Formats(芯片表达谱数据集文件)
    (1) GCT: Gene Cluster Text file format (*.gct)
    (2) RES: ExpRESsion (with P and A calls) file format (*.res)
    (3) PCL: Stanford cDNA file format (*.pcl)
    (4) TXT: Text file format for expression dataset (*.txt)
    2. Phenotype Data Formats(表型数据文件)
    ( 1) CLS: Categorical (e.g tumor vs normal) class file format (*.cls)
    ( 2) CLS: Continuous (e.g time-series or gene profile) file format (*.cls)
    3. Gene Set Database Formats ( 功能基因集文件)
    ( 1) GMX: Gene MatriX file format (*.gmx)
    ( 2) GMT: Gene Matrix Transposed file format (*.gmt)
    ( 3) GRP: Gene set file format (*.grp)
    ( 4) XML: Molecular signature database file format (msigdb_*.xml)
    4. Microarray Chip Annotation Formats ( 芯片注释文件)
    ( 1) CHIP: Chip file format (*.chip)
    5. Ranked Gene Lists
    ( 1) RNK: Ranked list file format (*.rnk)
    3.2 主要的文件格式简介:
    1. GCT: Gene Cluster Text 文件格式 (*.gct) 芯片表达谱数据集文件
     

    第一行:版本信息(一般是#1.2);
    第二行:第二行为两个数值,第一个数值是基因数目,第二个数值是样本数目,中间以制表符( tab)隔开;
    第三行:数据表头 NAME:基因名称, Description: 基因相关描述或注释,其他列:各样本名称,中间以制表
    符( tab)隔开;
    第四行及以后:数据内容为:基因名称,基因相关描述或注释和基因在各样本表达值数据,中间以制表符
    ( tab)隔开。
    2. CLS: Categorical (e.g tumor vs normal) class file format (*.cls) 表型数据文件

    Cls 文件包含三行内容
    第一行为样本数目,分类数目和 1, 中间以空格隔开,也可以 tab 隔开;
    第二行为#, 类别一的名字, 类别二的名字……, 中间以空格隔开, 也可以 tab 隔开;
    第三行为各样本的类别名,顺序为 gct 文件中样本顺序, 中间以空格隔开,也可以 tab 隔开。
    前面为处理,后面为对照
    3. 功能基因集文件
    功能基因集文件即是位于 MSigDB 中的各个信息文件,可有.gmx 和.gmt 两种格式。当功能基因集文件较多
    时,使用.gmt 格式更便于保存,但当功能基因集文件数<256 时,应用.gmx 格式更利于发挥 EXCEL 的编辑优势。通
    常在 GSEA 网站可直接下载.gmt 格式的功能基因集文件使用,而无需自行创建。
    4. 芯片注释文件
    芯片注释文件格式为"*.chip",常用的芯片类型其注释文件都可在 GSEA 网站下载。
    虽然此文件不直接在 GSEA 算法中使用,但它用于注释输出结果,也可用于将表达数据集中的每个探针集
    合折叠为单个基因载体。
    其他文档有相应的别的要求,如果文档格式出错,会有完整的提示,注意文档格式的要求很严格,稍有不对就无法
    完整导入数据。各个文档格式的详细要求:
    http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
    4. RUN
    4.1 数据导入(点击左侧的 load data)

    4.2 RUN 环节
    成功导入数据后,可以进入 run 环节,点击 Run GSEA 进入下图所示页面:

    各个文件选择好后,可以命名分析名称,选择输出文件路径, cpu 使用情况等。一切就绪后,点击下方的 Run 按钮
    即可开始运行。
    注意: 1.输出路径默认是本地用户组用户文件目录下的 gsea_home/output/X 月 X 文件夹内(可更改)
    2.cpu 使用情况有 Low、 high 和 Normal 三个选项。选择代表 cpu 的使用情况, 高、中、低三种情况。
    运行具体过程可点击下图所示位置,可弹出相应应用信息

    4.3 结果生成
    观察左下角的 GSEA reports 框可观察结果产生情况,状态有三种: 1.Running (正在运行) 2.Error (出错) Success (成
    功)。

    点击 Running 后可选择调整 cpu 使用情况(low、 normal 或者 high)和暂停正在运行的任务;
    点击 Error 后会出现详细的报错信息,可以查看报错的详细信息或者帮助等;
    点击 Success 后跳转到一个 GSEA 报告页面。具体情况如下文所述;
    点击下方的 Show results folder 会跳转到默认输出目录中,如果更改了输出目录,需要自己手动寻找结果目录。
    5. 结果解读
    1) Error 也会在结果文件夹内产生一个 error 文件夹,里面存在两个文件, index.html 和 xtools.css。基本没什么用
    处,可直接删除。
    2) Success 会在结果文件夹中产生以 Analysis name.Gsea.数字代码命名的一个文件夹,里面有各种结果文件,直接
    打开其中的 index.html 文件可打开结果报告网页。或者直接点击软件左下角, GSEA reports 框中的相应分析结
    果后的 Success 也可直接跳转至结果报告页面。
    结果报告例子如下图所示:

    结果报告一般分为七个部分(cls 文件中包含两个类别)
    (1) Enrichment in phenotype: class 1 (sample num)

     类别 1 中,在研究的 m 个基因子集中,有 n 个基因子集的表达上调(n/m);
    75%的正确率支持其中的 x 个基因子集中的基因表达上调;
    有 y 条基因子集中的基因表达上调的名义 p 值小于 0.01;
    有 z 条基因子集中的基因表达上调的名义 p 值小于 0.05。
    (2) Enrichment in phenotype: class 2 (sample num)

    类别 2 中,在研究的 m 个基因子集中,有 n 个基因子集的表达上调(n/m);
    75%的正确率支持其中的 x 个基因子集中的基因表达上调;
    有 y 条基因子集中的基因表达上调的名义 p 值小于 0.01;
    有 z 条基因子集中的基因表达上调的名义 p 值小于 0.05。
    (3) Dataset details

    原始芯片数据包含探针数目以及这些探针代表的基因数目。
    (4) Gene set details

    在选择的数据库中,基因集中包含基因子集数目,在定义的选择标准下(功能基因子集中所含基因数大于 15小于 500),有多少基因子集被过滤,有多少基因子集用于分析。
    (5) Gene markers for the class 1 versus class 2 comparison

    (6) Global statistics and plots

    (7) Other
     

    分析所用的各个参数
    3) 部分结果详解
    点击结果报告中分组富集信息中的 Snapshot 超链接
    其结果中有四个关键的统计量值,分别是富集得分(enrichment score,ES),标准化富集得分(normalized enrichmentscore, NES),错误发现率(false discovery rate, FDR)和名义 P 值(nominal P value)。


     富集得分(enrichment score,ES)
    ES 是 GSEA 分析的原始结果, ES 反映 S 在秩列表 L 顶部过表达的程度。 设总基因个数为 N, L=[g1,g2,……gN] 是依据各基因与表型间相关性 r 排序的基因列表。 基因集 S 富集得分 ES(S)的计算基于基因列表 L,从 ES(S)=0 开始,遇
    到 S 中的基因时增加 ES(S),相反则减小,增加或减少的幅度依赖于基因与表型的相关性 r(𝑔𝑗) =𝑟𝑗、 基因集的大小M、总基因个数 N。
    根据基因列表 L,从 i=1 到 N,若 i∈S,则 ES(S)增加|𝑟𝑖|𝑝/ ∑𝑔∈𝑆 |𝑟𝑖|𝑝; 若 i∉S,则 ES(S)减小 1/N-M。 当 ES 值为正,表示某一功能基因集富集在排序序列的前方,当 ES 值为负,表示某一功能基因集富集在排序序列的后方。 富集得
    分取绝对值的最大的数值。

    图上部的曲线表示动态 ES 值,最高点表示此通路的 ES 值。图的中间部分表示杂交数据的排序序列, 黑色竖线表示在此通路中出现的基因。此外还有一个领头亚集(leading edge subset),领头亚集中的基因是指对 ES 值贡献最大的
    基因集合。当 ES 为正值时,领头亚集位于 ES 值对应排序序列之前,反之,则位于 ES 值对应排序序列之后。显然,领头亚集的出现说明一方面这些基因在通路中有富集,非散在分布,另一方面,说明这些基因在通路中有共同的表
    达趋势。显然,在 ES 图中出现领头亚集的形状的,表明这个功能基因集在定义的实验条件下具有更显著的生物学意义。此外,图下部曲线表示的是排序值。排序值沿芯片数据的排序序列由大到小分布,在分类表型数据文件,正
    值表示该基因与第一个表型相关,负值表示与第二个表型相关。在连续表型数据文件,如时间系列,正值表示相关而负值表示负相关或没有相关性。 在所有的基因中,如果出现一个基因属于这个组合并且表达量在对应组里面表达
    高于另一组,富集分数就增加,反之就下降。
    标准化富集得分(normalized enrichment score, NES)
    由于 ES 是根据分析的芯片数据集中的基因是否在一个功能基因集中出现来计算的,但各个功能基因集中所包含的基因数不同,而且不同的功能基因集与芯片数据间的相关性也不同,因此在比较芯片数据集在不同功能基因集中的
    富集程度时,需要对 ES 进行标准化处理。 GSEA 定义 NES 为: NES = (某一功能基因集的 ES)/(数据文件集所有随机组合得到的 ES 的平均值)NES 是功能基因集富集分析结果的主要统计量。NES 是建立在数据文件集所有随机组合得到的 ES 平均值的基础上,因此,数据随机组合方法,随机组合次数,或表达数据文件集大小的改变都会影响 NES。
    错误发现率(false discovery rate, FDR)
    FDR 描述的是一个估计的可能性,即当一个功能基因集的 NES 值确定后,判断其中可能包含的错误的阳性发现率。
    例如, FDR=25%意味着对此 NES 值的确定, 4 次中可能有 1 次是错误的。在 GSEA 的结果报告中,高亮显示了 FDR值小于 25%的富集功能基因集,因为从这些富集功能基因集中最可能产生有意义的科学假设及促进进一步深入研究。
    在大多数情况下,选择 FDR 值为 25%来判定是否是富集的功能基因集是合适的,因为通常用于分析的芯片表达数据之间,大部分都缺乏一致性,而且每次分析的功能基因集数目不多。但是,当分析的芯片数据集较小,分析时选择
    的是探针间的随机组合(gene-set permutation)而不是表型间的随机组合(phenotype permutation), P 值采用的严格度又不高时,应该选择更加严格的 FDR 界值,如 FDR=5%。一般而言, NES 的绝对值越大, FDR 的值就越小, 说
    明富集度越高的功能基因集,分析结果的可信度就越高。
    名义 P 值(nominal P value)
    名义 P 值描述的是针对某一个功能基因子集得到的富集得分的统计显著性,很显然, P 值越小说明基因的富集性越好。但在 GSEA 分析结果的四个参数中,只有 FDR 值进行了功能基因子集大小和多重假设检验的校正,而 P 值没
    有。因此,当分析结果中出现一个高度富集的功能基因子集,具有很小的名义 P 值但却有很大的 FDR 值时,往往意味着其实和其它功能基因子集相比较,它的富集并不是很显著。原因可能是用于分析的芯片数据样本量较少,或
    者是杂交信号微弱,又或者是选择的功能基因子集并没有很好地反映样本的生物学意义。在 GSEA 的报告中, P=0.0 意味着实际的 P 值小于 1/ 随机组合的次数。例如,当选择的随机组合数为 100 时,报
    告中的 P=0.0 即是说实际的 P 值小于 0.001。所以提高随机组合数,就能得到更精确的 P 值。一般运行 GSEA 时随机组合次数都选择 1000,但不要超过 1000,因为一旦超过 1000,运行 GSEA 程序可能会耗尽内存。
    结果详解(http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html?_Interpreting_GSEA_Results )

    6. 实例
    6.1 转录组项目
    1. 数据准备
    (1) 各样本表达量文件
    (2) 各样本的分类情况
    (3) *.gct 文件的构建
    a) AllSamples.GeneExpression.FPKM.xls 文件,其中包含各个样本的表达量值。 其中包括:
    gene_id、transcript_id(s)、sample1_FPKM、sample2_FPKM、sample4_FPKM… Symbol 、Cellular Component、
    Molecular Function 、 Biological Process、 Kegg Orthology、 Nr Description
    b) 从中获取*.gct 文件需求的各样品的表达量(sample_FPKM) 以及基因名(symbol), 可获取相关注释描述(如 Nr Description)。
    例如:
    提取前:

    gene_idtranscript_id(s)sample1_FPKMsample2_FPKMSymbolDescription
    XXXXXXXXXXXXXXXXXXXXXXXX/NA

    提取后:(注意中间为制表符或者空格)

    SymbolDescriptionsample1_FPKMsample2_FPKM
    XXXXXX/NAXXXXXXXXX

    c) 根据 gct 文件格式的要求,需要统计基因数目和样本数目(注意中间为制表符或者空格)

    #1.2 (版本信息一般是#1.2)
    Num(基因总数)Num(样本总数)
    SymbolDescriptionsample1sample2
    XXXXXX/NAXXXXXXXXX

    (4) *.cls
    根据 cls 文件格式的要求配置文件
    例如:
    说明 a、 中间以制表符或者空格隔开;
    b、第一行依次表示样本数、类别数、 1;
    c、第二行依次表示#、类别 1 的名字、类别 2 的名字;
    d、第三行为 sample1 的类别、 sample2 的类别、 sample3 的类别(注意顺序和 gct 文件中样本顺序一致。
    2. 运行
    (1) Load data
    (2) 进入 RUN 阶段


    3. 结果


    左下角的 GSEA reports 中出现 Success 5 时意味成功,点击 success 5 会出现结果的网络版报告。出现 Error!意味出错,点击 Error!可展示详细报错内容,根据报错信息,调整后跑出成功结果。
    最终结果保存在 RUN 步骤中第 6 步选择的 Save results in this folder 的文件夹中。
    6.2 RNA-Seq
    类似转录组重测序,从GeneExp 提取 all.gene.FPKM.xls
    文件,后构建 gct 文件和 cls 文件。后续步骤类比上述。
    6.3 注意事项
    1. 建议使用 symbol,根据 chip platform 选择 Expression dataset(*.gct)。
    2. Gene sets database 请根据需求选择。
    3. Collapse dataset to gene symbol
    如果折叠数据集(将数据集折叠到基因符号参数= True),则折叠数据集中的基因通过基因符号识别,因此基因集中的基因标识符必须是基因符号;
    如果不折叠数据集(将数据集折叠到基因符号参数= False),则基因集中的基因标识符必须与表达式数据集中的基因标识符相同。
    4. 左下角的 GSEA reports 中出现 Error! 后请点击 Error! 后查看详细的报错信息和由软件提供的解决方案,更改错误后重新运行!


    5.根据样本数目选择对应的 Metric for ranking genes
    7. 附录
    7.1 Molecular Signatures Database v6.1(http://software.broadinstitute.org/gsea/msigdb/index.jsp)
    分子标记数据库(MSigDB)中的 17786 基因集分为 8 个主要集合和几个子集。
    1) H: hallmark gene sets(browse 50 gene sets)
    Hallmark 基因集合总结和代表特定的明确定义的生物状态或过程,并显示相关的表达。这些基因组基于鉴定基因组重叠和保留显示坐标表达的基因的计算方法产生。 Hallmark 减少噪声和冗余,并为 GSEA 提供更好的划定的生物空间。
    To cite your use of the collection, and for further information, please refer to Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell
    Syst. 2015 Dec 23;1(6):417-425.
    2) C1: positional gene sets (browse 326 gene sets)
    对应于每个人染色体和具有至少一个基因的每个细胞遗传带的基因集。 (细胞遗传位置从 HUGO 和 Unigene解析。当存在冲突时,使用 Unigene 条目。)这些基因组有助于鉴定与染色体缺失或扩增,剂量补偿,表观遗传沉默, 和其他区域效应。
    3) C2: curated gene sets (browse 4738 gene sets)
    此集合中的基因集来自以下方面:
    在线途径数据库:代表代谢和信号通路的基因集从此处列出的在线途径数据库导入。
    生物医学文献:在过去几年中,微阵列研究已经鉴定了几种重要的生物和临床状态(例如癌症转移,干细胞特征,耐药性)的标志。许多具有这些标志的基因集,最初是作为表格出版在论文中的,可提取出来作为 GSEA 分析的基因集。因此,其编辑了具有公开的基因表达标志的微阵列文章的列表,并且从每篇文章中,从主文本或补充信息的表中提取一个或多个基因组。目前,该集合包括来自超过 340 篇 PubMed 文章的基因集。其正在努力创建一个更自动化的方法来从文献中制定基因集。
    L2L:从公开的哺乳动物微阵列研究中汇编的基因集MYC 靶基因数据库:来自约翰霍普金斯大学医学院 MYC 靶基因数据库的 Chi Dang 博士策划的基因组。
    ( a) CGP: chemical and genetic perturbations(browse 3409 gene sets)
    代表遗传和化学变化表达标志的基因集。 许多这些基因集成对存在,例如: xxx_UP( xxx_DN)基因集代表由变化诱导(抑制)的基因。 每个基因集的基因组页面列出了它所基于的 PubMed 引文。
    ( b) CP: Canonical pathways(browse 1329 gene sets)
    来自通路数据库的基因集。 通常,这些基因集是由领域专家编译的生物过程的规范表示。
    ( c) CP:BIOCARTA: BioCarta gene sets(browse 217 gene sets)
    来自 BioCarta 途径数据库的基因集 http://www.genecarta.com/
    ( d) CP:KEGG: KEGG gene sets(browse 186 gene sets)
    来自 KEGG 途径数据库的基因集 http://www.genome.jp/kegg/pathway.html
    ( e) CP:REACTOME: Reactome gene sets(browse 674 gene sets)
    来自 Reactome 途径数据库的基因集 http://www.reactome.org/
    4) C3: motif gene sets(browse 836 gene sets)
    包含共享在人,小鼠,大鼠和狗基因组中保守的顺式调节元件的基因集。 这些基序被编目,同时这些基序表示启动子和 3'-UTR 中的已知或可能的调节元件。 这些基因组使得可能与微阵列实验中的变化与保守的推定的顺式调节元件相联系。
    ( a) MIR: microRNA targets(browse 221 gene sets)
    含有共享 3'-UTR 微小 RNA 结合基序的基因的基因集
    ( b) TFT: transcription factor targets(browse 615 gene sets)
    含有共享在 TRANSFAC (version 7.4, http://www.gene-regulation.com/ )数据库中定义的转录因子结合位点的基因的基因集。 这些基因组中的每一个由 TRANSFAC 记录注释。
    5) C4: computational gene sets(browse 858 gene sets)
    计算基因集是通过挖掘癌症导向微阵列数据结果的集合。
    ( a) CGN: cancer gene neighborhoods(browse 427 gene sets)
    由以 380 个癌症相关基因为中心的表达邻域定义的基因集合(Brentani, Caballero et al. 2003)。
    ( b) CM: cancer modules(browse 431 gene sets)
    由 Segal 等人定义的基因集合(Segal et al. 2004)。 作者从各种资源(如 KEGG, GO 等)编辑了基因集( '模块')。 通过挖掘癌症相关微阵列数据的大纲,他们确定了 456 个这样的模块,在各种癌症条件下显著变化。
    6) C5: GO gene sets(browse 5917 gene sets)
    该集合中的基因集是从 Gene Ontology 项目( www.geneontology.org )的中导出的。基因集是基于 GO 术语( go-basic.obo, 2016 年 5 月 3 日下载)及其与人类基因的相关性的( gene2go,从 2016 年 5 月 3 日从 NCBIFTP 服务器下载)。分为三部分:分子功能( MF),细胞组分(CC)和生物过程( BP)。基因产物可能与一个或多个分子功能相关或位于其中。
    GO 注释由描述特定 GO 项和基因产物之间的关联所基于的工作或分析的特定参考相关联的 GO 项组成。每个注释还必须包括证据代码,以指示如何支持对特定术语的注释( http://geneontology.org/page/guide-goevidence-codes )。GO 基因集省略了非常广泛的类别,如生物过程,同时也省略了具有少于 10 个基因的 GO 集(NCBI Entrez Gene ID)。如果 Jaccard 的系数大于 0.85,其将集合定义为“高度相似”。对于每一对高度相似的集合,其保持最大的集合并重复该过程,直到所有这样的对被解决。
    (a) BP: GO biological process(browse 4436 gene sets)
    生物过程(BP) http://www.geneontology.org/GO.process.guidelines.shtml
    (b) CC: GO cellular component(browse 580 gene sets)
    细胞组分(CC) http://www.geneontology.org/GO.component.guidelines.shtml
    (c) MF: GO molecular function(browse 901 gene sets)
    分子功能(MF) http://www.geneontology.org/GO.function.guidelines.shtml
    7) C6: oncogenic signatures(browse 189 gene sets)
    该基因集代表在癌症中经常失调的细胞途径。 大多数是直接从 NCBI GEO 的微阵列数据或者涉及已知癌基因的变化的内部未公开的分析实验产生的。 此外,少数致癌基因是从科学出版物找出的。
    8) C7: immunologic signatures(browse 4872 gene sets)
    免疫学标记集合(也称为 ImmuneSigDB)由代表细胞类型,状态和免疫系统内的变化的基因组组成。从人和人免疫学中已发表的研究中找寻的。其首先捕获在免疫学文献中公布的相关微阵列数据集的具有保存到基因表达综合征(GEO)的原始数据。 对
    于每个公开的研究,鉴定相关比较(例如 WT 对 KO;治疗前和治疗后等),并且创建简短的,生物学上有意义的描述。以相同的方式处理和归一化所有数据后用来鉴定基因组。
    7.2 Chip2Chip mapping
    GSEA 软件提供 Gene sets database 与 target chip 比对的功能,最后结果为两者的交集,最后产生的结果为两者交集的 gmt 文件,可用于 GSEA 分析。


    8. FAQ
    8.1 Error: Tool execution error
    运行开始后即跳出如下错误: Connection reset by peer: socket write error.


    A:此类错误一般由网络连接问题导致。 由于选择了在线的 GSEA 基因集文档和 chip 文档,因此软件在运行时需要连接相关在美国的 ftp 服务器,这可能会被一些学校、公司的网络所禁止。 为了避免在线分析, 我们只需要从官网页面 http://software.broadinstitute.org/gsea/downloads.jsp 下载我们所需要的 gmt、 chip 文件, 在 run 之前选择基因集等数据库时,选择 local 文件即可。
    8.2 GSEA 支持什么芯片平台和物种?
    A: 只要所需要分析的表达量数据有基因标识符(gene identifier),且能对应到 GSEA 的基因集中(gene sets), GSEA就能对此表达量数据做分析。GSEA 所利用的 MSigDB 中的基因都含有人类基因标识(human gene symbols), 其通过特殊的 chip file 将其他基因的标识转化为人类基因标识。 因此,需要 chip files 提供相关的 mapping 关系数据。所以,如果我们分析的数据是非人样本,我们需要确定是否利用 MSigDB 来进行分析。 具体情况如下:
    1) 利用非人物种作为模式物种来研究人类。 MSigDB 是正确的选择来进行分析,只需提供合适的 chip file 即可;
    2) 非人物种是研究的主要对象,没有需要去联合研究人类。MSigDB依然是可以选择的,只需要提供合适的chip file;
    3) 非人物种是研究的主要对象, 且不希望利用 MSigDB。 客户需要自己提供 gene sets 文件(GMT 或 GMX 格式),
    同样需要确认在 gene sets 中的基因标识能和测序结果的基因标识匹配得上, 如果不能匹配,则需要提供合适的 chip files 来进行匹配。
    PS: GSEA 中自带的 chip files 都只提供了匹配到人类基因标识的关系,不支持匹配到其他物种。
    8.3 RNA-Seq 数据和 Ensembl CHIP files
    A: GSEA 利用基因表达数据文件中的 gene identifiers 跟 MSigDB 的 HUGO human gene symbols 比对, 比对关系则依赖于 chip files 的信息。 但在 RNA-seq 中, 针对参考基因组来源于 Ensembl 数据库的分析结果, GSEA 提供了相关 chip
    file, 将人和鼠的 Ensembl IDs 转化为 HUGO human gene symbols。 目前支持的 Ensembl id 前缀是:人(ENSG 和 ENST),鼠(ENSMUSG 和 ENSMUST)。
    在设置参数的过程中,一般建议 gene sets 利用 Hallmarks collection 做分析,对于选择哪种 chip file, 规则如下:
    ENSEMBL_human_gene.chip => Ensembl ID prefix ENSG
    ENSEMBL_human_transcript.chip => Ensembl ID prefix ENST
    ENSEMBL_mouse_gene.chip => Ensembl ID prefix ENSMUSG
    ENSEMBL_mouse_transcript.chip => Ensembl ID prefix ENSMUST

    展开全文
  • IEEE-GRSL-GSEA-Code-源码

    2021-03-09 09:55:16
    代码共享
  • GSEA在全基因组表达谱芯片数据分析中的应用,比较系统的介绍了GSEA的用法,有利于初学者的学习
  • 如果之前看过生信宝典的一文掌握GSEA,超详细教程,一定会特别熟悉GSEA的原理和操作流程。当然越是理解,越是想不明白单个基因怎么做GSEA。当然如果您不熟悉GSEA,建议先看上一篇文章。 后来群友点拨理解了,不是对...

    生物信息学习的正确姿势

    NGS系列文章包括NGS基础、在线绘图、转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这)、ChIP-seq分析 (ChIP-seq基本分析流程)、单细胞测序分析 (重磅综述:三万字长文读懂单细胞RNA测序分析的最佳实践教程)、DNA甲基化分析、重测序分析、GEO数据挖掘(典型医学设计实验GEO数据分析 (step-by-step))、批次效应处理等内容。

    今天在讨论群看到有群友提问 单基因GSEA怎么做?。之前也看到过这个概念,但一直不清楚这个是什么含义,一直以为是用单个基因做GSEA。如果之前看过生信宝典的一文掌握GSEA,超详细教程,一定会特别熟悉GSEA的原理和操作流程。当然越是理解,越是想不明白个基因怎么做GSEA。当然如果您不熟悉GSEA,建议先看上一篇文章。

    后来群友点拨理解了,不是对单个基因做GSEA,是拿单个基因 (一般是感兴趣的基因)作为分组方式,探索与给定的单个基因相关的 (可以是表达相关,也可以是其它相关)基因富集在哪些调控通路和分子功能。

    分组方法有两种,一种是定性分组,一种是定量相关

    定性分组

    根据感兴趣的单个基因的度量值如表达值进行分组,比如按照该基因表达中位数分组,该基因表达值高于中位数的样品为一组,低于中位数的样品为一组,构建一个cls文件。如下,假如有8个样品,其中4个样品中该基因表达高 (samp1, samp3, samp7, samp8),4个样品中该基因表达低(samp 2 4 5 6),则cls文件如下 (一定注意样本顺序要一致):

    8 2 1
    #gene_high gene_low
    gene_high gene_high gene_high gene_high gene_low gene_low gene_low gene_low

    调整后的表达矩阵格式如下 (注意列的对应,high对高的样品。)

    Gene samp1 samp3 samp7 samp8 samp2 samp4 samp5 samp6
    A 4 4 4 4 1 1 1 1
    B . . . . . . . .
    C . . . . . . . .

    后续的操作就不说了,还是看生信宝典的一文掌握GSEA,超详细教程,看完就都会了。

    :也可以按照该基因表达的第一和三四分位数分组,小于第一四分位数的为一组,大于第三四分位数的为另一组。

    相关性排序

    与前面把样本分组不同,这里样本不进行分组了,而是把感兴趣基因的表达做为样本的一个属性。在做GSEA分析时,其它基因按照与感兴趣基因的表达相关性排序进行后续分析。

    这时应该怎么准备cls文件呢?

    咱们先以一个时间序列样本的cls文件为例:

    • #numeric为固定写法,第一行,不需要修改

    • #Time名字随便取,这里是时间序列,取名Time#是必须的。

    • 3行是每个样品的处理时间,00小时11小时;每个时间3个重复,所以写了3遍;总共5个时间点,15个样品。

    #numeric
    #Time
    0 0 0 1 1 1 6 6 6 24 24 24 48 48 48

    回到我们这个例子,还是8个样品,分别为samp 1 2 3 4 5 6 7 8,假如感兴趣基因是A,表达矩阵如下:

    Gene samp1 samp2 samp3 samp4 samp5 samp6 samp7 samp8
    A 9 8 7 6 3 4 1 2
    B . . . . . . . .
    C . . . . . . . .

    这时对应的cls文件这么写(注意一一对应关系)。Aexpr随便起的一个名字,代表A基因的表达。

    #numeric
    #Aexpr
    9 8 7 6 3 4 1 2

    然后导入GSEA就可以分析了。需要注意的是选择合适的Ranking metric,如pearson相关性CosineManhattanEuclidean

    基于相关性的GSEA操作展示

    直接看动画,数据格式也有展示,GMT文件是自己整理的。这是1我们单细胞和群体转录组课程的一个小环节 (回头把这部分视频拆出来放到腾讯课堂供访问)。

    公众号看不了动画,截图两张,点击阅读原文去查看吧。

    718d89dddb45c0542ac9b3de4929e5a5.png1baa5f69b3fd2c13aafcf1fdef65814e.png

    8fe77c0949eae828f0a1cbd92b747f15.png

    讨论学习是个很好的方式,欢迎大家有问题发到train@ehbio.com,信息的,问题可重现的,或有意思的开放问题我们都会给予解决,写个推文发出,既方便自己,又方便他人。

    当然如果类似转录组怎么分析, 宏基因组怎么分析,这样大的问题还是参加我们的线下培训班或购买网课吧,都在www.ehbio.com/Training

    • 这个只需一步就可做富集分析的网站还未发表就被CNS等引用超过350次

    • 什么,你算出的P-value看上去像齐天大圣变的庙?

    • DESeq2差异基因分析和批次效应移除

    • GO、GSEA富集分析一网打进

    • GSEA富集分析 - 界面操作

    • 无需写代码的高颜值富集分析神器

    • 去东方,最好用的在线GO富集分析工具

    • 没钱买KEGG怎么办?REACTOME开源通路更强大

    • 超简便的国产lncRNA预测工具LGC

    • 一文掌握GSEA,超详细教程

    • UCSC XENA - 集大成者(TCGA, ICGC)

    • ICGC数据库使用

    • TCGA数据库在线使用

    • BROAD开发的TCGA分析平台,强大的下载功能

    • cBioPortal功能强大的TCGA再分析平台

    • 这是数据更新最实时的TCGA网站,功能强大

    • 不懂R,如何进行GEO数据库表达谱的差异分析、富集分析、蛋白互作、可视化?

    • 典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集

    • 典型医学设计实验GEO数据分析 (step-by-step) - 数据获取到标准化

    • 为什么GEO2R/GEOquery的结果可能是错的?

    • WGCNA分析,简单全面的最新教程

    • psych +igraph:共表达网络构建

    • 一文学会网络分析——Co-occurrence网络图在R中的实现

    • 一文看懂PCA主成分分析

    往期精品(点击图片直达文字对应教程)

    740a6ef2b107b0f43cfd10fe8fd39b6b.png

    c7772ce5350c798363a755ced8ac227d.png

    f0ca9daba7bd27693f8839db8f754f69.png

    79c0b8c37d0e324d17a4ecc8d151a73f.png

    215b937a458b543c41276e485fa512ec.png

    c882954384aaacb1921225ae6191ca13.png

    8f72e0aa00c577b79d0a08ce5fd9920c.png

    e5404ab47b38572d03aee0babb968839.png

    9b91220a8a5c71782f536755637e6193.png

    f8d7b87d1ea862a5075b3a2bb998082e.png

    34e670d816b14088cc74f3eb9d18e230.png

    763b025f957723d7181bffb00371669a.png

    88d24507cd97fc7dfd0711416310a8c6.png

    def3b7ab5f33ac7cf6cb05bd3718a06c.png

    575318a74aaabef7f19521d871b89a12.png

    3c6ea772a6f648e5a46ee83836cf8004.png

    0d8c25a2a7be83fc242ea2c2c2f2bd33.png

    3eefe832ddf00311c4c7af3c2c46cd7f.png

    14d250cbfbad403d52c42233507dc447.png

    9f0fea14c6cd5818e093fada7b4d29b8.png

    c6b24c01386237b163e2d0d84c2e9c80.png

    da3a2925feb165d5222aed8dbae27e48.png

    c161a08fbccbc7e5e4b311136f3e421d.png

    0a951ad287e4110532c371e0cb692a85.png

    d9aac57dad7929e6ed11b84141fd70b9.png

    85dd7ec4b2a72ef3f225e9ffabc47f0a.png

    aefc3385f720f1e366b23791304e7a54.png

    f1d5e868840cb953c46cb90a1d5a3dc4.png

    后台回复“生信宝典福利第一波”或点击

    a017ada33038e30c66f9ffe9d47d9c3e.png

    477a679f8dab7e49d18b2ac1337d9d02.png

    979e33978c421d34f7d04b77ac83a5e8.png

    展开全文
  • 一文掌握GSEA富集分析-最详细教程

    万次阅读 多人点赞 2019-04-26 17:02:46
    之前总结了一篇关于GSEA富集分析的推文——《GSEA富集分析 - 界面操作》,大略介绍了GSEA的定义、GSEA原理、GSEA分析、Leading-edge分析等,不太了解的朋友可以点击阅读先理解下概念。 最近用自己数据实战分析时...

    生信宝典之前总结了一篇关于GSEA富集分析的推文——《GSEA富集分析 - 界面操作》,介绍了GSEA的定义GSEA原理GSEA分析Leading-edge分析等,不太了解的朋友可以点击阅读先理解下概念 (下面摘录一部分)。

    GSEA案例解析

    介绍GSEA分析之前,我们先看一篇Cell文章(https://sci-hub.tw/10.1016/j.cell.2016.11.033)的一个插图。

    img

    以下是文章原文对图的注解:GSEA analyses of genesets for cardiac (top) and endothelial/endocardial (bottom) development. NES, normalized enrichment score. FDR, false discovery rate. Positive and negative NES indicate higher and lower expression in iwt, respectively.

    关于文章中使用的GSEA分析方法和参数,我们截取对应原文:Gene Set Enrichment Analysis was performed using the GSEA software (https://www.broadinstitute.org/gsea/) with permutation = geneset, metric = Diff_of_classes, metric = weighted, #permutation = 2500.

    根据以上信息可知,上图是研究者使用GSEA软件所做的分析结果。文章通过GSEA分析,发现

    与心脏发育有关的基因集 (影响心脏的收缩力、钙离子调控和新陈代谢活力等)在iwt组 (GATA基因野生型)中普遍表达更高,而在G296S组 (GATA基因的一种突变体)中表达更低;

    而对于参与内皮或内膜发育的基因集,在iwt组中表达更低,在G296S组中表达更高。

    作者根据这个图和其它证据推测iwt组的心脏发育更加完善,而G296S组更倾向于心脏内皮或内膜的发育,即GATA基因的这种突变可能导致心脏内皮或内膜的过度发育而导致心脏相关疾病的产生。

    那么GSEA分析是什么?

    参考GSEA官网主页的描述:Gene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states (e.g. phenotypes). 在上述Cell文章中,作者更加关心参与心脏发育的基因集 (即a priori defined set of genes)与两个状态(突变体和野生型,状态的度量方式是基因表达)的关系,因此利用GSEA对其进行分析后发现,参与心脏发育 (收缩力、钙调控和新陈代谢)的基因集的表达模式更接近于iwt组的表型,而不是G296S组; 而参与心脏内皮或内膜发育的这些基因的表达模式更接近于G296S组的表型而不是iwt组的表型。

    这就是GSEA分析所适用的主要场景之一。它能帮助生物学家在两种不同的生物学状态 (biological states)中,判断某一组有特定意义的基因集合的表达模式更接近于其中哪一种。因此GSEA是一种非常常见且实用的分析方法,可以将数个基因组成的基因集与整个转录组、修饰组等做出简单而清晰的关联分析。

    除了对特定gene set的分析,反过来GSEA也可以用于发现两组样本从表达或其它度量水平分别与哪些特定生物学意义的基因集有显著关联,或者发现哪些基因集的表达模式或其他模式更接近于表型A、哪些更接近于表型B。这些特定的基因集合可以从GO、KEGG、Reactome、hallmark或MSigDB等基因集中获取,其中MSigDB数据库整合了上述所有基因集。研究者也可自定义gene set (即新发现的基因集或其它感兴趣的基因的集合)。

    GSEA分析似乎与GO分析类似但又有所不同。GO分析更加依赖差异基因,实则是对一部分基因的分析 (忽略差异不显著的基因),而GSEA是从全体基因的表达矩阵中找出具有协同差异 (concordant differences)的基因集,故能兼顾差异较小的基因。因此二者的应用场景略有区别。另外GO富集是定性的分析,GSEA考虑到了表达或其它度量水平的值的影响。所以,对于时间序列数据或样品有定量属性时,GSEA的优势会更明显,不需要每个分组分别进行富集,直接对整体进行处理。可以类比于之前的WGCNA分析

    GSEA定义

    Gene Set Enrichment Analysis (基因集富集分析)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势,从而判断其对表型的贡献。其输入数据包含两部分,一是已知功能的基因集 (可以是GO注释、MsigDB的注释或其它符合格式的基因集定义),一是表达矩阵 (也可以是排序好的列表),软件会对基因根据其与表型的关联度(可以理解为表达值的变化)从大到小排序,然后判断基因集内每条注释下的基因是否富集于表型相关度排序后基因表的上部或下部,从而判断此基因集内基因的协同变化对表型变化的影响。

    (The gene sets are defined based on prior biological knowledge, e.g., published information about biochemical pathways or coexpression in previous experiments. The goal of GSEA is to determine whether members of a gene set S tend to occur toward the top (or bottom) of the list L, in which case the gene set is correlated with the phenotypic class distinction.)

    这与之前讲述的GO富集分析不同。GO富集分析是先筛选差异基因,再判断差异基因在哪些注释的通路存在富集;这涉及到阈值的设定,存在一定主观性并且只能用于表达变化较大的基因,即我们定义的显著差异基因。而GSEA则不局限于差异基因,从基因集的富集角度出发,理论上更容易囊括细微但协调性的变化对生物通路的影响,尤其是差异倍数不太大的基因集。

    GSEA原理

    给定一个排序的基因表L和一个预先定义的基因集S (比如编码某个代谢通路的产物的基因, 基因组上物理位置相近的基因,或同一GO注释下的基因),GSEA的目的是判断S里面的成员sL里面是随机分布还是主要聚集在L的顶部或底部。这些基因排序的依据是其在不同表型状态下的表达差异,若研究的基因集S的成员显著聚集在L的顶部或底部,则说明此基因集成员对表型的差异有贡献,也是我们关注的基因集。

    GSEA计算中几个关键概念:

    1. 计算富集得分 (ES, enrichment score). ES反应基因集成员s在排序列表L的两端富集的程度。计算方式是,从基因集L的第一个基因开始,计算一个累计统计值。当遇到一个落在s里面的基因,则增加统计值。遇到一个不在s里面的基因,则降低统计值。

      每一步统计值增加或减少的幅度与基因的表达变化程度(更严格的是与基因和表型的关联度,可能是fold-change,也可能是pearson corelation值,后面有介绍几种不同的计算方式)是相关的,可以是线性相关,也可以是指数相关 (具体见后面参数选择)。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集,负值ES表示基因集在列表的底部富集。

    2. 评估富集得分(ES)的显著性。通过基于表型而不改变基因之间关系的排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少,也可基于基因集做排列检验 (permutation test),计算p-value

    3. 多重假设检验校正。首先对每个基因子集s计算得到的ES根据基因集的大小进行标准化得到Normalized Enrichment Score (NES)。随后针对NES计算假阳性率。(计算NES也有另外一种方法,是计算出的ES除以排列检验得到的所有ES的平均值)

    4. Leading-edge subset,对富集得分贡献最大的基因成员。

    本文通过总结多人学习使用过程中遇到的问题进一步记录软件操作过程和结果解读,力求讲清每个需要注意的细节点。

    从前文中我们了解到GSEA分析的目的是要判断S集基因(基于先验知识的基因注释信息,某个关注的基因集合)中的基因是随机分布还是聚集在排序好的L基因集的顶部或底部(这便是富集分析)。

    与GO富集分析的差异在于GSEA分析不需要指定阈值(p值或FDR)来筛选差异基因,我们可以在没有经验存在的情况下分析我们感兴趣的基因集,而这个基因集不一定是显著差异表达的基因。GSEA分析可以将那些GO/KEGG富集分信息中容易遗漏掉的差异表达不显著却有着重要生物学意义的基因包含在内。

    下面来看看软件具体操作和结果解读。

    一、软件安装

    软件下载地址:http://software.broadinstitute.org/gsea/downloads.jsp

    使用官方推荐的第一个软件javaGSEA Desktop Application,根据分析数据的大小和电脑内存多少可以选择下载不同内存版本的软件。该软件是基于java环境运行的,而且需要联网。若会出现打不开的现象(小编就是就碰到了),要么是没有安装java,要么是java版本太低了,安装或更新下java就能打开。也可能是网速太慢,或Java安全性问题,这时选择官网提供的第二个软件javaGSEA Java Jar file,同样依赖java运行,但不需联网,启动快。

    软件启动界面如下:

    img

    二、数据准备

    所有矩阵的列以tab键分割,不同类型的数据格式和后缀要求见下表。

    Data FileContentFormatSource
    Expression datasetContains features (genes or probes), samples, and an expression value for each feature in each sample. Expression data can come from any source (Affymetrix, Stanford cDNA, and so on).res, gct, pcl, or txtYou create the file. 一般的基因表达矩阵整理下格式就可以。如果是其它类型数据或自己计算rank也可以,后面有更多示例。(如果后缀为txt格式,传统的基因表达矩阵就可以,第一列为基因名字,名字与待分析的功能注释数据集一致,同为GeneSymbol或EntrezID或其它自定义名字,第一行为标题行,含样品信息。gct文件需要符合下面的格式要求。)
    Phenotype labelsContains phenotype labels and associates each sample with a phenotype.clsYou create the file or have GSEA create it for you. 一般是样品分组信息或样品属性度量值或时间序列信息。
    Gene setsContains one or more gene sets. For each gene set, gives the gene set name and list of features (genes or probes) in that gene set.gmx or gmtYou use the files on the Broad ftp site, export gene sets from the Molecular Signature Database (MSigDb) or create your own gene sets file. 欲检测是否富集的基因集列表。注意基因ID与表达矩阵基因ID一致。自己准备的基因集注意格式与官网提供的gmt格式一致。
    Chip annotationsLists each probe on a DNA chip and its matching HUGO gene symbol. Optional for the gene set enrichment analysis.ChipYou use the files on the Broad ftp site, download the files from the GSEA web site, or create your own chip file. 主要是为芯片探针设计的转换文件。如果表达矩阵的基因名与注释集基因名一致,不需要这个文件。

    1. 表达数据集文件

    GESA提供有Example Datasets,下载地址:http://software.broadinstitute.org/gsea/datasets.jsp

    在这里可以下载表达矩阵Expression dataset(gct文件,常见txt格式也可以)和样品分组信息Phenotype labelscls文件)

    img

    数据示例中两个gct文件都是表达矩阵,其中*hgu133a.gct文件第一列是探针名字,*collapsed.gct文件的第一列是gene symbol。

    • 第一行:#1.2,表示版本号,自己准备文件时照抄就行;
    • 第二行:两个数分别表示gene NAME的数量和样本数量(矩阵列数-2);
    • 矩阵:第一列是NAME;第二列Description,没有的话可以全用na或任意字符串填充;后面的就是基因在不同样本中标准化后的表达数据了 (部分统计量metrics for ranking genes计算需要log转换后的数据,后面会有提及。其它情况是否为log转换的数据都可用,GSEA关注的是差异,只要可比即可)。

    img

    2. 样品分组信息

    • 第一行:三个数分别表示:34个样品,2个分组,最后一个数字1是固定的;
    • 第二行:以#开始,tab键分割,分组信息(有几个分组便写几个,多个分组在比较分析时,后面需要选择待比较的任意2组);(样品分组中NGT表示正常耐糖者,DMT表示糖尿病患者,自己使用时替换为自己的分组名字)
    • 第三行:样本对应的组名。样本分组信息的第三行,同一组内的不同重复一定要命名为相同的名字,可以是分组的名字。例如相同处理的不同重复在自己试验记录里一般是Treat6h_1、Treat6h_2、Treat6h_3,但是在这里一定都要写成一样的值Treat6h。与表达矩阵的样品列按位置一一对应,名字相同的代表样品属于同一组。如果是样本分组信息,上图中的01也可以对应的写成NGTDMT,更直观。但是,如果想把分组信息作为连续表型值对待,这里就只能提供数字

    img

    3. 功能基因集文件(gene sets)

    GSEA官网提供了8种基因分类数据库,都是关于人类的数据,包括Marker基因,位置临近基因,矫正过的基因集,调控motif基因集,GO注释,癌基因,免疫基因,最新一次更新是在2018年7月,下载地址:http://software.broadinstitute.org/gsea/downloads.jsp#msigdb

    img

    官网提供的gmt文件有两种类型,*.symbols.gmt中基因以symbols号命名,*.entrez.gmt中基因以entrez id命名。注意根据表达矩阵的基因名字命名方式选择合适的基因集。表达数据和通路数据能关联在一起依赖的是基因名字相同,所以一定保证基因命名方式的统一。

    img
    gmt格式是多列注释文件,第一列是基因所属基因集的名字,可以是通路名字,也可以是自己定义的任何名字。第二列,官方提供的格式是URL,可以是任意字符串。后面是基因集内基因的名字,有几个写几列。列与列之间都是TAB分割。

    Pathway_description	Anystring	Gene1	Gene2	Gene3
    Pathway_description2	Anystring	Gene4	Gene2	Gene3	Gene5
    

    GSEA官网只提供了人类的数据,但是掌握了官网中基因表达矩阵和注释文件的数据格式,就可以根据自己研究的物种,在公共数据库下载对应物种的注释数据,自己制作格式一致的功能基因集文件,这样便就可以做各种物种的GSEA富集分析了。

    4. 芯片注释文件

    如果分析的表达数据是芯片探针数据就需要用到芯片注释文件(chip),用来做ID转换,把探针名字转换为基因名字。如果我们的表达数据文件中已经是基因名了就不再需要这个文件了。

    三、分析参数设置和软件运行

    演示使用的数据来自GSEA官网:

    • 表达矩阵:Diabetes_collapsed_symbols.gct
    • 样品分组信息:Diabetes.cls
    • 基因功能分类数据选择GO数据库:c5.all.v6.2.symbols.gmt
    • 因为表达矩阵没有选择芯片数据,则第四个文件不用选

    1. 数据导入

    img
    按照上图步骤依次点击Load data——Browse for file——在弹出文件框中找到待导入的文件,选中点击打开即可;

    若文件格式没问题会弹出一个提示There were no error的框,证明文件上传成功,并且会显示在5所示的位置;若出错,请仔细核对文件格式。

    注意:
    1)本地文件存放路径不要有中文、空格(用_代替空格)和其他特殊字符;
    2)所有用到的文件都需要通过上述方式先上传至软件;
    3)数据上传错误后可以通过点击工具栏file——clear recent file history进行清除。

    2. 指定参数

    点击软件左侧Run GSEA,将跳出参数选择栏。参数设置分为三个部分Require fields(必须设置的参数项)、Basic fields(基本参数设置栏)和Advanced fields(高级参数设置栏),后面两栏的参数一般不做修改,使用默认的就行。后面两部分参数设置,如果涉及到需要根据实验数据做调整的地方,会在后面的分析中会提到。

    1)Require fields

    img

    • Expression dataset: 导入表达数据集文件,点击后自动显示上一步中从本地导入软件内的文件,所以一定要确认上一步导入数据是否成功;
    • Gene sets database: 基因功能集数据库,可以从本地导入(上一步);在联网的情况下软件也可以为自动下载GSEA官网中的gene sets文件;
    • Number of permutations: 置换检验的次数,数字越大结果越准确,但是太大会占用太多内存,软件默认检验1000次。软件分析时会得到一个基因富集的评分(ES),但是富集评分是否具有统计学意义,软件就会采用随机模拟的方法,根据指定参数随机打乱1000次,得到1000个富集评分,然后判断得到的ES是否在这1000个随机产生的得分中有统计学意义。测试使用时建议填一个很小的数如10,先让程序跑通。真正分析时再换为1000。
    • Phenotype labels: 选择比较方式,如果文件只有2个组别的话就比较方便了,任意选一个就行,哪个在前在后全在自己怎么解释方便;如果数据有多组的话,GSEA会提供两两间比较的组合选项或者某一组与剩下所有组的比较。选择好后,GSEA会在分析过程中根据组别信息自动到表达数据集文件中提取对应的数据作比较。
    • Collapse dataset to gene symbols: 如果表达数据集文件中NAME已经与gene sets database中名字一致,选择FALSE,反之选择TRUE
    • Permutation type: 选择置换类型,phenotype或者gene sets每组样本数目大于7个时 ,建议选择phenotype,否则选择gene sets。
    • Chip platform: 表达数据集为芯片数据时才需要,目的是对ID进行注释转换,如果已经转换好了就不需要了。应该也适用于其它需要转换ID的情况,不过事先转换最方便。

    2)Basic fields

    img

    通常选择默认参数即可,在此简单介绍一下

    • Analysis name: 取名需要注意不能有空格,需要用_代替空格。如果做的分析多,最好选择一个有意义的名字,比如shengxinbaodian (生信宝典全拼),方便查找。

    • Enrichment statistic: 基因集富集分析(PNAS)的最后一部分给出了GSEA中所用方法的数学描述,感兴趣的可以查看一下论文。在此给出每种富集分析不同算法的参数情况:
      ▪ classic: p=0 若基因存在,则ES值加1;若基因不存在,则ES值减1
      ▪ weighted (default): p=1 若基因存在,则ES加rank值;若基因不存在,则ES减rank值
      ▪ weighted_p2: p=2 基因存在,ES加rank值的平方,不存在则减rank值的平方
      ▪ weighted_p1.5: p=1.5 基因存在,ES加rank值得1.5次方,不存在则减rank值得1.5次方

    备注:如果想用其它加权,就自己计算rank值,使用preranked mode

    • Metric for ranking genes: 基因排序的度量

      • 下面提到的均值也可以是中位数
      • 如果表型是分组信息,GSEA在计算分组间的差异值时支持5种统计方式,分别是signal2noiset-Testratio_of_classdiff_of_class(log2转换后的值计算倍数)和log2_ratio_of_class。下面公式很清楚。
      • 如果表型是连续数值信息(定量表型): GSEA通过表型文件(cls)和表达数据集文件(gct),使用pearson相关性CosineManhattanEuclidean指标之一计算两个配置文件之间的相关性。(注意:若是分组表型文件想转换为定量表型,cls文件中分类标签应该指定为数字)
    • Gene list sorting mode: 对表达数据集中的基因进行排序,按照排序度量的真实值(默认)或者绝对值排序;

    • Gene list ordering mode: 使用此参数确定表达数据集中基因是按照降序(默认)或者升序排列;

    • Max size & Min size: 从功能基因集中筛选出不属于表达数据集中的基因后,剩下基因总数在此范围内则保留下来做后续的分析,否则将此基因集排除;一般太多或太少都没有分析意义。

    • Save results in this folder: 在此可以选择分析文件在本地电脑的存储地址。

    3)Advanced fields

    img

    • Collapsing mode for probe sets => 1 gene: 多个探针对应一个基因时的处理方式。
    • Normalized mode: 富集得分的标准化方式。
    • Randomization mode:只用于phenotype permutation。
    • Median for class metrics: 计算metrics ranking时用中值而不是平均值。
    • Number of markers:红蝶图中展示的Gene Marker数目。
    • Plot graphs for the top sets of each phenotype:绘制多少GSEA plot,默认top 20,其它不绘制。一般会把这个值调高。
    • Seed for permutation:随机数种子,如果想让每次结果一致,这里需要设置同样的一个整数。

    以上参数都设置好后点击参数设置栏下方的一个绿色按钮Run,若软件左下方GSEA reports处的状态显示Running的话则表示运行成功,此过程大概需要十分钟左右,视数据大小而定。

    img

    • Command: 显示运行这个分析的命令行,以后就可以批量运行类似分析了。

    四、结果解读

    数据分析完后的结果会保存到我们设置的路径下,点开文件夹中的index.html就可以查看网页版结果,更加方便。

    结果报告分为多个子项目,其中最重要的是前面两部分,基因富集结果就在这里。从第三部分开始其实是软件在分析数据的过程产生的中间文件, 也很重要,读懂后可以加深对GSEA分析的认识,理解我们是如何从最初的基因表达矩阵得到最终的结果(即报告的前两个项目)。建议先从Dataset details看起,然后再返回看第一部分的结果报告。

    1. Enrichment in phenotype

    以正常人组NGT的17个样本数据为例解析最终结果。

    img

    报告首页文字总结信息表示:

    • 经过条件筛选后还剩下3953个GO条目,其中1697个GO条目在NGT组中富集;
    • 有36个GO基因条目在FDR<25%的条件下显著富集,这部分基因最有可能用于推进后续实验;
    • 在统计检验p<0.01, p<0.05的条件下分别有19和114个GO条目显著富集;
    • 结果有多种显示方式:图片快照(snapshot)、网页(html)和表格(Excel)形式;
    • 点击Guide to可以查看官方帮助解读结果的文档。

    1) 点击enrichment results in html,在网页查看富集结果,如下:
    img

    • GS:基因集的名字,GO条目的名字
    • SIZE:GO条目中包含表达数据集文中的基因数目(经过条件筛选后的值);
    • ES:富集评分;
    • NES:校正后的归一化的ES值。由于不同用户输入的基因数据库文件中的基因集数目可能不同,富集评分的标准化考虑了基因集个数和大小。其绝对值大于1为一条富集标准。计算公式如下:
      img
    • NOM p-val:即p-value,是对富集得分ES的统计学分析,用来表征富集结果的可信度;
    • FDR q-val:即q-value,是多重假设检验校正之后的p-value,即对NES可能存在的假阳性结果的概率估计,因此FDR越小说明富集越显著;
    • RANK AT MAX:当ES值最大时,对应基因所在排序好的基因列表中所处的位置;
      (注:GSEA采用p-value<5%,q-value<25%进行数据过滤)
    • LEADING EDGE:该处有3个统计值,tags=59%表示核心基因占该基因集中基因总数的百分比;list=21%表示核心基因占所有基因的百分比;signal=74%,将前两项统计数据结合在一起计算出的富集信号强度,计算公式如下:
      img
      其中n是列表中的基因数目,nh是基因集中的基因数目

    点击Details跳转至对应的详情结果。只有前20个GO富集详情可以查看,想要生成的结果报告可以查看更多的富集信息,可以通过在Advanced fields处设置参数Plot graphs for the top sets of each phenotype

    2)Details for gene set
    img
    首先是一个选定GOset下的汇总信息表,每一部分意思在上面已做解释,其中Upregulated in class表示该基因集在哪个组别中高表达,这个主要看富集分析后的leading edge分布位置。

    img
    接下来是富集分析的图示,该图示分为三部分,在图中已做标记:

    • 第一部分是Enrichment score折线图:显示了当分析沿着排名列表按排序计算时,ES值在计算到每个位置时的展示。最高峰处的得分 (垂直距离0.0最远)便是基因集的ES值。
    • 第二部分,用线条标记了基因集合中成员出现在基因排序列表中的位置,黑线代表排序基因表中的基因存在于当前分析的功能注释基因集。leading edge subset 就是(0,0)到绿色曲线峰值ES出现对应的这部分基因。
    • 第三部分是排序后所有基因rank值得分布,热图红色部分对应的基因在NGT中高表达,蓝色部分对应的基因在DMT中高表达,每个基因对应的信噪比(Signal2noise,前面选择的排序值计算方式)以灰色面积图显展示。

    在上图中,我们一般关注ES值,峰出现在排序基因集的前端还是后端(ES值大于0在前端,小于0在后端)以及Leading edge subset(即对富集贡献最大的部分,领头亚集);在ES图中出现领头亚集的形状,表明这个功能基因集在某处理条件下具有更显著的生物学意义;对于分析结果中,我们一般认为|NES|>1,NOM p-val<0.05,FDR q-val<0.25的通路是显著富集的。
    img
    最后还有一个该GO基因集下每个基因的详细统计信息表,RANK IN GENE LIST表示在排序好的基因集中所处的位置;RANK METRIC SCORE是基因排序评分,我们这里是Signal2noiseRUNNING ES是分析过程中动态的ES值;CORE ENRICHMENT是对ES值有主要贡献的基因,即Leading edge subset,在表中以绿色标记。

    2. Dataset details

    芯片原始数据和去重后的数据;如果分析的时候没有用到芯片数据或没涉及到名字转换则前后基因数目一样。
    img

    3. Gene set details

    我们分析提供的gmt文件中有多个GO条目,每个GO条目里又有多个基因;GSEA分析软件会在每个GO条目中搜索表达数据集gct文件中的基因,并判断有多少个在GO条目中;若经过筛选后保留在GO条目中的基因在15-500(闭区间)时该GO条目才被保留下来进行后续的分析。
    img

    此结果显示我们从5917个GO条目中淘汰了了1964个GO,剩下3953个GO条目用作后续分析。

    点击gene sets used and their sizes可以下载详细Excel表。

    Excel第一列是GO名称,第二列是GO条目中包含的基因数目,第三列是筛选后每个GO中还有多少基因属于表达数据集文件中的基因,不满足参数(15-500)的条目被抛弃,显示为Rejected不纳入后续分析。

    备注: 此处的筛选范围15-500是可调参数,在软件的参数basic fields处的Max sizeMin size处更改。
    img

    4. Gene markers for the NGT versus DMT comparison

    这部分展示的是我们提供的表达数据集文件中的基因在两个组别中的表达情况。

    输入的文件中总共有15056个基因,其中有7993个基因在正常人(NGT)中表达更高,占总基因数的53.1%;有7063个基因在糖尿病患者(DMT)中表达更高,占总基因数的46.9%。后面一个面积百分比,稍后看图的时候再做解释。

    img

    点击rank ordered gene list可以下载一个排序好的基因集Excel表,排序原则是根据Basic fields参数设置处的Metric for ranking genes决定的。我们选的是信噪比(signal2noise),显示在表格中的最后一列。根据NGT_vs_DMT评分得到一个降序排列的基因集,之后便可以做基因的富集分析了。

    GSEA基因富集分析的原理就是基于该排列好的基因集,从第一个基因开始判断该基因是否存在于经过筛选的GO功能基因集中,如果存在则加分,反之减分。所以评分过程是一个动态的过程,最终我们会得到一个评分峰值,那就是GO功能富集的评分。加分规则通过Basic fields参数设置处的Enrichment statistic决定的。

    img

    接着有一个分析的结果的热图和gene list相关性的图。

    热图中展示了分别在两组处理中高表达的前50个基因,总共100个基因的表达情况。
    img

    gene list相关性图如下。横坐标是已经排序好的基因,纵坐标是signal2noise的值。虚线左侧的基因是在NGT中高表达,右侧的基因在DMT中高表达。这部分结果报告中的面积比就是基于该图计算的,可以看出面积百分比和基因数目百分比有一定的差异,面积百分比可以从整体上反映组间信噪比的大小。
    img

    Butterfly plot显示了基因等级与排名指标评分之间的正相关(左侧)和负相关性(右侧)。左侧蓝色虚线和右侧红色虚线是真实的信噪比结果,其他颜色的线是软件对数据做了随机重排后的结果。默认情况下,图形只显示前100个基因,也就是排名第一和最后100个基因。可以使用运行GSEA页面上Advanced fields处的Number of markers来更改显示的基因数量。

    img

    5. Global statistics and plots

    这部分包含两个图:1) p值与归一化富集分数(NES)的对比图,这提供了一种快速、直观的方法来掌握有意义的富集基因集的数量。 2) 通过基因集的富集分数统计图,提供了一种快速、直观的方法来掌握富集的基因集的数量。

    img
    img

    理解了上面各个部分的结果后,再回过头看这张GSEA分析原理图就简单了。
    img

    Cytoscape富集网路可视化

    在GSEA软件的左侧提供了Enrichment Map Visualization的功能,点击后GSEA软件会自动调用Cytoscape,建议等待Cytoscape启动后再进行接下来的操作,且保证在分析过程中Cytoscape是处于开启状态

    选择一个GSEA分析结果,点击Load GSEA Results,其他项为默认值就行,点击Build Enrichment Map以展示基因富集结果的网络图。
    (备注:GSEA分析结果用的是和上面演示数据不同的文件,可自行更改)

    img
    运行成功之后会弹出下面的提示框,结果直接展现在了Cytoscape中,如下图所示:
    img
    img

    Graphpad作图比较多个ES

    GSEA富集分析可视化结果是给每个功能基因集富集情况单独出一张图,有的时候我们想要比较基因集在两个不同的GO中的富集情况,利用GSEA软件分析得到的Excel结果表,提取有用的数据结果,在graphpad里进行加工再出图,可以达到我们想要的结果!

    效果图如下:
    img

    《Graphpad,经典绘图工具初学初探》一文中介绍了graphpad入门的基础知识,基本操作可以单击回看。最近使用graphpad发现其多图排版功能十分强大,不仅可以实现多个图形排版还能实现图层叠加。上面这个图的作图思路也就是把该图拆分为两部分,Enrichment score和基因位置分布条带图。

    在GSEA分析结果文件夹里随便找一个感兴趣的GO条目分析结果Excel表,作图需要提取的信息即图中标黄的部分,RANK IN GENE LISTRUNNING ES
    img

    加工一下已有数据,添加一列high取值都为0.1,设置高度,黄色部分的数据就是用来绘制基因位置分布条带图的;绿色部分用来绘制动态的ES评分曲线。
    img

    打开graphpad之后,我们在XY类图下选择Enter and plot a single Y value for each point,将两部分数据分开粘贴到软件不同数据表格中(如下图左侧所示),下图中间展示两个图选择的不同绘图方式,调整参数后最终得到右侧的结果。
    img
    在左侧目录树处点击layout创建一个图形排版界面,将Graphs下的图形复制粘贴到layout1下,拖拉移动位置很快就能将两部分图对齐。

    img
    之后用同样地方式画另外一个富集结果,粘贴到layout1中便得到最开始展示的图。

    注意:设置X轴的范围是1到总排序基因数,Y轴是0到多个富集分析得分的最大值。

    GSEA分析自定义功能注释集

    软件安装和数据格式准备不清晰的请见前文。

    下面以拟南芥的转录组数据为例介绍GSEA的使用。待分析的数据是拟南芥的两组实验 (mut_vs_wt),每组各三个生物重复,故一共6个样本。各个样本有唯一的名称:mut_1、mut_2、mut_3、wt_1、wt_2、wt_3。

    准备数据

    在分析之前需要准备三个文件:样本分组文件 (Phenotype labels, cls格式)、基因表达文件 (Expression dataset)和基因集数据库文件(gene sets database)

    Phenotype labels file

    本质上是一个文本文件,需按照以下格式书写:

    img

    其中每两个字符串之间均用空格隔开

    第一行6 表示共有6个样本;2 表示是两个分组;1 是固定写法,不要改。

    第二行的#后面没有空格,mutwt为2个分组的名字,可以自行更改为自己试验设定的分组名字,名字内部不只允许有字母、数字、下划线

    第三行的书写方式需要格外注意,这里写的是6个字符串(代表6个样本),样本重复必须赋予相同的名字,以便GSEA判断样品所属的组,并且其排列的相对顺序必须与稍后要准备的基因表达文件(或基因表达矩阵)中对应的样本书写顺序相同。

    • 如果你使用windows来书写这个文件,建议使用NotePadd++UltraEdit等。不建议使用Excel或记事本(空格、换行及Tab键等格式不易把握)。写完之后注意需要将文件后缀名保存为.cls

    • 如果你使用bash环境的vi或vim来写这个文件(使用windows的请忽略此处),那么保存为.cls文件后使用cat -A查看这个文件,输出到屏幕上的应是 ($代表换行符,不需要自己书写):

    6 2 1$
    #mut wt$
    mut mut mut wt wt wt$
    

    基因表达文件(Expression dataset)

    img

    • 该文件建议使用Excel文件生成(测序公司一般都会提供这个文件)。只是需要按照下面几项稍作修改:

      • 第一行是固定写法:#1.2

      • 第二行25000表示此表格有25000个基因,注意25000不是这个表的总行数,而是从第四行开始数的基因数,并且无重复的基因名称6表是有6个样本(与样本设置文件中的6必须对应)。

      • 第三行NAMEDESCRIPTION为固定写法。DESCRIPTION列下面可以写na或任意字符串。

      • 第四行及以后的行中,基因的表达值最好是没有取过对数的。比如TPM、RPKM/FPKM或DESeq2标准化后的值均可。

      • 将此文件另存为文本文件(制表符分隔)(*.txt),然后再将后缀名称改为.gct即可

    • 若在bash环境中生成该文件(使用windows的请忽略此处),使用cat -A查看这个文件,输出到屏幕上的应是:

    #1.2$
    25000^I6$
    NAME^IDESCRIPTION^Imut_1^Imut_2^Imut_3^Iwt_1^Iwt_2^Iwt_3$
    AT3G45264^Ina^I296472.8029^I298638.6905^I280473.9967^I385382.3187^I355358.3136^I377348.3648$
    AT5G38382^Ina^I164223.8048^I169543.267^I164432.3066^I124346.784^I120357.6296^I123873.0918$
    

    符号^I表示此处是一个Tab键。生成文件后将后缀名称改为.gct

    基因集数据库文件(gene sets database)

    如果研究的物种是人,在MSigDB数据库中下载感兴趣的基因集即可。地址:http://software.broadinstitute.org/gsea/downloads.jsp

    其中含有GO、KEGG、Reactome、hallmark、BioCarta及MSigDB汇总的基因集,其中GO基因集又分为mf、bp和cc,且基因名称包含Entrez IDs和gene symbols等。因此在选择的时候需要仔细甄别一下再下载,且注意保证注释集中的基因名称与表达矩阵中基因名称类型一致

    基因集数据库文件的格式:

    img

    • 第一列是某一个生物过程、细胞组分、分子功能或信号通路等等,是一个基因集(gene set)的名称
    • 第二列是该基因集的编号或其它描述信息均可
    • 第三列及以后的每列是从属于该基因集的基因,有的是一个,有的是两个或多个。你也可以输入自己感兴趣的基因,只需满足以上格式即可。
    • 列与列之间以Tab键隔开,基因与基因之间也以Tab键隔开
    • 文件后缀名写为.gmt

    本文分析的对象是拟南芥,MSigDB并无提供,因此需要自己去GO数据库或拟南芥信息资源网站(The Arabidopsis Information Resource, TAIR)上下载。再将其加工修改为一个gmt文件。以bash环境为例:

    wget https://www.arabidopsis.org/download_files/GO_and_PO_Annotations/Gene_Ontology_Annotations/ATH_GO_GOSLIM.txt
    
    # 提取biological process分类
    # 注意按照自己的文件筛选合适的列,修改下面的$8,$1,$5,$6.
    awk 'BEGIN{OFS=FS="\t"}{if($8=="P") print $1,$5,$6}' ATH_GO_GOSLIM.txt | sort -u | awk 'BEGIN{OFS=FS="\t"}{anno=$2"\t"$3; a[anno]=a[anno]==""?$1:a[anno]"\t"$1}END{for(i in a) print i,a[i]}'  >  ATH_GO_GOSLIM_LocusName_bp_useme.gmt
    

    这个ATH_GO_GOSLIM_LocusName_bp_useme.gmt就是一个新建的基因集数据库文件。注意其中的基因名称需要与已经做好的基因表达文件中的基因名称对应

    使用javaGSEA软件开始分析

    Java安装或更新后即可双击gsea_4096m.jnlp启动GSEA,如下图:

    img

    • 点击1处的Load data后,点击2处的Browse for files;将已经准备好的样本分组文件(.cls)、基因表达文件(.gct)和基因集数据库文件(.gmt)一并上传。注意:如果分析的是人的基因数据,且调用GSEA数据库中的gmt文件,软件会自动下载,此处可以无需上传本地的gmt文件。

    • 上传后点击3处的Run GSEA,弹出下图。在4处的Expression dataset选择上传的基因表达文件(.gct文件)。在gene sets database右侧点击复选框,选择Gene Matrix(local gmx/gmt),选择已上传的基因集数据库文件(.gmt)。然后回来,在Phenotype labels右侧点击复选框,选择已上传的样本分组文件(.cls)并选择一个phenotype

    图6

    img

    • Collapse dataset to gene symbols处选择false,这是因为在基因表达文件(.gct)和基因集数据库文件(.gmt)中已经有对应好的基因名称,不需要再做转换。

    • Permutation type处选择gene_set

    • 6处即可点击Run来运行GSEA。

    • 7处会显示运行的状态running,当出现success后点点击 (等同于点击Show result folder并打开文件夹中的index.html),即可弹出分析报告,如下图。

    img

    • 点击图中Detailed enrichment results in html format可查看排序好的基因集富集结果的表格(依照NES值排序),并可点击Details..查看详情。
    • 点击Snapshot即可看到GSEA分析的结果图 (后面解释)。
    • Snapshot页面点击图片即可进入对应基因集的详细页面,其中的信息如GSEA Results SummaryEnrichment plotGSEA details(可看到该基因集的基因数目)等。

    img

    img

    图中1说明描绘的是对核糖体生物合成(ribosome biogenesis)基因集的GSEA分析结果。

    • 首先GSEA将所有基因按照与表型的相关度排序,获得一个排序后的基因列表,排序值以热图形式展示 (图中的4)。
    • 然后在排序好的基因列表从前向后挨个查看每个基因是否出现在此图对应的基因集中,若初选用竖线表示(图中的3(共119个基因)),并且增加一个统计值( a running-sum statistic),反之则减少一个统计值,最终每个位置的基因各自获得一个累计值(图中的绿色曲线)。
    • 每个基因对应的累计值就叫做富集得分 (Enrichment score, ES) (图8中的2),而这个基因集的富集得分 (ES)则定义为遍历基因列表时遇到的离零的最大偏差,即峰值 (此图中为0.597)。峰值为正值表示基因集富集在列表的顶部(mut),负值表示富集在底部(wt)。
    • 标准化富集评分(NES)是检验基因集富集结果的主要统计指标。GSEA富集分析时,不同的用户所输入的基因集数据库文件中的基因集数目可能不同。富集得分的标准化考虑了基因集个数和大小的差异。因此NES可以用来比较不同基因集之间的分析结果。NES计算方法如下:

    img

    因此,GSEA对基因集富集分析的排序 (如Snapshot页面或Detailed enrichment results in html format页面)是依据NES的值,而不是ES或p-val等值。这也就是为什么文章开头的那篇cell文章需要在两个基因集富集分析结果中(图1)明确标明NES值。

    • NOM p-value,即Nominal p value,是对富集得分(ES)的统计学分析,但未根据基因集的大小或多重假设检验来校正,因此在比较不同基因集的作用上有限。

    • FDR q-value,即假阳性率(False discovery rate)。是对标准化富集得分(NES)可能存在的假阳性结果的概率估计。因此FDR越小说明富集越可靠。FDR比NOM p-value更有参考意义 (图1也只标了FDR)。

    参考文献

    [1] Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005 Oct 25;102(43):15545-50. Epub 2005 Sep 30.

    [2] Ang YS, Rivas RN, Ribeiro AJS, Srivas R, Rivera J, Stone NR, Pratt K, Mohamed TMA, Fu JD4, Spencer CI, Tippens ND, Li M, Narasimha A, Radzinsky E, Moon-Grady AJ, Yu H, Pruitt BL, Snyder MP, Srivastava D. Disease Model of GATA4 Mutation Reveals Transcription Factor Cooperativity in Human Cardiogenesis. Cell. 2016 Dec 15;167(7):1734-1749.e22. doi: 10.1016/j.cell.2016.11.033.

    展开全文
  • 参考文章:GSEA学习笔记 老菜鸟做了几次GSEA分析,现在还没完全搞明白,之前写的记录有些杂乱,这里重新整理一下,以便以后学习。 1.使用GSEA软件进行绘图需要准备的文件。 该分析需要使用三个输入文件,分别是: ...
  • 一文掌握GSEA,超详细教程!

    千次阅读 2020-05-24 18:32:00
    生信宝典之前总结了一篇关于GSEA富集分析的推文——GSEA富集分析:从概念理解到界面实操,介绍了GSEA的定义、GSEA原理、GSEA分析、Leading-edge分析等,是全网最流行...
  • GSEA富集分析

    2021-05-05 10:53:14
    GSEA定义 Gene Set Enrichment Analysis (基因集富集分析)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势,从而判断其对表型的贡献。其输入数据包含两部分,一是已知功能的基因集 ...
  • R绘图实战|GSEA富集分析图

    千次阅读 2021-03-18 20:26:17
    后台难得有读者私信,请教了下图中文章的GSEA图能不能用R来画,今天就来简单写个教学。 GSEA(Gene Set EnrichmentAnalysis),即基因集富集分析,它的基本思想是使用预定义的基因,将基因按照在两类样本中的差异...
  • GSEA原理的通俗易懂的理解

    千次阅读 2020-07-01 15:13:59
    GSEA分析GSEA介绍GSEA原理数据矩阵文件GSEA计算中几个关键概念典型富集结果解读GSEA与传统的比较 GSEA介绍 我们先提出问题:在解读传统的富集分析(基于超几何分布或Fisher检验)结果时,经富集分析筛选的功能通路中...
  • omp算法matlab代码GSEA 用于LINCS数据的基因集富集分析(GSEA)工具以并行方式 内容: 1. Description. 2. Benchmark. 3. Compilation and installation. 4. Description of Tools. 4.1 Matlab Tools: matlab_for_...
  • 基因富集分析 GSEA for time-course

    千次阅读 2019-09-17 21:37:44
    基因富集分析(Gene Set Enrichment Analysis,GSEA)是一种针对全基因组表达谱芯片数据的分析方法,将基因与预定义的基因集进行比较。即综合现有的对基因的定位、性质、功能、生物学意义等信息基础,构建一个分子...
  • 今天在讨论群看到有群友提问 单基因GSEA怎么做?。之前也看到过这个概念,但一直不清楚这个单是什么含义,一直以为是用单个基因做GSEA。如果之前看过生信宝典的一文掌握GSEA,超详细教程,一定会特别熟悉GSEA的原理...
  • 今天和大家分享的是2020年一月份发表在Front Bioeng Biotechnol(IF:5.122)的一篇文章,作者在主要针对头颈癌,从表达量入手,对PGRMC1高低表达组中的差异基因近行GO,GSEA等分析,并与头颈癌患者的临床信息相联系,...
  • GSEA分析结果详细解读

    万次阅读 多人点赞 2018-11-05 10:49:41
    还是这张原理图,GSEA的输入是一个基因表达量矩阵,其中的样本分成了A和B两组,首先对所有基因进行排序,在之前的文章中也有提到排序的标准,这里简单理解就是foldchange, 用来表示基因在两组间表达量的变化趋势。...
  • 常规的高通量数据分析思路是一种趋势分析。即从几万个基因中,通过逻辑的思路,一步一步的缩小范围,最终找到与疾病关联的关键基因/通路。这也是我们生信分析中...在这里给大家推荐一个软件GSEA,可以非常好的解决...
  • 我在使用R语言进行GSEA分析时提示这个请教大家是哪里出了问题呢? <p style="text-align:center"><img alt="" src="https://img-ask.csdnimg.cn/upload/1611716066341.png" /></p> ...
  • GSEA需要的gene set是现成的没有问题,但是genelist没有,这里我们可以把所有基因跟单个基因的相关性系数当做LogFC,有正有负,就解决了geneList的问题。这个想法不是我的,是我的一个学员的,不过他要解决的是...
  • 软件包旨在更好地概述GSEA。 目录: 安装 从git仓库安装软件包: devtools :: install_github( " nicolash2/ggbrace " ) 浓缩GSEA condenseGsea采用data.frame,该数据框必须具有一个带有基因的列。 这些可能是...
  • 5.DNp63与TAp63的相关信号通路预测 分别做BLCA、BRCA、LUSC人群中高表达DNp63和高表达TAp63的队列的GSEA基因富集分析,(用Wilcoxon-P值来对GSEA的基因集进行排序,并识别统计富集的基因集)。发现DNp63高表达的肿瘤...
  • GSEA的使用

    2019-06-18 17:40:00
    java -Xmx2048m -cp /biosoftware/gsea-3.0.jar xtools.gsea.Gsea -res Flox.HFHC.16W_vs_CKO.HFHC.16W.fpkm.txt -cls Flox.HFHC.16W_vs_CKO.HFHC.16W.class.cls -gmx mmu.kegg_optimal.symbol.gmt -out Flox.HFHC....
  • 现在都跑GSEA!” 我睁开我蒙昧的小眼睛:“老师,啥叫GSEA啊?” 老师愣了一下,“这么简单都不会,自己查查去”。 。。。 我。。。。我的老板应该不知道GSEA是什么。。。 GSEA定义 Gene Set Enrichment Analysis ...
  • Single sample gene set enrichment analysis(ssGSEA) ssGSEA(single sample GSEA)主要针对单样本无法做GSEA而提出的一种实现方法,原理上与GSEA是类似的。 ssGSEA根据表达谱文件计算每个基因基因表达值的rank值,...
  • 方法:在本研究中,从人类自噬数据库(HADb)和基因集富集分析(GSEA)网站下载了自噬相关基因。采用最小绝对收缩和选择算子(LASSO)回归和多变量Cox回归分析来确定构建风险特征的基因。将风险特征与临床病理因素...
  • 基于ssGSEA(单样本GSEA)的免疫基因集文章套路 --生信自学网光俊今天我们给大家介绍下生信自学网的” 基于ssGSEA的免疫基因集文章套路”课程,该课程根据最新发表的5.6分纯生信文章录制。该课程的核心思想是通过...
  • 然后,从已发表的数据中得到了HDR评分(想要的小伙伴去文章的参考文献找哈)和DDR通路分布,并通过R包GSVA采用单样本GSEA计算了一组调节DNA损伤修复(DDR)通路的基因集(共276个基因)的富集分数,然后R包pheatmap被用来...
  • 作者使用GSEA来进行KEGG和GO基因集富集分析,对比高风险组与低风险组之间的功能表型,发现5个KEGG通路(neuroactive ligand receptor interaction, calcium signaling pathway, long term potentiation, type II ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 705
精华内容 282
关键字:

gsea