生信实操丨一个生信菜鸟的上道经验分享-转录组测序(富集分析绘图篇)

上篇文章小编为各位小伙伴介绍了转录组分析的第八步——富集分析。通过富集分析可了解各个基因行使的功能。那么这次小编就以2020年8月发表在Viruses上的文章” Transcriptome Analysis of Rice Reveals the lncRNA–mRNA Regulatory Network in Response to Rice Black-Streaked Dwarf Virus Infection ”为例介绍该文章中有关富集分析的图片的绘制。该文章主要研究水稻感染RBSDV后mRNA和lncRNA的相互作用。文章中对感染前后的水稻样本进行RNA测序分析,其中富集分析的可视化结果包括柱状图,气泡图和有向无环图。

1.柱状图

柱状图用于展示注释到各个GO terms和KEGG通路差异表达基因的分布。横坐标为GO terms,纵坐标为富集到各个GO terms的差异表达基因数量。

绘制柱状图的数据文件是我们上期富集分析的结果。

1. go=read.table('./Desktop/go.txt',header=T,sep='\t')

2. head(go,1)

3.            ONTOLOGY         ID                       Description GeneRatio   BgRatio       pvalue     p.adjust       qvalue

4. GO:0030198       BP GO:0030198 extracellular matrix organization   69/1215 368/18670 9.924215e-16 3.207017e-12 2.543615e-12

5.                                                                                                                                                                                                                                                                                                                                                                                                                                       geneID

6. GO:0030198 MMP1/MMP12/ITGB7/MMP9/CTSV/MMP7/PRDX4/FERMT1/KLK7/PRSS2/COL9A3/MMP15/EGFL6…

7.            Count

8. GO:0030198    69

我们挑选分别注释到BP,CC和MF分类下富集最显著的前10位GO terms绘制柱状图,使用dplyr包对输入文件的qvalue值进行排序并取出前10位GO terms,最后使用rbind函数进行合并。

1. BP=(go %>% arrange(qvalue) %>% filter(ONTOLOGY=='BP'))[1:10,]

2. CC=(go %>% arrange(qvalue) %>% filter(ONTOLOGY=='CC'))[1:10,]

3. MF=(go %>% arrange(qvalue) %>% filter(ONTOLOGY=='MF'))[1:10,]

4. GO=rbind(BP,CC,MF)

5. dim(GO)

6. [1] 30 10

使用ggplot2绘制柱状图。

1. p1=ggplot(GO,aes(x=Description,y=Count,fill=ONTOLOGY))+geom_bar(stat='identity',width=0.8)

2. p1

我们可以看到绘制的柱状图横坐标的GO terms的名称都重叠到一起且名称过长,同一GO terms没有放到一起,图片美观性不强,横坐标、纵坐标和图例标签都需要修改。因此,需要对柱状图进一步优化。

1. f1=function(x){

2.   if(nchar(x)>30){

3.     x=substr(x,1,30)

4.     x=paste(x,'...',sep='')

5.     return(x)

6.   }else{

7.     return(x)

8.   }

9. }

10. GO$label=sapply(GO$Description,f1)

11. d=factor(as.integer(rownames(GO)),labels=GO$label)

12. p1=ggplot(GO,aes(x=d,y=Count,fill=ONTOLOGY))+geom_bar(stat='identity',width=0.8)+

13.   theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+

14.   theme(axis.text.x = element_text(color='gray50',size=6,angle=70,vjust=1,hjust=1))+

15.   ylab('Number of genes')+xlab('GO terms')+

16.   scale_fill_discrete(name='Ontology',breaks=c('BP','CC','MF'),labels=c('biological process','cellular component','molecular function'))

17. p1

还可以互换横坐标和纵坐标,将柱状图变为横向。

1. p2=ggplot(GO,aes(x=Count,y=d,fill=ONTOLOGY))+geom_bar(stat='identity',width=0.8)+

2.   theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+

3.   theme(axis.text.x = element_text(color='gray50',size=6,vjust=1,hjust=1))+

4.   xlab('Number of genes')+ylab('GO terms')+

5.   scale_fill_discrete(name='Ontology',breaks=c('BP','CC','MF'),labels=c('biological process','cellular component','molecular function'))

6. p2

2.气泡图

气泡图可展示富集到各个GO term或KEGG通路的基因数量及其显著程度。横坐标为RichFactor,纵坐标为GO term或KEGG通路,点的大小为富集到该功能下的基因数量而点的颜色为富集的显著程度。

输入文件同样为富集分析结果,本次以kegg富集分析结果为例。

1. kegg=read.table('./Desktop/kegg.txt',header=T,sep='\t')

2. head(kegg,1)

3.                ID                    Description GeneRatio  BgRatio

4. hsa05165 hsa05165 Human papillomavirus infection  303/5917 331/8081

