分类：实验技术

TRIzol法同时提取RNA、DNA、蛋白质_实验方法_丁⾹通

蛋白质的提取


准备试剂：异丙醇含0.3M盐酸胍的95%乙醇无水乙醇1%SDS。

操作步骤：

1.取沉淀DNA后剩余的上清，用异丙醇沉淀蛋白质。每使用1ml TRIzol加1.5ml异丙醇，室温放置10分钟，2～8℃12000×g离心10分钟弃上清。

2.用含0.3M盐酸胍的95%乙醇洗涤蛋白质沉淀。每使用1mlTRIzol加2ml洗涤液，室温放置20分钟，2～8℃7500×g离心5分钟，弃上清，重复两次。用2ml无水乙醇同样方法再洗一次。

3.真空抽干蛋白质沉淀5～10分钟，用1%SDS溶解蛋白质，反复吸打，50℃温浴使其完全溶解，不溶物2～8℃10000×g离心10分钟除去。分离得到的蛋白质样品可用于Western Blot或-5至-20℃保存备用。

注意事项：

　　1.蛋白质沉淀可保存在含0.3M盐酸胍的95%乙醇或无水乙醇中2～8℃一个月以上或-5至-20℃一年以上。

　　2.用0.1% SDS在2～8℃透析三次，10000×g离心10分钟取上清即可用于Western Blot。

点击数：0

科学合理节约样品：Trizol提RNA后提取蛋白的方法_min

1. 以1mL Trizol提取RNA为例，将水相完全吸掉后，加入300 μL无水乙醇，颠倒混匀涡旋混匀30s，使中间层快状物体完全溶解，置冰上室温5 min，离心4℃ 4000 rpm（2000g） 5min（图2）

图2

2. 取出上清（0.5 mL）置新EP管中，加入3倍体积1.5 mL异丙醇（也有方法学写丙酮、甲醇、氯仿），颠倒混匀30s，置冰上10 min，离心4℃12000 rpm 10 min 。（将有机溶液例如丙酮和乙醇加入蛋白质溶液时，和高浓度盐类似产生沉淀，这是因为它们能够降低蛋白质的溶解度。在温度为10 ℃左右时蛋白质在有机溶剂中容易变性，所以在沉淀时必须特别注意温度0 ℃-4 ℃）离心后就能看见图1的沉淀，我是把上清分成两管加3倍体积有机溶剂沉淀蛋白质的，因为如果蛋白过多，后期可能不易于溶解。

3. 弃上清，加入1ml 0.3M盐酸胍含2.5%甘油的95%乙醇（首先，甘油可以降低体系的极性，体系的极性降低了，有助于蛋白结构稳定；第二，甘油可以降低体系的凝固点，这样利于蛋白在低温下保存，如果不加甘油的话，蛋白反复冻融数次就会失活，甚至变性），将沉淀用枪吹散，室温放置20分钟，4℃12000 rpm 5 min，弃上清，重复两次。加入甘油的试剂尽量1个月用完。

4. 重复步骤3。

5. 加入无水乙醇洗涤沉淀，静置10-20 min后4℃12000 rpm 5 min，弃上清。

6. 室温放置10 min, 使乙醇挥发掉，不要采用真空抽干，避免过度挥干，降低蛋白的溶解度，加入200-300 μL1%SDS置50℃温浴使蛋白完全溶解。浓度大约是1-2 mg。如果较难溶解，可以水浴30分钟就放在超声清洗机里10 min。最后完全不溶解的蛋白可以短暂离心去除。（网上还有网友推荐使用8MUrea+2MThiourea+4%CHAPS溶解沉淀所得蛋白质。大家可以尝试。）

7. 准备上样即可

官方说明书链接：https://www.thermofisher.com/document-connect/document-connect.html?url=https://as//://assets.thermofisher.com/TFS-Assets/LSG/manuals/trizol_reagent.pdf

二、常见问题

1. 得率低：样品裂解或匀浆处理不彻底。一定要充分裂解样品，注意加入Trizol后静置的时间。

