R语言可视化通路富集网络图
2021-12-01 本文已影响0人
单细胞空间交响乐
今天我们来点简单的，实现通路网络图，如下
图片.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 转换功能，这些都值得我们进一步探索。
生活很好，有你更好
