#! /usr/bin/env python 
###################################
##    Davi Ortega - 05/18/2011   ##
###################################
import bitk
import sys
import copy
import random
import numpy
import os

if '-h' in sys.argv:
	print 'This program will produce a csv file to be loaded to excel table with several information about the alignemnt.\n\n \
	Requirements:\n \
	\t1) All sequences must be the same size in all alignments (including gaps) - Ver 1.1 might not needed \n \
	\tOptions:\n \
	\t\t -g2 group1.fa group2.fa\t Include comparison tests between the groups. \n \
	\t\t -in sequence_tag\tSequence Tag of reference organism for a.a. coordinates\n \
	\t\t -ss\tSecondary Structure: The fasta file with the ss information must be aligned with the other alignments.\n \
	\t\t -stride pdb_file.pdb sequence_tag\tRuns STRIDE from VMD and calculate the secondary structure of each residue using the structure file (pdb). The sequence_tag is a tag that must be present in one of the MSA \n \
	\t\t -ns\tNo sequence: won\'t display the sequences in the table\n \
	\t\t -c\tConsensus at 100%, 90%, 80%, 70%, 60% and 50% \n \
	\t\t -cg\tConsensus at 100%, 90%, 80%, 70%, 60% and 50% class dependent \n \
	\t\t -chi\tPerforme Chi_square test of two samples.\n \
	\t\t -lda\tPerforme Linear Discriminant Alalysis.\n \
	\t\t -aa_d\tBuild amino acid distribution table to be used in excel model form, aa_distr (import in A1)\n \
	\t\t -f\t Force lenght fix of alignments by building a profile alignemtn with mafft-profile and separating the sequences\n \
	\t\t -sasa pdb_file.pdb sequence_tag\t Include Solvent Accessible Surface Area score (using sasa and VMD). The sequence_tag is a tag that must be present in one of the MSA\n \
	\n\n\tSintaxe example 1:\tpython make_table2.py Fasta_alignment_class1.fa Fasta_alignment_class2.fa -ss SS_alignment.fa -in input_file\n \
	\t Sintaxe example 2:\tpython make_table2.py fasta_alignment.fa\n'
	sys.exit()


directory = os.getcwd()

# reading the alignements
if '-g2' in sys.argv:
	g2 = sys.argv.index('-g2')
	class1, class1_list = bitk.fastareader(sys.argv[g2+1])
	class2, class2_list = bitk.fastareader(sys.argv[g2+2])
else:
	class1, class1_list = bitk.fastareader(sys.argv[1])
	g2 = 0
if '-ss' in sys.argv:
	ss = sys.argv.index('-ss')
	SS, SS_list = bitk.fastareader(sys.argv[ss + 1])
else:
	ss = 0

if '-chi' in sys.argv:
	chi = 1
else:
	chi = 0

if '-lda' in sys.argv:
	lda = 1
else:
	lda = 0


# reading the input file
flag = 1
if '-in' not in sys.argv:
	ref_bug = class1.keys()[0]
	bias = 1
else:
	input = sys.argv.index('-in')
	ref_bug = sys.argv[input + 1]
	try:
		if int(sys.argv[input + 2]) == 0:
			bias = 1
		else:
			bias = int(sys.argv[input + 2])
	except ValueError:
		bias = 1
	except IndexError:
		bias = 1

# Check if the alignemnts have the same number of residues
len_old = 0
ALL = {}
all_seq = {}
for bug,seq in class1.iteritems():
	if len_old == 0 :
		len_old = len(seq)
	else:
		len_new = len(seq)
		if len_new != len_old:
			print 'Correct the alignements'
			print bug
			sys.exit()
		else:
			len_old = len_new
	ALL[bug]=seq
	all_seq[bug] = seq
	example = seq
print 'Class 1 check!'
print len_old
len1 = len_old
len_old = 0
if g2 != 0:
	for bug, seq in class2.iteritems():
		if len_old == 0 :
	                len_old = len(seq)
		else:
			len_new = len(seq)
			if len_new != len_old:
				print 'Correct the alignements'
				print bug
				print len_new
				sys.exit()
			else:
				len_old = len_new
		ALL[bug] = seq
		all_seq[bug] = seq
		example = seq
	print('Class 2 check')
	len2 = len_old
	if len1 != len2:
		print "MSAs not aligned"
		if "-f" in sys.argv:
			print "Force lenght fix with mafft profile by adj_msa_len"
			os.system("adj_msa_len " + sys.argv[g2+1] + " " + sys.argv[g2+2])
			class1, class1_list = bitk.fastareader(sys.argv[g2+1])
		        class2, class2_list = bitk.fastareader(sys.argv[g2+2])
			new_command = ''
			for arg in sys.argv:
				new_command += arg + " "
			print "Starting over..."
			print new_command
			os.system(new_command)
			sys.exit()
		
			

