osprai/trunk/converter.py

"""
Converter: Perform the datafile reading, pre-processing and writing.
Yuhang Wan, Christopher Lausted
Last modified on 100427 (yymmdd) by YW

Typical Pipeline:
>import converter as cv
>dat = cv.datafile_read("example.txt") 
>dat = cv.outlier_remove(dat)
>dat = cv.calibrate(dat)
>dat = cv.background_subtract(dat)
>dat = cv.curve_split(dat, 100, 800, 200)
>dat = cv.data_output_clamp(dat, "/dir", 1)
>dat = cv.data_output_biosensor(dat, "/dir", 1,2,3,4)
>dat = cv.obj_save(dat, "saved_dataobj")
"""
__version__ = "100427"

## Import libraries.
import numpy as np
import pylab as plt
##import sys
import xlrd
import time
import os
import copy
import pickle
import SPRdataclass as spr

## Global constants and variables.
FMT_UNKNOWN = 0
FMT_SPRIT = 1        ## Old Plexera SPRit format.
FMT_ICMTXT = 2       ## Plexera Instrument Control Module text.
FMT_CLAMP34 = 3      ## Morton & Myszka CLAMP Version 3.4+
FMT_DAMXLS = 4       ## Plexera Data Analysis Module Excel export.
FMT_DATAOBJ = 5      ## Packed spr Data Object

def check_format_input(fname):
    """Examine the format of an input .txt data file."""
    ## Modified 100407 CGL.
    
    ## Open text file. Exit if file not found.
    try:
        fp = open(fname, "r")
    except IOError:
        print "Error: Unable to open file %s for reading." % fname
        return
    
    ## Examine first line of text file to guess format.
    fmt = FMT_UNKNOWN   
    sprit_line1="Elapsed Time (Seconds)\tAverage Intensity (Pixel Intensity)" 
    clamp_line1="Vers 3.4"
    line1 = fp.readline()
    if (sprit_line1 in line1):
        fmt = FMT_SPRIT      
    elif (line1[2]=="/" and line1[5]=="/"):
        ## This format has date mm/dd/yyy in first column.
        fmt=FMT_ICMTXT
    elif (clamp_line1 in line1):
        fmt = FMT_CLAMP34
    
    fp.close()
    return fmt, line1


def sprit_format(fname):
    """Read the SPRit formatted data and reshape it.
    Return the data in the form of list.
    """
    # open a file, if file doesn't exit, then exit
    try:
        fp = open(fname, "r")
    except IOError:
        print 'Cannot open file %s for reading' % fname
        return
    status = {'DatainputType':'SPRit data'}
    # skip the first line of the input file
    Tmpstr = fp.readline()
    # check if the second line is BEGIN
    TmpStr = fp.readline()
    if ("BEGIN" not in TmpStr):
        print "Warning: Second line of data file is not 'BEGIN'"
        ##sys.exit(0)
    # count the lines of each data spot
    TmpStr = fp.readline() # skip first "0.000000e+000"        
    num_line = 1 # so we start to count from 1
        # count until read the "0.000000e+000"
    while True:
        TmpStr = fp.readline()
        if TmpStr.split('\t')[0] == "0.000000e+000":
            break
        else:
            num_line = num_line + 1
    # count the spots in each file
    num_spot = 0
        # skip the first two lines
    fp.seek(0)    
    for i in range(0,2):
        TmpStr = fp.readline()
    text = fp.readlines()
    for i in text:
        num_spot += i.count('0.000000e+000')
    fp.close()
    #---------------------- Reshaping the data --------------------
    # reshape the data in the string form    num_line*num_spot
    Shaped_str = np.reshape(text,(num_spot,-1))
    tmpshape=np.shape(Shaped_str)
    # reshape the data in the number form:    num_line*(num_spot+1) 
    # only the first column is time: time spotl spot2 ... spot n
    time = np.zeros(num_line)
    Tmpdata = Shaped_str[0,:]
    for n,i in enumerate(Tmpdata):
        time[n]=float(i.split('\t')[0])
    Shaped_data_1 = np.zeros((tmpshape[0]+1,tmpshape[1]))
    Shaped_data_1[0,:]=time
    for n1,i in enumerate(Shaped_str):
        for n2,j in enumerate(i):
            Shaped_data_1[n1+1][n2]=float(Shaped_str[n1][n2].split('\t')[1])
    # reshape the data in the number form:    
    # time1 spotl time2 spot2 ... timen spot n
    Shaped_data_2 = np.zeros((tmpshape[0]*2,tmpshape[1]))
    for i in range(tmpshape[0]):
        Shaped_data_2[i*2]=Shaped_data_1[0]
        Shaped_data_2[i*2+1] = Shaped_data_1[i+1]
        
