1 files changed, 254 insertions, 0 deletions
diff --git a/scripts/maintenance/readProbeSetSE_v7.py b/scripts/maintenance/readProbeSetSE_v7.py
new file mode 100755
index 00000000..fd6f0bb8
--- /dev/null
+++ b/scripts/maintenance/readProbeSetSE_v7.py
@@ -0,0 +1,254 @@
+#!/usr/bin/python2
+"""This script use the nearest marker to the transcript as control, increasing permutation rounds according to the p-value"""
+########################################################################
+# Last Updated Sep 27, 2011 by Xiaodong
+# This version fix the bug that incorrectly exclude the first 2 probesetIDs
+########################################################################
+
+import string
+import sys
+import MySQLdb
+import getpass
+import time
+
+
+def translateAlias(str):
+	if str == "B6":
+		return "C57BL/6J"
+	elif str == "D2":
+		return "DBA/2J"
+	else:
+		return str
+
+########################################################################
+#
+#  Indicate Data Start Position, ProbeFreezeId, GeneChipId, DataFile
+#
+########################################################################
+
+dataStart = 1
+
+GeneChipId = int( raw_input("Enter GeneChipId:") )
+ProbeSetFreezeId = int( raw_input("Enter ProbeSetFreezeId:") )
+input_file_name = raw_input("Enter file name with suffix:")
+
+fp = open("%s" % input_file_name, 'rb')
+
+
+try:
+	passwd = getpass.getpass('Please enter mysql password here : ')
+        con = MySQLdb.Connect(db='db_webqtl',host='localhost', user='username',passwd=passwd)
+
+	db = con.cursor()
+        print "You have successfully connected to mysql.\n"
+except:
+	print "You entered incorrect password.\n"
+        sys.exit(0)
+	
+time0 = time.time()
+########################################################################
+#
+#  Indicate Data Start Position, ProbeFreezeId, GeneChipId, DataFile
+#
+########################################################################
+
+#GeneChipId = 4 
+#dataStart = 1 	
+#ProbeSetFreezeId = 359 #JAX Liver 6C Affy M430 2.0 (Jul11) MDP
+#fp = open("GSE10493_AllSamples_6C_Z_AvgSE.txt", 'rb')
+
+
+#########################################################################
+#
+#  Check if each line have same number of members
+#  generate the gene list of expression data here
+#
+#########################################################################
+print 'Checking if each line have same number of members'
+
+GeneList = []
+isCont = 1
+header = fp.readline()
+header = string.split(string.strip(header),'\t')
+header = map(string.strip, header)
+nfield = len(header)
+line = fp.readline()
+
+kj=0
+while line:
+	line2 = string.split(string.strip(line),'\t')
+	line2 = map(string.strip, line2)
+	if len(line2) != nfield:
+		print "Error : " + line
+		isCont = 0
+
+	GeneList.append(line2[0])
+	line = fp.readline()
+
+	kj+=1
+	if kj%100000 == 0:
+		print 'checked ',kj,' lines'
+
+GeneList = map(string.lower, GeneList)
+GeneList.sort()
+	
+if isCont==0:
+	sys.exit(0)
+
+
+print 'used ',time.time()-time0,' seconds'
+#########################################################################
+#
+#  Check if each strain exist in database
+#  generate the string id list of expression data here
+#
+#########################################################################
+print 'Checking if each strain exist in database'
+
+isCont = 1
+fp.seek(0)
+header = fp.readline()
+header = string.split(string.strip(header),'\t')
+header = map(string.strip, header)
+header = map(translateAlias, header)
+header = header[dataStart:]
+Ids = []
+for item in header:
+	try:
+		db.execute('select Id from Strain where Name = "%s"' % item)
+		Ids.append(db.fetchall()[0][0])
+	except:
+		print item,'does not exist, check the if the strain name is correct'
+		isCont=0
+
+if isCont==0:
+	sys.exit(0)
+
+
+print 'used ',time.time()-time0,' seconds'
+########################################################################
+#
+# Check if each ProbeSet exist in database
+#
+########################################################################
+print 'Check if each ProbeSet exist in database'
+
+##---- find PID is name or target ----##
+line = fp.readline()
+line = fp.readline()
+line2 = string.split(string.strip(line),'\t')
+line2 = map(string.strip, line2)
+PId = line2[0]
+
+db.execute('select Id from ProbeSet where Name="%s" and ChipId=%d' % (PId, GeneChipId))
+results = db.fetchall()
+IdStr = 'TargetId'
+if len(results)>0:
+	IdStr = 'Name'
+
+
+##---- get Name/TargetId list from database ----##
+db.execute('select distinct(%s) from ProbeSet where ChipId=%d order by %s' % (IdStr, GeneChipId, IdStr))
+results = db.fetchall()
+	
+Names = []
+for item in results:
+	Names.append(item[0])
+Names = map(string.lower, Names)
+Names.sort() # -- Fixed the lower case problem of ProbeSets affx-mur_b2_at  doesn't exist --#
+
+##---- compare genelist with names ----##
+x=y=0
+x1=-1
+GeneList2=[]
+while x<len(GeneList) and y<len(Names):
+	if GeneList[x]==Names[y]:
+		x += 1
+		y += 1
+	elif GeneList[x]<Names[y]:
+		if x!=x1:
+			GeneList2.append(GeneList[x])
+			x1 = x
+		x += 1
+	elif GeneList[x]>Names[y]:
+		y += 1
+
+	if x%100000==0:
+		print 'check Name, checked %d lines'%x
+
+while x<len(GeneList):
+	GeneList2.append(GeneList[x])
+	x += 1
+
+isCont=1
+ferror = open("ProbeSetError.txt", "wb")
+for item in GeneList2:
+	ferror.write(item + " doesn't exist \n")
+	print item, " doesn't exist"
+	isCont = 0
+        
+if isCont==0:
+	sys.exit(0)
+
+
+print 'used ',time.time()-time0,' seconds'
+#############################
+#Insert new Data into SE
+############################
+db.execute("""
+	select ProbeSet.%s, ProbeSetXRef.DataId from ProbeSet, ProbeSetXRef 
+	where ProbeSet.Id=ProbeSetXRef.ProbeSetId and ProbeSetXRef.ProbeSetFreezeId=%d"""
+	 % (IdStr, ProbeSetFreezeId))
+results = db.fetchall()
+
+ProbeNameId = {}
+for Name, Id in results:
+	ProbeNameId[Name] = Id
+
+ferror = open("ProbeError.txt", "wb")
+
+DataValues = []
+
+fp.seek(0) #XZ add this line
+line = fp.readline() #XZ add this line
+line = fp.readline()
+
+kj = 0
+while line:
+	line2 = string.split(string.strip(line),'\t')
+	line2 = map(string.strip, line2)
+
+	CellId = line2[0]
+	if not ProbeNameId.has_key(CellId):
+		ferror.write(CellId + " doesn't exist\n")
+		print CellId, " doesn't exist"
+	else:
+		DataId = ProbeNameId[CellId]
+		datasorig = line2[dataStart:]
+
+		i = 0
+		for item in datasorig:
+			if item != '':
+				value = '('+str(DataId)+','+str(Ids[i])+','+str(item)+')'
+				DataValues.append(value)
+			i += 1
+
+	kj += 1
+	if kj % 100 == 0:
+		Dataitems = ','.join(DataValues)
+		cmd = 'insert ProbeSetSE values %s' % Dataitems
+		db.execute(cmd)
+
+		DataValues = []
+		print 'inserted ',kj,' lines'
+		print 'used ',time.time()-time0,' seconds'
+	line = fp.readline()
+
+if len(DataValues)>0:
+	DataValues = ','.join(DataValues)
+	cmd = 'insert ProbeSetSE values %s' % DataValues
+	db.execute(cmd)
+
+con.close()
+
+