该文献的题目是:Identifying disease genes and module biomarkers by differential interactions
本文受益于R的开源,受限于本人的能力和水平,也没有在github等平台上找到作者的代码,故通过综合各种信息写了这些代码,至少我自己用R最新版是能跑出来的。
第四步以早期对照为例计算相关性系数矩阵,获得三元组即筛选出需要的达到一定相关性系数的基因
load(file = "step5output.Rdata")
e_e_c <- t(expr_early_contol) #把数据框转置
colnames(e_e_c) <- e_e_c[1,] #把第一行设置为行名
e_e_c <- e_e_c[-1,] #去除第一行
write.table(e_e_c,file = "e_e_c.csv",quote = F,sep = ",") #保存为csv文件,否则算cov时会报错(要求x为数值)
e_e_c_n <- read.csv("e_e_c.csv",header = TRUE,stringsAsFactors = F) #读取csv文件
cor_early_contol <- cor(e_e_c_n,method = "spearman") #计算相关性系数矩阵
cor_early_contol[lower.tri(cor_early_contol)] <- 0 #删除重相互作用,也就是只取半三角矩阵
cor_early_contol[lower.tri(cor_early_contol,diag = TRUE)] <- 0 #删除自相互作用,也就是去除对角线元素
cor_early_contol <- abs(cor_early_contol) #对矩阵内元素取绝对值
selec_cor_early_control <- (cor_early_contol>0.85)*cor_early_contol
#将相关系数>0.75的设置为逻辑值TRUE,矩阵相乘,TRUE的就会获得其原本的相关系数值,FALSE的就会变为0
install.packages("tidyverse")
library(tidyverse)
write.table(selec_cor_early_control,file = "selec_cor_early_control.csv")
library(Matrix) #用于获得三元组
Tlist_early_control <- as.data.frame(summary(Matrix(selec_cor_early_control)))
#把相关矩阵的行列基因名作为一个数据框,供匹配三元组用
genename_early_control <- data.frame(genename1=colnames(selec_cor_early_control),
genename2=rownames(selec_cor_early_control),
number=1:ncol(selec_cor_early_control))
#将Tlist_early_control里的i,j列替换为genename中基因名对应的序号
Tlist_early_control$i <- genename_early_control[match(Tlist_early_control$i,genename_early_control$number),2]
Tlist_early_control$j <- genename_early_control[match(Tlist_early_control$j,genename_early_control$number),2]
early_control_gene_net <- Tlist_early_control[,-3] #删除表示相关系数的第三列
colnames(early_control_gene_net) <- c("protein1","protein2")
第五步处理蛋白质PPIN数据,用于将基因相互作用与蛋白相互作用对应
ppi <- read.table("9606.protein.links.v12.0.txt",header = T,na.strings = c("NA"))
range(ppi$combined_score) #查看ppi第三列的值范围
ppi_cor <- subset(ppi,ppi$combined_score>750) #按照string官网上帮助的说法取结合数>0.75的方法
#ENS表示这是Ensembl ID,P表示这是蛋白质
ppi1 <- data.frame(protein1=str_replace(ppi_cor$protein1,"9606.",""),
protein2=str_replace(ppi_cor$protein2,"9606.",""))
#把结合得分去掉,再把每一个字符前的人类代码“9606.”去掉
#去除蛋白质重相互作用和自相互作用
ppi_unique <- unique(as.data.frame(t(apply(ppi1,1,sort))))
#先把蛋白质每行按列排序(就是把a-b和b-a都调成a-b),然后去除重相互作用(包含自相互作用)
rownames(ppi_unique) <- 1:nrow(ppi_unique)
colnames(ppi_unique) <- c("protein1","protein2")
library(biomaRt)
library(tidyverse)
listEnsembl() #列出可使用的数据库,基因ID转换选genes
listDatasets() #列出数据集
listFilters(ensembl)
listAttributes(ensembl)
#设置数据集和数据库,相当于biomart网站上第一步dataset栏
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
protein_1 <- getBM(attributes = c("ensembl_peptide_id","external_gene_name"),
filters = "ensembl_peptide_id",
values = ppi_unique$protein1,mart=ensembl)
#将ppi_unique里的protein1列的替换为protein_1里ensembl_peptide_id对应的gene name
ppi_unique$protein1 <- protein_1[match(ppi_unique$protein1,protein_1$ensembl_peptide_id),2]
ppi_unique$protein2 <- protein_1[match(ppi_unique$protein2,protein_1$ensembl_peptide_id),2]
#删除一切带有NA的行
library(tidyr)
ppi_unique[ppi_unique==""] <- NA #把空格填上NA
ppi_NA <- ppi_unique[complete.cases(ppi_unique),]
rownames(ppi_NA) <- 1:nrow(ppi_NA)
write.table(ppi_NA,file = "ppi_NA.csv",quote = F,sep = ",")
save(ppi_NA,ppi_unique,file = "step8output_ppi.Rdata")
第六步将蛋白质相互作用与基因相互作用取交集,获得在蛋白质相互作用里存在的筛选出的基因相互作用
load(file = "step8output_ppi.Rdata")
early_control_gene <- paste(early_control_gene_net$protein1,early_control_gene_net$protein2,sep = " ")
PPI_PASTE <- paste(ppi_NA$protein1,ppi_NA$protein2)
PPIN_IN_EARLY_control <- intersect(PPI_PASTE,early_control_gene)
PPIN_IN_early_control <- dplyr::intersect(early_control_gene_net,ppi_NA)
write.table(PPIN_IN_early_control,file = "PPIN_IN_early_control.csv")
save(cor_early_contol,PPIN_IN_early_control,PPIN_IN_EARLY_control,file = "early_control_PPIN.Rdata")