R语言复现一篇2012年生物网络标志物的文献（4）

该文献的题目是：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")