今天我们来点简单的，实现通路网络图，如下

图片.png

这里我们需要调整很多想要的内容，比如颜色，线路等等，我们来学习以下

1、安装以及参数介绍

#安装
BiocManager::install("pathview")
library(pathview)

简单介绍主要参数

pathview(gene.data = NULL,  #是我们上传的数据,可以是单样本或多样本,包含基因ID以及向量,
         cpd.data = NULL,   #化合物的数据,与gene.data类似。ID要与KEGG compound ID相对应
         pathway.id,        #是要指定的KEGG通路ID,通常是5个字符
         species = "hsa",   #物种信息,默认为人类,也可以用"Homo sapiens"或"human",字符型
         kegg.dir = ".",    #输出的结果所保存的地址,默认为当前路径
         cpd.idtype = "kegg",     #化合物数据的ID类型,默认为标准的KEGG基因ID,字符型
         gene.idtype ="entrez",   #基因数据的ID类型,默认为"entrez",字符型
                                  #data(gene.idtype.list); gene.idtype.list,可用来查看ID list
         gene.annotpkg = NULL,    #其他的基因注释包,用于ID转换
         min.nnodes = 3,          #通路最小的节点数,默认为3,整数型
         kegg.native = TRUE,      #当为T时,结果图是KEGG原始图,png格式；当为F时,结果图是由Graphviz引擎绘制,pdf格式
         map.null = TRUE,         #是否在通路图上绘制NULL的基因数据或化合物数据,
                                  #当kegg.native=T时,map.null默认为T
         expand.node = FALSE,     #当kegg.native=F时,是否将多个节点合并为一个节点,合并的节点继承原来节点的所有关系
         split.group = FALSE,     #当kegg.native=F时,是否将合并的节点分裂为多个节点
         map.symbol =TRUE,        #在输出的结果中,是否将map基因IDs转化为symbols,
                                  #此参数只有在kegg.native=F,或者kegg.native=T,same.layer=F时有用
         map.cpdname = TRUE,      #在输出的结果中,是否将map化合物IDs转化为常用名称,只有在kegg.native=F时有用
         node.sum = "sum",        #当多个基因或者化合物对应一个节点时,计算节点数据的方法,
                                  #选项有"sum","mean", "median", "max", "max.abs"和"random"
         discrete=list(gene=FALSE,cpd=FALSE),        #基因或化合物中的数据是连续型变量还是离散型变量
                                                     #dsicrete=list(gene=FALSE, cpd=FALSE),都是连续型变量
         limit = list(gene = 1, cpd = 1),            #基因或化合物的数据转换成颜色的限值
         bins = list(gene = 10, cpd= 10),            #基因或化合物的数据转换成颜色的组数
         both.dirs = list(gene = T, cpd = T),        #基因或化合物的数据是单方向性（0~1）还是双方向性（-1~1）
         trans.fun = list(gene =NULL, cpd = NULL),   #如何处理基因和化合物数据,例如"log"或"abs"
         low = list(gene = "green", cpd = "blue"),   #基因数据和化合物数据图例中最低值的颜色
         mid = list(gene = "gray", cpd = "gray"),    #基因数据和化合物数据图例中中间值的颜色
         high = list(gene = "red", cpd ="yellow"),   #基因数据和化合物数据图例中最高值的颜色
         na.col = "transparent", ...)                #基因数据和化合物数据中NA数据的填充颜色

2、数据

加载基因数据

我们输入的数据包含 gene ID 和 vector(单样本)部分，这里的 gene ID 是一个通用概念，可以是基因、转录本、酶或蛋白质。这里的 vector 可以是样本的表达量、倍数变化, p-value, 组蛋白修饰数据等可测量的属性。下面我们以一个 RNA-seq 差异分析后的数据为例，来学习 pathview 的用法。

#加载基因数据
data <- read.delim('D:/lately/weixin/pathview/data/gene_data.txt',row.names=1,header = T)
head(data)  #这里使用的gene ID是"SYMBOL"类型的,不是标准的KEGG基因ID"entrez"类型,所以在使用时需要转换ID

          log2FoldChange
Plekha6       -0.2108345
Pgk1-ps3      -1.6656907
Catsperg2      0.8888933
P3h3           0.7348740
Gm4859        -1.7972257
Def6           4.9005827

