Identify cluster specific network

From ZhangLabWiki
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)
# }
# }
#
# }
#}