Posted on
############################################################
# Producer=Jin Wang (https://www.jingege.wang)
############################################################
# 安装STRINGdb 软件包
# source("https://bioconductor.org/biocLite.R")
# biocLite("STRINGdb")
# install.packages("rlist")
############################################################
library(STRINGdb)
library("rlist")
# 设置程序参数
work_dir <- "/Users/zhangqiuxue/Documents/Train/TCGA/lab/PPI" 
deg_file <- "/Users/zhangqiuxue/Documents/Train/TCGA/lab/DEG/Gene/DE_genes.txt"

setwd(work_dir)


# 获取物种的分类编号
# get_STRING_species(version="10", species_name=NULL) 
# 9606 代表人类
string_db <- STRINGdb$new(version="10", species=9606,
                           score_threshold=700, input_directory= work_dir)

# 读取差异表达的文件,获得差异表达基因列表
degs = read.table(deg_file,header=T,comment.char = "",check.names=F)
degs$gene <- rownames(degs)
head(degs)

# 查看有多少差异表达的基因需要分析 
cat("Total deg genes:", dim(degs)[1])

# 将基因的ID map 到string 数据库中, 不一定每个基因都能map上
deg_mapped <- string_db$map( degs, "gene", removeUnmappedRows = TRUE )

# 查看有多少ID map 上了 
cat("Total String id mapped :", dim(deg_mapped)[1])


# 设置绘图相关的参数
#options(network_file=list(fig=function()
#  par(mar=c(2.1, 0.1, 4.1, 2.1))))

# 筛选出一部分结果,进行绘图
hits <- deg_mapped$STRING_id[1:200]

options(par(mar=c(2.1, 0.1, 4.1, 2.1)))
# 绘图
network_file = "network_plot.png"
png(file=network_file, height = 900, width = 1600)
string_db$plot_network( hits,required_score = 400)
dev.off()


# 将所有的结果输出到文件,后面采用cytoscape 进行网络分析
info <- string_db$get_interactions(deg_mapped$STRING_id)
write.table(info, file = "STRING_info.txt",sep="\t", row.names =F, quote = F)

# 采用igraph 进行聚类分析
clustersList <- string_db$get_clusters(deg_mapped$STRING_id[1:1000])
list.save(clustersList,"clustersList.json")
# 设置绘图参数
#options(SweaveHooks=list(fig=function()
#  par(mar=c(2.1, 0.1, 4.1, 2.1))))

# 绘制前4个聚类图
cluster_plot_file = "cluster_plot.png"
png(file=cluster_plot_file, height = 900, width = 1600)
par(mfrow=c(2,2))
for(i in seq(1:4)){
  string_db$plot_network(clustersList[[i]])
  
}
dev.off()



发表评论

邮箱地址不会被公开。 必填项已用*标注