用limma包进行多组差异表达分析_庄无因的博客-CSDN博客_limma 差异表达分析

用limma包进行多组差异表达分析

置顶
庄无因
2020-03-15 00:20:16

4026

收藏

16

分类专栏:
生物信息学基础知识

写在前面:最近在使用limma包进行差异表达分析,参考了网上许多教程都觉得说的云里雾里,很不清楚。经过我自己一段时间非常痛苦的钻研,弄明白了,解决了我的实际需求。于是决定将我的分析经验写下来,分享给需要的人。

 首先加载前期预处理好的表达矩阵。(我的原始表达矩阵文件已附在文后,大家有需要可以下载实践)

  1. 1
    setwd("E://Rworkplace")
  2. 2
    data<-read.table("GSE98793.txt")

表达矩阵格式如下图所示,行名为基因名,列名为样本名。

创建样本分类信息表。

所谓样本分类信息表,实际上就是值你的样本名一一对应的属性是什么(譬如case还是control)。

我自己根据GEO数据库中该芯片的分类信息先自己在excel中整理了一张表。因为这张芯片原作者在上传数据的时候,样本分类数据很乱,实验组和对照组的样本有时候交错分布,譬如ID为GSM2612287的样本是第192个样本,是对照组即健康人群。而ID为GSM2612128的样本是第33个样本,取样自重症抑郁症患者。我们在对样本进行差异表达分析之前,一定要做的工作是弄明白每一个样本的属性,我做这个工作,实际上是对样本的属性进行了手工注释。

然后我将上图中的第二第三列,复制单独保存在一个txt文档中。(即gene_list.txt)

第一列(1,2,3,……)只是为了辅助排序的东西,因为这个1,2,3,就是样本在总的样本中的位置,这个表格就是告诉我们第几个样本,它是对照组,第几个样本,取样自重症抑郁症患者。

1 HC
2 HC
3 HC
4 HC
…… ……

加载上述分类表gene_list.txt。并按照数据框的第一列升序排序。

取第二列,这样处理得到的就是每一个位置所对应的样本属性信息。我们最终想要得到的也就是向量group。(实际情况大家可能不会像我这样麻烦,具体问题具体分析,只要得到一个类如group的向量)

  1. 1
    group_list<-read.table("gene_list.txt")
  2. 2
    new_group<- group_list[order(group_list[,1]),]
  3. 3
    group<-new_group[,2]

注:其中HC代表的是健康对照组,MDD指重症抑郁症;anxiety,即焦虑症。

 

接着,加载limma包,构造如下矩阵。(这一段可以照抄。根据自己的实际情况替换变量)

  1. 1
    suppressMessages(library(limma))
  2. 2
    design <- model.matrix(~0+factor(group))
  3. 3
    colnames(design)=levels(factor(group))
  4. 4
    rownames(design)=colnames(data)

 

然后,我们规定哪一组数据与哪一组数据比较。

这里也是原来我比较困惑的地方,因为我的数据是有MDD,anxiety两个实验组,HC对照组。而我实际只需要比较MDD与HC的差异表达情况,而一般的教程中并未涉及或者介绍模糊。

下面这段的代码的意思就是比较design矩阵中HC与MDD level。(关于多组分析的其它疑惑,可以参考知乎文章

contrast.matrix<-makeContrasts("HC-MDD",levels=design)

完成以上两个表格之后,就可以直接跑下面的代码,得到差异表达数据了。(有时间具体分析什么意思)

  1. 1
    ##step1
  2. 2
    fit <- lmFit(data,design)
  3. 3
    ##step2
  4. 4
    fit2 <- contrasts.fit(fit, contrast.matrix)
  5. 5
    fit2 <- eBayes(fit2)
  6. 6
    ##step3
  7. 7
    tempOutput = topTable(fit2, coef=1, n=Inf)
  8. 8
    nrDEG = na.omit(tempOutput)

分析得到的nrDEG如下图所示:

 保存分析好的差异表达数据。

write.table(nrDEG,"limma_GSE98793.txt")

 全部代码:

  1. 1
    setwd("E://Rworkplace")
  2. 2
    data<-read.table("GSE98793.txt")
  3. 3
    group_list<-read.table("gene_list.txt")
  4. 4
    new_group<- group_list[order(group_list[,1]),]
  5. 5
    group<-new_group[,2]
  6. 6
  7. 7
    suppressMessages(library(limma))
  8. 8
  9. 9
  10. 10
    design <- model.matrix(~0+factor(group))
  11. 11
    colnames(design)=levels(factor(group))
  12. 12
    rownames(design)=colnames(data)
  13. 13
  14. 14
    contrast.matrix<-makeContrasts("HC-MDD",levels=design)
  15. 15
  16. 16
    ##step1
  17. 17
    fit <- lmFit(data,design)
  18. 18
    ##step2
  19. 19
    fit2 <- contrasts.fit(fit, contrast.matrix)
  20. 20
    fit2 <- eBayes(fit2)
  21. 21
    ##step3
  22. 22
    tempOutput = topTable(fit2, coef=1, n=Inf)
  23. 23
    nrDEG = na.omit(tempOutput)
  24. 24
    write.table(nrDEG,"limma_GSE98793.txt")

 熬夜写到这儿,一定会有些疏漏和错误,欢迎指正,共同进步!

【数据链接】https://me.csdn.net/download/weixin_40640700

重新修改了上传数据所需的积分,设置为零。可是系统还是自动调为了1。

如果参考数据不能顺利下载的话,可以邮箱2456392738@qq.com,我邮箱发送。

点击数:0

R包clusterProfiler: 转换ID、GO/KEGG富集分析

R包clusterProfiler: 转换ID、GO/KEGG富集分析

12020.04.23 16:21:40字数 1,707阅读 3,773

参考这篇
ID转换不用怕,clusterProfiler帮你忙
这篇详细讲了无参考基因组该怎么办,值得一看

简单介绍一下几种常用的ID:
Ensemble id:一般以ENSG开头,后边跟11位数字。如TP53基因:ENSG00000141510
Entrez id:通常为纯数字。如TP53基因:7157
Symbol id:为我们常在文献中报道的基因名称。如TP53基因的symbol id为TP53
Refseq id:NCBI提供的参考序列数据库:可以是NG、NM、NP开头,代表基因,转录本和蛋白质。如TP53基因的某个转录本信息可为NM_000546

1 转换ID

1.1 载入包library(clusterProfiler)

library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)
246

org.Hs.eg.db是人类基因组注释,如果需要其他物种的注释信息,可以去bioconductor这里获得

1.2. keytypes()查看注释包中支持的ID转换类型,基本包括了所以常用的ID类型

keytypes()

1.3. bitr()函数转换ID

函数全部内容如下

bitr(geneID, fromType, toType, OrgDb, drop = TRUE
2

geneID:输入待转换的geneID
fromType:输入的ID类型
toType:希望输出的ID类型
OrgDb:注释对象的信息
drop:drop = FALSE 保留空值

TP53为例,希望输出’ENTREZID’,’ENSEMBL’,’REFSEQ’

my_id <- c("TP53")

output <- bitr(my_id,
               fromType = 'SYMBOL',
               toType = c('ENTREZID','ENSEMBL','REFSEQ'),
               OrgDb = 'org.Hs.eg.db')
24681012
转换后的结果

1.4. bitr_kegg()函数进行基因ID与蛋白质ID的转换

函数全部内容如下,与bitr()类似

bitr_kegg(geneID, fromType, toType, organism, drop = TRUE)
2

在参数这里与bitr()有些区别!
geneID:输入待转换的geneID
fromType:输入的ID类型只能是kegg(与entrez id相同), ncbi-geneid, ncbi-proteinid or uniprot中的一个
toType:希望输出的ID类型,也只能是kegg(与entrez id相同), ncbi-geneid, ncbi-proteinid or uniprot中的一个
organismhsa,代表人类,其他的物种名称可以参考kegg的网站
drop:去除空值与否

TP53entrez id7157

kegg <- bitr_kegg(7157,
                   fromType = 'kegg',
                   toType = c('ncbi-geneid', 'ncbi-proteinid', 'uniprot'),
                   organism = 'hsa')
2468

报错了,所以一次应该只能进行一种转换!

改成这样↓

kegg <- bitr_kegg(7157,
                  fromType = 'kegg',
                  toType = 'ncbi-proteinid',
                  organism = 'hsa')

246810

转换为uniprot

kegg <- bitr_kegg(7157,
                  fromType = 'kegg',
                  toType = 'uniprot',
                  organism = 'hsa')
2468

出现了三个uniprot

UniProt网站上查询一下为什么会有三个:
P04637显示的是TP53基因的蛋白质表达水平,级别是Reviewed,就是其来源为UniProtKB/Swiss-Prot。

P04637

K7PPA8显示的是转录本水平的表达,级别是Unreviewed,就是其来源为UniProtKB/TrEMBL

K7PPA8

Q53GA5显示的是转录本水平的表达,级别是Unreviewed,就是其来源为UniProtKB/TrEMBL

Q53GA5

一般ID转换仅仅为开始的准备工作,将自己的数剧转换好之后可以进行后续的分析。
另外,利用clusterProfiler包可以进行许多丰富的下游分析,比如GO分析KEGG分析等等

2 富集分析

参考这篇
最常用的基因注释工具是GOKEGG注释,这基本上是和差异基因分析一样,必做的事,GO可在BP(生物过程),MF(分子功能),CC(细胞组分)三方面进行注释

GO中的注释信息

KEGG数据库全部的物种及其简写(三个字母)

有很多在线(DAVIDGatherGOrilla,revigo)或者客户端版的工具(IPA(IPA不是用的GO和KEGG数据库)和FUNRICH)可以用来分析差异基因,我主要实践一下clusterProfiler

clusterProfiler提供有两种富集方式:

① ORA(Over-Representation Analysis)差异基因富集分析(不需要表达值,只需要gene name)
② FCS,它的代表就是GSEA,即基因集(gene set)富集分析(不管有无差异,需要全部genes表达值)

① ORA 差异基因富集分析(不需要表达值,只需要gene name)

1) 先读进来差异基因
sig.gene <- read.delim("你的路径/final_result.txt(之前计算出的差异基因文件)",  sep = 't', stringsAsFactors = FALSE, check.names = FALSE)
gene<-sig.gene[,1]
head(gene)