5.                pvalue     p.adjust       qvalue

6. hsa05165 1.334503e-17 2.558173e-15 8.673174e-16

7. geneID

8. hsa05165 890/3695/1978/9636/898/9134/4599/6772/8326/1299/28227/1869/22798/50617/8638…

9.          Count

10. hsa05165   303

我们挑选富集程度最显著的前10位pathway绘制气泡图,首先计算富集因子并使用ggplot2绘图,我们参考柱状图的绘制代码对气泡图进行优化。

1. kegg$RichFactor=kegg$Count/as.integer(strsplit(kegg$GeneRatio,'/')[[1]][2])

2. kegg_top10=(kegg %>% arrange(qvalue))[1:10,]

3. p3=ggplot(kegg_top10,aes(x=RichFactor,y=Description))+

4.   geom_point(aes(size=Count,color=qvalue))+

5.   theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+

6.   theme(axis.text.y = element_text(size=9))+

7.   scale_color_gradient(low='green',high='red')+

8.   labs(color='Qvalue',size='Gene_number')+

9.   labs(x='RichFactor',y='Pathway')+

10.   theme(text=element_text(size=9))

11. p3

图中Qvalue的颜色渐变不是十分明显,我们可以通过设置变量scale_color_gradientn来更改。

1. p3=ggplot(kegg_top10,aes(x=RichFactor,y=Description))+

2.   geom_point(aes(size=Count,color=qvalue))+

3.   theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+

4.   theme(axis.text.y = element_text(size=9))+

5.   scale_color_gradientn(colours=rainbow(10))+

6.   labs(color='Qvalue',size='Gene_number')+

7.   labs(x='RichFactor',y='Pathway')+

8.   theme(text=element_text(size=9))

9. p3

这样图片是不是更好看些呢~

3.有向无环图

有向无环图展示GO富集term的关系,由上至下代表GO terms的包含关系,各个方框的颜色代表该term的显著程度。使用R包topGO绘制有向无环图。

输入文件包含两列,第一列为差异基因id,第二列为该基因注释到的GO terms的GO id并用逗号分隔。该文件可用clusterProfile的GO注释结果整理。

1. if(!requireNamespace('BiocManager',quietly=TRUE)){

2.     install.packages('BiocManager')

3. }

4. BiocManager::install('topGO')

5. library(topGO)

6. geneID2GO=readMappings(file='./Desktop/gene.map')

7. str(head(geneID2GO))

8. List of 6

9.  $ MMP1 : chr [1:10] 'GO:0030198' 'GO:0043062' 'GO:0050900' 'GO:0032963' ...

10.  $ MMP12: chr [1:15] 'GO:0030198' 'GO:0043062' 'GO:0032963' 'GO:0022617' ...

11.  $ ITGB7: chr [1:7] 'GO:0030198' 'GO:0043062' 'GO:0050900' 'GO:0031589' ...

12.  $ MMP9 : chr [1:19] 'GO:0030198' 'GO:0043062' 'GO:0032963' 'GO:0022617' ...

13.  $ CTSV : chr [1:22] 'GO:0030198' 'GO:0043062' 'GO:0048608' 'GO:0061458' ...

14.  $ MMP7 : chr [1:16] 'GO:0030198' 'GO:0043062' 'GO:0032963' 'GO:0022617' ...

15. GO2geneID=inverseList(geneID2GO)

16. str(head(GO2geneID))

17. List of 6

18.  $ GO:0000070: chr [1:31] 'CDCA8' 'CDC20' 'KIF23' 'CENPE' ...

19.  $ GO:0000075: chr [1:12] 'CDC20' 'NDC80' 'MAD2L1' 'CDT1' ...

20.  $ GO:0000076: chr [1:2] 'CDT1' 'CDC6'

21.  $ GO:0000079: chr [1:3] 'CCNB1' 'CDC6' 'PLK1'

22.  $ GO:0000082: chr [1:4] 'CDT1' 'CCNB1' 'KIF14' 'CDC6'

23.  $ GO:0000083: chr [1:2] 'CDT1' 'CDC6'

24. geneNames=names(geneID2GO)

25. head(geneNames)

26. [1] 'MMP1'  'MMP12' 'ITGB7' 'MMP9'  'CTSV'  'MMP7'

27. myIntersetingGenes=sample(geneNames,length(geneNames)/10)

28. head(myIntersetingGenes)

29. [1] 'ZWINT'   'BGN'     'COL17A1' 'MMP7'    'HPN'     'TRIP13'

30. geneList=factor(as.integer(geneNames %in% myIntersetingGenes))

31. names(geneList)=geneNames

32. str(geneList)

33. Factor w/ 2 levels '0','1': 1 1 1 1 1 2 1 1 1 2 ...

