Ns126:Calendar/NOTES/2017-4-30
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)),])) COR[[i]]<-cor(t(data),method="spearman") print(i) } gene<-c() for(i in 1:nrow(COR[[1]])){ for(j in 1:nrow(COR[[1]])){ for(z in 1:26){ SD=sd(c(COR[[1]][i][j],COR[[2]][i][j],COR[[3]][i][j], COR[[4]][i][j],COR[[5]][i][j],COR[[6]][i][j], COR[[7]][i][j],COR[[8]][i][j],COR[[9]][i][j], COR[[10]][i][j],COR[[11]][i][j],COR[[12]][i][j], COR[[13]][i][j],COR[[14]][i][j],COR[[15]][i][j], COR[[16]][i][j],COR[[17]][i][j],COR[[18]][i][j], COR[[19]][i][j],COR[[20]][i][j],COR[[21]][i][j], COR[[22]][i][j],COR[[23]][i][j],COR[[24]][i][j], COR[[25]][i][j],COR[[26]][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(COR[[1]])){ # print(i) # for(j in 1:nrow(COR[[1]])){ # for(z in 1:4){ # SD=try(sd(c(COR[[1]][i,j],COR[[2]][i,j],COR[[3]][i,j],COR[[4]][i,j]),na.rm=T)) # if(!is.na(SD) && SD>0.6){ # gene<-c(gene,i,j) # } # } # # } #}