WGCNA衍生：超几何分布识别模块Pivot regulators

realNet_TF<-read.table(file.choose(),header=T,sep="\t")
head(realNet_TF);dim(realNet_TF)

realNet_LNC<-read.table(file.choose(),header=T,sep="\t")
head(realNet_LNC);dim(realNet_LNC)

realNet_miR<-read.table(file.choose(),header=T,sep="\t")
head(realNet_miR);dim(realNet_miR)

#resOrdered<-as.data.frame(realNet[order(realNet$Score),])
#head(resOrdered)
#resOrdered_f<-resOrdered[resOrdered$Methods!='Prediction',] #miR
#head(resOrdered_f)
#write.table(resOrdered_f,'RAID_miRNA_mRNA_score-0.5.txt',sep = '\t',row.names = F,quote = F,col.names = colnames(realNet))
#realNet<-resOrdered_f[,c(6,8)]

realNet_TF<-realNet_TF[,c(2,3)]; head(realNet_TF)
realNet_LNC<-realNet_LNC[,c(6,8)]; head(realNet_LNC) 
realNet_miR<-realNet_miR[,c(4,6)]; head(realNet_miR)

realNet_TF<-realNet_TF[,c(1,3)]; head(realNet_TF)
realNet_LNC<-realNet_LNC[,c(2,8)]; head(realNet_LNC) 
realNet_miR<-realNet_miR[,c(2,6)]; head(realNet_miR)

net_mRNA_TF<-unique(as.vector(realNet_TF[,2]))
length(unique(as.vector(realNet_TF[,1])))#543(animal-trrust) 785(TRRUST-TFs) 1371(RAID-lnc) 2596(RAID-miR) 0.5置信node
length(unique(as.vector(realNet_TF[,2])))#2319(animal-trrust) 7809(RAID-lnc) 15802(RAID-miR) edges

net_mRNA_LNC<-unique(as.vector(realNet_LNC[,2]))
length(unique(as.vector(realNet_LNC[,1])))
length(unique(as.vector(realNet_LNC[,2])))

net_mRNA_miR<-unique(as.vector(realNet_miR[,2]))
length(unique(as.vector(realNet_miR[,1])))
length(unique(as.vector(realNet_miR[,2])))

library(igraph)
g_TF<-graph.data.frame(realNet_TF,directed = F)
g_TF<-simplify(g_TF) 
net_TF<-get.adjacency(g_TF)

g_LNC<-graph.data.frame(realNet_LNC,directed = F)
g_LNC<-simplify(g_LNC) 
net_LNC<-get.adjacency(g_LNC)

g_miR<-graph.data.frame(realNet_miR,directed = F)
g_miR<-simplify(g_miR) 
net_miR<-get.adjacency(g_miR)

getNeighbor<-function(gene,net){
  neighbor=c()
  for (i in gene){
    neighbor=c(neighbor,which(net[i,]>0))
  }
  neighbor=colnames(net)[unique(neighbor)]
}

cc<-c('4','6','17','29')
pfc<-c('2','6','8')

write.table(paste('sample','module','Gene_id/Mirbase_id', 'Pivot_type','module_links','p_value','target_genes', sep = '\t'),"Pivot regulators of crosstalk modules.txt",append=T,row.names=F,col.names = F,quote = F);
write.table(paste('sample','module','pivot','m_gene','Pivot_type',sep = '\t'),"all_pivots_edges.txt",append=T,row.names=F,col.names = F,quote = F);

All_module_gene<-read.table(file.choose(),stringsAsFactors = F,header = T) #模块号+基因
number<-sort(unique(All_module_gene[,1]));number
head(All_module_gene)