34.  - attr(*, 'names')= chr [1:100] 'MMP1' 'MMP12' 'ITGB7' 'MMP9' ...

35. head(geneList)

36. MMP1 MMP12 ITGB7  MMP9  CTSV  MMP7

37.     0     0     0     0     0     1

38. Levels: 0 1

39.

40. ##构建MF topGOdata

41. MF_GOdata=new('topGOdata',ontology='MF',

42.               allGenes=geneList,

43.               annot=annFUN.gene2GO,

44.               gene2GO=geneID2GO)

45. Building most specific GOs .....

46.     ( 25 GO terms found. )

47.

48. Build GO DAG topology ..........

49.     ( 58 GO terms and 70 relations. )

50.

51. Annotating nodes ...............

52.     ( 74 genes annotated to the GO terms. )

53. ##使用Kolmogorov-Smirnov算法做检验

54. result_KS.elim=runTest(MF_GOdata,

55.                        algorithm='elim',

56.                        statistic = 'ks')

57.  -- Elim Algorithm --

58.

59.          the algorithm is scoring 58 nontrivial nodes

60.          parameters:

61.              test statistic: ks

62.              cutOff: 0.01

63.              score order: increasing

64.

65.      Level 9:   1 nodes to be scored    (0 eliminated genes)

66.

67.      Level 8:   3 nodes to be scored    (0 eliminated genes)

68.

69.      Level 7:   3 nodes to be scored    (0 eliminated genes)

70.

71.      Level 6:   8 nodes to be scored    (0 eliminated genes)

72.

73.      Level 5:   10 nodes to be scored   (0 eliminated genes)

74.

75.      Level 4:   16 nodes to be scored   (0 eliminated genes)

76.

77.      Level 3:   12 nodes to be scored   (0 eliminated genes)

78.

79.      Level 2:   4 nodes to be scored    (35 eliminated genes)

80.

81.      Level 1:   1 nodes to be scored    (35 eliminated genes)

82. ##绘图

83. showSigOfNodes(MF_GOdata,score(result_KS.elim),

84.                firstSigNodes = 10,

85.                useInfo = 'all')

在构建topGOdata时,修改ontology可绘制CC和BP的有向无环图。

1. ##构建CC topGOdata

2. CC_GOdata=new('topGOdata',ontology='CC',

3.               allGenes=geneList,

4.               annot=annFUN.gene2GO,

5.               gene2GO=geneID2GO)

6. ##使用Kolmogorov-Smirnov算法做检验

7. CC_result_KS.elim=runTest(CC_GOdata,

8.                        algorithm='elim',

9.                        statistic = 'ks')

10. ##绘图

11. showSigOfNodes(CC_GOdata,score(CC_result_KS.elim),

12.                firstSigNodes = 10,

13.                useInfo = 'all')

1. ##构建BP topGOdata

2. BP_GOdata=new('topGOdata',ontology='BP',

3.               allGenes=geneList,

4.               annot=annFUN.gene2GO,

5.               gene2GO=geneID2GO)

6. ##使用Kolmogorov-Smirnov算法做检验

7. BP_result_KS.elim=runTest(BP_GOdata,

8.                           algorithm='elim',

9.                           statistic = 'ks')

10. ##绘图

11. showSigOfNodes(BP_GOdata,score(BP_result_KS.elim),

12.                   firstSigNodes = 10,

13.                   useInfo = 'all')

这样有关GO的三大分类的有向无环图就绘制完成啦,是不是十分简单呢~

4.经验总结

本文介绍了常见的富集分析的可视化图片。绘制图形主要使用R进行绘制,R的下载和安装方法小伙伴可参考https://www.r-project.org/。建议安装R4.0版本以上。绘制图形时需要安装一些所需的R包,可使用install.packages(“R包名称”)进行安装。至此,转录组的分析我们就介绍完了,希望对小伙伴有过帮助。我们下期再见吧~

(0)

