搜索
查看: 4751|回复: 5

[software] R的bioconductor里面的stringDB包来做PPI分析

[复制链接]

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-2-28 08:05:19 | 显示全部楼层 |阅读模式
本帖最后由 学习最快乐 于 2017-2-28 08:09 编辑

R的bioconductor里面的stringDB包来做PPI分析
功能:可以找出基因对应的STRING的蛋白的ID,可以用来建立PPI网络及绘制图,可以进行KEGG富集分析和GO富集分析,也可以找到在STRING中与某个蛋白相互作用的蛋白。
主要是下面的网址里面的内容:http://www.bioconductor.org/pack ... inst/doc/STRINGdb.R
在下面的实例中,我们默认人类蛋白相互作用网络。人类的分类表示符为9606,如果需要查询其他物种的分类表示符,请访问(http://www.ncbi.nlm.nih.gov/taxonomy)或者本包中的物中分类标识符表。好啦,接下来程序滚滚而来:
## 整个包不是用roxygen2来写帮助文档的,而且自己把所有函数放在了string_db对象里面,用$符合来调用各个函数,也可以查看函数的帮助文档!
## 首先选定物种及数据库的版本!由于高亮代码没有R语言,所以选择了python,事实是R语言哦
[Python] 纯文本查看 复制代码
string_db <- STRINGdb$new( version="10", species=9606,
score_threshold=0, input_directory="" )
###################################################
### code chunk number 3: help
###################################################
STRINGdb$methods() # To list all the methods available.
STRINGdb$help("get_graph") # To visualize their documentation.
## 列出该包所包含的所有函数,并且可以具体查看某个函数的帮助文档。
###################################################
### 下载数据
###################################################
data(diff_exp_example1)#下载自带数据
head(diff_exp_example1)#查看数据,只显示前6行哦
##一个测试数据,三列,如下:
# pvalue logFC gene
# 0.0001018 3.333461 VSTM2L
# 0.0001392 3.822383 TBC1D2
# 通常就是差异分析的结果
###################################################
### code chunk number 5: map
###################################################
example1_mapped <- string_db$map( diff_exp_example1, "gene", removeUnmappedRows = TRUE )
## 因为我们的差异分析是以基因来标识的,需要map到string数据库的蛋白ID
STRINGdb$help("map")
# 查看帮助文档,明白map函数如何使用,以及该函数返回的是什么!
# 本质上就是根据输入一个data.frame,返回的data.frame多了一列!多了的这一列就是蛋白ID,例如9606.ENSP00000270202这种形式
###################################################
### code chunk number 6: STRINGdb.Rnw:118-121
###################################################
options(SweaveHooks=list(fig=function()
par(mar=c(2.1, 0.1, 4.1, 2.1))))
#par(mar=c(1.1, 0.1, 4.1, 2.1))))
## 设置画图的属性
###################################################
### code chunk number 7: get_hits
###################################################
hits <- example1_mapped$STRING_id[1:200]
# 这里简单的挑选了前面的200个蛋白来进行下一步的分析!
## 请记住,这个例子是在随机挑选,事实上我们应该挑选自定义的差异基因
###################################################
### code chunk number 8: plot_network
###################################################
string_db$plot_network( hits )
## 只有有蛋白ID就可以进行画网络图,ID越多,耗时越长!最后会出来一个彩色蛋白蛋白相互作用网络图,很好看。
## 函数会根据输入的ID列表在string数据库里面找到所有的PPI数据,然后画网络图
## STRINGdb$help("plot_network")#查看plot的帮助文件
###################################################
### code chunk number 9: add_diff_exp_color
###################################################
# filter by p-value and add a color column
# (i.e. green down-regulated gened and red for up-regulated genes)
example1_mapped_pval05 <- string_db$add_diff_exp_color( subset(example1_mapped, pvalue<0.05),
logFcColStr="logFC" )
## 上面简单的网络图一般不满足需求,比如我们需要定位基因的上下调关系,还有联系的紧密与否,可以用红绿色的深浅来刻画。
## 用add_diff_exp_color函数得到的对象还是data.frame,但是增加了一列是color
STRINGdb$help("add_diff_exp_color")
###################################################
### code chunk number 10: post_payload
###################################################
# post payload information to the STRING server
payload_id <- string_db$post_payload( example1_mapped_pval05$STRING_id,
colors=example1_mapped_pval05$color )
## 前面add_diff_exp_color函数为我们的data.frame增加了一列是color,还需要用post_payload函数来把string的蛋白ID跟color对应成功,返回一个payload_id对象给画图函数。
STRINGdb$help("post_payload")
###################################################
### code chunk number 11: plot_halo_network
###################################################
# display a STRING network png with the "halo"
string_db$plot_network( hits, payload_id=payload_id )
## 同样是画网络图,但是增加了一个color的属性。
## 可以看出来,基因太多了,画的图其实很拥挤
###################################################
### code chunk number 13: plot_ppi_enrichment
###################################################
# plot the enrichment for the best 1000 genes 画一个PPI富集的图
string_db$plot_ppi_enrichment( example1_mapped$STRING_id[1:1000], quiet=TRUE )
STRINGdb$help("plot_ppi_enrichment")
你可以看到结果列表中基因的贡献。该图可以帮助你定义阈值,如果实验是成功的,你可以看到基因列表开头基因的富集大于结尾的基因的富集。如果知道的可以在下面留言。谢谢。困惑很久了。
## 这个代码我没有看懂在干吗,上面的只是瞎猜测
###################################################
### code chunk number 14: enrichment
###################################################
enrichmentGO <- string_db$get_enrichment( hits, category = "Process", methodMT = "fdr", iea = TRUE )
enrichmentKEGG <- string_db$get_enrichment( hits, category = "KEGG", methodMT = "fdr", iea = TRUE )
head(enrichmentGO, n=7)
head(enrichmentKEGG, n=7)
### 直接根据 string 数据库的蛋白ID来做富集分析,此函数会自动下载一些数据。默认是以人类的蛋白库作为背景,但是大部分情况下是需要改变的,否则P值就算的不准确啦
#################################################
# code chunk number 15: background (eval = FALSE)
#################################################
# 这里修改背景值,人类本来有两万多个基因,这里变成只有2000个了
backgroundV <- example1_mapped$STRING_id[1:2000] # as an example, we use the first 2000 genes
string_db$set_background(backgroundV)
## string_db 是一个全局变量,之前是直接选择人类的V10.0版本,现在被修改了,只是做一个测试,一定要记得改回去!!!
###################################################
### code chunk number 16: new_background_inst (eval = FALSE)
###################################################
string_db <- STRINGdb$new( score_threshold=0, backgroundV = backgroundV )
###################################################
### code chunk number 17: enrichmentHeatmap (eval = FALSE)
###################################################
eh <- string_db$enrichment_heatmap( list( hits[1:100], hits[101:200]),
list("list1","list2"), title="My Lists" )##画富集热图,比较两个或多个基因列表
## 我们还是把 string_db 修改回来吧!
string_db <- STRINGdb$new( version="10", species=9606,
score_threshold=0, input_directory="" )
###################################################
### code chunk number 18: clustering1
###################################################
# get clusters
clustersList <- string_db$get_clusters(example1_mapped$STRING_id[1:600])
###################################################
### code chunk number 19: STRINGdb.Rnw:254-256
###################################################
options(SweaveHooks=list(fig=function()
par(mar=c(2.1, 0.1, 4.1, 2.1))))
###################################################
### code chunk number 20: clustering2
###################################################
# plot first 4 clusters
par(mfrow=c(2,2))
for(i in seq(1:4)){
string_db$plot_network(clustersList[[i]])
}
## 把4个cluster画在同一个画布上面!
###################################################
### code chunk number 21: proteins
###################################################
string_proteins <- string_db$get_proteins()
## 下面是一下其它小工具,比如找两个蛋白的interaction呀,还有两个蛋白直接相互作用的paper呀,还有找某个蛋白在其它物种的同源蛋白呀!
###################################################
### code chunk number 22: atmtp
###################################################
tp53 = string_db$mp( "tp53" )
atm = string_db$mp( "atm" )##只定位一个或少数蛋白的STRING ID
###################################################
### code chunk number 23: neighbors (eval = FALSE)
###################################################
## string_db$get_neighbors( c(tp53, atm) )
###################################################
### code chunk number 24: interactions
###################################################
string_db$get_interactions( c(tp53, atm) )
#找到两蛋白的相互作用以及score
###################################################
### code chunk number 25: pubmedInteractions (eval = FALSE)
支持他两相互作用的文献
###################################################
## string_db$get_pubmed_interaction( tp53, atm )
###################################################
### code chunk number 26: homologs (eval = FALSE)
得到同源蛋白
###################################################
## # get the reciprocal best hits of the following protein in all the STRING species
## string_db$get_homologs_besthits(tp53, symbets = TRUE)
###################################################
### code chunk number 27: homologs2 (eval = FALSE)
###################################################
## # get the homologs of the following two proteins in the mouse (i.e. species_id=10090)
##限定种类为老鼠
## string_db$get_homologs(c(tp53, atm), target_species_id=10090, bitscore_threshold=60 )
###################################################
### code chunk number 28: benchmark1,
##首先提供一个两列分别为交互蛋白以及一列为信心分数
##在此中KEGG数据库作为金标准
###################################################
data(interactions_example)
interactions_benchmark = string_db$benchmark_ppi(interactions_example, pathwayType = "KEGG",
max_homology_bitscore = 60, precision_window = 400, exclude_pathways = "blacklist")
###################################################
### code chunk number 29: STRINGdb.Rnw:391-393
###################################################
options(SweaveHooks=list(fig=function()
par(mar=c(4.1, 4.1, 4.1, 2.1))))
###################################################
### code chunk number 30: benchmark2
###################################################
plot(interactions_benchmark$precision, ylim=c(0,1), type="l", xlim=c(0,700),
xlab="interactions", ylab="precision")
###################################################
### code chunk number 31: benchmark3
###################################################
interactions_pathway_view = string_db$benchmark_ppi_pathway_view(interactions_benchmark, precision_threshold=0.2, pathwayType = "KEGG")
head(interactions_pathway_view)

##列出各种TP相互作用所属的途径以及统计学分数

回复

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-3-4 19:34:17 | 显示全部楼层
这帖子不错
人生若只如初见!
回复 支持 反对

使用道具 举报

2

主题

17

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-3-6 20:46:45 | 显示全部楼层
越学习  越快乐  不过初学者好痛苦啊
回复 支持 反对

使用道具 举报

13

主题

33

帖子

243

积分

版主

Rank: 7Rank: 7Rank: 7

积分
243
发表于 2017-3-17 12:30:09 | 显示全部楼层
写的好认真,点赞!

不过看到后面就看不懂了!
回复 支持 反对

使用道具 举报

0

主题

5

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-9-26 11:26:16 | 显示全部楼层
楼主,请教一下,我在mapping的时候报错了
example1_mapped <- string_db$map( diff_exp_example1, "gene", removeUnmappedRows = TRUE )
试开URL’http://string.uzh.ch/permanent/s ... n_aliases_tf.tsv.gz
Content type 'application/x-gzip' length 13939907 bytes (13.3 MB)
downloaded 13.3 MB

试开URL’http://string.uzh.ch/permanent/s ... 06__proteins.tsv.gz
Content type 'application/x-gzip' length 1815122 bytes (1.7 MB)
downloaded 1.7 MB

Warning:  we couldn't map to STRING 14% of your identifiers>

不知道什么原因?
还有一个问题,如果没有目标物种,有什么别的办法吗?
回复 支持 反对

使用道具 举报

16

主题

32

帖子

273

积分

中级会员

Rank: 3Rank: 3

积分
273
发表于 2018-5-20 16:10:53 | 显示全部楼层
问问大家,我也是用string_db$map下载数据只能下载两个,其他的下不下来,是网络的问题吗?可是前两个倒是挺快的,由什么方法下载吗?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-26 10:17 , Processed in 0.033104 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.