cmp_cr

#!/usr/bin/python

#this program is for calculating the clusters overlapping between to PRIDE Cluster results
#usage: compare clusterresult1 clusterresult2
#author mingze
#2016.10


from __future__ import division 
import re
import sys 
import json

fileName1 = sys.argv[1]
fileName2 = sys.argv[2]

resultFile1 = open(fileName1, 'r')
resultFile2 = open(fileName2, 'r')

similarSpectraNoinF1 = 0;
spectraNoinF1 = 0;
similarSpectraNoinF2 = 0;
spectraNoinF2 = 0;

nodes = [];
edges = [];


def readAClusterFromLines(lineList):
	specList = []
	
	i = 0
	line = lineList[i] 
	idMatch = re.match(r'id=(.*?)\n',line)
	while(not idMatch): #to the id line
		i = i +1
		if i >= len(lineList):
			return (None, specList)
		line = lineList[i]
		idMatch = re.match(r'id=(.*?)\n',line)
 
	if idMatch:
		id = idMatch.group(1)
	else:
		print "error or EndofFile, it is not begin with 'id=*'"
		print idLine
		return (None, specList)

	i = i + 6
	line = lineList[i]  #jump to the SPEC line
	
	SPECMatch = re.match(r'SPEC\t(.*?)\t',line)
	while(SPECMatch):	
		specList.append(SPECMatch.group(1))
		i = i + 1
		if i >= len(lineList) : break
		line = lineList[i]
		SPECMatch = re.match(r'SPEC\t(.*?)\t',line)

	return (id,specList)	


def readLinesForACluster(resultFile):
	lines = []

	line = resultFile.readline()
	if line == None or line == '':
		return lines	

	while(line != '=Cluster=\n' ):
		lines.append(line)
		line = resultFile.readline()
		if line == None or line == '':
			break

	return lines


def isClusterSame(cluster1, cluster2):
	(id1, len1, specList1) = cluster1
	(id2, len2, specList2) = cluster2
#	print "comparing: " + id1 + "---" + id2	
	if len1 != len2 :
		return False
	for i in range(0,len1):
		if cmp(specList1[i], specList2[i]) != 0:
			return False
	return True

def getClusterSimilarity(cluster1, cluster2):
	global similarSpectraNoinF1
	global similarSpectraNoinF2
	(id1, len1, specList1) = cluster1
	(id2, len2, specList2) = cluster2
#	print "comparing: " + id1 + "---" + id2	
	n = 0
	i = 0
	j = 0
	while (i < len1 and j < len2):
		cmp_score = cmp(specList1[i], specList2[j])
		if  cmp_score == 0:
			i = i + 1
			j = j + 1
			n = n+1	
		if cmp_score < 0:
			i = i + 1
		if cmp_score > 0:
			j = j + 1
	
	similarSpectraNoinF1 = similarSpectraNoinF1 + n;
	similarSpectraNoinF2 = similarSpectraNoinF2 + n;
	return (n,2*n/(len1+len2))



clusterList1 = []
clusterList2 = []


distribute1 = {}
distribute2 = {}

#read from files
print "start reading file 1: " + fileName1 
i = 1
while(i>=1):
	bufferLines = readLinesForACluster(resultFile1)
	if len(bufferLines) < 1:
		break
	(id1,specList1) = readAClusterFromLines(bufferLines)
	
	if(id1 == None):
		continue	
	if len(specList1 ) > 1: 
		clusterList1.append((id1, len(specList1), sorted(specList1)))
		node_data = {
			"name":"Grp1_" + str(i),
			"id":id1,
			"groupNo":1,
			"length":len(specList1)
		}
		node = {
			"data":node_data
		}
		nodes.append(node)
		distribute1[len(specList1)] = distribute1.get(len(specList1),0) + 1
		spectraNoinF1 = spectraNoinF1 + len(specList1)	
	i = i + 1
print "end of file 1"
print "start reading file 2: " +fileName2 

i = 1
while(i>=1):
	bufferLines = readLinesForACluster(resultFile2)
	if len(bufferLines) < 1:
		break

	(id2,specList2) = readAClusterFromLines(bufferLines)
	if(id2 == None):
		continue	
	if len(specList2 ) > 1: 
		clusterList2.append((id2, len(specList2), sorted(specList2)))
		node_data = {
			"name":"Grp2_" + str(i),
			"id":id2,
			"groupNo":2,
			"length":len(specList2)
		}
		node = {
			"data":node_data
		}
		nodes.append(node)
		distribute2[len(specList2)] = distribute2.get(len(specList2),0) + 1
		spectraNoinF2 = spectraNoinF2 + len(specList2)	
	i = i + 1
