Pathview包教程|3.数据整合

今天我们来跟进一步学习Pthview包的功能，设计的内容主要是Pathview包的数据整合能力。

1. 数据整合功能

Pathview提供为数据整合提供了强大的支持。包括基因表达量、蛋白质表达量、遗传学关联、代谢物、基因组数据以及其它可以映射到基因通路上的数据。值得注意的是，当数据映射到KEGG同源路径时，它可以直接用于宏基因组、微生物组或未知物种数据。使用者提供不同类型的ID类型，都可以准确的映射到KEGG通路上。

Pathview可以处理多种类型的数据，包括连续变量和离散变量，矩阵形式或向量形式，单样本或多样本。此外，Pathview可以在一副图内整合并比较多样本数据或状态。

2. 药物以及基因数据整合

上面的例子我们用到了一个典型的信号通路，现在我们想看看代谢通路。在代谢通路当中，除了基因节点，还多出了化合物的节点。因此，我们可以在上面尝试整合并可视化基因数据和化合物数据。

此处的基因数据是一个广泛的概念，包括基因、转录组、蛋白质、酶及其表达、修饰和任何可测量的属性。同样化合物数据可以指代代谢物、药物、它们的测量值和属性。

下面我们以6个乳腺癌样本的基因数据和一个模拟的化合物数据来做展示。

2.1 sim.mol.data()函数

在开始之前，首先要说一下sim.mol.data()函数，这个函数可以模拟自定义样本数量的基因数据或化合物数据，这些基因或化合物的ID有多种选择，并且可以设置变量是连续性的还是离散性的。

可以用data(gene.idtype.bods)或data(cpd.simtypes) 来查询该模拟函数可用的基因ID以及化合物ID的类型

sim.mol.data(mol.type = c("gene", "gene.ko", "cpd")[1], #选择数据类型
             id.type = NULL,    #基因或化合物的ID
             species="hsa",     #物种
             discrete = FALSE, #变量的类型
             nmol = 1000,     #样本数
             nexp = 1, 
             rand.seed=100)


data(gene.idtype.bods) #可用的基因ID类型
data(cpd.simtypes)     #可用的化合物ID类型

2.2 绘制整合后的通路图

我们这里创建一个3000个样本的cpd数据，选择一个代谢通路，来看看基因数据加上化合物数据的效果，只需要在此前的基础上设置cpd.data参数即可

#创建一个3000个样本的cpd数据
sim.cpd.data=sim.mol.data(mol.type="cpd", nmol=3000)


#绘制整合了cpd数据的通路图
i <- 3
print(demo.paths$sel.paths[i])#选择一个代谢通路
[1] "00640"


pv.out <- pathview(gene.data = gse16873.d[, 1],
                   cpd.data = sim.cpd.data,   #添加化合物数据
                   pathway.id = demo.paths$sel.paths[i],
                   species = "hsa", 
                   out.suffix = "gse16873.cpd",
                   keys.align = "y",
                   kegg.native = T,            #用原始KEGG图片形式
                   key.pos = demo.paths$kpos1[i])

此时我们先看一下pv.out的一个结构，我们会发现多了一个plot.data.cpd的数据。

> str(pv.out)
List of 2
 $ plot.data.gene:'data.frame':  24 obs. of  10 variables:
  ..$ kegg.names: chr [1:24] "4329" "31" "23417" "18" ...
  ..$ labels    : chr [1:24] "ALDH6A1" "ACACA" "MLYCD" "ABAT" ...
  ..$ all.mapped: chr [1:24] "4329" "31,32" "23417" "18" ...
  ..$ 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 ...
  ..$ mol.data  : num [1:24] 0.747 -0.483 -0.251 2.785 0.77 ...
  ..$ mol.col   : chr [1:24] "#EF3030" "#5FDF5F" "#8FCE8F" "#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" ...

结果图如下，可以看到包括基因和化合物对应的节点，都根据自身的值有了不同的颜色：

2.3 绘制Graphviz形式的通路图

我们同样可以产生Graphviz形式的通路图，设置kegg.native = FALSE即可