#转换ID
gene.df<-bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Hs.eg.db,
              drop = FALSE)

head(gene.df)
246810121416182022
2) GO富集
ego_cc<-enrichGO(gene = gene.df$ENSEMBL,
                 OrgDb  = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont  = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)

ego_bp<-enrichGO(gene = gene.df$ENSEMBL,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont  = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
24681012141618202224262830

参数意义:
keyType 指定基因ID的类型,默认为ENTREZID
OrgDb 指定该物种对应的org包的名字
ont 代表GO的3大类别,BP, CC, MF,也可是全部ALL

pAdjustMethod 指定多重假设检验矫正的方法,有“ holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一种
cufoff 指定对应的阈值
readable=TRUE 代表将基因ID转换为gene symbol。

barplot() #富集柱形图
dotplot() #富集气泡图
cnetplot() #网络图展示富集功能和基因的包含关系
emapplot() #网络图展示各富集功能之间共有基因关系
heatplot() #热图展示富集功能和基因的包含关系

3) barplot()作图显示,依赖于ggplot2包
library(ggplot2)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")  
24

ggplot2参数意义具体可以看这篇

barplot()

dotplot()

emapplot()
4) 或者是KEGG分析

输入的gene要为vector格式的entrezid,可以用is.vector()判断是否是vector格式

library(stringr)

kk<-enrichKEGG(gene = gene.df$ENTREZID,
               organism = "hsa",
               keyType = "kegg",          
               pAdjustMethod = "BH",
               pvalueCutoff = 0.01)
2468101214

KEGG

但是这样产生的kk几乎都是Na,更别提画图了
于是我去看了Y叔的手册,发现需要使用DOSE构建,改成

library(DOSE)

data(geneList, package = "DOSE")
gene <- names(geneList)[abs(geneList) > 2]

kk<-enrichKEGG(gene = gene,
               organism = "hsa",
               keyType = "kegg",
               pAdjustMethod = "BH",
               pvalueCutoff = 0.01)
2468101214161820

就可以了,继续往下进行

#bar图
barplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
#点图
dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
246810121416
barplot()

dotplot()

cnetplot()

heatplot()

② GSEA 基因集(gene set)富集分析(不管有无差异,需要全部genes表达值)

这篇为网页版操作
好处:可以发现被差异基因舍弃的genes可能参与了某重要生理过程或信号通路
要点:用所有差异基因来跑GSEA,而不是按筛选条件筛选出来的,否则GSEA就失去了意义

1)读入全部差异基因数据
library(dplyr)
all_deseq <- read.delim("你的路径/DESeq2-res.txt",  sep = 't', stringsAsFactors = FALSE, check.names = FALSE)
24
部分内容
2) 提取出gene_idlog2FoldChange
ensemblid_log2fdc <-dplyr::select(all_deseq, gene_id, log2FoldChange)
2
3) 可以看到这时gene_id还不是整数形式,所以对这列gsub()取整,命名为symbol列添加到ensemblid_log2fdc
ensemblid_log2fdc$symbol <- gsub("\.\d*", "", all_deseq[,2])
2
4)转换id:将ENSEMBLID转为ENTREZID
entrezid <-bitr(ensemblid_log2fdc$symbol, 
              fromType = "ENSEMBL", 
              toType = c("ENTREZID"),
              OrgDb = org.Hs.eg.db,
              drop = TRUE)
246810
5)因为6)中merge需要相同列名,所以将symbol列列名改为ENSEMBL
names(ensemblid_log2fdc) <- c("gene_id", "log2FoldChange", "ENSEMBL")
2
6)ensemblid_log2fdcentrezid拼接到一起,依靠的是它们具有相同列名的列——ENSEMBL
final <- merge(ensemblid_log2fdc, entrezid)
2
7) 提取entrezidlog2fdc
library(dplyr)    
entrezid_log2fdc <- dplyr::select(final, ENTREZID, log2FoldChange)
24
8) 用arrange()log2fdc排序,desc()表示降序
final.sort <- arrange(entrezid_log2fdc, desc(log2FoldChange))
2
9) 构建一下要GSEA分析的数据
genelist <- final.sort$log2FoldChange
names(genelist) <- final.sort[,1]
24
10) 进行GSEA分析,需要entrezid,因此我前面才会进行ID转化
gsemf <- gseGO(genelist,
           OrgDb = org.Hs.eg.db,
           keyType = "ENTREZID",
           ont="BP",
           pvalueCutoff = 1)
head(gsemf)
24681012
11)画出GSEA图
gseaplot(gsemf, 1)
2
GSEA图

3 可视化

1.GO富集分析结果可视化

barplot

barplot(ego, showCategory = 10) #默认展示显著富集的top10个,即p.adjust最小的10个

dotplot

dotplot(ego, showCategory = 10)

DAG有向无环图

plotGOgraph(ego) #矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色,对应p值从大到小。

igraph布局的DAG

goplot(ego)

GO terms关系网络图(通过差异基因关联)

emapplot(ego, showCategory = 30)

GO term与差异基因关系网络图

cnetplot(ego, showCategory = 5)

2.Pathway富集分析结果可视化

barplot

barplot(kk, showCategory = 10)

dotplot

dotplot(kk, showCategory = 10)

pathway关系网络图(通过差异基因关联)

emapplot(kk, showCategory = 30)

pathway与差异基因关系网络图

cnetplot(kk, showCategory = 5)

pathway映射

browseKEGG(kk, “hsa04934”) #在pathway通路图上标记富集到的基因,会链接到KEGG官网

禁止转载,如需转载请通过简信或评论联系作者。

25人点赞

更多精彩内容下载简书APP
“小礼物走一走,来简书关注我”

还没有人赞赏,支持一下

总资产2 (约0.15元)共写了1.3W字获得98个赞共87个粉丝

点击数:0

R语言:两组数据关联分析,pheatmap热图和cytoscape网络图 – 简书

R语言:两组数据关联分析,pheatmap热图和cytoscape网络图

0.9642020.05.28 17:56:18字数 543阅读 2,781

导读

举例展示R语言组学关联分析的方法。宏基因组数据以KO-样品丰度表为例。代谢组数据以metabolite-样品丰度表为例。基本方法是用R语言psych包corr.test函数进行两组数据的相关分析,结果经格式化后用pheatmap可视化得热图。

一、模拟输入

1. KO丰度表

ko_abun = as.data.frame(matrix(abs(round(rnorm(200, 100, 10))), 10, 20))
colnames(ko_abun) = paste("KO", 1:20, sep="_")
rownames(ko_abun) = paste("sample", 1:10, sep="_")
ko_abun
1234

2. metabolite丰度表

metabo_abun = as.data.frame(matrix(abs(round(rnorm(200, 200, 10))), 10, 20))
colnames(metabo_abun) = paste("met", 1:20, sep="_")
rownames(metabo_abun) = paste("sample", 1:10, sep="_")
metabo_abun
1234

二、相关分析函数

思路

  1. 参数一:other -> KO或其他组学丰度表
  2. 参数二:metabo -> 代谢物丰度表
  3. 参数三:route -> 输出目录【提前创建】
  4. corr.test进行两组数据相关分析
  5. 用stringr split将ko-metabolite结果列拆成两列
  6. 结果保留r_value p_value
library(psych)
library(stringr)

correlate = function(other, metabo, route)
{
    #  读取方式:check.name=F, row.names=1, header=T
    # 计算相关性:
    #other = data
    #metabo = env
    #route="gut"
    result=data.frame(print(corr.test(other, metabo, use="pairwise", method="spearman", adjust="fdr", alpha=.05, ci=TRUE, minlength=100), short=FALSE, digits=5))
    # FDR矫正
    result_raw=data.frame(print(corr.test(other, metabo, use="pairwise", method="spearman", adjust="none", alpha=.05, ci=TRUE, minlength=100), short=FALSE, digits=5))
    # 原始P value

    # 整理结果
    pair=rownames(result)  # 行名
    result2=data.frame(pair, result[, c(2, 4)])  # 提取信息

    # P值排序
    # result3=data.frame(result2[order(result2[,"raw.p"], decreasing=F),])

    # 格式化结果【将细菌代谢物拆成两列】
    result4=data.frame(str_split_fixed(result2$pair, "-", 2), result2[, c(2, 3)], p_value=result_raw[, 4])
    colnames(result4)=c("feature_1", "feature_2", "r_value", "fdr_p_value", "raw_p_value")

    # 保存提取的结果
    write.table(result4, file=paste(route, "Correlation_result.txt", sep="/"), sep="t", row.names=F, quote=F)
}
1234567891011121314151617181920212223242526272829