if ref_bug not in ALL.keys():
	print('referenced organism not found in the alignment: please check the sequence name')
	sys.exit()

print ' Pass lenght check test '


if '-sasa' in sys.argv:
	os.system('rm sasa.dat')
        sasa = sys.argv.index('-sasa')

	if sys.argv[sasa + 2] not in ALL.keys():
		print(' sasa referenced organism not found in the alignment: please check the sequence name or include the pdb sequence in the alignment')
	        sys.exit()
	ref_bug_sasa = sys.argv[sasa + 2]
        pdb_file = sys.argv[sasa + 1]
        tcl_script = '#This script runs on VMD and it is part of the program MSAAT1\n \
cd ' + directory + '\n \
set filesasa [open sasa.dat w]\n \
mol new {' + pdb_file +'} type {pdb} first 0 last 1 step 1 waitfor 1\n \
set protein [atomselect top protein]\n \
set chains [lsort -unique [$protein get chain]]\n \
foreach chain $chains {\n \
\tputs "I am here with chain $chain" \n \
\tset mol_res [[atomselect top "name CA and protein and chain $chain"] get resid]\n \
\tforeach res $mol_res {\n \
\t\tset sasa_res [measure sasa 1.4 $protein -restrict [atomselect top "chain $chain and resid $res"] ]\n \
\t\tputs $filesasa "$res;$chain;$sasa_res"\n \
\t\t}\n \
\t} \n \
puts "I am done... closing file"\n \
close $filesasa\n \
exit'
	file_sasa = open('sasa_script_tmp.tcl', 'w')
	file_sasa.write(tcl_script)
	file_sasa.close()
	print 'File writen'
	os.system("vmd -dispdev none -e " + directory + "/sasa_script_tmp.tcl")
#	os.system('ls -ll')
	file_data_sasa = open("sasa.dat", 'r')
	sasa = {}
	for line in file_data_sasa:
		field = line.split(';')
		if field[1] not in sasa.keys():
			sasa[field[1]] = {}
		sasa[field[1]][field[0]] = field[2].replace('\n','') #sasa {chain : residue : sasa score}
	file_data_sasa.close()
	
		
if '-stride' in sys.argv:
	stride = sys.argv.index('-stride')
	if sys.argv[stride + 2] not in ALL.keys():
                print(' sasa referenced organism not found in the alignment: please check the sequence name or include the pdb sequence in the alignment')
                sys.exit()
	ref_bug_stride = sys.argv[stride + 2]		

	pdb_file = sys.argv[stride + 1]
	tcl_script = '#This script runs on VMD and it is part of the program MSAAT1\n \
cd ' + directory + '\n \
global env\n \
set arch_name [vmdinfo arch]\n \
$env(VMDDIR)/stride_$arch_name ' + directory + '/' + pdb_file + ' > ' + directory + '/stride.dat\n \
exit'
	file_stride = open('stride_script_tmp.tcl', 'w')
        file_stride.write(tcl_script)
        file_stride.close()
	os.system("vmd -dispdev none -e " + directory + "/stride_script_tmp.tcl")
	file_data_stride = open("stride.dat", 'r')
	stride = {}
	for line in file_data_stride:
		if 'ASG' in line:
			field = line.split(' ')
			while '' in field:
				field.remove('')
			if field[2] not in stride.keys():
        	                stride[field[2]] = {}
                	stride[field[2]][field[3]] = field[6].replace('\n','') #sasa {chain : residue : sasa score}
        file_data_stride.close()



if ss != 0:
	for struct, seq in SS.iteritems():
		len_new = len(seq)
		if len_new != len_old:
			print 'Correct the alignments'
			print struct
			print len_new
			print seq
			print example
			sys.exit()
		else:
			len_old = len_new
		ALL[struct] = seq
		example = seq
	



# aa_dist

if g2 == 0 and "-aa_d" in sys.argv:
	os.system("aa_distr " + sys.argv[g2+1] + " all")