##    status = {'Datainput Type':'SPRit data'}
    dataformat = 0
    ROInum = num_spot
    dataobj=spr.DataPass(Shaped_data_1, ROInum, dataformat)
    dataobj.status = {}
    dataobj.updateStatus(**status)
    print "There are %d ROIs, and %d lines of data" % (ROInum, num_line)  
    
    return dataobj


def plexera_icm_format(fname):
    """Read the txt data exported from Plexera ICM software and reshape it.
    Return the data in the form of list.
    """
    fp=file(fname,'r')
    Tmpdata,tmptime=[],[]
    status = {'DatainputType':'Plexera ICM data'}
    for eachline in fp:
        Tmpdata.append(eachline.split('\t'))
    num_line,num_spot=np.shape(Tmpdata)
    # initialize the Shaped_data
    Shaped_data = np.zeros((num_spot,num_line))
    num_spot=num_spot-1         # the real number of spot
    Tmpdata=np.transpose(Tmpdata)
    # delete the last '\r\n' of the data in last column of the datafile
    for n,i in enumerate(Tmpdata[-1]):
        Tmpdata[-1][n]=i.rstrip()
    # time converting: convert the first column of the datafile(date and time information) into numerical time
    for i in Tmpdata[0]:
        tmptime_str=i.split(' ')[1]
        timelist=tmptime_str.split(':')
        now=float(timelist[0])*60*60+float(timelist[1])*60+float(timelist[2])
        tmptime.append(now)
    Shaped_data[0]=tmptime
    Shaped_data[0]=Shaped_data[0]-Shaped_data[0][0]
    # assign the signal data
    for i in range(num_spot):
        Shaped_data[i+1]=Tmpdata[i+1]
 
    dataformat = 0
    ROInum = num_spot
    dataobj=spr.DataPass(Shaped_data, ROInum, dataformat)
    dataobj.status = {}
    dataobj.updateStatus(**status)   
    
    print "There are %d ROIs, and %d lines of data" % (ROInum, num_line)  
    return dataobj   


def dam_single_sheet(sh,DataList,ROIinfo,start,end, colsize):
    """Retrieve and shape data form one single sheet from the Spreadsheet.
    Return the data, ROI infomation in the form of list.
    """
    ## and return the minimal spot Id (start_spot)
    ## retrieve the head of the Spreadsheet, SpotList is used to store the spot Id of each sheet
    
    Head = sh.row_values(1)
    SpotList = []
    num_spot = sh.row_values(0).count('Relative Time')  
    
    for i in range(num_spot):
        TmpHead=unicode.encode(Head[colsize*i])
        # now there are 5 columns for each ROI for Raw Intensity
        # and 4 columns for each ROI for Satellite subtracted Intensity
        DataList.append(sh.col_values(2+colsize*i,start,end))
        
        key, val = [], []
        tmpinfo = TmpHead.split('\n')
        for j in tmpinfo:
            key.append(j.split(': ')[0])
            val.append(j.split(': ')[1])
        tmpdic = dict(zip(key,val))
        try:                                    # spot id
            tmpdic['ID'] = int(tmpdic['ID'])    # either in form of "spot1" or "1"
        except ValueError:
            if str.upper(tmpdic['ID']).startswith('SPOT'):
                tmpdic['ID'] = int(tmpdic['ID'][4:])
            else:
                print "Warning: Illegal galfile format."
                print "The ID is neither 'spot1' form nor '1' form. \n", "'ID:'", tmp
        
        tmpdic['Position'] = (int(tmpdic['Block'])+1, int(tmpdic['Row'])+1, int(tmpdic['Column'])+1)
        ROIinfo.append(tmpdic)

    return DataList,ROIinfo
    
def dam_spreadsheet_format(book):
    """Read the spreadsheet exported from Plexera DAM software and reshape it.
    Return the shaped data and ROI information.
    """
    ROIinfo, newROIinfo, n_sheet = [], [], 0
    DataList=[]
    status = {'DatainputType':'Plexera DAM data'}
    sh0=book.sheet_by_index(0)
    if sh0.row_values(0).count('Raw Intensity') == 0 and sh0.row_values(0).count('Subtracted Intensity') > 0:
        print 'This is the satellite subtracted data.'
        colsize = 4
        status ['DAM spreadsheet'] = 'Satellite subtracted Intensity'
    elif sh0.row_values(0).count('Raw Intensity') > 0 and sh0.row_values(0).count('Subtracted Intensity') == 0:
        print 'This is the raw intensity data.'
        colsize = 5
        status ['DAM spreadsheet'] = 'Raw Intensity'
        
    for i in range(2,sh0.nrows):
        if sum(sh0.col_values(2,1,i))>0:
            break
    start,end = i-1,sh0.nrows
    num_line = end - start    
