Identify cluster specific network
Jump to navigation
Jump to search
library(Seurat) setwd("/media/Home_Raid1/zhl002/NAS1/RNA_seq/hiseq_020617/seurat_analysis") load("./ipsnt_33k_latest3.RData") cluster<-as.matrix(pbmc33k.merged@data[,which(pbmc33k.merged@ident==1)]) COR<-list() for(i in 1:26) { cluster<-as.matrix(pbmc33k.merged@data[,which(pbmc33k.merged@ident==i)]) data<-(na.omit(cluster[match(as.character(pbmc33k.merged@var.genes),rownames(cluster)),])) CORi<-cor(t(data),method="spearman") print(i) } gene<-c() for(i in 1:nrow(COR1)){ for(j in 1:nrow(COR1)){ for(z in 1:26){ SD=sd(c(COR1[i][j],COR2[i][j],COR3[i][j], COR4[i][j],COR5[i][j],COR6[i][j], COR7[i][j],COR8[i][j],COR9[i][j], COR10[i][j],COR11[i][j],COR12[i][j], COR13[i][j],COR14[i][j],COR15[i][j], COR16[i][j],COR17[i][j],COR18[i][j], COR19[i][j],COR20[i][j],COR21[i][j], COR22[i][j],COR23[i][j],COR24[i][j], COR25[i][j],COR26[i][j],na.rm=T),na.rm=T) if(SD>0.6){ gene<-c(gene,i,j) } } } } ### example #gene<-c() #for(i in 1:nrow(COR1)){ # print(i) # for(j in 1:nrow(COR1)){ # for(z in 1:4){ # SD=try(sd(c(COR1[i,j],COR2[i,j],COR3[i,j],COR4[i,j]),na.rm=T)) # if(!is.na(SD) && SD>0.6){ # gene<-c(gene,i,j) # } # } # # } #}