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