三、相关性分析

dir.create("Result")  # 创建结果目录
correlate(ko_abun, metabo_abun, "Result")
12
  • 结果Correlation_result.txt如下,行数为20X20

四、pheatmap可视化函数

思路:

  1. 参数一:infile -> Correlation_result.txt相关性分析结果
  2. 参数二:route -> Route输出路径
  3. 用reshape2 dcast函数把feature_1 feature_2 p_value做成matrix,作为pheatmap输入文件
  4. 用reshape2 dcast函数把feature_1 feature_2 r_value做成matrix,作为pheatmap display
  5. 显著相关标记* 0.05>=p>0.01; ** 0.01>=p>0.001; *** 0.001>=p
library(reshape2)
library(pheatmap)

correlate_pheatmap = function(infile, route)
{
    data=read.table(paste(route, infile, sep="/"), sep="t", header=T)

    data_r=dcast(data, feature_1 ~ feature_2, value.var="r_value")
    data_p=dcast(data, feature_1 ~ feature_2, value.var="raw_p_value")
    rownames(data_r)=data_r[,1]
    data_r=data_r[,-1]
    rownames(data_p)=data_p[,1]
    data_p=data_p[,-1]
    
    # 剔除不显著的行
    del_row = c()
    for(i in 1:length(data_p[, 1]))
    {
        if(all(data_p[i, ] > 0.05))
        {
            del_row = c(del_row, i)
        }
    }

    # 剔除不显著的列
    del_col = c()
    for(j in 1:length(data_p[1, ]))
    {
        if(all(data_p[, j] > 0.05))
        {
            del_col = c(del_col, j)
        }
    }
    
    # null值处理
    if(is.null(del_row) && !(is.null(del_col)))
    {
        data_p = data_p[, -del_col]
        data_r = data_r[, -del_col]
    }else if(is.null(del_col) && !(is.null(del_row)))
    {
        data_p = data_p[-del_row,]
        data_r = data_r[-del_row,]
    }else if(is.null(del_row) && is.null(del_col))
    {
        print("delete none")
    }else if(!(is.null(del_row)) && !(is.null(del_col)))
    {
        data_p = data_p[-del_row, -del_col]
        data_r = data_r[-del_row, -del_col]
    }
    
    # data_p = data_p[-del_row, -del_col]
    # data_r = data_r[-del_row, -del_col]
    write.csv(data_p, file=paste(route, "data_p.csv", sep="/"))
    write.csv(data_r, file=paste(route, "data_r.csv", sep="/"))

    # 用"*"代替<=0.05的p值,用""代替>0.05的相对丰度
    data_mark=data_p
    for(i in 1:length(data_p[,1])){
        for(j in 1:length(data_p[1,])){
            #data_mark[i,j]=ifelse(data_p[i,j] <= 0.05, "*", "")
            if(data_p[i,j] <= 0.001)
            {
                data_mark[i,j]="***"
            }
            else if(data_p[i,j] <= 0.01 && data_p[i,j] > 0.001)
            {
                data_mark[i,j]="**"
            }
            else if(data_p[i,j] <= 0.05 && data_p[i,j] > 0.01)
            {
                data_mark[i,j]="*"
            }
            else
            {
                data_mark[i,j]=""
            }
        }
    }
    write.csv(data_mark, file=paste(route, "data_mark.csv", sep="/"))

    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.pdf", sep="/"))
    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.png", sep="/"))
}
12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485

五、pheatmap绘制热图

correlate_pheatmap("Correlation_result.txt", "Result")
1
  • 结果目录,新增结果图(一个png一个pdf):
  • 打开pdf,如下:
  • 随机模拟的数据,没有显著的不奇怪。如果是做项目一般两组数据的相关分析都可以得到一些显著相关的结果,如下:

两组数据完全相同的话:

六、绘图展示 |r| > 0.8的相关

1 函数

select_pheatmap = function(infile, route)
{
    data=read.table(paste(route, infile, sep="/"), sep="t", header=T)

    data_r=dcast(data, feature_1 ~ feature_2, value.var="r_value")
    data_p=dcast(data, feature_1 ~ feature_2, value.var="raw_p_value")
    rownames(data_r)=data_r[,1]
    data_r=data_r[,-1]
    rownames(data_p)=data_p[,1]
    data_p=data_p[,-1]
    
    # 记下0.8 > r > -0.8 的行
    del_row = c()
    for(i in 1:length(data_p[, 1]))
    {
        if(all(data_r[i, ] < 0.8 & data_r[i, ] > -0.8))
        {
            del_row = c(del_row, i)
        }
    }

    # 记下0.8 > r > -0.8 的列
    del_col = c()
    for(j in 1:length(data_p[1, ]))
    {
        if(all(data_r[, j] < 0.8 & data_r[, j] > -0.8))
        {
            del_col = c(del_col, j)
        }
    }
    
    # null值处理
    if(is.null(del_row) && !(is.null(del_col)))
    {
        data_p = data_p[, -del_col]
        data_r = data_r[, -del_col]
    }else if(is.null(del_col) && !(is.null(del_row)))
    {
        data_p = data_p[-del_row,]
        data_r = data_r[-del_row,]
    }else if(is.null(del_row) && is.null(del_col))
    {
        print("delete none")
    }else if(!(is.null(del_row)) && !(is.null(del_col)))
    {
        data_p = data_p[-del_row, -del_col]
        data_r = data_r[-del_row, -del_col]
    }
    
    # data_p = data_p[-del_row, -del_col]
    # data_r = data_r[-del_row, -del_col]
    write.csv(data_p, file=paste(route, "data_p.csv", sep="/"))
    write.csv(data_r, file=paste(route, "data_r.csv", sep="/"))

    # 用"*"代替<=0.05的p值,用""代替>0.05的相对丰度
    data_mark=data_p
    for(i in 1:length(data_p[,1])){
        for(j in 1:length(data_p[1,])){
            if(data_p[i,j] <= 0.001)
            {
                data_mark[i,j]="***"
            }
            else if(data_p[i,j] <= 0.01 && data_p[i,j] > 0.001)
            {
                data_mark[i,j]="**"
            }
            else if(data_p[i,j] <= 0.05 && data_p[i,j] > 0.01)
            {
                data_mark[i,j]="*"
            }
            else
            {
                data_mark[i,j]=""
            }
        }
    }
    # r值mark 
    data_mark_r = data_r
    for(i in 1:length(data_r[,1])){
        for(j in 1:length(data_r[1,])){
            if(data_r[i,j] >= 0.8 | data_r[i,j] <= -0.8)
            {
                data_mark_r[i,j]=round(data_r[i,j], 2)
            }
            else
            {
                data_mark_r[i,j]=""
            }
        }
    }
    write.csv(data_mark, file=paste(route, "data_mark_p.csv", sep="/"))
    write.csv(data_mark_r, file=paste(route, "data_mark_r.csv", sep="/"))

    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=17, filename=paste(route, "Correlation_result_p.pdf", sep="/"))
    pheatmap(data_r, display_numbers=data_mark_r, cellwidth=20, cellheight=20, fontsize_number=7, filename=paste(route, "Correlation_result_r.pdf", sep="/"))
        
    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=17, filename=paste(route, "Correlation_result_p.png", sep="/"))
    pheatmap(data_r, display_numbers=data_mark_r, cellwidth=20, cellheight=20, fontsize_number=7, filename=paste(route, "Correlation_result_r.png", sep="/"))
}
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899

2 应用

dir.create("metabo_species_0.8r")
select_pheatmap("Correlation_result.txt", "metabo_species_0.8r")
12

七、cytoscape绘图文件

1 准备文件

# 读取关联分析结果
data = read.table("Result/Correlation_result.txt", sep="t", header=T)
# 画图文件
data = data[data$raw_p_value <= 0.05,]
r_label = c()
for(i in 1:nrow(data))
{
    if(data[i, 3] < 0)
    {
        r_label = c(r_label, "neg")
    }
    else
    {
        r_label = c(r_label, "pos")
    }
}

data$r_label = r_label
write.table(data, file="Result/input_pre.txt", sep="t", row.names=F, quote=F)
12345678910111213141516171819
data$r_value = abs(data$r_value)
input_network = data[, c(1,2,3,6)]
write.table(input_network, file="Result/input_network.txt", sep="t", row.names=F, quote=F)
123
g1 = data.frame(feature=input_network[,1], group=rep("tax", nrow(input_network)))
g2 = data.frame(feature=input_network[,2], group=rep("ko", nrow(input_network)))
input_group = unique(rbind(g1, g2))
write.table(input_group, file="Result/input_group.txt", sep="t", row.names=F, quote=F)
1234