library(tidyr)
for(j in cc){ #字符时可以in本身，数字时用length
  #j=4
  m<-as.vector(All_module_gene[All_module_gene$module==j,][,2])# 模块基因
  mn<-intersect(net_mRNA_TF,m);length(mn) # 多少在网络里
  m_mi<-getNeighbor(mn,net_TF); length(m_mi) #第一次正向邻居（不同mRNA可能有相同邻居，但R包取邻居不用去重）
  
  nodes=1:nrow(net_TF); #1 to pivot的全部节点数
  names(nodes)=rownames(net_TF);#nodes有两行，上边是名字，下边是本身（序号）
  #mgene=as.data.frame(m);#列名为‘m’
  #modNode=merge(realNet,mgene,by.x ='gene',by.y = 'm' )
  
  for(k in 1:length(m_mi)){
    nodeNei=getNeighbor(c(m_mi[k]),net_TF);# 第二次某pivot反向邻居（调控的全部网络基因）
    sig=intersect(nodeNei,m);#pivot调控的模块基因(一个模块内不用去重)
    neiAndM=length(sig);
    
    if(neiAndM > 1){
      p=1-phyper(neiAndM-1, length(nodeNei),nrow(net_TF)-length(nodeNei),length(mn));
      #（pivot调控模块内基因数-1，pivot反向总调控数，网络总节点数-pivot反向总调控数，模块基因在网数）
      if(p<0.05 & p>0){
        print(neiAndM/length(m))
        
        neiAndM2<-as.data.frame(t(sig)); colnames(neiAndM2)<-c(1:dim(neiAndM2)[2])
        neiAndM3<-unite(neiAndM2,"combine",colnames(neiAndM2), sep="/", remove = F)[,1]
        
        write.table(paste('CC',j,names(nodes[m_mi[k]]),'TF', neiAndM,p,neiAndM3, sep = '\t'),"Pivot regulators of crosstalk modules.txt",append=T,row.names=F,col.names = F,quote = F);
        write.table(paste('CC',j,names(nodes[m_mi[k]]),sig,'TF',sep = '\t'),"all_pivots_edges.txt",append=T,row.names=F,col.names = F,quote = F);
      }
    }
  }
}

for(j in cc){ #字符时可以in本身，数字时用length
  #j=4
  m<-as.vector(All_module_gene[All_module_gene$module==j,][,2])# 模块基因
  mn<-intersect(net_mRNA_LNC,m);length(mn) # 只是看看多少在网络里
  m_mi<-getNeighbor(mn,net_LNC); length(m_mi) #第一次正向邻居（不同mRNA可能有相同邻居，但R包取邻居不用去重）
  
  nodes=1:nrow(net_LNC);
  names(nodes)=rownames(net_LNC);#nodes有两行，上边是名字，下边是本身（序号）
  #mgene=as.data.frame(m);#列名为‘m’
  #modNode=merge(realNet,mgene,by.x ='gene',by.y = 'm' )
  
  for(k in 1:length(m_mi)){
    nodeNei=getNeighbor(c(m_mi[k]),net_LNC);# 某pivot反向邻居（调控的全部网络基因）
    sig=intersect(nodeNei,m);#调控的模块基因(一个模块内不用去重)
    neiAndM=length(sig);
    
    if(neiAndM > 1){
      p=1-phyper(neiAndM-1, length(nodeNei),nrow(net_LNC)-length(nodeNei),length(mn));
      if(p<0.05 & p>0){
        print(neiAndM/length(m))
        
        neiAndM2<-as.data.frame(t(sig)); colnames(neiAndM2)<-c(1:dim(neiAndM2)[2])
        neiAndM3<-unite(neiAndM2,"combine",colnames(neiAndM2), sep="/", remove = F)[,1]
        
        write.table(paste('CC',j,names(nodes[m_mi[k]]),'LncRNA', neiAndM,p,neiAndM3, sep = '\t'),"Pivot regulators of crosstalk modules.txt",append=T,row.names=F,col.names = F,quote = F);
        write.table(paste('CC',j,names(nodes[m_mi[k]]),sig,'LncRNA',sep = '\t'),"all_pivots_edges.txt",append=T,row.names=F,col.names = F,quote = F);
      }
    }
  }
}