data(gene.idtype.list)
gene.idtype.list  #查看可转化的ID列表

 [1] "SYMBOL"       "GENENAME"     "ENSEMBL"      "ENSEMBLPROT"  "UNIGENE"      "UNIPROT"
 [7] "ACCNUM"       "ENSEMBLTRANS" "REFSEQ"       "ENZYME"       "TAIR"         "PROSITE"
[13] "ORF"` </pre>

获取通路 ID

在 KEGG PATHWAY Database 查询,例如查询小鼠的"Cell Cycle"这条通路:

图片

图片

得到通路 ID 为"04110",物种为"mmu"

画图

我们通过指定 gene.data 和 pathway.id 来观察我们数据里的基因在信号通路“Pathways in cancer”上的表达变化:

pv.out <- pathview(gene.data = gene_data,
                   pathway.id = "04110",
                   species = "mmu",
                   out.suffix ="data",        #输出文件的后缀
                   kegg.native = TRUE,        #输出的结果图为原始的KEGG pathway图,图片格式为".png"
                   gene.idtype = "SYMBOL")    #指明我们的ID类型,以便转换ID

list.files(pattern="mmu04110", full.names=T)  #查看我们输出了哪些结果
[1] "./mmu04110.data.png" "./mmu04110.png"      "./mmu04110.xml"

str(pv.out)   #查看pv.out的内容

List of 2     #由两部分内容组成，基因数据目前如下，化合物数据为空
 $ plot.data.gene:'data.frame': 92 obs. of  10 variables:
  ..$ kegg.names: chr [1:92] "12578" "56371" "17215" "18392" ...
  ..$ labels    : chr [1:92] "Cdkn2a" "Fzr1" "Mcm3" "Orc1" ...
  ..$ all.mapped: chr [1:92] "12578" "56371" "17215,17216,17217,17218,17219,17220" "18392,18393,26428,26429,50793,56452" ...
  ..$ type      : chr [1:92] "gene" "gene" "gene" "gene" ...
  ..$ x         : num [1:92] 532 919 553 494 919 919 188 432 123 77 ...
  ..$ y         : num [1:92] 124 536 556 556 297 519 519 191 704 687 ...
  ..$ width     : num [1:92] 46 46 46 46 46 46 46 46 46 46 ...
  ..$ height    : num [1:92] 17 17 17 17 17 17 17 17 17 17 ...
  ..$ ge1       : num [1:92] -0.644 -0.054 4.753 1.211 -0.464 ...
  ..$ mol.col   : chr [1:92] "#30EF30" "#BEBEBE" "#FF0000" "#FF0000" ...
 $ plot.data.cpd : NULL

以下是输出结果图

图片

相比于原始的 KEGG 图，我们可以使用 graphviz 产生一个新的布局，并且输出 PDF 格式的文件:

pv.out <- pathview(gene.data = gene_data,
                   pathway.id = "04110",
                   species = "mmu",
                   out.suffix ="data",
                   kegg.native = FALSE,
                   gene.idtype = "SYMBOL",
                   sign.pos = "bottomleft" )

以下是输出结果图

图片

3、结果图优化

增加一个图层

如果我们想要运行的更快一点，并且不介意输出图片的大小，我们可以分图层，用 same.layer = F 将节点颜色和标签添加到另一个图层中，并且原来的 KEGG 基因标签会变成官方的 gene symbols:

pv.out <- pathview(gene.data = gene_data,
                   pathway.id = "04110",
                   species = "mmu",
                   out.suffix ="data.2layer",
                   kegg.native = TRUE,
                   gene.idtype = "SYMBOL",
                   same.layer = F)            #将节点颜色和标签添加到另一个图层中

图片

在此基础上，修改 kegg.native = FALSE，我们就可以得到一个主图与图例分成两个页面的 PDF 文件

节点拆分与汇总

在原始的 KEGG 视图中，一个基因节点可能代表具有相似或者冗余功能的基因/蛋白质，我们可以将这种包含多个基因的节点拆分成独立的节点，这样可以更好的从基因层面而不是节点层面来查看数据。同时也可以通过汇总基因数据来可视化节点数据:

#split
pv.out <- pathview(gene.data = gene_data,
                   pathway.id = "04110",
                   species = "mmu",
                   out.suffix ="data.split",
                   kegg.native = F,
                   gene.idtype = "SYMBOL",
                   split.group = T)

#expand
pv.out <- pathview(gene.data = gene_data,
                   pathway.id = "04110",
                   species = "mmu",
                   out.suffix ="data.split.expand",
                   kegg.native = F,
                   gene.idtype = "SYMBOL",
                   split.group = T,
                   expand.node = T)

节点拆分结果如下

图片

汇总结果如下

图片

为了画面有更好的清晰度和可读性，默认不分裂节点，也不单独标记每个成员基因。

4、化合物数据

代谢途径中，除了基因节点还有化合物节点，我们可以尝试利用代谢途径( Propanoate metabolism)整合基因数据和化合物数据。这里的化合物数据包括代谢物、药物，对它们的测量和它们的属性。在这里我们仍然使用之前 RNA-seq 差异分析的数据作为 gene data,然后，我们生成模拟化合物或代谢组数据，并加载适当的化合物 ID 类型以进行演示:

#数据准备
sim.cpd.data=sim.mol.data(mol.type="cpd", nmol=3000)   
#生成模拟化合物数据
head(sim.cpd.data)

     C00232      C01881      C02424      C07418      C13756      C07378
-1.09759772  0.12891537  2.07851027  0.93299245 -0.00786048 -0.09023300

#作图
pv.out <- pathview(gene.data = gene_data,
                   cpd.data = sim.cpd.data,            #模拟的化合物的数据
                   pathway.id = "00640",               #查询丙酸代谢通路的KEGG ID为00640
                   species = "mmu",
                   out.suffix ="gene_cpd",
                   keys.align = y,
                   kegg.native = T,
                   gene.idtype = "SYMBOL")

str(pv.out)        #查看pv.out的内容，仍然由两部分内容组成，但是有了化合物数据

List of 2
 $ plot.data.gene:'data.frame': 24 obs. of  10 variables:
  ..$ kegg.names: chr [1:24] "104776" "100705" "56690" "268860" ...
  ..$ labels    : chr [1:24] "Aldh6a1" "Acacb" "Mlycd" "Abat" ...
  ..$ all.mapped: chr [1:24] "104776" "100705,107476" "56690" "268860" ...
  ..$ type      : chr [1:24] "gene" "gene" "gene" "gene" ...
  ..$ x         : num [1:24] 159 159 159 276 377 ...
  ..$ y         : num [1:24] 325 252 204 378 327 409 494 390 390 229 ...
  ..$ width     : num [1:24] 46 46 46 46 46 46 46 46 46 46 ...
  ..$ height    : num [1:24] 17 17 17 17 17 17 17 17 17 17 ...
  ..$ ge1       : num [1:24] -0.85 -0.249 0.184 2.367 -0.476 ...
  ..$ mol.col   : chr [1:24] "#00FF00" "#8FCE8F" "#BEBEBE" "#FF0000" ...
 $ plot.data.cpd :'data.frame': 40 obs. of  10 variables:
  ..$ kegg.names: chr [1:40] "C00222" "C00804" "C01013" "C00099" ...
  ..$ labels    : chr [1:40] "C00222" "C00804" "C01013" "C00099" ...
  ..$ all.mapped: chr [1:40] "" "C00804" "" "C00099" ...
  ..$ type      : chr [1:40] "compound" "compound" "compound" "compound" ...
  ..$ x         : num [1:40] 225 225 324 325 222 ...
  ..$ y         : num [1:40] 327 451 327 388 228 105 105 105 157 228 ...
  ..$ width     : num [1:40] 8 8 8 8 8 8 8 8 8 8 ...
  ..$ height    : num [1:40] 8 8 8 8 8 8 8 8 8 8 ...
  ..$ mol.data  : num [1:40] NA 0.564 NA 0.871 NA ...
  ..$ mol.col   : chr [1:40] "#FFFFFF" "#DFDF5F" "#FFFFFF" "#FFFF00" ...

#生成Graphviz视图，对于代谢途径，我们需要解析xml文件中的反应条目，并将其转换为基因和化合物节点之间的关系。我们用椭圆表示复合节点
pv.out <- pathview(gene.data = gene_data,
                   cpd.data = sim.cpd.data,
                   pathway.id = "00640",
                   species = "mmu",
                   out.suffix ="gene_cpd",
                   keys.align = y,
                   kegg.native = F,
                   gene.idtype = "SYMBOL",
                   cpd.lab.offset = -0.8)

结果如下

图片

5、多样本多状态的数据整合

pathview 可以集成并将多个样本或状态绘制成一个图，我们可以使用多个重复样本模拟化合物数据:

#载入多样本基因数据
multiple_data <- read.delim('D:/lately/weixin/pathview/data/multiple_data.txt',row.names=1,header = T)
head(multiple_data)

                 CD3         CD5        CDX
Plekha6   -0.2108345 -0.45169573  1.4176037
Catsperg2  0.8888933 -3.93830440  1.6417349
P3h3       0.7348740  0.71500376  0.9836691
Gm4859    -1.7972257 -0.03267321 -1.5157978
Def6       4.9005827  3.41884068 -1.9512619
Ttc7      -0.9966976  0.07041968  0.5702768

#模拟化合物数据
sim.cpd.data2 = matrix(sample(sim.cpd.data, 18000,replace = T), ncol = 6)
rownames(sim.cpd.data2) = names(sim.cpd.data)
colnames(sim.cpd.data2) = paste("exp", 1:6, sep = "")
head(sim.cpd.data2, 3)

            exp1       exp2       exp3       exp4       exp5        exp6
C00232 0.5127650  1.3607983  1.3061348 -0.9911892  0.8583468 -0.01232217
C01881 0.8933992 -2.5994239 -0.6302143 -0.6226316 -0.1347750  1.19753281
C02424 0.1837953 -0.5941693 -0.6696285  1.5044386 -1.1817240 -0.41272992

#画图
pv.out <- pathview(gene.data = multiple_data,
                   cpd.data = sim.cpd.data2[,1:2],
                   pathway.id = "00640",
                   species = "mmu",
                   out.suffix ="gene_3_cpd_2",
                   keys.align = y,
                   kegg.native = T,
                   gene.idtype = "SYMBOL",
                   match.data = F,                      #如果基因数据与化合物数据样本是成对的，可以选择T
                   multi.state = T)                     #多样本数据要设置

结果如下，可以看到基因节点和化合物节点被分成多份，对应不同的样本:

图片

6、处理离散型数据

我们可以根据将化合物数据分为绝对值大于 5 和小于 5 两类，构成一组离散型数据:

sel.cpds <- names(sim.cpd.data)[abs(sim.cpd.data) > 0.5]  #只保留绝对值大于0.5的化合物

pv.out <- pathview(gene.data = data,
                   cpd.data = sel.cpds,
                   pathway.id = "00640",
                   species = "mmu",
                   out.suffix ="gene_sel_cpd",
                   keys.align = y,
                   kegg.native = T,
                   gene.idtype = "SYMBOL",
                   limit = list(gene = 1, cpd = 1),         #设置转换为颜色后的假值分几个区间。
                   bins = list(gene = 10, cpd = 10),        #设置转换为颜色后的颜色分几个区间。
                   na.col = "gray",                         #设置NULL值的节点为灰色
                   discrete = list(gene = F, cpd = T))

结果如下:

图片

7、节点颜色设置

pv.out <- pathview(gene.data = gene_data,
                   cpd.data = sim.cpd.data,
                   pathway.id = "00640",
                   species = "mmu",
                   out.suffix ="gene_cpd",
                   keys.align = y,
                   kegg.native = T,
                   gene.idtype = "SYMBOL",
                   low = list(gene = "#003cff", cpd ="#9900ff"),   #基因数据和化合物数据图例中最低值的颜色
                   mid = list(gene = "gray", cpd = "gray"),        #基因数据和化合物数据图例中中间值的颜色
                   high = list(gene = "#cc3333", cpd ="#ffff00"),  #基因数据和化合物数据图例中最高值的颜色
                   na.col = "transparent")

图片

3小结

Pathview 包中的主函数是 pathview()，有着各种参数，是我们用到最多的函数。在这篇文章中，我们介绍了 pathview()的比较常见的用法,包括包安装，数据准备，以及其他有用的特性。我们也可以使用 pathxiew 的网页版，地址是 https://pathview.uncc.edu/。此外，Pathview 在数据整合方面有很强大的功能，包含 4800 个物种，能处理的数据属性和格式包括 连续/离散数据、矩阵/矢量、单个/多个样本数据，包中还具有强大的 ID 转换功能，这些都值得我们进一步探索。

生活很好，有你更好