2 cytoscape绘图

过程见:Cytoscape绘制相关网络图

3 几种常用的layout

layout -> group attributes layout -> group

layout -> group attributes layout -> degree

layout -> attribute circle layout -> degree

layout -> attribute circle layout -> group

于2020.9.21更新:fdr_p raw_p的对应关系
于2020.9.22更新:输出绘图文件data_p data_r data_mark
于2020.10.14更新:增加了根据R值筛选相关的select_pheatmap

17人点赞

R

更多精彩内容下载简书APP
“小礼物走一走,来简书关注我”

还没有人赞赏,支持一下

小白菜学生信简书记录成长<br>欢迎简信打扰
总资产78 (约4.30元)共写了10.7W字获得992个赞共981个粉丝

点击数:0

R语言:两组数据关联分析,pheatmap热图和cytoscape网络图 – 简书

R语言:两组数据关联分析,pheatmap热图和cytoscape网络图

0.9642020.05.28 17:56:18字数 543阅读 2,775

导读

举例展示R语言组学关联分析的方法。宏基因组数据以KO-样品丰度表为例。代谢组数据以metabolite-样品丰度表为例。基本方法是用R语言psych包corr.test函数进行两组数据的相关分析,结果经格式化后用pheatmap可视化得热图。

一、模拟输入

1. KO丰度表

ko_abun = as.data.frame(matrix(abs(round(rnorm(200, 100, 10))), 10, 20))
colnames(ko_abun) = paste("KO", 1:20, sep="_")
rownames(ko_abun) = paste("sample", 1:10, sep="_")
ko_abun
1234

2. metabolite丰度表

metabo_abun = as.data.frame(matrix(abs(round(rnorm(200, 200, 10))), 10, 20))
colnames(metabo_abun) = paste("met", 1:20, sep="_")
rownames(metabo_abun) = paste("sample", 1:10, sep="_")
metabo_abun
1234

二、相关分析函数

思路

  1. 参数一:other -> KO或其他组学丰度表
  2. 参数二:metabo -> 代谢物丰度表
  3. 参数三:route -> 输出目录【提前创建】
  4. corr.test进行两组数据相关分析
  5. 用stringr split将ko-metabolite结果列拆成两列
  6. 结果保留r_value p_value
library(psych)
library(stringr)

correlate = function(other, metabo, route)
{
    #  读取方式:check.name=F, row.names=1, header=T
    # 计算相关性:
    #other = data
    #metabo = env
    #route="gut"
    result=data.frame(print(corr.test(other, metabo, use="pairwise", method="spearman", adjust="fdr", alpha=.05, ci=TRUE, minlength=100), short=FALSE, digits=5))
    # FDR矫正
    result_raw=data.frame(print(corr.test(other, metabo, use="pairwise", method="spearman", adjust="none", alpha=.05, ci=TRUE, minlength=100), short=FALSE, digits=5))
    # 原始P value

    # 整理结果
    pair=rownames(result)  # 行名
    result2=data.frame(pair, result[, c(2, 4)])  # 提取信息

    # P值排序
    # result3=data.frame(result2[order(result2[,"raw.p"], decreasing=F),])

    # 格式化结果【将细菌代谢物拆成两列】
    result4=data.frame(str_split_fixed(result2$pair, "-", 2), result2[, c(2, 3)], p_value=result_raw[, 4])
    colnames(result4)=c("feature_1", "feature_2", "r_value", "fdr_p_value", "raw_p_value")

    # 保存提取的结果
    write.table(result4, file=paste(route, "Correlation_result.txt", sep="/"), sep="t", row.names=F, quote=F)
}
1234567891011121314151617181920212223242526272829

三、相关性分析

dir.create("Result")  # 创建结果目录
correlate(ko_abun, metabo_abun, "Result")
12
  • 结果Correlation_result.txt如下,行数为20X20

四、pheatmap可视化函数

思路:

  1. 参数一:infile -> Correlation_result.txt相关性分析结果
  2. 参数二:route -> Route输出路径
  3. 用reshape2 dcast函数把feature_1 feature_2 p_value做成matrix,作为pheatmap输入文件
  4. 用reshape2 dcast函数把feature_1 feature_2 r_value做成matrix,作为pheatmap display
  5. 显著相关标记* 0.05>=p>0.01; ** 0.01>=p>0.001; *** 0.001>=p
library(reshape2)
library(pheatmap)

correlate_pheatmap = function(infile, route)
{
    data=read.table(paste(route, infile, sep="/"), sep="t", header=T)

    data_r=dcast(data, feature_1 ~ feature_2, value.var="r_value")
    data_p=dcast(data, feature_1 ~ feature_2, value.var="raw_p_value")
    rownames(data_r)=data_r[,1]
    data_r=data_r[,-1]
    rownames(data_p)=data_p[,1]
    data_p=data_p[,-1]
    
    # 剔除不显著的行
    del_row = c()
    for(i in 1:length(data_p[, 1]))
    {
        if(all(data_p[i, ] > 0.05))
        {
            del_row = c(del_row, i)
        }
    }

    # 剔除不显著的列
    del_col = c()
    for(j in 1:length(data_p[1, ]))
    {
        if(all(data_p[, j] > 0.05))
        {
            del_col = c(del_col, j)
        }
    }
    
    # null值处理
    if(is.null(del_row) && !(is.null(del_col)))
    {
        data_p = data_p[, -del_col]
        data_r = data_r[, -del_col]
    }else if(is.null(del_col) && !(is.null(del_row)))
    {
        data_p = data_p[-del_row,]
        data_r = data_r[-del_row,]
    }else if(is.null(del_row) && is.null(del_col))
    {
        print("delete none")
    }else if(!(is.null(del_row)) && !(is.null(del_col)))
    {
        data_p = data_p[-del_row, -del_col]
        data_r = data_r[-del_row, -del_col]
    }
    
    # data_p = data_p[-del_row, -del_col]
    # data_r = data_r[-del_row, -del_col]
    write.csv(data_p, file=paste(route, "data_p.csv", sep="/"))
    write.csv(data_r, file=paste(route, "data_r.csv", sep="/"))

    # 用"*"代替<=0.05的p值,用""代替>0.05的相对丰度
    data_mark=data_p
    for(i in 1:length(data_p[,1])){
        for(j in 1:length(data_p[1,])){
            #data_mark[i,j]=ifelse(data_p[i,j] <= 0.05, "*", "")
            if(data_p[i,j] <= 0.001)
            {
                data_mark[i,j]="***"
            }
            else if(data_p[i,j] <= 0.01 && data_p[i,j] > 0.001)
            {
                data_mark[i,j]="**"
            }
            else if(data_p[i,j] <= 0.05 && data_p[i,j] > 0.01)
            {
                data_mark[i,j]="*"
            }
            else
            {
                data_mark[i,j]=""
            }
        }
    }
    write.csv(data_mark, file=paste(route, "data_mark.csv", sep="/"))

    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.pdf", sep="/"))
    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.png", sep="/"))
}
12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485

五、pheatmap绘制热图

correlate_pheatmap("Correlation_result.txt", "Result")
1
  • 结果目录,新增结果图(一个png一个pdf):
  • 打开pdf,如下:
  • 随机模拟的数据,没有显著的不奇怪。如果是做项目一般两组数据的相关分析都可以得到一些显著相关的结果,如下:

两组数据完全相同的话:

六、绘图展示 |r| > 0.8的相关

1 函数

