可以用網(wǎng)頁版做，但有上限2000個基因的限制。所以今天開發(fā)一下怎么用R飛一波。
1. 下載STRING數(shù)據(jù)庫中蛋白質相互作用網(wǎng)絡
2. 下載Uniprot ID轉換文件
打開terminal
wget -c ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/
得到的是網(wǎng)頁形式，亂碼了
繼續(xù)剛才的步驟，拷貝鏈接
wget -c ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/by_organism/
重復上述步驟，拷貝鏈接，wget
3. 萬事俱備，我們現(xiàn)在手里有三個文件，我整理一下他們之間的關系
1）input gene：我們關注的基因，symbol格式，可用Y叔的包轉化成Uniprot ID
2）id文件：uniprot和ENSG之間的對應關系
3）PPI互作文件：ENSG格式
接下來就開始愉快的數(shù)據(jù)清洗啦～～
- 準備工作：
library(dynamicTreeCut)
library(openxlsx)
library(stringr)
library(Matrix)
library(WGCNA)
Sys.setenv(LANGUAGE = "en") #顯示英文報錯信息
options(stringsAsFactors = FALSE) #禁止chr轉成factor
setwd("/Users/baiyunfan/desktop")
- 讀取這三個文件
idmapping<-read.table("HUMAN_9606_idmapping_selected.tab",header = F,as.is=T,sep="\t")
ppi <- read.table("9606.protein.actions.v11.0.txt",header=T,sep = "\t")
gene<-read.table("turquoise.txt",sep=",")
idmapping
gene
ppi
- 將我們的輸入基因SYMBOL轉化成UNIPROT ID
library(clusterProfiler)
m<-bitr(gene[,2],fromType = "SYMBOL",toType = "UNIPROT",OrgDb = "org.Hs.eg.db")
colnames(idmapping)[1]<-"UNIPROT"
- 通過idmapping文件，將UNIPROT，SYMBOL，ENSP三種ID聯(lián)系到一起
n<-merge(m,idmapping[,c(1,21)],by="UNIPROT",all.x=T)
n<-n[-which(n[,3]==""),]
- 上圖可看出，第三列有多個ENSP擠在一個格里，按照分號給拆分一下
prots<-str_split(n[,3],"[;]")
names(prots)<-n[,1]
prots_tmp<-unlist(lapply(1:length(prots), function(x){paste(names(prots)[x], prots[[x]],sep=";")}))
prots_mat <- str_split(prots_tmp,"[;]",2,simplify = T)
colnames(prots_mat) <- c("uniprot","ensemblprot")
prots_mat
- ppi文件前面多個9606.，需要清洗掉
ppi$item_id_a <- str_replace(ppi$item_id_a,"9606.","")
ppi$item_id_b <- str_replace(ppi$item_id_b,"9606.","")
- 將PPI中的目標基因留下，其余的刪掉，并刪去重復的
ppi<-ppi[which(ppi[,1] %in% prots_mat[,2] & ppi[,2] %in% prots_mat[,2]),]
ppi$identical<-paste0(ppi[,1],ppi[,2])
ppi<-ppi[!duplicated(ppi$identical),]
ppi$identical<-paste0(ppi[,2],ppi[,1])
ppi<-ppi[!duplicated(ppi$identical),]
ppi1
最后兩行就是我們的目標蛋白互作啦～