##    print 'start: ', start
    # the time axis initialization and assignment
    Time=sh0.col_values(1,start,end)
    # in case we lose the time axis during the DAM processing
    if sum(Time)==0:    Time=np.arange(num_line)
    DataList.append(Time)
        
    # list DataList to store the Shaped_data of each ROI in each sheet, 
    # list ROIinfo to store the Spot id and the ligand name of each ROI in each sheet
    for i in range(book.nsheets):
        sh=book.sheet_by_index(i)      
        if sh.ncols!=0:        
            DataList,ROIinfo = dam_single_sheet(sh,DataList,ROIinfo,start,end, colsize)
            n_sheet=n_sheet+1
        else:
            break
    # calculate the dim of the Shaped_data, and initialize the matrix        
    num_spot=len(DataList)-1
    num_line=end-start
    Shaped_data=np.zeros((num_spot+1,num_line))    
    Shaped_data[0]=DataList[0]
    
    # in case that the ROIs doesn't count from spot1
    newROIinfo = copy.deepcopy(ROIinfo)
    StartspotList=np.array([i['ID'] for i in ROIinfo])
    start_spot=min(StartspotList)

    for i in range(num_spot):
        tmp = ROIinfo[i]['ID']
        j=tmp-(start_spot-1)             
        Shaped_data[j]=DataList[i+1]
        newROIinfo[j-1]=ROIinfo[i]   
        
    # pack the date and relative information into SPRdata obj
    dataformat = 0
    ROInum = num_spot
    dataobj=spr.DataPass(Shaped_data, ROInum, dataformat, newROIinfo)
    dataobj.status = {}
    dataobj.updateStatus(**status)
    print "There are %d ROIs, and %d lines of data" % (ROInum, num_line)  
    
    return dataobj


def read_clamp_data(fname):
    """Read the Clamp format data 
    Return the data in the form of list.
    """
    fp=open(fname,'r')
    Tmpstr=fp.readline()
    # examine the file head
    # e.g. 'Vers 3.41 cTBA-1 Conc:250 nM\n'
    ROInum, dataformat, Conc, ROIdic = 1, 4, [], {}
    try:
        name = Tmpstr.split(' ')[2].strip()
        ROIdic['ROI name'] = name
        stmp = Tmpstr.split(' ')[3]
        rest = Tmpstr[Tmpstr.find(stmp):].strip()
        ROIdic['concentration'] = rest.split(':')
    except IndexError:
        pass
        
    for Tmpstr in fp:
        try:
            float(Tmpstr.split('\t')[0])
            n_conc=(Tmpstr.count('\t')+1)/2
            ##print Tmpstr, 'the beginning of data'
            break
        except ValueError:
            if Tmpstr.startswith('Conc1'):
                n_conc=Tmpstr.count('\t')
                Conc=Tmpstr.strip().split('\t')[1:]
            elif Tmpstr.startswith('Start1'):
                Start=Tmpstr.strip().split('\t')[1:]
                n_conc=Tmpstr.count('\t')
            elif Tmpstr.startswith('Stop1'):
                Stop=Tmpstr.strip().split('\t')[1:]
                n_conc=Tmpstr.count('\t')
            
        
    # write the data into matrix
    Tmpdata=[]
    for eachline in fp:
        a=eachline.strip().split('\t')
        Tmpdata.append(np.array(a,dtype='f'))
    Data=np.transpose(Tmpdata)
    
    # write the infomation
    ROIinfo = [ROIdic]
    status = {'DatainputType':'Clamp data'}
    sampleinfo, injinfo = [], []    
    injpara = [map(float,i) for i in [Start,Stop]]
    
    if len(Conc) >= 1:
        for i,c in enumerate(Conc):
            sampledic = {}
            name = 'Sample '+str(i+1)
            sampledic = dict(zip(['Name','Concentration'],[name,c]))
##             sampledic['Name'] = name
##             sampledic['Concentration'] = c[injpara[0][i],injpara[1][i]]
            sampleinfo.append(sampledic)
            injdic = dict(zip(['ton','toff'],[float(Start[i]),float(Stop[i])]))
            injinfo.append(injdic)
    dataobj=spr.DataPass(Data, ROInum, dataformat, ROIinfo, sampleinfo, injinfo)
    dataobj.status = {}
    dataobj.updateStatus(**status)
    
    return dataobj
    

def keyfile_read(fkey):
    # The Key File contains 
    """Function: Read the key file for old SPR instrument.
    Input: the filename.
    Return: the ROI information.    
    """
    try:
        fp = open(fkey, "r")
    except IOError:
        print "Cannot open file %s for reading" % fname
        return
    ROIinfo = []
    firstline=fp.readline()     # skip the first line