select_pheatmap = function(infile, route)
{
    data=read.table(paste(route, infile, sep="/"), sep="t", header=T)

    data_r=dcast(data, feature_1 ~ feature_2, value.var="r_value")
    data_p=dcast(data, feature_1 ~ feature_2, value.var="raw_p_value")
    rownames(data_r)=data_r[,1]
    data_r=data_r[,-1]
    rownames(data_p)=data_p[,1]
    data_p=data_p[,-1]
    
    # 记下0.8 > r > -0.8 的行
    del_row = c()
    for(i in 1:length(data_p[, 1]))
    {
        if(all(data_r[i, ] < 0.8 & data_r[i, ] > -0.8))
        {
            del_row = c(del_row, i)
        }
    }

    # 记下0.8 > r > -0.8 的列
    del_col = c()
    for(j in 1:length(data_p[1, ]))
    {
        if(all(data_r[, j] < 0.8 & data_r[, j] > -0.8))
        {
            del_col = c(del_col, j)
        }
    }
    
    # null值处理
    if(is.null(del_row) && !(is.null(del_col)))
    {
        data_p = data_p[, -del_col]
        data_r = data_r[, -del_col]
    }else if(is.null(del_col) && !(is.null(del_row)))
    {
        data_p = data_p[-del_row,]
        data_r = data_r[-del_row,]
    }else if(is.null(del_row) && is.null(del_col))
    {
        print("delete none")
    }else if(!(is.null(del_row)) && !(is.null(del_col)))
    {
        data_p = data_p[-del_row, -del_col]
        data_r = data_r[-del_row, -del_col]
    }
    
    # data_p = data_p[-del_row, -del_col]
    # data_r = data_r[-del_row, -del_col]
    write.csv(data_p, file=paste(route, "data_p.csv", sep="/"))
    write.csv(data_r, file=paste(route, "data_r.csv", sep="/"))

    # 用"*"代替<=0.05的p值,用""代替>0.05的相对丰度
    data_mark=data_p
    for(i in 1:length(data_p[,1])){
        for(j in 1:length(data_p[1,])){
            if(data_p[i,j] <= 0.001)
            {
                data_mark[i,j]="***"
            }
            else if(data_p[i,j] <= 0.01 && data_p[i,j] > 0.001)
            {
                data_mark[i,j]="**"
            }
            else if(data_p[i,j] <= 0.05 && data_p[i,j] > 0.01)
            {
                data_mark[i,j]="*"
            }
            else
            {
                data_mark[i,j]=""
            }
        }
    }
    # r值mark 
    data_mark_r = data_r
    for(i in 1:length(data_r[,1])){
        for(j in 1:length(data_r[1,])){
            if(data_r[i,j] >= 0.8 | data_r[i,j] <= -0.8)
            {
                data_mark_r[i,j]=round(data_r[i,j], 2)
            }
            else
            {
                data_mark_r[i,j]=""
            }
        }
    }
    write.csv(data_mark, file=paste(route, "data_mark_p.csv", sep="/"))
    write.csv(data_mark_r, file=paste(route, "data_mark_r.csv", sep="/"))

    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=17, filename=paste(route, "Correlation_result_p.pdf", sep="/"))
    pheatmap(data_r, display_numbers=data_mark_r, cellwidth=20, cellheight=20, fontsize_number=7, filename=paste(route, "Correlation_result_r.pdf", sep="/"))
        
    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=17, filename=paste(route, "Correlation_result_p.png", sep="/"))
    pheatmap(data_r, display_numbers=data_mark_r, cellwidth=20, cellheight=20, fontsize_number=7, filename=paste(route, "Correlation_result_r.png", sep="/"))
}
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899

2 应用

dir.create("metabo_species_0.8r")
select_pheatmap("Correlation_result.txt", "metabo_species_0.8r")
12

七、cytoscape绘图文件

1 准备文件

# 读取关联分析结果
data = read.table("Result/Correlation_result.txt", sep="t", header=T)
# 画图文件
data = data[data$raw_p_value <= 0.05,]
r_label = c()
for(i in 1:nrow(data))
{
    if(data[i, 3] < 0)
    {
        r_label = c(r_label, "neg")
    }
    else
    {
        r_label = c(r_label, "pos")
    }
}

data$r_label = r_label
write.table(data, file="Result/input_pre.txt", sep="t", row.names=F, quote=F)
12345678910111213141516171819
data$r_value = abs(data$r_value)
input_network = data[, c(1,2,3,6)]
write.table(input_network, file="Result/input_network.txt", sep="t", row.names=F, quote=F)
123
g1 = data.frame(feature=input_network[,1], group=rep("tax", nrow(input_network)))
g2 = data.frame(feature=input_network[,2], group=rep("ko", nrow(input_network)))
input_group = unique(rbind(g1, g2))
write.table(input_group, file="Result/input_group.txt", sep="t", row.names=F, quote=F)
1234

2 cytoscape绘图

过程见:Cytoscape绘制相关网络图

3 几种常用的layout

layout -> group attributes layout -> group

layout -> group attributes layout -> degree

layout -> attribute circle layout -> degree

layout -> attribute circle layout -> group

于2020.9.21更新:fdr_p raw_p的对应关系
于2020.9.22更新:输出绘图文件data_p data_r data_mark
于2020.10.14更新:增加了根据R值筛选相关的select_pheatmap

17人点赞

R

更多精彩内容下载简书APP
“小礼物走一走,来简书关注我”

还没有人赞赏,支持一下

小白菜学生信简书记录成长<br>欢迎简信打扰
总资产78 (约4.21元)共写了10.7W字获得992个赞共981个粉丝

点击数:0

使用pheatmap包绘制热图 – 简书

使用pheatmap包绘制热图

152018.10.03 19:58:50字数 33阅读 94,335

加载所需R包

library(pheatmap)
1

设置工作路径

setwd("/Users/Davey/Desktop/VennDiagram/")
# 清除当前环境中的变量
rm(list=ls())
123

构建测试数据集

test = matrix(rnorm(200), 20, 10)
test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
colnames(test) = paste("Test", 1:10, sep = "")
rownames(test) = paste("Gene", 1:20, sep = "")
head(test)
1234567
##          Test1      Test2    Test3      Test4    Test5       Test6
## Gene1 4.064973  0.7535271 3.024070 -2.1294440 4.407945 -0.35677097
## Gene2 2.360043  1.6974946 3.273425 -2.3341406 3.839523  0.16982944
## Gene3 3.253465 -0.9011582 1.716257 -0.2294471 4.636610 -0.24520382
## Gene4 4.070226 -0.6191941 3.734437  1.9348314 4.426825 -0.17730957
## Gene5 3.821414  0.5584876 1.871479 -0.2784607 2.633761  0.01332901
## Gene6 3.012469  0.1738285 3.652423 -2.0083435 4.124951 -0.67899611
##          Test7      Test8    Test9        Test10
## Gene1 3.602764  1.2903843 2.044119  1.826159e+00
## Gene2 3.083160  0.2642755 2.855381  1.988289e-01
## Gene3 3.417809 -0.1362079 3.858884 -8.390304e-01
## Gene4 2.911934  0.4299550 4.128398 -3.011521e+00
## Gene5 2.651758 -1.6884728 3.001079  1.861780e+00
## Gene6 1.934270  0.5811059 2.297763  6.878644e-05
1234567891011121314
# 默认绘图
pheatmap(test)
12
image.png
# scale = "row"参数对行进行归一化
pheatmap(test, scale = "row")
12
image.png
# clustering_method参数设定不同聚类方法,默认为"complete",可以设定为'ward', 'ward.D', 'ward.D2', 'single', 'complete', 'average', 'mcquitty', 'median' or 'centroid'
pheatmap(test,scale = "row", clustering_method = "average")
12
image.png
# clustering_distance_rows = "correlation"参数设定行聚类距离方法为Pearson corralation,默认为欧氏距离"euclidean"
pheatmap(test, scale = "row", clustering_distance_rows = "correlation")
12
image.png
# color参数自定义颜色
pheatmap(test, color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
12
image.png
# cluster_row = FALSE参数设定不对行进行聚类
pheatmap(test, cluster_row = FALSE)
12
image.png
# legend_breaks参数设定图例显示范围,legend_labels参数添加图例标签
pheatmap(test, legend_breaks = c(1:5), legend_labels = c("1.0","2.0","3.0","4.0","5.0"))
12
image.png
# legend = FALSE参数去掉图例
pheatmap(test, legend = FALSE)
12
image.png
# border_color参数设定每个热图格子的边框色
pheatmap(test, border_color = "red")
12
image.png
# border=FALSE参数去掉边框线
pheatmap(test, border=FALSE)
12
image.png
# show_rownames和show_colnames参数设定是否显示行名和列名
pheatmap(test,show_rownames=F,show_colnames=F)
12
image.png
# treeheight_row和treeheight_col参数设定行和列聚类树的高度,默认为50
pheatmap(test, treeheight_row = 30, treeheight_col = 50)
12
image.png
# display_numbers = TRUE参数设定在每个热图格子中显示相应的数值,number_color参数设置数值字体的颜色
pheatmap(test, display_numbers = TRUE,number_color = "blue")
12
image.png
# number_format = "%.1e"参数设定数值的显示格式
pheatmap(test, display_numbers = TRUE, number_format = "%.1e")
12
image.png
# 自定义数值的显示方式
pheatmap(test, display_numbers = matrix(ifelse(test > 5, "*", ""), nrow(test)))
12
image.png
# cellwidth和cellheight参数设定每个热图格子的宽度和高度,main参数添加主标题
pheatmap(test, cellwidth = 15, cellheight = 12, main = "Example heatmap")
12
image.png
# 构建列注释信息
annotation_col = data.frame(
  CellType = factor(rep(c("CT1", "CT2"), 5)), 
  Time = 1:5
)
rownames(annotation_col) = paste("Test", 1:10, sep = "")
head(annotation_col)
1234567
##       CellType Time
## Test1      CT1    1
## Test2      CT2    2
## Test3      CT1    3
## Test4      CT2    4
## Test5      CT1    5
## Test6      CT2    1
1234567
# 构建行注释信息
annotation_row = data.frame(
  GeneClass = factor(rep(c("Path1", "Path2", "Path3"), c(10, 4, 6)))
)
rownames(annotation_row) = paste("Gene", 1:20, sep = "")
head(annotation_row)
123456
##       GeneClass
## Gene1     Path1
## Gene2     Path1
## Gene3     Path1
## Gene4     Path1
## Gene5     Path1
## Gene6     Path1
1234567
# annotation_col参数添加列注释信息
pheatmap(test, annotation_col = annotation_col)
12
image.png
# annotation_legend = FALSE参数去掉注释图例
pheatmap(test, annotation_col = annotation_col, annotation_legend = FALSE)
12
image.png
# annotation_col和annotation_row参数同时添加行和列的注释信息
pheatmap(test, annotation_row = annotation_row, annotation_col = annotation_col)
12
image.png
# 自定注释信息的颜色列表
ann_colors = list(
  Time = c("white", "firebrick"),
  CellType = c(CT1 = "#1B9E77", CT2 = "#D95F02"),
  GeneClass = c(Path1 = "#7570B3", Path2 = "#E7298A", Path3 = "#66A61E")
)
head(ann_colors)
1234567
## $Time
## [1] "white"     "firebrick"
## 
## $CellType
##       CT1       CT2 
## "#1B9E77" "#D95F02" 
## 
## $GeneClass
##     Path1     Path2     Path3 
## "#7570B3" "#E7298A" "#66A61E"
12345678910
# annotation_colors设定注释信息的颜色
pheatmap(test, annotation_col = annotation_col, annotation_colors = ann_colors, main = "Title")
12
image.png
pheatmap(test, annotation_col = annotation_col, annotation_row = annotation_row, 
         annotation_colors = ann_colors)
12
image.png
pheatmap(test, annotation_col = annotation_col, annotation_colors = ann_colors[2]) 
1
image.png
# gaps_row = c(10, 14)参数在第10和14行处添加gap, 要求对行不进行聚类
pheatmap(test, annotation_col = annotation_col, cluster_rows = FALSE, gaps_row = c(10, 14))
12
image.png
# cutree_col = 2参数将列按聚类树的结果分成两部分, 要求对列进行聚类
pheatmap(test, annotation_col = annotation_col, cluster_rows = FALSE, gaps_row = c(10, 14),
         cutree_col = 2)
123
image.png
# 对行和列都不聚类,自定义划分行和列的gap
pheatmap(test, annotation_col = annotation_col, cluster_rows = FALSE, cluster_cols = FALSE, 
         gaps_row = c(6, 10, 14), gaps_col = c(2, 5, 8))
123
image.png
# 自定义行的标签名
labels_row = c("", "", "", "", "", "", "", "", "", "", "", "", "", "", "", 
               "", "", "Il10", "Il15", "Il1b")
# labels_row参数添加行标签
pheatmap(test, annotation_col = annotation_col, labels_row = labels_row)
12345
image.png
# 自定义聚类的距离方法
drows = dist(test, method = "minkowski")
dcols = dist(t(test), method = "minkowski")
# clustering_distance_rows和clustering_distance_cols参数设定行和列的聚类距离方法
pheatmap(test, clustering_distance_rows = drows, clustering_distance_cols = dcols)
12345
image.png
# fontsize参数设定标签字体大小,filename参数设定图片保存名称
pheatmap(test, cellwidth = 15, cellheight = 12, fontsize = 8, filename = "test.pdf")
12

将热图结果按聚类后的顺序输出

aa=pheatmap(test,scale="row")  #热图,归一化,并聚类
1
image.png
# 简要查看热图对象的信息
summary(aa)
12
##          Length Class  Mode   
## tree_row 7      hclust list   
## tree_col 7      hclust list   
## kmeans   1      -none- logical
## gtable   6      gtable list
12345
order_row = aa$tree_row$order  #记录热图的行排序
order_col = aa$tree_col$order    #记录热图的列排序
datat = data.frame(test[order_row,order_col])   # 按照热图的顺序,重新排原始数据
datat = data.frame(rownames(datat),datat,check.names =F)  # 将行名加到表格数据中
colnames(datat)[1] = "geneid" 
write.table(datat,file="reorder.txt",row.names=FALSE,quote = FALSE,sep='t')  #输出结果,按照热图中的顺序
123456
sessionInfo()
1
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.3
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] pheatmap_1.0.10
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.18       digest_0.6.16      rprojroot_1.3-2   
##  [4] grid_3.5.1         gtable_0.2.0       backports_1.1.2   
##  [7] magrittr_1.5       scales_1.0.0       evaluate_0.11     
## [10] stringi_1.2.4      rmarkdown_1.10     RColorBrewer_1.1-2
## [13] tools_3.5.1        stringr_1.3.1      munsell_0.5.0     
## [16] yaml_2.2.0         compiler_3.5.1     colorspace_1.3-2  
## [19] htmltools_0.3.6    knitr_1.20
12345678910111213141516171819202122232425

