diff options
Diffstat (limited to 'scripts/maintenance')
-rwxr-xr-x | scripts/maintenance/readProbeSetSE_v7.py | 512 |
1 files changed, 256 insertions, 256 deletions
diff --git a/scripts/maintenance/readProbeSetSE_v7.py b/scripts/maintenance/readProbeSetSE_v7.py index 2a1d44ff..0b15ce09 100755 --- a/scripts/maintenance/readProbeSetSE_v7.py +++ b/scripts/maintenance/readProbeSetSE_v7.py @@ -1,256 +1,256 @@ -#!/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:
- isCont = 0
- print("Error : " + line)
-
- 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:
- isCont = 0
- print(item, 'does not exist, check the if the strain name is correct')
-
-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")
- isCont = 0
-
- print(item, " doesn't exist")
-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")
- 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 = []
- line = fp.readline()
- print(CellId, " doesn't exist")
- print('inserted ', kj, ' lines')
- print('used ', time.time()-time0, ' seconds')
-
-if len(DataValues) > 0:
- DataValues = ','.join(DataValues)
- cmd = 'insert ProbeSetSE values %s' % DataValues
- db.execute(cmd)
-
-con.close()
+#!/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: + isCont = 0 + print("Error : " + line) + + 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: + isCont = 0 + print(item, 'does not exist, check the if the strain name is correct') + +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") + isCont = 0 + + print(item, " doesn't exist") +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") + 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 = [] + line = fp.readline() + print(CellId, " doesn't exist") + print('inserted ', kj, ' lines') + print('used ', time.time()-time0, ' seconds') + +if len(DataValues) > 0: + DataValues = ','.join(DataValues) + cmd = 'insert ProbeSetSE values %s' % DataValues + db.execute(cmd) + +con.close() |