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

上篇文章小编为各位小伙伴介绍了转录组分析的第八步——富集分析。通过富集分析可了解各个基因行使的功能。那么这次小编就以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包名称”)进行安装。至此，转录组的分析我们就介绍完了，希望对小伙伴有过帮助。我们下期再见吧~