2. 最后得到的蛋白质沉淀未完全溶解：蛋白过度变性不易于溶解或蛋白量大。不宜采用强有机溶剂变性，量大可分装后加入有机溶剂。

3. 电泳时条带变形：蛋白质沉淀洗涤不充分，有有机溶剂残留，可增加步骤5静置时间。

三、注意事项

1.以上各试剂的用量都是以1mLTrizol为准，Trizol体积变化，各试剂用量随之等比例变化。

2.蛋白质沉淀可保存在含0.3M盐酸胍的95%乙醇或无水乙醇中2-8℃一个月以上或-5至-20℃一年以上。

点击数：0

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

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

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

参考这篇
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)

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

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

keytypes()

1.3. bitr()函数转换ID

函数全部内容如下

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

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')
转换后的结果

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

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

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

在参数这里与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中的一个
organism：hsa，代表人类，其他的物种名称可以参考kegg的网站
drop：去除空值与否

TP53的entrez id为7157

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

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

改成这样↓

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

转换为uniprot

kegg <- bitr_kegg(7157,
                  fromType = 'kegg',
                  toType = 'uniprot',
                  organism = 'hsa')
出现了三个uniprot

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

P04637

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

K7PPA8

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

Q53GA5

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

2 富集分析

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

GO中的注释信息

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

有很多在线（DAVID、Gather、GOrilla,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)
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)

参数意义：
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 ")

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)
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)

就可以了，继续往下进行

#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))
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)
部分内容
2） 提取出gene_id 和 log2FoldChange列
ensemblid_log2fdc <-dplyr::select(all_deseq, gene_id, log2FoldChange)
3) 可以看到这时gene_id还不是整数形式，所以对这列gsub()取整，命名为symbol列添加到ensemblid_log2fdc中
ensemblid_log2fdc$symbol <- gsub("\.\d*", "", all_deseq[,2])
4）转换id：将ENSEMBLID转为ENTREZID
entrezid <-bitr(ensemblid_log2fdc$symbol, 
              fromType = "ENSEMBL", 
              toType = c("ENTREZID"),
              OrgDb = org.Hs.eg.db,
              drop = TRUE)
5）因为6）中merge需要相同列名，所以将symbol列列名改为ENSEMBL
names(ensemblid_log2fdc) <- c("gene_id", "log2FoldChange", "ENSEMBL")
6）ensemblid_log2fdc 和 entrezid拼接到一起，依靠的是它们具有相同列名的列——ENSEMBL
final <- merge(ensemblid_log2fdc, entrezid)
7) 提取entrezid 和log2fdc列
library(dplyr)    
entrezid_log2fdc <- dplyr::select(final, ENTREZID, log2FoldChange)
8) 用arrange()按log2fdc排序，desc()表示降序
final.sort <- arrange(entrezid_log2fdc, desc(log2FoldChange))
9) 构建一下要GSEA分析的数据
genelist <- final.sort$log2FoldChange
names(genelist) <- final.sort[,1]
10) 进行GSEA分析，需要entrezid，因此我前面才会进行ID转化
gsemf <- gseGO(genelist,
           OrgDb = org.Hs.eg.db,
           keyType = "ENTREZID",
           ont="BP",
           pvalueCutoff = 1)
head(gsemf)
11）画出GSEA图
gseaplot(gsemf, 1)
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官网

点击数：0

用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中的一个
organism：hsa，代表人类，其他的物种名称可以参考kegg的网站
drop：去除空值与否

TP53的entrez id为7157

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 富集分析

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

GO中的注释信息

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

有很多在线（DAVID、Gather、GOrilla,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_id 和 log2FoldChange列
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_log2fdc 和 entrezid拼接到一起，依靠的是它们具有相同列名的列——ENSEMBL
final <- merge(ensemblid_log2fdc, entrezid)
2
7) 提取entrezid 和log2fdc列
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)：设置x、y展示区间
  • pCutoff = 0.001、FCcutoff = 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 （需要先注册，不过一个邮箱地址即可）
    包含GO，KEGG，REACTOME，HALLMARKS等等，需要什么下载什么

    图片.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