#clear self correlations
library(igraph)
library(ggplot2)
igraphme<-function(csv,censor="RESTBP_SDCensor")
{	print("Starting graph analysis, Good Luck!")
	subj=strsplit(csv,"_")[[1]][1]
	cen=read.delim(paste(subj,censor,sep=""),header=F)
	trimcsv=read.table(csv,header=T)[cen$V1==1,c(-1:-2)]
#	trimcsv=read.delim(csv)[cen$V1==1,c(-1:-3)]
	connMat=cor(trimcsv)
	diag(connMat) = rep(0, length(diag(connMat)))
	write.table(connMat,paste("corrmat_",subj,".mat",sep=""),row.names=F, col.names=F,sep="\t")
	 density = 0.1
	 nEdges = length(connMat[upper.tri(connMat)])*density
	 thresh = sort( connMat[upper.tri(connMat)], decreasing=T)[nEdges]
	 adj = 1*(connMat >= thresh)
	 bingraph = graph.adjacency(adj, mode="undirected", weighted=NULL)
	 components = clusters(bingraph)
	 maxID = which(components$csize == max(components$csize))[1]
	 adj[components$membership!=maxID,] = 0
	 adj[,components$membership!=maxID] = 0
	 graph = graph.adjacency( adj, mode="directed", weighted=NULL )
	 deg = degree(graph)
	 deg[deg==0] = NA
	 pathsmat = shortest.paths(graph, weights=NA)
	 pathsmat[!is.finite(pathsmat)] = NA
	 paths = rowMeans(pathsmat, na.rm=TRUE)
	 paths[paths==0] = NA
	 clust = transitivity(graph, type="local")
	 clust[deg < 2] = NA
	 pager = page.rank(graph)$vector
	 pager[deg < 2] = NA
	 leff <- numeric(length(deg))
	 goodnodes <- which(deg > 1)
	 leff[goodnodes] <- sapply(goodnodes, function(x) {
	  neighbs <- neighbors(graph, v=x)
	  g.sub <- induced.subgraph(graph, neighbs)
	  Nv <- vcount(g.sub)
	  lpaths <- shortest.paths(g.sub, weights=NA)
	  lpaths <- paths[upper.tri(lpaths)]
	  pathsup <- lpaths[upper.tri(lpaths)]
	  2 / Nv / (Nv - 1) * sum(1 / lpaths[which(is.na(lpaths)==FALSE)]) })
	 leff[deg < 2] = NA
	 leff[which(is.na(deg)==TRUE)] = NA
	 nNodes = length(deg)
	 cnode.dat = data.frame(Node=rep(1:nNodes,5))
	 cnode.dat$Value = c( deg, paths, leff, clust, pager )
	 cnode.dat$Metric = c( rep("Degree", nNodes), rep("Shortest Path", nNodes), rep("Local Efficiency", nNodes), rep("Clustering Coefficient", nNodes), rep("Page-Rank", nNodes) )
	 cnodePlot = ggplot(cnode.dat, aes(x=Node, y=Value, group=Metric, fill=Metric, colour=Metric))
	 cnodePlot = cnodePlot + geom_point()
	 cnodePlot = cnodePlot + theme(text=element_text(size=10), legend.position="none")
	 cnodePlot = cnodePlot + ggtitle("Node metrics")
	 cnodePlot = cnodePlot + facet_grid(Metric~ ., scales="free")
	 cnode.csv = data.frame(cbind(deg, paths, leff, clust, pager ))
	 geff<-1/(shortest.paths(graph))
	 geff[!is.finite(geff)]<-NA
	 geff<-mean(geff,na.rm=TRUE)
	 cc = transitivity(graph)
	 write.table(cnode.csv,paste("gMetrics_",subj,"igraph.txt",sep=""),row.names=T,sep="\t") #row.names=F & col.names=F for connmat save
	 pdf(paste("plot_",subj,"igraph.pdf",sep=""))
	  heatmap(connMat)
	  print(cnodePlot)
	 print("Output goes to NBS for further graph and network analysis.Hopefully, you can publish this in a high-impact journal and secure more funding!")
	 dev.off()
}

csvs=Sys.glob('*trim.csv')
for(csv in csvs){
	try(igraphme(csv))
}