paper-rank-cells.py

#!/usr/bin/env python2
from __future__ import print_function
import sys
sys.path.append('../lib')
import os
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

isNorm = True

import protocols
import model_ikr as m

savepath = './figs'
if not os.path.isdir(savepath):
    os.makedirs(savepath)

#
# Protocols
#
protocol_funcs = {
    'staircaseramp': protocols.leak_staircase,
    'pharma': protocols.pharma,  # during drug application
    'apab': 'protocol-apab.csv',
    'apabv3': 'protocol-apabv3.csv',
    'ap05hz': 'protocol-ap05hz.csv',
    'ap1hz': 'protocol-ap1hz.csv',
    'ap2hz': 'protocol-ap2hz.csv',
    'sactiv': None,
    'sinactiv': None,
}
protocol_dir = '../protocol-time-series'
protocol = 'pharma'

data_dir = '../data'
file_dir = './out'
file_list = [
        'herg25oc1',
        ]
temperatures = np.array([25.0])
temperatures += 273.15  # in K
fit_seed = '542811797'
withfcap = False

#
# Load up all selected cells
#
RANK = {}
for i_temperature, (file_name, temperature) in enumerate(zip(file_list,
    temperatures)):

    # Get selected cells
    files_dir = os.path.realpath(os.path.join(file_dir, file_name))
    searchwfcap = '-fcap' if withfcap else ''
    selectedfile = './manualv2selected-%s.txt' % (file_name)
    selectedwell = []
    with open(selectedfile, 'r') as f:
        for l in f:
            if not l.startswith('#'):
                selectedwell.append(l.split()[0])

    # Model
    protocol_def = protocol_funcs[protocol]
    if type(protocol_def) is str:
        protocol_def = '%s/%s' % (protocol_dir, protocol_def)
    model = m.Model('../mmt-model-files/kylie-2017-IKr.mmt',
                    protocol_def=protocol_def,
                    temperature=temperature,  # K
                    transform=None,
                    useFilterCap=False)  # ignore capacitive spike

    for cell in selectedwell:
        # Fitted parameters
        param_file = '%s/%s-staircaseramp-%s-solution%s-%s.txt' % \
                (files_dir, file_name, cell, searchwfcap, fit_seed)
        obtained_parameters = np.loadtxt(param_file)

        # Data
        data = np.loadtxt('%s/%s-%s-%s.csv' % (data_dir, file_name, protocol,
            cell), delimiter=',', skiprows=1)
        times = np.loadtxt('%s/%s-%s-times.csv' % (data_dir, file_name,
            protocol), delimiter=',', skiprows=1)
        assert(data.shape == times.shape)

        # Simulation
        simulation = model.simulate(obtained_parameters, times)
        if False:
            for _ in range(5):
                assert(all(simulation == 
                    model.simulate(obtained_parameters, times)))
        voltage = model.voltage(times)

        norm = np.max(simulation) if isNorm else 1.

        idx = (np.abs(times - 1.1)).argmin()
        idx_r = 100  # dt=0.2ms -> 20 ms
        rank_value = np.mean(simulation[idx - idx_r:idx + idx_r]) / norm
        RANK[file_name + cell] = rank_value

#
# Rank and output
#
RANKED_CELLS = [key for key, value in sorted(RANK.iteritems(), key=lambda (k,v): (v,k))]

with open('./paper-rank-cells.txt', 'w') as f:
    f.write('# experiment-id + cell-id, ordered using %s\n' % __file__)
    for c in RANKED_CELLS:
        f.write(c + '\n')

'''
# Can use the following for colouring
import seaborn as sns
# color_list = sns.cubehelix_palette(len(SORTED_CELLS))
color_list = sns.color_palette('Blues', n_colors=len(SORTED_CELLS))
color_list.as_hex()
'''

print('Done')