#    Table.append([firstline.split('\t')[0],firstline.split('\t')[1],firstline.split('\t')[2],firstline.split('\t')[3]])
    for eachline in fp:        
        ROIdic = {}
        tmplist = eachline.strip().split('\t')
        
        try:                                    # spot id
            ROIdic['ID'] = int(tmplist[0])      # either in form of "spot1" or "1"
        except ValueError:
            if str.upper(tmplist[0]).startswith('SPOT'):
                ROIdic['ID'] = int(tmplist[0][4:])
            else:
                print "Error: Unable to read Spot IDs"
                
        ROIdic['Name'] = tmplist[1]     # name of the protein immobilized
        try:
            ROIdic['Conc'] = tmplist[2]
            ROIdic['Background Spot'] = map(int, tmplist[3].split(';')) # in case more than one spot set as background with ";" as SEPARATOR
        except IndexError:
            pass
        ROIinfo.append(ROIdic)
 
    fp.close()    
    print 'There are %d ROIs in the Key file.' % len(ROIinfo)
    return ROIinfo

def keyfile_read_fake(obj):
    ROIinfo = []
    spot_num = obj.ROInum
    for i in range(spot_num):
        ROIdic = {}
        ROIdic['ID'] = i+1
        ROIdic['Name'] = "Ligand"+str(i+1)
        ROIinfo.append(ROIdic)
    return ROIinfo
        

def galfile_read(fname):
    """Function: Read the GAL file.
    Input: the filename.
    Return: the ROI information.    
    """
    fp=file(fname,'r')
    ROIinfo = []
    Tmpstr = fp.readline()      # skip the first line
    Tmpstr = fp.readline()      
    Tmpn = int(Tmpstr.split('\t')[0])
    for i in range(Tmpn):
        labelline = fp.readline()

        if labelline.startswith('"Block1'):
            a = labelline.strip().split(',')
            columnsize, rowsize = int(a[3]), int(a[5])
##            print labelline
    def IDcal(p):
        id = (p[0]-1)*columnsize*rowsize + (p[1]-1)*columnsize + p[2]
        return id
        
    for eachline in fp:
        ROIdic = {}
        tmplist = eachline.strip().split('\t')
        pos = tuple(map(int,tmplist[:3]))
        ROIdic['Position'] = pos
        ROIdic['IDcal'] =IDcal(pos)         
        ##(pos[0]-1)*columnsize*rowsize + (pos[1]-1)*columnsize + pos[2]

        ##
        ROIdic['Name'] = tmplist[3]
##        ROIdic['ID'] = tmplist[4]
        try:                                    # spot id
            ROIdic['ID'] = int(tmplist[4])      # either in form of "spot1" or "1"
        except ValueError:
            if str.upper(tmplist[4]).startswith('SPOT'):
                ROIdic['ID'] = int(tmplist[4][4:])
            else:
                print "Error: Unable to read Spot IDs"
        
        try:
            ROIdic['Conc'] = tmplist[5]
            ROIdic['Family'] = tmplist[6]
            blist = tmplist[7][1:-1].split(';')     # always something like '"1,1,2; 1,2,2"'
            Backn = len(blist)
            Bid = []

            for i in blist:
                if '' in i.split(','):
                    break 
                Bpos = tuple(map(int, i.split(',')))
                Bid.append(IDcal(Bpos))
                ##
            ROIdic['Background Spot'] = Bid     ##tuple(Bid)
        except IndexError:
            pass
        
        if ROIdic['ID'] != ROIdic['IDcal']:
            print 'The spot ID should be consistent with the position. Please check the ID.' 
            break
        else:
            ROIdic['ID'] = ROIdic['IDcal']
        ROIinfo.append(ROIdic)
    fp.close()
    print 'There are %d ROIs in the Gal file.' % len(ROIinfo)
    
    return ROIinfo
    
    
def method_read(fname):
    """Function: Read the analyte table.
    Input: the filename
    Return: the sample information
    """
    # create a concentration dictionary for the sample injected into the flowcell  
    # export a table containing the concentration of the sample injected, and the duration
    fp=file(fname,'r')
    sampleinfo = []
    ##TmpTable,TmpTable2={},{}
    labelline=fp.readline()
    for n,eachline in enumerate(fp):
        sampledic = {}
        sampledic['Location'] = eachline.split('\t')[0]
        sampledic['Name'] = eachline.split('\t')[1]
        sampledic['Concentration'] = float(eachline.split('\t')[2])
        sampledic['Duration'] = float(eachline.split('\t')[4])
        sampledic['Flow Rate'] = float(eachline.split('\t')[3])
        sampledic['Analyte Series'] = eachline.split('\t')[6]
        sampledic['Buffer Blank Series'] = eachline.split('\t')[7]
        sampleinfo.append(sampledic)

    print 'Name','\t','Concentration'
    for i in sampleinfo:
        print i['Name'],'\t', i['Concentration']
    
    return sampleinfo

def method_read_fake():
    sampleinfo = []
    sampleinfo.append({'Location':0, 'Name':0, 'Concentration':0, 'Duration':0, 'Flow Rate':0, 'Analyte Series':0, 'Buffer Blank Series':0})
    return sampleinfo