elif g2 != 0 and "-aa_d" in sys.argv:
	os.system("aa_distr " + sys.argv[g2+1] + " all")
	os.system("mv aa_distr.dat " + sys.argv[g2+1][:-2] + "dat")
	os.system("aa_distr " + sys.argv[g2+2] + " all")
        os.system("mv aa_distr.dat " + sys.argv[g2+2][:-2] + "dat")
	os.system("concat_aa_tables " + sys.argv[g2+1][:-2] + "dat " + sys.argv[g2+2][:-2] + "dat")


#calculations:




# Indentity (conservations score)

id_rank_all = bitk.identity_rank(all_seq)
if g2 != 0:
	id_rank_class1 = bitk.identity_rank(class1)
	id_rank_class2 = bitk.identity_rank(class2)

#Concensus
if '-c' in sys.argv:
	cons_all = bitk.consensus(all_seq, [100,90,80,70,60,50])
else:
	list_cons = [100,90,80,70,60,50]

if '-cg' in sys.argv:
	cons_class1 = bitk.consensus(class1, [100,90,80,70,60,50])
	cons_class2 = bitk.consensus(class2, [100,90,80,70,60,50])
	
	
#Rule number 1: select the residues 100% conserved - function aa_hist
 
hist_all = bitk.aa_hist(all_seq)
if g2 != 0:
	hist_class1 = bitk.aa_hist(class1)
	hist_class2 = bitk.aa_hist(class2)

# organizing data
output = 'ID;ref: ' + ref_bug + ';'
#seq_output = 'ID;ref: ' + ref_bug + ';'
seq_output = ''

if ss != 0:
	for struct in SS.keys():
		seq_output += struct + ';'
for bug in class1.keys():
	seq_output += bug + ';'
if g2 != 0:
	for bug in class2.keys():
		seq_output += bug + ';'

#if '-ns' not in sys.argv:
#	output = seq_output
	
list_bugs = seq_output[:-1].split(';')
#print list_bugs

if '-c' in sys.argv:
	list_cons = cons_all.keys()
	list_cons.sort()
	list_cons.reverse()
	for cons in list_cons:
		output += 'Consensus ' + str(cons) + '%;'
else:
	for cons in list_cons:
                output += 'Consensus ' + str(cons) + '%;'

if '-cg' in sys.argv:
	list_cons_1 = cons_class1.keys()
	list_cons_1.sort()
	list_cons_1.reverse()
	for cons in list_cons_1:
		output += 'Consensus ' + str(cons) + '% Group 1;'
	list_cons_2 = cons_class2.keys()
	list_cons_2.sort()
	list_cons_2.reverse()
	for cons in list_cons_2:
		output += 'Consensus ' + str(cons) + '% Group 2;'

output += 'Identity Rank ALL;'
if g2 != 0:
	output += 'Identity Rank Group 1;'
	output += 'Identity Rank Group 2;'
output += 'Rule 1: 100% conserved;'
output += 'Rule 2: 100% - 1 conserved;'
if g2 != 0:
	output += 'Rule 3: 100% conserved within groups;'
	output += 'Rule 4: 100% - 1 conserved within groups;'
	if chi != 0:
		output += 'Chi-square test between the two groups;'
		output += 'P-value of the Chi-square;'
	if lda != 0:
		output += 'Signal to noise Ratio;'
		output += 'Class Specific (1 yes,0 no);'

if '-sasa' in sys.argv:
	for chain in sasa.keys():
		output += 'SASA Score chain ' + chain + ';'

if '-stride' in sys.argv:
        for chain in stride.keys():
                output += 'STRIDE Assignment chain ' + chain + ';'


#if g2 != 0:
#	for bug in class2.keys():
#		seq_output += bug + ';'

if '-ns' not in sys.argv:
	output += seq_output


output += '\n'



J = {} #That's the dictionary of alignment coordinates and the reference coordinates
j = 1
j_sasa = 1
j_stride = 1



print 'Initiating Requested Analysis'

for i in range(len_new):
#	print '\nPosition ' + str(i) + ' being analyzed'
	if ALL[ref_bug][i] == '-':
		line = str(i) + ';X;'
	else:
		line = str(i) + ';' + str(j + bias - 1) + ';'
		j += 1
	J[j-1] = i

	if '-sasa' in sys.argv:
		if ALL[ref_bug_sasa][i] != '-':
			j_sasa += 1
	
	if '-stride' in sys.argv:
		if ALL[ref_bug_stride][i] != '-':
	                j_stride += 1	
        
	J[j-1] = i

	seq_line = ''
	for bug in list_bugs:
		seq_line += ALL[bug][i] + ';'
	