pv.out <- pathview(gene.data = gse16873.d[, 1], 
                    cpd.data = sim.cpd.data,
                    pathway.id = demo.paths$sel.paths[i],
                    species = "hsa", 
                    out.suffix = "gse16873.cpd",
                    keys.align = "y",
                    kegg.native = F, 
                    key.pos = demo.paths$kpos2[i],
                    sign.pos = demo.paths$spos[i], 
                    cpd.lab.offset = demo.paths$offs[i])

样式如下，Graphviz视图更好地显示了层次结构。对于代谢途径，pathview需要解析xml文件并将其转换为基因和化合物节点之间的关系。这里使用椭圆作为化合物节点。化合物的名称是从CHEMBL数据库检索的标准名称。

3. 多样本多状态的数据整合

在之前的例子当中，我们用到的是单样本数据，它可以是向量，也可以是只有一列的矩阵。Pathview还可以处理多样本数据，只需要我们提供数据矩阵即可。我么先来创建一个6*3000的重复样本的化合物数据，也就是把我们2.2中创建的cpd数据重复抽取六次。

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

想要展示多样本的整合数据，我们在函数中设置multi.state = TRUE即可：

pv.out <- pathview(gene.data = gse16873.d[, 1:3],
                   cpd.data = sim.cpd.data2[, 1:2],
                   pathway.id = demo.paths$sel.paths[i],
                   species = "hsa", 
                   out.suffix = "gse16873.cpd.3-2s",
                   keys.align = "y",
                   kegg.native = T, 
                   match.data = F,  #如果基因数据与化合物数据相匹配，可以选择T
                   multi.state = T,
                   same.layer = T)

结果如下，在这些图中，我们看到基因节点和化合物节点被分割成多个片段，对应于不同的状态或样本，而颜色的变化代表着样本在该基因表达量上的变化。这是一种类似于热图的形式。

3.1 匹配数据

此外由于gene.data和cpd.data的样本大小不同，如果两种数据类型中的样本实际上是成对的，我们可以选择匹配数据，设置match.data = TRUE即可。之后我们只要让我们相同的样本的基因数据和cpd数据是能匹配上的就行，比如在各自矩阵的排列顺序相同。

3.2 多样本分开绘制

如果我们想把多个样本的基因和化合物数据分开绘制，可以设置multi.state = FALSE，其他参数不变即可。比如下面我想分别绘制出三个样本整合后的基因数据和cpd数据。

pv.out <- pathview(gene.data = gse16873.d[, 1:3],
                   cpd.data = sim.cpd.data2[, 1:2], 
                   pathway.id = demo.paths$sel.paths[i],
                   species = "hsa", 
                   out.suffix = "gse16873.cpd.3-2s", 
                   keys.align = "y",
                   kegg.native = T, 
                   match.data = F,
                   multi.state = F, 
                   same.layer = T)

此时根据样本数分别绘制了各自的通路图，映射有整合后的基因数据和cpd数据：

> dir(pattern = "DCIS", full.names = T)
[1] "./hsa00640.gse16873.cpd.3-2s.DCIS_1.png" "./hsa00640.gse16873.cpd.3-2s.DCIS_2.png"
[3] "./hsa00640.gse16873.cpd.3-2s.DCIS_3.png"

4. 放在两个图层上来加速

之前我们讲过，绘制单一图层的KEGG视图需要更久的时间，但是体积小，使用的是原始的KEGG ID。如果我们不介意丢失原始的KEGG基因标签（或EC编号），我们可以选择使用两个图层来绘制通路图，这样可以以加快绘图的过程，设置same.layer = FALSE即可。

pv.out <- pathview(gene.data = gse16873.d[, 1:3],
                    cpd.data = sim.cpd.data2[, 1:2],
                   pathway.id = demo.paths$sel.paths[i],
                   species = "hsa", 
                   out.suffix = "gse16873.cpd.3-2s.2layer",
                    keys.align = "y", 
                   kegg.native = T, 
                   match.data = F, 
                   multi.state = T,
                   same.layer = F)

5. 小结

今天我们学习了Pathview包强大的数据整合能力，在此基础上我们可以将多个样本的基因数据和化合物数据整合到同一张通路图上，也可以用拆分的方法批量绘制多个样本的通路图。