def datafile_read(fname):
    """Function: Read the raw data from the Plexera instrument, 
        and pack the raw data (Shaped_data) into dataobject. 
        Therefore, the initial dataformat = 0
    Input: the filename
    Return: the packed SPR.DataPass type data object
    """
    #===================check the format of the input files-------
    if fname.split('.')[-1] == 'txt':
    
        Format_input, Firstline = check_format_input(fname)    
        if Format_input == FMT_SPRIT:   # for SPRit file format 
            dataobj = sprit_format(fname)  
            print '-'*10,' This is a SPRit data file. ','-'*10        
            
        elif Format_input == FMT_ICMTXT: # for Plexera ICM file format 
            dataobj = plexera_icm_format(fname)
            print '-'*10,' This is a Plexera ICM exported data file. ','-'*10

        elif Format_input == FMT_CLAMP34:         # for Clamp datafile format
            dataobj = read_clamp_data(fname)
            print '-'*10,' This is a Clamp file. ','-'*10
        
        else:
            print '-'*10,' unreadable file input! ','-'*10
            return
        
        if Format_input == FMT_SPRIT or Format_input == FMT_ICMTXT:
            flag = str.upper(raw_input('Load the key file? (y/n): '))
            if flag == 'Y':
                fkey = raw_input('Input the path of the key file: ')
                ROIinfo = keyfile_read(fkey)
                dataobj.updateROIinfo(ROIinfo)        
            else:
                dataobj.updateROIinfo(keyfile_read_fake(dataobj))   
    
    elif fname.split('.')[-1] == 'xls':
        # for the DAM Spreadsheet file format 
        print '-'*10,' This is a Plexera DAM exported data file. ','-'*10
        Format_input = FMT_DAMXLS
        book = xlrd.open_workbook(fname)
        #shape the data in the Spreadsheet, whether single sheet or multiple sheets 
        dataobj=dam_spreadsheet_format(book)
        flag = str.upper(raw_input('Load the gal file? (y/n): '))
        if flag == 'Y':
            fgal = raw_input('Input the path of the gal file: ')
            ROIinfo = galfile_read(fgal)
            dataobj.updateROIinfo(ROIinfo)
            
    elif fname.split('.')[-1] == 'obj':
        # for the saved data object file format
        print '-'*10,' This is a saved data object file. ','-'*10
        Format_input = FMT_DATAOBJ
        fp = file(fname, 'r+')
        dataobj = pickle.load(fp)   
        fp.close()        
    
    flag = str.upper(raw_input('Load the experimental analyte file? (y/n): '))
    if flag == 'Y':
        fprotocol = raw_input('Input the path of the analyte table: ')
        sampleinfo = method_read(fprotocol)
        dataobj.updateSampleinfo(sampleinfo)
    elif Format_input != FMT_CLAMP34 and Format_input != FMT_DATAOBJ:     
        # clamp data and obj data may contain the experimental protocol already
        dataobj.updateSampleinfo(method_read_fake())

    return dataobj
    

def outlier_remove(obj):
    ##dataobj.viewAll()
    """Function: Remove unwanted data points from noisy periods.
    Return: the data object with clean data
    """
    dataobj = copy.deepcopy(obj)
    print ('Select the time area that you want to remove.')
    tmpstr = raw_input('in the format of "200:240, 500:510" :')
    areatodel = tmpstr.split(',')
    data = dataobj.data
    for i in areatodel:
        ##if 
        t1, t2 = float(i.split(':')[0]), float(i.split(':')[1])
        t = data[0]     # for dataformat=0 Shaped_data
        tmpind = plt.find((t>t1) & (t<t2))  # find the outliar indice
##        data[1:,tmpind] = 0
        tmplist = data.tolist()
        tmplist2 = []
        for j in tmplist:
            del j[min(tmpind):max(tmpind)]
            tmplist2.append(j)
        data = np.array(tmplist2)
##        tmplist = tmplist
    dataobj.data = data
#    dataobj.status = 'Outliar removed'
    dataobj.updateStatus(outlier_remove=True)
    return dataobj

def calibrate(obj):
    """Function: calibrate the Intensity response.
    Return: the data object with calibrated data
    Note: at present, this function is valid only when 
        the whole procedure includes calibration precedure. 
    """
    dataobj = copy.deepcopy(obj)
    data = dataobj.data
    caldata = copy.deepcopy(data)
    t = data[0]
    if dataobj.dataformat != 0:
        print 'We only calibrate the raw data.'
        return
    else:
        t1 = input('Enter the time position of the TOP response during calibration:')
        s1 = input('The refractive index of the TOP response is: ')
        t2 = input('Enter the time position of the BOTTOM response during calibration:')
        s2 = input('The refractive index of the BOTTOM response is: ')
        ind1 = plt.find(abs(t-t1) == min(abs(t-t1)))[0]
        ind2 = plt.find(abs(t-t2) == min(abs(t-t2)))[0]
        
        for i,y in enumerate(data[1:]):
            y1, y2 = y[ind1], y[ind2]
            slope = (s1-s2)/(y1-y2)
            offset = s1- slope*y1
            caldata[i+1] = slope*y+offset
    dataobj.data = caldata    
    dataobj.updateStatus(calibrate=True)
    return dataobj