#	if '-ns' not in sys.argv:
#		line += seq_line
		
	if seq_line != 'X;' + '-;'*(len(list_bugs)-1):
		# Consensus
		if '-c' in sys.argv:
			for cons in list_cons:
				line += cons_all[cons][i] + ';'
		if '-cg' in sys.argv:
			for cons in list_cons_1:
				line += cons_class1[cons][i] + ';'
			for cons in list_cons_2:
				line += cons_class2[cons][i] + ';'
		# Conservation score
		if g2 != 0:
			line += str(id_rank_all[i]) + ';' + str(id_rank_class1[i]) + ';' + str(id_rank_class2[i]) + ';'
		else:
			line += str(id_rank_all[i]) + ';'
		# Rule number 1: 100% conserved positions
		if len(all_seq.keys()) in hist_all[i].values():
			if g2 != 0:
				line += '1;0;0;0;'
			else:
				line += '1;0;'
		else:
			line += '0;'
			#Rule number 2: 100% - 1 conserved position
			if len(all_seq.keys()) - 1  in hist_all[i].values():
				if g2 != 0:
					line += '1;0;0;'
				else:
					line += '1;'
			else:
				line += '0;'
				if g2 != 0:
					#Rule number 3 and 4 : 100% conserved within class
					if len(class1.keys()) in hist_class1[i].values():
						line += '1;'
					elif len(class2.keys()) in hist_class2[i].values():
						line += '1;'
					else:
						line += '0;'
					
					#Rule number 4: 100% - 1 conserved position within class
					if len(all_seq.keys()) - 1  in hist_all[i].values():
						line += '1;'
					elif len(class2.keys()) - 1 in hist_class2[i].values():
						line += '1;'
					else:
						line += '0;'
		if g2 != 0 and chi != 0:
		
			# Calculation of significant difference between the samples of the two classes.
			aa_list = hist_class1[i].keys()
			aa_list.sort()
			P1 = []
			P2 = []
			P3 = []
			cont_table = []
			sumP1 = 0
			sumP2 = 0
			for aa in aa_list:
				if hist_class1[i][aa] +  hist_class2[i][aa] != 0:
					sumP1 += hist_class1[i][aa]
					sumP2 += hist_class2[i][aa]
					P1.append(hist_class1[i][aa])
					P2.append(hist_class2[i][aa])
					P3.append(hist_class1[i][aa] + hist_class2[i][aa])
	
			P1.append(sumP1)
			P2.append(sumP2)
			P3.append(sumP1 + sumP2)
			cont_table.append(P1)
			cont_table.append(P2)
			cont_table.append(P3)
		

			score = bitk.Xsq2sample(cont_table)
			line += str(score[0]) + ';' + str(score[1]) + ';'
		
		if g2 != 0 and -lda != 0:
			
################ TRANINING -----------------------------------------------------------------#####################
#			
			AA_list = ['A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','Y','W','-']
			#Dimensionality reduction

			allclass_hist = bitk.aa_distr(all_seq, i)
			red_dim_aa_list = []

			for AA_LDA in allclass_hist.keys():
				if allclass_hist[AA_LDA] == 0:
					AA_list.remove(AA_LDA)
#			print 'New list ' + str(AA_list)
					
			#print 'Initiating Linear Discriminant Analysis'


			#Number of member of the smallest group
			N1 = len(class1.keys())
			N2 = len(class2.keys())
			Nmin = min(N1,N2)

			if i == 0:
				print 'Nmin = ' + str(Nmin)
						
			Nsam = 7000 #Number of samples

			#building vectors

			# Matrix of the sample vectors
			class1_X = []
			class2_X = []

			# Mean vectors
			class1_M = numpy.zeros(len(AA_list))
			class2_M = numpy.zeros(len(AA_list))
			#E_c1c2 = numpy.zeros(len(AA_list))

#			'Starting training'

			for nsam in range(Nsam):
				random.seed()
				class1_sample_list = random.sample(class1.keys(), int(Nmin*0.8))
				class2_sample_list = random.sample(class2.keys(), int(Nmin*0.8))
				
#				print class1_sample_list

				class1_selected = {}
				class2_selected = {}

				for seq_name in class1_sample_list:
					class1_selected[seq_name] = class1[seq_name]

#				print class1_selected

				for seq_name in class2_sample_list:
					class2_selected[seq_name] = class2[seq_name]

				class1_X_dic = bitk.aa_distr(class1_selected, i)
				class2_X_dic = bitk.aa_distr(class2_selected, i)