189人点赞

更多精彩内容下载简书APP
“小礼物走一走,来简书关注我”
共1人赞赏

Davey1220中山大学中山医学院-基础医学博士在读<br>欢迎关注个人公众号<br><br>公众号:bioi…
总资产167 (约8.94元)共写了9.4W字获得1,240个赞共1,620个粉丝

点击数:0

R语言_生成带基因名的火山图_EnhancedVolcano_[空城·崛起]的博客-CSDN博客

R语言_生成带基因名的火山图_EnhancedVolcano

[空城·崛起]
2020-10-10 21:46:03

2139

收藏

3

文章标签:
生物学
r语言

环境:R_x64_4.0.2 & RStudio_1.2.1335
相比网站生成火山图,使用R语言生成火山图可以满足更多的要求,但相关文章不甚清晰,遂记录一下生成 带标签火山图 的过程,留与媛媛查阅。

效果预览

一、EnhancedVolcano安装方法

1.安装

RStudio中输入

install.packages('devtools');devtools::install_github('kevinblighe/EnhancedVolcano');
  • 1

这里代码意义是下载Github上的EnhancedVolcano包,如果选择节点记得选择国内节点

等待安装直到出现如下提示:

即已经完成安装。

2.检验

检验是否成功安装,输入

library(EnhancedVolcano)
  • 1

声明包,若无ERROR字样即成功安装。

二、使用EnhancedVolcano绘制火山图

1.使用library(EnhancedVolcano)声明EnhancedVolcano包。
2.导入数据,以.csv文件为例

data=read.csv(file="test.csv",header=T,row.names=1,sep=',')
  • 1

参数说明:

  • file="" : 引号内为文件名称,文件要放在默认文件夹下,可打开右下Fils并将文件拖入,见下图
  • header=T:列名默认设置。
  • row.names=1:将第一列作为行名,就是我们要在图中展示的标签。
  • sep=',':以,为分隔符,如果txt输入将参数换为" "(空格)。

示例数据:

3.使用EnhancedVolcano()语句生成火山图,举例如下:

EnhancedVolcano(data, lab = rownames(data), x = 'Foldchange', y = 'Pvalue',xlim = c(-17, 13),ylim=c(1,6),pCutoff = 0.001,FCcutoff = 2)
  • 1

参数说明:

  • data:导入的数据文件。
  • lab = rownames(data):设置标签,方便展示(此参数须在2.步中设置row.names才能使用)
  • x = 'Foldchange'y = 'Pvalue':设置横纵坐标轴与data中列的对应
  • xlim = c(-17, 13)ylim=c(1,6):设置xy展示区间
  • pCutoff = 0.001FCcutoff = 2:自定义阈值线。

在右侧即可看到生成的火山图

点击数:0

R做GSEA富集分析

R做GSEA富集分析

22020.06.30 16:30:37字数 577阅读 7,023

首先感谢Y叔的clusterprofiler神包,做富集分析优点是在线爬取数据,结果很可信,但是缺点也是网络问题,网络差点就要等很久,不过GSEA有自带GMT文件,因此下载好离线数据,这些就可以摆脱在线的问题,单机就可以操作GSEA了
GSEA也就是基因集富集分析,不需要分析是上调基因还是下调基因,分析的是所有基因,所以结果应该更可靠,根据官方操作说明, 需要有两组数据,第一组是表达矩阵,第二组是分组,然后用软件操作,但是至始至终出的图奇丑无比,而且还只能是png格式,远达不到300dpi的发表级,虽然目前有很多再次作图方法,但依然比较曲折,Y叔的包解决了这个问题,下面分享一下教程

  • 第一步:数据格式
    表格需要两列, 第一列是基因名,第二列是Foldchange或者Log2FC
    以前我一直不明显第二列是如何来的,明明是矩阵和分组,怎么算Fc,但是也怪自己学艺不精,其实有列矩阵和分组就可以算出来了,懂R的可以用各种差异分析包如“DESq”,“limma”和“edgeR”算出来,不懂R的用Excel也可以,可以看我之前的教程,其实Foldchange就是两组数据的均值之比,也就是先算出各组的平均值,然后拿比较组除以对照组就算出来了 (没事还是要多读书啊)
  • 第二步,如基因名是symbol,需要基因ID转换为entrezid格式
gene<-read.csv("你的文件.csv")  #csv可读性较好,带表头,需要有一列是symbol
library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
gene<-str_trim(d$symbol,"both") #定义gene
#开始ID转换
gene=bitr(gene,fromType="SYMBOL",toType="ENTREZID",OrgDb="org.Hs.eg.db") #会有部分基因数据丢失,或者ENSEMBL
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)
gene_df <- data.frame(logFC=gene$logFC, #可以是foldchange
SYMBOL = gene$symbol) #记住你的基因表头名字
gene_df <- merge(gene_df,gene,by="SYMBOL")
24681012141618202224
  • 第三步,定义基因列表,对logFC进行从高到低排序