for(j in cc){ #字符时可以in本身，数字时用length
  #j=4
  m<-as.vector(All_module_gene[All_module_gene$module==j,][,2])# 模块基因
  mn<-intersect(net_mRNA_miR,m);length(mn) # 只是看看多少在网络里
  m_mi<-getNeighbor(mn,net_miR); length(m_mi) #第一次正向邻居（不同mRNA可能有相同邻居，但R包取邻居不用去重）
  
  nodes=1:nrow(net_miR);
  names(nodes)=rownames(net_miR);#nodes有两行，上边是名字，下边是本身（序号）
  #mgene=as.data.frame(m);#列名为‘m’
  #modNode=merge(realNet,mgene,by.x ='gene',by.y = 'm' )
  
  for(k in 1:length(m_mi)){
    nodeNei=getNeighbor(c(m_mi[k]),net_miR);# 某pivot反向邻居（调控的全部网络基因）
    sig=intersect(nodeNei,m);#调控的模块基因(一个模块内不用去重)
    neiAndM=length(sig);
    
    if(neiAndM > 1){
      p=1-phyper(neiAndM-1, length(nodeNei),nrow(net_miR)-length(nodeNei),length(mn));
      if(p<0.05 & p>0){
        print(neiAndM/length(m))
        
        neiAndM2<-as.data.frame(t(sig)); colnames(neiAndM2)<-c(1:dim(neiAndM2)[2])
        neiAndM3<-unite(neiAndM2,"combine",colnames(neiAndM2), sep="/", remove = F)[,1]
        
        write.table(paste('CC',j,names(nodes[m_mi[k]]),'miRNA', neiAndM,p,neiAndM3, sep = '\t'),"Pivot regulators of crosstalk modules.txt",append=T,row.names=F,col.names = F,quote = F);
        write.table(paste('CC',j,names(nodes[m_mi[k]]),sig,'miRNA',sep = '\t'),"all_pivots_edges.txt",append=T,row.names=F,col.names = F,quote = F);
      }
    }
  }
}

#=================
All_module_gene<-read.table(file.choose(),stringsAsFactors = F,header = T) #基因+模块号
number<-sort(unique(All_module_gene[,1]));number
head(All_module_gene)

for(j in pfc){ #字符时可以in本身，数字时用length
  #j=4
  m<-as.vector(All_module_gene[All_module_gene$module==j,][,2])# 模块基因
  mn<-intersect(net_mRNA_TF,m);length(mn) # 只是看看多少在网络里
  m_mi<-getNeighbor(mn,net_TF); length(m_mi) #第一次正向邻居（不同mRNA可能有相同邻居，但R包取邻居不用去重）
  
  nodes=1:nrow(net_TF);
  names(nodes)=rownames(net_TF);#nodes有两行，上边是名字，下边是本身（序号）
  #mgene=as.data.frame(m);#列名为‘m’
  #modNode=merge(realNet,mgene,by.x ='gene',by.y = 'm' )
  
  for(k in 1:length(m_mi)){
    nodeNei=getNeighbor(c(m_mi[k]),net_TF);# 某pivot反向邻居（调控的全部网络基因）
    sig=intersect(nodeNei,m);#调控的模块基因(一个模块内不用去重)
    neiAndM=length(sig);
    
    if(neiAndM > 1){
      p=1-phyper(neiAndM-1, length(nodeNei),nrow(net_TF)-length(nodeNei),length(mn));
      if(p<0.05 & p>0){
        print(neiAndM/length(m))
        
        neiAndM2<-as.data.frame(t(sig)); colnames(neiAndM2)<-c(1:dim(neiAndM2)[2])
        neiAndM3<-unite(neiAndM2,"combine",colnames(neiAndM2), sep="/", remove = F)[,1]
        
        write.table(paste('PFC',j,names(nodes[m_mi[k]]),'TF', neiAndM,p,neiAndM3, sep = '\t'),"Pivot regulators of crosstalk modules.txt",append=T,row.names=F,col.names = F,quote = F);
        write.table(paste('PFC',j,names(nodes[m_mi[k]]),sig,'TF',sep = '\t'),"all_pivots_edges.txt",append=T,row.names=F,col.names = F,quote = F);
      }
    }
  }
}

