读取数据

library(pathview)
expr_data <- read.csv("gse16873.d3.txt", sep = "\t", header=T)
head(expr_data)
##      GeneID      DCIS_1      DCIS_2       DCIS_3
## 1     10000 -0.30764480 -0.14722769 -0.023784808
## 2     10001  0.41586805 -0.33477259 -0.513136907
## 3     10002  0.19854925  0.03789588  0.341865341
## 4     10003 -0.23155297 -0.09659311 -0.104727283
## 5 100048912 -0.04490724 -0.05203146  0.036390376
## 6     10004 -0.08756237 -0.05027725  0.001821133
dim(expr_data)
## [1] 11979     4
entrez <- as.numeric(expr_data$GeneID)
expr_data <- expr_data[!is.na(entrez),]
rownames(expr_data) <- expr_data$GeneID
head(expr_data)
##              GeneID      DCIS_1      DCIS_2       DCIS_3
## 10000         10000 -0.30764480 -0.14722769 -0.023784808
## 10001         10001  0.41586805 -0.33477259 -0.513136907
## 10002         10002  0.19854925  0.03789588  0.341865341
## 10003         10003 -0.23155297 -0.09659311 -0.104727283
## 100048912 100048912 -0.04490724 -0.05203146  0.036390376
## 10004         10004 -0.08756237 -0.05027725  0.001821133
# The expression may be log(Intensity)
expr_data$FC1 <- expr_data$DCIS_2 - expr_data$DCIS_1
expr_data$FC2 <- expr_data$DCIS_3 - expr_data$DCIS_1
expr_data$FC3 <- expr_data$DCIS_3 - expr_data$DCIS_2
head(expr_data)
##              GeneID      DCIS_1      DCIS_2       DCIS_3          FC1
## 10000         10000 -0.30764480 -0.14722769 -0.023784808  0.160417109
## 10001         10001  0.41586805 -0.33477259 -0.513136907 -0.750640633
## 10002         10002  0.19854925  0.03789588  0.341865341 -0.160653377
## 10003         10003 -0.23155297 -0.09659311 -0.104727283  0.134959859
## 100048912 100048912 -0.04490724 -0.05203146  0.036390376 -0.007124218
## 10004         10004 -0.08756237 -0.05027725  0.001821133  0.037285120
##                   FC2          FC3
## 10000      0.28386000  0.123442886
## 10001     -0.92900495 -0.178364320
## 10002      0.14331609  0.303969466
## 10003      0.12682569 -0.008134168
## 100048912  0.08129762  0.088421835
## 10004      0.08938350  0.052098384

搜索感兴趣的生物通路

https://www.kegg.jp/pathway/hsa05171 (COVID-19) 通路 如果R联网有问题，可以点击Download链接，预先下载生物通路对应的KGML和Image png (1X)文件

Generate pathway view

pv.out <- pathview(gene.data = expr_data[, 5, drop = FALSE], cpd.data = NULL,
                   pathway.id = "hsa05171",
                   species = "hsa",
                   gene.idtype = "entrez",
                   out.suffix = "gse16873_d5",
                   kegg.native = T)

Result

hsa05171.gse16873_d5
hsa05171.gse16873_d5

Multiple comparisons

pv.out <- pathview(gene.data = expr_data[, c(5, 6, 7), drop = FALSE], 
                   cpd.data = NULL,
                   pathway.id = "hsa05171",
                   species = "hsa",
                   gene.idtype = "entrez",
                   out.suffix = "gse16873_d5_d6", 
                   kegg.native = TRUE)

Result

hsa05171.gse16873_d5_d6
hsa05171.gse16873_d5_d6

Default ender pathway graph as native KEGG graph (.png) now use graphviz layout engine

pv.out <- pathview(gene.data = expr_data[, 5, drop = FALSE], cpd.data = NULL,
                   pathway.id = "hsa04110",
                   species = "hsa",
                   gene.idtype = "entrez",
                   out.suffix = "gse16873_04110_d5",
                   kegg.native = F)
##      [,1] [,2] 
## [1,] "9"  "300"
## [2,] "9"  "306"

Result

hsa04110.gse16873_04110_d5
hsa04110.gse16873_04110_d5