geneList<-gene_df $logFC #第二列可以是folodchange,也可以是logFC
names(geneList)=gene_df $ENTREZID #使用转换好的ID
geneList=sort(geneList,decreasing = T) #从高到低排序
246
  • 第四步,下载离线GMT文件
    https://www.gsea-msigdb.org/gsea/downloads.jsp (需要先注册,不过一个邮箱地址即可)
    包含GOKEGGREACTOMEHALLMARKS等等,需要什么下载什么

    图片.png

  • 第五步,读取GMT格式,直接GSEA分析并出图,以Kegg为例

 kegmt<-read.gmt("c2.cp.kegg.v7.1.entrez.gmt") #读gmt文件
 KEGG<-GSEA(geneList,TERM2GENE = kegmt) #GSEA分析
 library(ggplot2)
 dotplot(KEGG) #出点图 
dotplot(KEGG,color="pvalue")  #按p值出点图 
246810
默认p.ajust

p值
dotplot(KEGG,split=".sign")+facet_grid(~.sign) #出点图,并且分面激活和抑制
2
图片.png
 dotplot(KEGG,split=".sign")+facet_wrap(~.sign,scales = "free") #换个显示方式
2
图片.png
 library(enrichplot)
#特定通路作图
 gseaplot2(KEGG,1,color="red",pvalue_table = T) # 按第一个做二维码图,并显示p值
246
图片.png
 gseaplot2(KEGG,1:10,color="red") #按第一到第十个出图,不显示p值
2
图片.png

23人点赞

更多精彩内容下载简书APP
“小礼物走一走,来简书关注我”

还没有人赞赏,支持一下

欧阳松石河子大学医学院第一附属医院泌尿外科主治医师<br>华中科技大学同济医学院附属同济医院泌尿外科…
总资产7 (约0.54元)共写了7433字获得174个赞共186个粉丝

点击数:0

带你深度剖析基因集富集分析-GSEA_通路

带你深度剖析基因集富集分析-GSEA

2020-08-26 11:00
来源:医学方

各位医学方的朋友,大家好,我们你们的老朋友Flyman,时隔许久,今天我们又和大家见面了。今天我们和大家分享一下大家司空见惯的话题基因通路富集分析。谈到富集分析,大家首先会想到Y叔叔写的神包clusterprofiler,几行代码轻松搞定,目前主流的富集分析的算法寻找相关通路本质就是一个超几何分布,那我们今天剖析的GSEA算法则和超几何分布是略有不同。

GSEA全称是 gene-set enrichment analysis,即基因集富集分析,而不是简单地基因富集分析。上次我们和大家已经简述了GSEA的实现过程,那我们这次更深度的解析一下GSEA的算法本质,并更详细的实现相应的操作。

首先我们看一下上面的这张原理图,第一步我们需要准备一个基因表达量矩阵作为GSEA输入,在上图我们看到样本是可以分为A和B两组,类似我们常规的给药组和模型组以及模型组与正常组,接着对所有基因进行排序,这一步也是最关键的一步,对基因排序,一般我们是按照foldchange值从大到小排序,注意这里不是绝对值之后的排序,而是原始值的排序,用来表示基因在两组间表达量的变化趋势。比如AvsB排序之后的基因列表,大于0的基因列表表示在A组高表达,B组低表达,小于0的基因代表在A低表达,B高表达,那对于这样的基因列表,其顶部的基因可以看做是上调的差异基因,其底部是下调的差异基因。

GSEA分析的是一个基因集下的所有基因是否在这个排序列表的顶部或者底部富集,如果在顶部富集,我们可以说,从总体上看,该基因集是上调趋势,反之,如果在底部富集,则是下调趋势。

下面进入正题,基于R语言操作这里使用的是Y叔的R包clusterProfiler:

1

安装并加载包,关于包的安装,已经讲过多次,直接上代码:

2

构建基因排序信息,这里我们纳入采用内置的数据集geneList对象,该对象是包含每个基因的logFC值大小,由于该对象存在于DOSE中,所以大家也要安装好DOSE,另外该数据集如下:

每个数值的名字属性,代表该基因的Entrz ID。数据准备好之后,我们便可以完成GSEA分析了。

3

针对KEGG的GSEA分析,如下,直接借助gseKEGG函数实现。

4

结果解读,生成的结果如下:

我们需要关注两列,enrichmentScore和NES列,一般这两列方向具有同号性,一般大于0可以认为该通路被激活或者高表达,小于0认为被抑制或者低表达。

5

首先结果可视化一下一条通路:

很明显当NES<0,其主要在底部富集,当NES>0,主要在顶部富集,如下:

当然我们还可以进行更好看一些的颜色搭配,采用的是gseaplot2如下:

6

多条通路同时展示,如下:

Ok,今天的推文就到这,我们分享了什么是GSEA以及如何使用clusterProfiler工具进行GSEA分析,希望能对大家有所帮助,最后,欢迎大家留言,有不正确的地方,也请大家留言指正。

END

征 稿 启 事

“医学方”始终致力于服务“医学人”,将最前沿、最有价值的临床、科研原创文章推送给各位临床医师、科研人员。返回搜狐,查看更多

声明:该文观点仅代表作者本人,搜狐号系信息发布平台,搜狐仅提供信息存储空间服务。

Measure
Measure

点击数:0

一文掌握GSEA,超详细教程 – 云+社区 – 腾讯云

一文掌握GSEA,超详细教程

2019-05-152019-05-15 10:47:02阅读 23K0

生信宝典之前总结了一篇关于GSEA富集分析的推文——《GSEA富集分析 – 界面操作》,介绍了GSEA的定义GSEA原理GSEA分析Leading-edge分析等,是全网最流行的原理+操作兼备教程,不太了解的朋友可以点击阅读先理解下概念 (为了完整性,下面也会摘录一部分)。

GSEA案例解析

介绍GSEA分析之前,我们先看一篇Cell文章(https://sci-hub.tw/10.1016/j.cell.2016.11.033)的一个插图 (SCI-HUB客户端(文献神器V4.0)——下载文献如此简单)。

以下是文章原文对图的注解: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运行,但不需联网,启动快。

软件启动界面如下:

二、数据准备

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

Data File

Content

Format

Source

Expression dataset

Contains 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 txt

You create the file. 一般的基因表达矩阵整理下格式就可以。如果是其它类型数据或自己计算rank也可以,后面有更多示例。(如果后缀为txt格式,传统的基因表达矩阵就可以,第一列为基因名字,名字与待分析的功能注释数据集一致,同为GeneSymbol或EntrezID或其它自定义名字,第一行为标题行,含样品信息。gct文件需要符合下面的格式要求。)

Phenotype labels

Contains phenotype labels and associates each sample with a phenotype.

cls

You create the file or have GSEA create it for you. 一般是样品分组信息或样品属性度量值或时间序列信息。

Gene sets

Contains 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 gmt

You 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 annotations

Lists each probe on a DNA chip and its matching HUGO gene symbol. Optional for the gene set enrichment analysis.

Chip

You 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文件)

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

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

2. 样品分组信息

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

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

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

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

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. 数据导入

按照上图步骤依次点击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

  • 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

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

  • 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

  • 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的话则表示运行成功,此过程大概需要十分钟左右,视数据大小而定。

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

四、结果解读

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

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

1. Enrichment in phenotype

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

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

  • 经过条件筛选后还剩下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,在网页查看富集结果,如下:

  • GS:基因集的名字,GO条目的名字
  • SIZE:GO条目中包含表达数据集文中的基因数目(经过条件筛选后的值);
  • ES:富集评分;
  • NES:校正后的归一化的ES值。
    由于不同用户输入的基因数据库文件中的基因集数目可能不同,富集评分的标准化考虑了基因集个数和大小。
    其绝对值大于1为一条富集标准。
    计算公式如下:
  • 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%,将前两项统计数据结合在一起计算出的富集信号强度,计算公式如下:

其中n是列表中的基因数目,nh是基因集中的基因数目

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

2)Details for gene set

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

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

  • 第一部分是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的通路是显著富集的。

最后还有一个该GO基因集下每个基因的详细统计信息表,RANK IN GENE LIST表示在排序好的基因集中所处的位置;RANK METRIC SCORE是基因排序评分,我们这里是Signal2noiseRUNNING ES是分析过程中动态的ES值;CORE ENRICHMENT是对ES值有主要贡献的基因,即Leading edge subset,在表中以绿色标记。

2. Dataset details

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

3. Gene set details

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

此结果显示我们从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处更改。

4. Gene markers for the NGT versus DMT comparison

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

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

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

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

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

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

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

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

5. Global statistics and plots

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

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

Cytoscape富集网路可视化

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

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

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

Graphpad作图比较多个ES

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

效果图如下:

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

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

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

打开graphpad之后,我们在XY类图下选择Enter and plot a single Y value for each point,将两部分数据分开粘贴到软件不同数据表格中(如下图左侧所示),下图中间展示两个图选择的不同绘图方式,调整参数后最终得到右侧的结果。