for(j in pfc){ #字符时可以in本身，数字时用length
  #j=5
  m<-as.vector(All_module_gene[All_module_gene$module==j,][,2])# 模块基因
  mn<-intersect(net_mRNA_LNC,m);length(mn) # 只是看看多少在网络里
  m_mi<-getNeighbor(mn,net_LNC); length(m_mi) #第一次正向邻居（不同mRNA可能有相同邻居，但R包取邻居不用去重）
  
  nodes=1:nrow(net_LNC);
  names(nodes)=rownames(net_LNC);#nodes有两行，上边是名字，下边是本身（序号）
  #mgene=as.data.frame(m);#列名为‘m’
  #modNode=merge(realNet,mgene,by.x ='gene',by.y = 'm' )
  
  for(k in 1:length(m_mi)){
    nodeNei=getNeighbor(c(m_mi[k]),net_LNC);# 某pivot反向邻居（调控的全部网络基因）
    sig=intersect(nodeNei,m);#调控的模块基因(一个模块内不用去重)
    neiAndM=length(sig);
    
    if(neiAndM > 1){
      p=1-phyper(neiAndM-1, length(nodeNei),nrow(net_LNC)-length(nodeNei),length(mn));
      if(p<0.05 & p>0){
        print(neiAndM/length(m))
        
        neiAndM2<-as.data.frame(t(sig)); colnames(neiAndM2)<-c(1:dim(neiAndM2)[2])
        neiAndM3<-unite(neiAndM2,"combine",colnames(neiAndM2), sep="/", remove = F)[,1]
        
        write.table(paste('PFC',j,names(nodes[m_mi[k]]),'LncRNA', neiAndM,p,neiAndM3, sep = '\t'),"Pivot regulators of crosstalk modules.txt",append=T,row.names=F,col.names = F,quote = F);
        write.table(paste('PFC',j,names(nodes[m_mi[k]]),sig,'LncRNA',sep = '\t'),"all_pivots_edges.txt",append=T,row.names=F,col.names = F,quote = F);
      }
    }
  }
}

for(j in pfc){ #字符时可以in本身，数字时用length
  #j=2
  m<-as.vector(All_module_gene[All_module_gene$module==j,][,2])# 模块基因
  mn<-intersect(net_mRNA_miR,m);length(mn) # 只是看看多少在网络里
  m_mi<-getNeighbor(mn,net_miR); length(m_mi) #第一次正向邻居（不同mRNA可能有相同邻居，但R包取邻居不用去重）
  
  nodes=1:nrow(net_miR);
  names(nodes)=rownames(net_miR);#nodes有两行，上边是名字，下边是本身（序号）
  #mgene=as.data.frame(m);#列名为‘m’
  #modNode=merge(realNet,mgene,by.x ='gene',by.y = 'm' )
  
  for(k in 1:length(m_mi)){
    nodeNei=getNeighbor(c(m_mi[k]),net_miR);# 某pivot反向邻居（调控的全部网络基因）
    sig=intersect(nodeNei,m);#调控的模块基因(一个模块内不用去重)
    neiAndM=length(sig);
    
    if(neiAndM > 1){
      p=1-phyper(neiAndM-1, length(nodeNei),nrow(net_miR)-length(nodeNei),length(mn));
      if(p<0.05 & p>0){
        print(neiAndM/length(m))
        
        neiAndM2<-as.data.frame(t(sig)); colnames(neiAndM2)<-c(1:dim(neiAndM2)[2])
        neiAndM3<-unite(neiAndM2,"combine",colnames(neiAndM2), sep="/", remove = F)[,1]
        
        write.table(paste('PFC',j,names(nodes[m_mi[k]]),'miRNA', neiAndM,p,neiAndM3, sep = '\t'),"Pivot regulators of crosstalk modules.txt",append=T,row.names=F,col.names = F,quote = F);
        write.table(paste('PFC',j,names(nodes[m_mi[k]]),sig,'miRNA',sep = '\t'),"all_pivots_edges.txt",append=T,row.names=F,col.names = F,quote = F);
      }
    }
  }
}
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 194,457评论 5赞 459
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 81,837评论 2赞 371
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 141,696评论 0赞 319
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 52,183评论 1赞 263
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 61,057评论 4赞 355
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 46,105评论 1赞 272
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 36,520评论 3赞 381
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 35,211评论 0赞 253
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 39,482评论 1赞 290
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,574评论 2赞 309
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,353评论 1赞 326
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,213评论 3赞 312
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 37,576评论 3赞 298
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,897评论 0赞 17
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,174评论 1赞 250
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 41,489评论 2赞 341
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 40,683评论 2赞 335