def get_background_value(obj,bgids):
    """Get the averaged value of background spots for each ROI."""
    bgsignal = obj.data[0]*0
    for j in bgids:
        if j == int(obj.ROIinfo[j-1]['ID']):    # in case the ID in ROIinfo is string type
            bgsignal = bgsignal + obj.data[j]
                
        elif j != int(obj.ROIinfo[j-1]['ID']):
            for n,eachspot in enumerate(obj.ROIinfo):
                if j == int(eachspot['ID']):
                    bgsignal = bgsignal + obj.data[n+1]
                    
    bgsignal = bgsignal/len(bgids)
    return bgsignal
    

def background_subtract(obj, *bgspot):
    """Function: Perform the Background subtraction for the UNSPLIT curve.
    Input besides "obj" is the id of the spot taken as background.
        The following inputs are acceptable:
        1. background_subtract(obj): the default background in Galfile
            will be subtracted.
        2. background_subtract(obj, 1, 6): the average value of spot1 
            and spot6 will be subtracted.
    """
    dataobj = copy.deepcopy(obj)
    ROIinfo = obj.ROIinfo
    data = dataobj.data
        
    if dataobj.dataformat == 0:
        newdata = np.zeros(np.shape(data))
        newdata[0] = data[0]    # the time scale
        if bgspot == ():
        # The average of the default background spots listed in the galfile 
        #are to be subtracted.
            for i in range(1,dataobj.ROInum+1):
                bgids = ROIinfo[i-1]['Background Spot']
                bgsignal = get_background_value(dataobj,bgids)
                newdata[i] = data[i]-bgsignal
            # update the status of the data object.
            dataobj.updateStatus(BackgroundType='Default in Gal')
            print "The background spot embodied in gal is used." 
                
        else:
        # The average of the manually input background spots 
        #are to be subtracted.
            for i in range(1,dataobj.ROInum+1):
                bgsignal = get_background_value(dataobj,bgspot)
                newdata[i] = data[i]-bgsignal
            dataobj.updateStatus(BackgroundType='Manually chosen')
            print "The background spot is manually chosen." 
            
        dataobj.data = newdata
        dataobj.updateStatus(background_subtraction=True)
        return dataobj
        
    else: 
        print 'The Background Subtraction should be run at the beginning, with the UNsplit curve.'    
    return

def curve_split(obj, t_before, t_after, *t_inj):
    """Function: Split the whole curve that contains several injections 
        into pieces.
        The sample information will be checked during the split, 
            if the injection number in sampleinfo is not consistant
            with the input injection parameter, it will help you to 
            update the sample information.
    Input besides obj:
        t_before: time duration before the injection time point
        t_after: time duration after the injection time point
        t_inj: the exact start time for each injection
    Return: the data object with splitted data, and updated sample 
        information
    """
    dataobj = copy.deepcopy(obj)
    ROInum = dataobj.ROInum
    t=dataobj.data[0]
    ind=[]    
    # find how many lines of data each injection 
    for i in range(len(t_inj)):
        ind.append(np.shape(plt.find((t>(t_inj[i]-t_before))&(t<(t_inj[i]+t_after))))[0])
    m,j = max(ind),max(ind)-min(ind)    
    # the new dataset must be a m*((ROInum+1)*num_inj) matrix: 
    # T1,S_ROI1_c1,S_ROI2_c1,S_ROI3_c1...T2,S_ROI1_c2,S_ROI2_c2... 
    tmpdata = np.zeros(((ROInum+1)*len(t_inj),m))
    
    for i in range(len(t_inj)):
        TmpInd = plt.find((t>(t_inj[i]-t_before))&(t<(t_inj[i]+t_after)))
        Tmp1,Tmp2 = i*(ROInum+1),(i+1)*(ROInum+1)
        tmpdata[Tmp1:Tmp2,0:len(TmpInd)] = dataobj.data[:,TmpInd]
    if j!=0:
        Split_data = tmpdata[:,:-j]
    else:
        Split_data = tmpdata[:]
    # zero and align the split data
    for i in range(np.shape(Split_data)[0]):
        Split_data[i]=Split_data[i]-Split_data[i][0]
    
    newsampleinfo, injectinfo = check_sample_info(dataobj,t_before, t_after, t_inj)
    if newsampleinfo != None :        
        dataobj.data = Split_data
        dataobj.updateStatus(DataSplit=True)
        dataobj.updateSampleinfo(newsampleinfo)
        dataobj.dataformat = 2
        dataobj.injectinfo = injectinfo
        return dataobj
    else:
        dataobj.updateStatus(DataSplit=False)
        print 'Curve Splitting not successful!'
    
    return
    
    