在左侧目录树处点击layout创建一个图形排版界面,将Graphs下的图形复制粘贴到layout1下,拖拉移动位置很快就能将两部分图对齐。

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

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

本文分享自微信公众号 – 生信宝典(Bio_data),作者:生信宝典

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-04-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

Measure
Measure

点击数:0

GSEA-基因富集分析

GSEA-基因富集分析

Posted on

2018-05-08

Edited on
2019-11-09

In

Bioinformatics-Notes

,

Basic

Symbols count in article:
4.8k

Reading time ≈
4 mins.

说到富集,富集是将基因根据一些先验的知识(也就是常见的注释)进行分类的过程。我们一般会想到最常见的是GO/KEGG富集,其思路是先筛选差异基因,然后确定这些差异基因的GO/KEGG注释,然后通过超几何分布计算出哪些通路富集到了,通常会选择一个阈值来卡一下,比如p值和FDR等。因此这会涉及到人为的阈值选择,具有一定的主观性,而且只能用于差异较大的基因,所以结果可能有一定的局限性。 根据上述情况,有了GSEA(Gene Set Enrichment Analysis),其思路是发表于2005年的Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,主要是要有两个概念:预先定义的基因集S(基于先验知识的基因注释信息)和待测基因集L(一般是表达矩阵);然后GSEA目的就是为了判断S基因集中的基因是随机分布于L(排序后的数据集),还是聚集分布在L的顶部或者底部(这也就是富集)。如果待测基因集中的某些基因显著富集在L的顶部或者底部,这说明这些基因的表达(因为其是根据表达谱数据)对你定义的分组(预先分组)的差异有显著影响(一致性),从而找到我们关注的基因集;在富集分析的理论中,GSEA可以认为是第二代,即Functional Class Scoring (FCS) Approaches

GSEA-homegraphic.gif

GSEA的使用


这里不详细说算法,具体可看GSEA的文章,因为我也是一知半解。。。

下载地址http://software.broadinstitute.org/gsea/downloads.jsp,PS.会验证下你的邮箱,先注册下

第一次使用的话,而且数据不大的话,建议使用javaGSEA Desktop Application,Java的图形化界面,使用较为友好点,根据你的数据大小,选择不同内存的版本[1-8G不等],PS.记得看看系统的Java版本,2G内存开始的GSEA版本需要的是64位的Java 1.8版

打开后的界面如下:

gsea_software

数据准备

因为GSEA分析一般只作用于人物种的,所以我准备以TCGA的BRCA的mRNA数据作为测试数据,正好也试下UCSC xena 浏览器才是最简单的TCGA数据下载途径这个方法来下载TCGA数据(数据还蛮新的,2017年的)

数据下载地址:https://xenabrowser.net/datapages/?dataset=TCGA-BRCA%2FXena_Matrices%2FTCGA-BRCA.htseq_fpkm.tsv&host=https%3A%2F%2Fgdc.xenahubs.net

其还提供了ID/Gene Mapping的文件(整理好的),正好可以拿来用,因为虽然GSEA有EnsemblID转化的chip文件,但是感觉有些数据有点问题(可能是由于Ensembl的版本一直在更新的缘故),HUGO gene symbol最好

然后用R处理下,将癌组织和对应的癌旁组织的数据分别提取出来分别作为两组的表达矩阵(gct文件)以及或者分组文件(cls文件)

library(dplyr)
library(stringr)

data <- read.table(file = "TCGA-BRCA.htseq_fpkm.tsv", sep = "t", header = T, stringsAsFactors = F, check.names = F)

data_df <- data[, -1]

normal_df <- data_df[, str_sub(names(data_df), 14, 15) %in% 11:19]
length(names(normal_df))

tumor_df <- data_df[, grepl(paste(str_sub(names(normal_df), 1, 12), collapse = "|"), names(data_df)) & !names(data_df) %in% names(normal_df)]
length(names(tumor_df))

idmapping <- read.table(file = "gencode.v22.annotation.gene.probeMap", sep = "t", header = T, stringsAsFactors = F)
geneid <- data.frame(id = data$Ensembl_ID, stringsAsFactors = F)
geneid2symbol <- left_join(geneid, idmapping, by = "id")

all_df <- cbind(Name = geneid2symbol$gene, DESCRIPTION = "na", tumor_df, normal_df)

group <- c(rep("Tumor", 118), rep("Normal", 113))
group <- paste(group, collapse = " ")
group <- c(paste(c(118+113, 2, 1), collapse = " "), "# Tumor Normal", group)

write.table(file = "BRCA_fpkm.txt", all_df, sep = "t", col.names = T, row.names = F, quote = F)
write.table(file = "group.cls", group, col.names = F, row.names = F, quote = F)

从上述代码,我获得118个癌组织样本和对应的113个癌旁样本的表达谱数据,并且将Ensembl ID均转化为了Gene symbol(避免之后用GSEA时,再用chip做ID转化);然后可以直接将txt文件作为输入,也可以将BRCA_fpkm.txt文件做一些处理变成GSEA标准的gct文件,如下图所示:分别加上前2行内容,第一行照抄就行,第二行则是geneid数目和样本数目

BRCA_gct_file

接着是Phenotype labels文件(上述代码直接出了),即cls文件,格式如下图所示:第一行231代表样本数目,2代表分2组,空格间隔,1照抄;第二行井号注释说明分组信息;第三行为每个样本对应的组名,空格分隔

BRCA_cls_file

上述文件的详细格式可参照网站:http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

如果网络不佳的话,接下来最好将Gene sets file(也就是GSEA软件上需要输入的Gene sets database),作者将gene sets都储存在Signature Database (MSigDb)中,去官网下载即可http://software.broadinstitute.org/gsea/downloads.jsp,比如下载个c2.cp.kegg.v6.1.symbols.gmt文件作为Gene sets database

如果数据是芯片数据或者需要GSEA的chip文件做ID转化的话,则也可以先将chip文件下载下来,FTP地址:ftp://ftp.broadinstitute.org/pub/gsea/annotations

软件使用

因为是windows桌面式软件,使用就比较简单了。首先将BRCA_fpkm.txtgroup.cls文件导入,如果网速不好的话(比如我自己),那也把上述下载下来的c2.cp.kegg.v6.1.symbols.gmt文件也导入

GSEA_inputfile

接下来点击RUN GSEA,就是几个指定参数的选择了,如下图所示:

GSEA_paramter

  1. Expression dataset:输入的表达矩阵
  2. Gene sets database:gene sets文件,这里是c2.cp.kegg.v6.1.symbols.gmt文件;PS.这个文件也可以自己创建哦
  3. Number of permutations:置换检验的次数
  4. Phenotype labels:选择比较组,如果你输入的文件就只有2个组别的话,这个就很方便选一个就行了;如果你输入的有三个组别及以上的话,则这里就要跟你的需要选择两个组别的比较组,而且GSEA也会根据你的组别信息去表达矩阵中提取相对应的数据
  5. Collapse dataset to gene symbols: 如果你已经ID转化为HUGO gene symbol,那么这里选FALSE,否则选择TRUE
  6. Permutation type:选择置换的类型,是random phenotype还是random gene sets,一般每组样本数目大于7个时,建议选择phenotype,否则选择gene sets(这句话一直没在官网上找到。。。似乎在文章里?)
  7. Chip platform:早些时候主要都是芯片数据,那时必须要将芯片id转化为HUGO gene symbol,所以这个参数一直叫做chip(我猜的。。。),其实就是对ID进行注释,即ID转化,选择你ID对应的chip文件即可,如果已自行转化了ID的话,则这里空着就行(记得Collapse dataset to gene symbols选择否)

除了Required field参数外,下面还有Basic fields和Advanced fields,具体参见官网吧(注:或者鼠标悬浮在对应参数名称上,有简单的参数介绍哦)

最后点击RUN,等待左下角的Running变成Success,然后点击Success即可查看完整的结果,也可以点击Show results folder,GSEA将所有结果都放在一个文件夹中了!!!

分析结果

来看下文章里最常见的GSEA的结果图片,如下图所示:

GSEA_result

从图上,我们一般关注ES值,峰出现在前端还是后端(ES值大于0在前端,小于0在后端)以及Leading-edge subset(即对富集贡献最大的部分,领头亚集);在ES图中出现领头亚集的形状,表明这个功能基因集在某处理条件下具有更显著的生物学意义;对于分析结果中,我们一般认为|NES|>1,NOM p-val<0.05,FDR q-val<0.25的通路下的基因集合是有意义的

除了上述的结果外,GSEA还提供了Running the Leading Edge Analysis等操作,也可以看看

GSEA的结果解读我也不是太熟悉,还是得多看看文献中的解释说明啦

多于多个样本的批处理,GSEA也有服务器版本,通过命令行即可操作,适合批处理操作;其还提供了R脚本可供使用(但官网上说似乎并一定可行,需要自己调整?),反正我也正准备都试试看。。。

参考资料:

功能数据库专题-GSEA
GSEA富集分析 – 界面操作
http://software.broadinstitute.org/gsea/doc/desktop_tutorial.jsp
http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html
GSEA使用(初级)

本文出自于http://www.bioinfo-scrounger.com转载请注明出处

Measure
Measure

点击数:0