Merge branch 'testing' of https://github.com/genenetwork/genenetwork2 into testing

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()

+

+