def check_sample_info(dataobj, t_before, t_after, t_inj):
    """Check if the sample information consistent with the injection 
        parameter input. If not, help the user to update sample information
        for the splitted curve.
    """
    injectinfo = []
    if len(dataobj.sampleinfo) == len(t_inj):
        for i,j in enumerate(t_inj):                
                injectinfo.append([t_before, t_before+dataobj.sampleinfo[i]['Duration'], t_before+t_after])
        return dataobj.sampleinfo, injectinfo
            
    elif len(dataobj.sampleinfo) > len(t_inj):
        print """    The number of injection you just input does not match that in sample infomation.
    If you just want to split part of the curve, please follow the following direction to update the sampleinfo in the new object."""
        choose = str.upper(raw_input('Continue? (y/n)'))
        if choose == 'Y':
            newsampleinfo, injectinfo = [], []
            print '    Your input indicates %s injections out of %s times in experiment. ' % (len(t_inj), len(dataobj.sampleinfo))
            for n,i in enumerate(dataobj.sampleinfo):
                print '--',n,'.',i['Name'],':',i['Concentration']
            print '    Enter the ID number of the injections you just select, in the form of "2-10" or "4,5,6",'
            tmpstr = raw_input(':')
            if ',' in tmpstr:
                injID = map(int, tmpstr.split(','))
            elif '-' in tmpstr:
                injID = range(int(tmpstr.split('-')[0]), int(tmpstr.split('-')[1])+1)
            else:
                print '    Illegal input!'
                return
                
            if len(injID) != len(t_inj):
                print '    Please input %s injection ID numbers instead of %s You just input.' %(len(t_inj),len(injID))
                return
                
            for i in injID:
                newsampleinfo.append(dataobj.sampleinfo[i])
                injectinfo.append([t_before, t_before+dataobj.sampleinfo[i]['Duration'], t_before+t_after])
            return newsampleinfo, injectinfo
        else:
            return
    
    elif len(dataobj.sampleinfo) < len(t_inj):
        print '    The number of injection you just input does not match that in sample infomation.'
        print '    Your input indicates %s injections, where there is only %s times in experiment. ' % (len(t_inj), len(dataobj.sampleinfo))
        print '    Please check.'
        
        return
    
    
def write_data(Head,data_out,ROI_out,DataDir,output_format):
    """Write the data into a txt file.
    """
    # Modified 100405 CGL
    # create a newfolder with time stamp to save the data 
    if output_format=='B':
        ftitle='biosensor'
        Tmpstr = str(ROI_out).replace(', ','_') 
        Tmpstr = Tmpstr.replace('(','-')
        foutname = ftitle+Tmpstr[:-1]+'.txt'
    elif output_format=='C':
        foutname='clamp-'+str(ROI_out)+'.txt'
    else:
        ftitle=''
        
    # Change directories, if requested.
    if (DataDir != ""):
        os.chdir('..')        
        try:
            os.makedirs(DataDir)
        ##except WindowsError:
        except OSError:
            pass
        os.chdir(DataDir)
    
    np.savetxt('temp.txt',np.transpose(data_out),fmt='%f')
    fdt=file('temp.txt','r')
    
##    foutname=ftitle+Tmpstr[:-1]+'.txt'     ##raw_input('Enter the file name of the output: ')
    fout=file(foutname,'w+')
    # write the head of the data first
    fout.write(Head+'\n')
    count = len(fdt.readlines())    
    fdt.seek(0)
    for i in range(count):
        temp=fdt.readline()    
        fout.write(temp.replace(' ','\t'))
    fdt.close()
    fout.close()
    print '-------The path of the exported data is "', os.getcwd(),'"-------'
    return
    
    
def data_output_biosensor(obj, DataDir, *ROI_out):
    """Export the data into a biosensor format txt file, 
        which can be operated in Scrubber.
    Input: the data object, the output path, the ROIs selected to export
    """
    # initialize the data output matrix, which will contain 'len(spot_out)*2)*n_inj' columns
    dataobj = copy.deepcopy(obj)
    data = dataobj.data
    n_inj = len(dataobj.sampleinfo)
    ROIinfo = dataobj.ROIinfo
    ROInum = dataobj.ROInum
    
    data_out = np.zeros(((len(ROI_out)*2)*n_inj,np.shape(data)[1]))
    head1 = ""
    new_ROIinfo = []
    for n,i in enumerate(ROI_out):
        new_ROIinfo.append(ROIinfo[i-1])        ##the new ROI info
        if i != ROIinfo[i-1]['IDcal']:
            # check if the ID number matches position
            print n,i, ROIinfo[i-1]
            return
        else:
            headcomp = ROIinfo[i-1]['ID']+': '+ROIinfo[i-1]['Name']
            try:
                headcomp = headcomp + ', '+ROIinfo[i-1]['Conc']
            except KeyError:
                pass
            headcomp = headcomp+' Fc='+str(n+1)+' - '
        for j in range(n_inj):
            headtmp = headcomp+str(j+1)+'_'
            head1 = head1+headtmp+'X\t'+headtmp+'Y\t'
            data_out[2*n*n_inj+2*j]=data[j*(ROInum+1)]
            data_out[2*n*n_inj+2*j+1]=data[j*(ROInum+1)+i]
    head1=head1[:-1]
    
    dataobj.data = data_out
    dataobj.dataformat = 3      # biosensor dataformat
    dataobj.updateStatus(DataType='Biosensor data')
    dataobj.updateStatus(Head=head1)    
    dataobj.updateROIinfo(new_ROIinfo)
    
    write_data(head1,data_out,ROI_out,DataDir,'B')
    
    return dataobj
    
