Difference between revisions of "Python:DNA melting data analysis"
From Course Wiki
(create page) |
(added code so far) |
||
Line 1: | Line 1: | ||
− | + | I'm working on this. It is way not finished. | |
+ | |||
+ | <pre> | ||
+ | import numpy | ||
+ | import scipy | ||
+ | import scipy.signal as sig | ||
+ | import matplotlib.pyplot as plt | ||
+ | from math import pi | ||
+ | |||
+ | def loadData(expt): | ||
+ | """ | ||
+ | Loads a raw data file into python. The data file is actually not a | ||
+ | Matlab .m file, just two columns of numbers. (First column is RTD voltages | ||
+ | and second column is fluorescence signal voltages.) | ||
+ | """ | ||
+ | f = open(expt.filename) | ||
+ | rtdList = [] | ||
+ | fluorList = [] | ||
+ | for line in f: | ||
+ | rtd, fluor = line.split() | ||
+ | rtdList.append(rtd) | ||
+ | fluorList.append(fluor) | ||
+ | expt.vRTD = numpy.array(rtdList) | ||
+ | expt.fluorSignal = numpy.array(fluorList) | ||
+ | f.close() | ||
+ | |||
+ | def plotFilter(f): | ||
+ | """ | ||
+ | Takes a filter designed with scipy.signal.firwin() and plots | ||
+ | its frequency response. This is more complicated in Python than in | ||
+ | Matlab because you have to tell scipy and matplotlib how to talk to each | ||
+ | other--it's not integrated like in matlab--but it works the same. | ||
+ | """ | ||
+ | # BUT WHY IS THE MAGNITUDE DIFFERENT | ||
+ | # I AM CONFUSED | ||
+ | |||
+ | w, h = sig.freqz(f) | ||
+ | w = w / pi # normalize frequency | ||
+ | fig = plt.figure() | ||
+ | plt.title("Digital filter frequency response") | ||
+ | ax1 = fig.add_subplot(111) | ||
+ | plt.semilogy(w, numpy.abs(h), 'b') | ||
+ | plt.xlabel("Normalized frequency (* pi rad/sample)") | ||
+ | plt.ylabel("Amplitude (dB)", color='b') | ||
+ | plt.grid() | ||
+ | ax2=ax1.twinx() | ||
+ | angles = numpy.unwrap(numpy.angle(h)) | ||
+ | plt.plot(w, angles, 'g') | ||
+ | plt.ylabel("Angle (rad)", color='g') | ||
+ | plt.show() | ||
+ | |||
+ | def convertVoltages(expt): | ||
+ | """ | ||
+ | Accepts an Experiment object with RTD and Fluor voltages and converts them | ||
+ | to temperature and relative fluorescence values. | ||
+ | """ | ||
+ | # Convert RTD voltage to temperature. | ||
+ | # V_rtd = 15 * R_rtd/(R + R_rtd) | ||
+ | # --> R_rtd = R * (V_rtd - 15)/15 | ||
+ | # | ||
+ | # R_rtd = 1000 + 3.85(T-273) | ||
+ | # --> T = 273 + (R_rtd - 1000)/3.85 | ||
+ | |||
+ | expt.rRTD = expt.rOther * (expt.vRTD - 15.0) / 15.0 | ||
+ | expt.temp = 273.15 + (expt.rRTD - 1000.0)/3.85 | ||
+ | |||
+ | # Convert fluorescence signal voltage to relative fluorescence value. | ||
+ | # (relative fluorescence is the same as dsDNA fraction) | ||
+ | |||
+ | # Normalize to be between 1 and 0 after denoising | ||
+ | # NOT FINISHED | ||
+ | </pre> |
Latest revision as of 05:02, 13 April 2011
I'm working on this. It is way not finished.
import numpy import scipy import scipy.signal as sig import matplotlib.pyplot as plt from math import pi def loadData(expt): """ Loads a raw data file into python. The data file is actually not a Matlab .m file, just two columns of numbers. (First column is RTD voltages and second column is fluorescence signal voltages.) """ f = open(expt.filename) rtdList = [] fluorList = [] for line in f: rtd, fluor = line.split() rtdList.append(rtd) fluorList.append(fluor) expt.vRTD = numpy.array(rtdList) expt.fluorSignal = numpy.array(fluorList) f.close() def plotFilter(f): """ Takes a filter designed with scipy.signal.firwin() and plots its frequency response. This is more complicated in Python than in Matlab because you have to tell scipy and matplotlib how to talk to each other--it's not integrated like in matlab--but it works the same. """ # BUT WHY IS THE MAGNITUDE DIFFERENT # I AM CONFUSED w, h = sig.freqz(f) w = w / pi # normalize frequency fig = plt.figure() plt.title("Digital filter frequency response") ax1 = fig.add_subplot(111) plt.semilogy(w, numpy.abs(h), 'b') plt.xlabel("Normalized frequency (* pi rad/sample)") plt.ylabel("Amplitude (dB)", color='b') plt.grid() ax2=ax1.twinx() angles = numpy.unwrap(numpy.angle(h)) plt.plot(w, angles, 'g') plt.ylabel("Angle (rad)", color='g') plt.show() def convertVoltages(expt): """ Accepts an Experiment object with RTD and Fluor voltages and converts them to temperature and relative fluorescence values. """ # Convert RTD voltage to temperature. # V_rtd = 15 * R_rtd/(R + R_rtd) # --> R_rtd = R * (V_rtd - 15)/15 # # R_rtd = 1000 + 3.85(T-273) # --> T = 273 + (R_rtd - 1000)/3.85 expt.rRTD = expt.rOther * (expt.vRTD - 15.0) / 15.0 expt.temp = 273.15 + (expt.rRTD - 1000.0)/3.85 # Convert fluorescence signal voltage to relative fluorescence value. # (relative fluorescence is the same as dsDNA fraction) # Normalize to be between 1 and 0 after denoising # NOT FINISHED