print "end of reading file 2"


#which one is bigger?#
allClusters = {}
if len(clusterList1) >= len(clusterList2):
	smallList = 2
	bigList = 1
	NoSpectraInSmallList = spectraNoinF2 
	NoSpectraInBigList = spectraNoinF1 
	allClusters = {
		"bigClusterList" : clusterList1,	
		"smallClusterList" : clusterList2
	}
else:
	smallList = 1
	bigList = 2
	NoSpectraInSmallList = spectraNoinF1 
	NoSpectraInbigList = spectraNoinF2 
	allClusters = {
		"bigClusterList" : clusterList2,	
		"smallClusterList" : clusterList1
	}


#calculate and store the similarities
#prepare the network data
#we assume the arrow point from bigNode(small cluster List) to smallNode (big cluster List)
identical_num = 0
similarityDist = {}
for i in range(0,len(allClusters['smallClusterList'])):
	for j in range(0, len(allClusters['bigClusterList'])):
		(sharedSpectraNo,similarity) = getClusterSimilarity(allClusters['smallClusterList'][i], allClusters['bigClusterList'][j])
		if similarity == 1.0:
			identical_num = identical_num + 1
		if sharedSpectraNo > 0:
			edge_data = {
			"source":allClusters['smallClusterList'][i][0], #get id
			"target":allClusters['bigClusterList'][j][0],
			"sharedSpectraNo":sharedSpectraNo,
			"similarity":similarity
			}
			edge = {"data":edge_data}
			edges.append(edge)
		similarityDist[int(similarity*10)] = similarityDist.get(int(similarity*10),0) + 1
		

#output the statistics of the files
print "\n\n"
print "the statistics of the files"
print "----------------------------------------"
print str(len(clusterList1)) + "\t clusters in file 1: " + fileName1
print str(len(clusterList2)) + "\t clusters in file 2: " + fileName2
print str(identical_num) + "\t same clusters between them"
print str(similarSpectraNoinF1/(1.0*spectraNoinF1)) + "x" + str(spectraNoinF1) + "=" + str(similarSpectraNoinF1) + " shared spectra from " + fileName1
print str(similarSpectraNoinF2/(1.0*spectraNoinF2)) + "x" + str(spectraNoinF2) + "=" + str(similarSpectraNoinF2) + " shared spectra from file" + fileName2
print "----------------------------------------"


#output the distribution of similarities:
print "\n\n"
print "the distribution of cluster size in file 1 " + fileName1
print "----------------------------------------"
print "cluster size\tNo."
print "----------------------------------------"
for i in sorted(distribute1.keys()):
	print str(i) + "\t\t" + str(distribute1[i])
print "----------------------------------------"

print "\n\n"
print "the distribution of cluster size in file 2 " + fileName2
print "----------------------------------------"
print "cluster size\tNo."
print "----------------------------------------"
for i in sorted(distribute2.keys()):
	print str(i) + "\t\t" + str(distribute2[i])
print "----------------------------------------"

print "\n\n"
print "the distribution of similarity between them" 
print "----------------------------------------"
print "similarity\tNo."
print "----------------------------------------"
for i in sorted(similarityDist.keys()):
	print str(i) + "\t\t" + str(similarityDist[i])
print "----------------------------------------"


#calculate the stars#
starList = []
outterList = []
for i in range(0, len(nodes)):
	node = nodes[i]['data']
	if node['name'][3] == str(smallList) :
		outterIds = []
		star = {
			"starId":node['id'],
			"NoSpectra":node['length'],
			"outterIds": outterIds,
			"NoSharedSpectra": 0
		}
		starList.append(star)
	else :
		starIds = []	
		outter = {
			"outterId":node['id'],
			"NoSpectra":node['length'],
			"starIds": starIds,
			"NoSharedSpectra": 0
		}
		outterList.append(outter)

starConnectFlag = False
for i in range(0, len(edges)):
	edge = edges[i]['data']
	source = edge['source']
	target = edge['target']

	for j in range(0, len(starList)):
		if starList[j]['starId'] == source :	
			starList[j]['outterIds'].append(target)
			starList[j]['NoSharedSpectra'] += edge['sharedSpectraNo']
		
		if starList[j]['starId'] == target:	
			starList[j]['outterIds'].append(source)
			starList[j]['NoSharedSpectra'] += edge['sharedSpectraNo']
		
	for k in range(0, len(outterList)):
		if outterList[k]['outterId'] == source :	
			outterList[k]['starIds'].append(target)
			outterList[k]['NoSharedSpectra'] += edge['sharedSpectraNo']
			starConnectFlag = True
		
		if outterList[k]['outterId'] == target:	
			outterList[k]['starIds'].append(source)
			outterList[k]['NoSharedSpectra'] += edge['sharedSpectraNo']
			starConnectFlag = True


