-
Notifications
You must be signed in to change notification settings - Fork 0
/
MST_location.R
66 lines (46 loc) · 1.42 KB
/
MST_location.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#Args <- commandArgs(TRUE)
library('ConsensusClusterPlus')
argv <- commandArgs(TRUE)
work_dir=argv[1]
work_dir
#set the working dir
setwd(work_dir)
marker='Col10a1'
data=read.table("PCA_result.tsv",sep='\t',header=F,row.names=1);
#data=t(data[,1:3])
#distence=1-cor(data)
distence=dist(data[,1:3],upper=T,diag=T)
require(igraph)
g <- graph.adjacency(as.matrix(distence), weighted=T, mode = 'undirected')
g <- simplify(g)
#extract the MST
mst <- minimum.spanning.tree(g)
#get handtune structure
tkid<-tkplot(mst)
layid<-tkplot.getcoords(tkid)
gene_expression=read.csv('raw_data.csv',header=T,row.names=1);
genes=rownames(gene_expression)
cells=colnames(gene_expression)
k=0
for(i in 1:length(genes))
if (genes[i]==marker)
k=i
gene_expression=gene_expression[k,]
require(plotrix)
gene_expression=t(gene_expression)
mypalette <- color.scale(gene_expression,c(0,1,1),c(1,1,0),0,color.spec="rgb")
plot(mst, vertex.size=5,layout =layid,vertex.color=mypalette,vertex.label=NA,edge.width=3)
col.labels<-c("Low","Medium","High")
testcol<-color.gradient(c(0,1,1),c(1,1,0),0,nslices=100)
color.legend(-1.4,0.5,-1.3,1.3,rect.col=testcol,legend=col.labels,gradient="y")
#dev.off()
#calculating the longest path
#diameter(g)
#calculating the longes path and giving back the corresponding node IDs
#farthest.nodes(g)
st=cells[farthest.nodes(g)[1]]
#write.table(distence,'distence.txt')
write.table(st,'start.txt')
sink('edges.txt')
str(mst)
sink()