相关推荐

  • GO、KEGG富集分析(一)有参情况

    一.GO分析的理论知识 what is Gene Ontology(GO)? 基因'本体论' Gene Ontology中最基本的概念是 term.GO里面的每一个entry都有一个唯一的数字标记,形 ...

  • 关于功能富集分析的基础知识

    富集分析基因富集分析(gene set enrichment analysis)是在一组基因或蛋白中找到一类过表达的基因或蛋白.研究方法可分为三种:Over-Repressentation Analy ...

  • 不要怀疑,你的基因就是没办法富集到统计学显著的通路

    经常收到粉丝提问,明明是按照我课程视频操作,也是按照我的代码在处理他自己的数据,但是做kegg数据库富集的时候,就是返回值为空. 另外,插一个题外话,因为黑粉瞎举报,我们生信技能树已经被取消了半个月的 ...

  • 你只想做ID转换却不知道为什么要转换

    最近咱们<生信技能树>的VIP答疑群,有这样的提问:   我觉得很有代表性,很多人仅仅是学了个皮毛,知道是需要进行ID转换,也能够运行代码.但是却搞不懂,不理解自己为什么进行ID转换,以及 ...

  • GEO联合TCGA数据挖掘文献分享

    今天要介绍的这篇章是我们中国人写的,发表在Med Sci Monit上,这篇文章主要是通过下载GEO和TCGA的数据,通过差异表达分析,GO富集分析.KEGG富集分析,PPI分析,COX回归分析,筛选 ...

  • 手把手教你用R做GSEA分析

    GSEA是非常常见的富集分析方式,以前我们做GSEA需要用依赖java的GSEA软件,那个时候准备分析的文件可能要花上很长时间,报错还不知道如何处理.现在我们来学习一下R语言进行GSEA分析. 加载R ...

  • 转录组学习八(功能富集分析)

    任务 选择p<0.05而且abs(log2FC)大于1的基因为显著差异表达基因集,对这个基因集用R包做KEGG/GO超几何分布检验分析. 把表达矩阵和分组信息分别作出cls和gct文件,导入到G ...

  • 上下调基因各自独立进行GO数据库的3分类富集(求美图代码)

    火山图大家应该是也基本上都没有问题,下面的MA图其实跟火山图非常的类似,两者都是log2FC信息,不同的是火山图展现P值,而MA图展现的是表达量情况! 火山图是为了说明log2FC比较大的一般来说具有 ...

  • 生信实操|一个生信素人的上道经验分享-转录组测序(绘图篇)

    转录组测序技术(RNA-seq)作为目前二代测序领域最普遍的技术手段,自从转录组测序问世以来,已经开发了数百种分析工具.根据转录组分析内容可大致将其分析流程分为比对,转录本组装,差异表达分析和差异基因 ...

  • 生信实操丨高级转录组分析WGCNA应该这么做

    上期的看图说话小编以一篇WGCNA分析的文章为例,教大家如何解读WGCNA的分析结果.今天,小编就来教大家如何用R进行WGCNA分析,快来和我一起学习吧~ 1.R包安装 WGCNA分析使用R包'WGC ...

  • 生信实操丨带你复现单细胞转录组纯分析文章(一)

    生信实操 随着测序技术的进步开发了一种单细胞转录组测序(scRNA-seq)技术,单细胞转录组测序技术可以一次检测成千上万个细胞的转录水平,在单细胞水平上检测和定量基因表达水平变化,从而揭示bulk ...

  • 肿瘤学术盛典丨26位肺癌大咖的诊疗经验分享

    几十年前,肺癌几乎可以和死亡划上等号,而随着科技的发展.研究的进行,医学界对于肺癌机制和诊疗的探索认知也逐步加深. 特别是近年来,肺癌领域风起云涌,无论是外科微创技术和手术技巧,还是内科靶向治疗和免疫 ...

  • 一个简单转录组测序数据发两篇sci(你也可以!)

    你还认为普通转录组测序没有用吗? 最近看到群里有小伙伴在讨论一个数据集 GSE140275 ,我发现它这个简单转录组测序数据发两篇sci,是关于  acute ischemic stroke  这个疾 ...

  • 自然建筑 | 建造一个壁炉,温暖整个冬天(经验分享)

    邀请Amber写这篇分享,有这个想法是因为壁炉文化在国内还是个盲区.我发现很多哪怕国内高级别墅配备的砖砌壁炉,建造仍存在很多缺陷和不合理设计,费柴.取暖效率低下使得大多壁炉成了只在节假日烘托氛围的摆设 ...

  • 【干货】如何运作才能使企业参展效果价值最大化?一个十几年媒体人的经验分享,值得收藏

    现在展会繁多,企业不可能一一出现,如何运作准备?才能使企业参展投资的大额人力物力不打水漂,才不会变成企业负担和鸡肋,本人为<半导体行业采购指南>编辑部责任编辑,从事媒体行业十几年,每年参加 ...

  • protocol丨外泌体免疫胶体金电镜经验分享

    外泌体的电镜是鉴定外泌体形态大小的一个常规手段,但对于研究外泌体表面蛋白或者更确切的去鉴定外泌体则需要通过电镜来观察蛋白质在外泌体中的位置,就需要做外泌体免疫胶体金电镜了.比如我们可以检测一些外泌体特 ...

  • (20)一个WES实战-生信菜鸟团博客2周年精选文章集

    WES就是通常我们说的全外显子测序,已经逐步应用到医疗领域了!本次例子的测试数据就是一个自闭症家系的全外显子测序数据结果,我本来雄心勃勃的想帮助他解决病因的问题,后来发现我还是太嫩了,只是会一点数据分 ...