#				print class1_X_dic

				
				class1_Xn = numpy.zeros(len(AA_list))
				class2_Xn = numpy.zeros(len(AA_list))

				for AA_LDA in range(len(AA_list)):
					class1_Xn[AA_LDA] = class1_X_dic[AA_list[AA_LDA]]
					class1_M[AA_LDA] += float(class1_Xn[AA_LDA])/float(Nsam)

					class2_Xn[AA_LDA] = class2_X_dic[AA_list[AA_LDA]]
					class2_M[AA_LDA] += float(class2_Xn[AA_LDA])/float(Nsam)

					#E_c1c2 += float(class1_Xn[AA_LDA])*float(class2_Xn[AA_LDA])/float(Nsam)



				class1_X.append(class1_Xn)
				class2_X.append(class2_Xn)
			
#			for vector in class1_X:
#				print vector
	
			#sys.exit()

#			print 'Sampling complete'

			#Reduce dimensionality avoid singularity on covariance matrix

#			red_dim_list = []
			
#			for AA_LDA in range(len(AA_list)):
##				S = 0
#				for k_lda in range(Nsam):
#					S += class1_X[k_lda][AA_LDA] + class2_X[k_lda][AA_LDA]
#				
#				if S == 0:
#					red_dim_list.append(AA_LDA)

#			print red_dim_list

			#cleanning

#			for k_lda in range(Nsam):
#				class1_X[k_lda] = numpy.array([class1_X[k_lda][x] for x in range(len(AA_list)) if x not in red_dim_list])
#				class2_X[k_lda] = numpy.array([class2_X[k_lda][x] for x in range(len(AA_list)) if x not in red_dim_list])
				

#			class1_M = numpy.array([class1_M[x] for x in range(len(AA_list)) if x not in red_dim_list])
#			class2_M = numpy.array([class2_M[x] for x in range(len(AA_list)) if x not in red_dim_list])

#			AA_list = [AA_list[x] for x in range(len(AA_list)) if x not in red_dim_list]
			
#			print AA_list




#			print class1_M
#			print class2_M

			# Calculating the covariance:

			Cvar1 = numpy.mat(numpy.zeros([len(AA_list), len(AA_list)]))
			Cvar2 = numpy.mat(numpy.zeros([len(AA_list), len(AA_list)]))

			for i_lda in range(len(AA_list)):
				for j_lda in range(len(AA_list)):
					E_LDA1 = 0
					E_LDA2 = 0
					for k_lda in range(Nsam):
						E_LDA1 += class1_X[k_lda][i_lda] * class1_X[k_lda][j_lda] / float(Nsam)
						E_LDA2 += class2_X[k_lda][i_lda] * class2_X[k_lda][j_lda] / float(Nsam)
					
#				print E_LDA1
				#	if i_lda == 20 and j_lda == 3:
				#		print E_LDA1

					Cvar1[i_lda, j_lda] = E_LDA1 - (class1_M[i_lda] * class1_M[j_lda]) 
					Cvar2[i_lda, j_lda] = E_LDA2 - (class2_M[i_lda] * class2_M[j_lda])




					#Cvar1[i_lda, j_lda] += (class1_X[k_lda][i_lda] - class1_M[i_lda]) , (class1_X[k_lda][j_lda] - class1_M[j_lda])) / float(Nsam)
					#Cvar2[i_lda, j_lda] += numpy.dot((class2_X[k_lda][i_lda] - class2_M[i_lda]) , (class2_X[k_lda][j_lda] - class2_M[j_lda]) / float(Nsam)

			#float(len(AA_list)*len(AA_list)

			# Determining the dicriminant vector:

#			print Cvar1 + Cvar2


			#W_LDA = numpy.dot((Cvar1 + Cvar2).getI(), (class1_M - class2_M))

			try:
				W_LDA = numpy.array(numpy.dot((Cvar1 + Cvar2).getI(), (class1_M - class2_M)))
			except numpy.linalg.linalg.LinAlgError:
				print 'Cvar1+2 failed'
				W_LDA = numpy.zeros(len(AA_list))

			# Calculating the signal-to-noise ratio

			#S_LDA = numpy.dot((Cvar1 + Cvar2), W_LDA[:,numpy.newaxis])
			S_LDA = ((numpy.dot(W_LDA, (class1_M - class2_M))) ** 2) /float(numpy.dot(W_LDA, numpy.dot((Cvar1 + Cvar2), W_LDA.T)))

			line += str(S_LDA[0]) + ';'