#get the stand alone nodes, calculate the No of Lost Spectra and the distribution of dividing factors
standAloneStars = []
standAloneOutters = []
totalLostSpectraInStars = 0
totalLostSpectraInOutters = 0
dividFactorDist = {}
for i in range(0,len(starList)):
	if len(starList[i]['outterIds']) == 0 :
		standAloneStars.append(starList[i])

	divideFactor = len(starList[i]['outterIds']) + starList[i]['NoSpectra'] - starList[i]['NoSharedSpectra']
	dividFactorDist[divideFactor] = dividFactorDist.get(divideFactor,0) + 1

	if starList[i]['NoSpectra'] < starList[i]['NoSharedSpectra'] :
		print "!!!!!!!!!!!!Wrong here, total spectra No is less than Shared Spectra with outters"
		print "with star " + str(i) + str(starList[i]['NoSpectra']) + "is less than" + str(starList[i]['NoSharedSpectra'])
	
	if starList[i]['NoSpectra'] > starList[i]['NoSharedSpectra'] :
		totalLostSpectraInStars += starList[i]['NoSpectra'] - starList[i]['NoSharedSpectra']

for i in range(0,len(outterList)):
	if len(outterList[i]['starIds']) == 0 :
		standAloneOutters.append(outterList[i])

	if outterList[i]['NoSpectra'] < outterList[i]['NoSharedSpectra'] :
		print "!!!!!!!!!!!!Wrong here, total spectra No is less than Shared Spectra with stars"
		print "with star " + str(i) + str(outterList[i]['NoSpectra']) + "is less than" + str(outterList[i]['NoSharedSpectra'])

	if outterList[i]['NoSpectra'] > outterList[i]['NoSharedSpectra'] :
		totalLostSpectraInOutters += outterList[i]['NoSpectra'] - outterList[i]['NoSharedSpectra']

#output the No. of losted spectra 
print "\n\n----------------------------------------"
print "totally lost " + str(totalLostSpectraInStars) + "(" + str(totalLostSpectraInStars/(0.01 * NoSpectraInSmallList ))  +"%)" + "spectra in Stars"
print "totally lost " + str(totalLostSpectraInOutters) + "(" + str(totalLostSpectraInOutters/(0.01 * NoSpectraInBigList ))  +"%)" + "spectra in Outters"
print "----------------------------------------"


#output the average dividing factor
#each of a losted spectrum will be considered as one singleton cluster, which contribute equelly to the normal clusters
aveDvdFactor = (len(edges)+totalLostSpectraInStars)/(1.0 * len(starList))
print "\n\n"
print "average dividing factor ---(each of a losted spectrum will be considered as one singleton cluster, which contribute equelly to the normal clusters)"
print "----------------------------------------"
print aveDvdFactor
print "----------------------------------------"


#output the distribution of dividing factors
print "\n\n"
print "the distribution of dividing factors"
print "----------------------------------------"
print "Divide Factors \t No."
print "----------------------------------------"
for i in dividFactorDist.keys() :
	print str(i) + "\t\t" + str(dividFactorDist[i])
print "----------------------------------------"


#output stand alone nodes
print "\n\n"
print "No of stand alone clusters"
print "----------------------------------------"
print "small group\tbig group"
print "----------------------------------------"
print str(len(standAloneStars)) + "\t\t" + str(len(standAloneOutters))


#output star connecting info 
if (starConnectFlag == True) :
	print "\n\n"
	print "some stars are connected by these outters"
	print "----------------------------------------"
	NoConnectOutters = 0
	for i in range(0,len(outterList)):
		if len(outterList[i]['starIds']) > 1 :
			NoConnectOutters += 1
#			print "outter " + str(i) + "has " + str(len(outterList[i]['starIds'])) + " lines out"
		if len(outterList[i]['starIds']) > 2 :
			print "outter " + str(i) + "has " + str(len(outterList[i]['starIds'])) + " lines out(>1), should be careful "
	print "total No."
	print "----------------------------------------"
	print str(NoConnectOutters)
	print "----------------------------------------"



#output the Network data to cytoscape JSON file

networkDataInJson = {
	"elements":{
		"nodes" : nodes,
		"edges" : edges
	}      
}

networkString = json.dumps(networkDataInJson)

cytoscapeFile  = open("networkData.json", 'w')
cytoscapeFile.write(networkString)
cytoscapeFile.close()

resultFile1.close()
resultFile2.close()