def data_output_clamp(obj, DataDir, ROI_out):
    """Export the data into a Clamp format txt file, 
        which can be operated in Clamp.
    Input: the data object, the output path, the selected ROI
    """
    # initialize the data output matrix, which will contain '2*n_inj' columns
    dataobj = copy.deepcopy(obj)
    data = dataobj.data
    n_inj = len(dataobj.sampleinfo)
    n_inj = max(n_inj, 1)  # Must have at least one injection.
    ROIinfo = dataobj.ROIinfo
    ROInum = dataobj.ROInum
    sampleinfo = dataobj.sampleinfo
    
    data_out = np.zeros((2*n_inj,np.shape(data)[1]))
##  if ROI_out != ROIinfo[ROI_out-1]['IDcal']:
    if (ROI_out != int(ROIinfo[ROI_out-1]['ID'])):
        print 'Please check the ROI information.'
        print ROI_out, ROIinfo[ROI_out-1]
        return
    else:
        tmp = ROIinfo[ROI_out-1]['ID']+': '+ROIinfo[ROI_out-1]['Name']
    try:
        tmp = tmp +', Conc:'+ROIinfo[ROI_out-1]['Conc']
    except KeyError:
        pass
    Clamphead = 'Vers 3.41 '+tmp
    conc=''
    t_start,t_stop='',''
    head1='Time\tData\t'*n_inj
    head1=head1[:-1]
    
    for i in range(n_inj):
        data_out[2*i]=data[i*(ROInum+1)]
        data_out[2*i+1]=data[i*(ROInum+1)+ROI_out]
        name = sampleinfo[i]['Name']
        conc=conc+'\t'+str(sampleinfo[i]['Concentration'])
        t_start=t_start+'\t'+str(dataobj.injectinfo[i][0])
        t_stop=t_stop+'\t'+str(dataobj.injectinfo[i][1])
        
        hconc,h_start,h_stop='\nConc1'+conc,'\nStart1'+t_start,'\nStop1'+t_stop
    for i in (hconc,h_start,h_stop):
        if i[-1]!='':
            Clamphead=Clamphead+i
        
    Clamphead=Clamphead+'\n'+head1
    
    dataobj.data = data_out
    dataobj.dataformat = 4      # clamp dataformat
    dataobj.updateStatus(DataType='Clamp data')
    dataobj.updateStatus(Head=Clamphead)    
    dataobj.updateROIinfo(ROIinfo[ROI_out-1])
    
    write_data(Clamphead,data_out,ROI_out,DataDir,'C')
    
    return dataobj
##    return data_out,Clamphead

def obj_save(obj, foutname):
    foutname = foutname + '.obj'        # the file extension of the saved data object is .obj
    fout = file(foutname, 'w+')
    pickle.dump(obj, fout)
    fout.close()
    print '-------The path of the exported obj is "', os.getcwd(),'"-------'
    return

if (__name__ == '__main__'):
    # Print a short message if this module is run like an app.
    print "This is the OSPRAI file type conversion tool."

Revision:	25
Committed:	Wed Apr 28 20:22:24 2010 UTC (9 years, 2 months ago) by rjaynes
File size:	35996 byte(s)
Log Message:	Add py and obj files to allow modeling of more SPR experiments with converter and curvefitting modules. This is the work of Yuhang Wan and Rui Hou. 1. In "converter.py": Add the saving and reading function for the sprclass data object. Also add function "keyfile_read_fake" to provide default information for SPRit and ICM formats in case of the bug when do background_subtract. Fix the bugs in "background_subtract". Tested by DAM and ICM formats. 2. In model modules: "modelclass.py" is the parent class for all the other model classes that performs the theoretical simulating, loading and saving of the parameter or simulated data. Rui and I also add some other model modules like competing model, twostate model, parallel model, and the time variable concentrated models, where the simulated result is compared with Clamp's simulation to make sure the equations are correct. The basicmodel and basicmodel_varyC class are tested. 3. In "curvefitting.py": Add typical pipeline for operation. The examples are packed with the file. Add function to show the Elapsed time for each fitting.