#			print 'Mean vector class1: ' + str(class1_M)
#			print 'Mean vector class2: ' + str(class2_M)
			print 'Signal to noise ratio: ' + str(S_LDA)
#			print 'Discriminant Vector: '
#			print W_LDA
#
#			print 'W.m1 =' + str( numpy.dot(W_LDA, class1_M))
#			print 'W.m2 =' + str( numpy.dot(W_LDA, class2_M))

#			print 'Distriminant determined\n Tranining complete\n Initiating test'
		
################### Testing ###########################################################################

			Ntes = 30000 # number of tests
			co_lda = 0.999

			class1p = 0
			class1n = 0
			class2p = 0
			class2n = 0

			for ntes in range(Ntes):
				random.seed()
				class1_test_list = random.sample(class1.keys(), int(Nmin*0.2))
				class2_test_list = random.sample(class2.keys(), int(Nmin*0.2))

				class1_selected = {}
				class2_selected = {}
				
				for seq_name in class1_test_list:
					class1_selected[seq_name] = class1[seq_name]
				for seq_name in class2_test_list:
					class2_selected[seq_name] = class2[seq_name]
				
				class1_X_dic = bitk.aa_distr(class1_selected, i)
				class2_X_dic = bitk.aa_distr(class2_selected, i)
				class1_Xn = numpy.zeros(len(AA_list))
				class2_Xn = numpy.zeros(len(AA_list))

				for AA_LDA in range(len(AA_list)):
					class1_Xn[AA_LDA] = class1_X_dic[AA_list[AA_LDA]]
					class2_Xn[AA_LDA] = class2_X_dic[AA_list[AA_LDA]]

				# projection
#				print '\nTest class 1' + str(class1_Xn)
#				print 'Test class 2' + str(class2_Xn)

#				print 'W.X1_' +  str(ntes) + ' = ' + str (numpy.dot(W_LDA, class1_Xn))

				


				if numpy.dot(W_LDA, class1_Xn) > 0:
					class1p += 1
				else:
					class1n += 1
				
				if numpy.dot(W_LDA, class2_Xn) > 0:
					class2p += 1
				else:
					class2n += 1

				if (class1p > Ntes*(1-co_lda) and class1n > Ntes*(1-co_lda)) or (class2p > Ntes*(1-co_lda) and class2n > Ntes*(1-co_lda)):
					print str(i) + ' no good after ' + str(ntes)
					break
					

			# stats with results

#			print 'class1: + ' + str(class1p) + '\t - ' + str(class1n)
#			print 'class2: + ' + str(class2p) + '\t - ' + str(class2n)


			if abs(class1p + class2n) >= 2*Ntes*co_lda or abs(class2p + class1n) >= 2*Ntes*co_lda:
				line += '1;'
				print '\nPosition ' + str(i) + ' is class specific (LDA)\n'
			else:
				line += '0;'
			



#			class1_SCR = class1_+ - class1_-
#			class1_sig = class1_SCR/abs(class1_+ - class1_-)
#			class2_SCR = class2_+ - class2_-
#			class1_sig = class2_SCR/abs(class2_+ - class2_-)
#			
#			if class1_sig == class2_sig:
#				output += '0;'
#			elif class1_sig > 0:
#				if class1_+ + class2_- <
#				
##				
#				
#				max(class1_+, class1_-) < Ntes*0.95 or max(class2_+, class2_-) < Ntes*0.95:
#				output += '0.5;'
#			else:
#				output += '1;'

			
			
		#else:
		#	line += '0;'

		res_sasa_old = 0

		if '-sasa' in sys.argv:
			res_sasa = j_sasa + bias - 1
			for chain in sasa.keys():
				if str(res_sasa) in sasa[chain].keys() and res_sasa != res_sasa_old:
					line += str(sasa[chain][str(res_sasa)]) + ';'
				else:
					line += '-;'
			res_sasa_old = res_sasa
		
		res_stride_old = 0
		if '-stride' in sys.argv:
                        res_stride = j_stride + bias - 1
                        for chain in stride.keys():
                                if str(res_stride) in stride[chain].keys() and res_stride != res_stride_old:
                                        line += str(stride[chain][str(res_stride)]) + ';'
                                else:
                                        line += '-;'
                        res_stride_old = res_stride


		if '-ns' not in sys.argv:
			line += seq_line


			
		output += line + '\n'




#writing data

#print sasa["A"].keys()

outfile = open('table.csv','w')
outfile.write(output)
outfile.close()






