Back again. The script I'm working on for lower limit of detection (LLoD) for laboratory assays ran fine yesterday, but now that I'm trying to direct the output to a .txt, I'm running into all kinds of errors. Have tried the sys.stdout = open("test.txt", "w") with sys.stdout.close() at the end, to no avail (further, this throws a ValueError: I/O operation on closed file. error if I run the script again, even if I removed the sys.stdout stuff). My most recent attempt was to use with open(outName, 'w') as f: print(hitTable, '\n', file=f) for the first block to be printed to a .txt (note: before this block, I have the following:
ts = time.time() dt = datetime.datetime.fromtimestamp(ts).strftime('%Y%m%d') tm = datetime.datetime.fromtimestamp(ts).strftime('%H:%M') outName = dt + " Probit Analysis " + tm +".txt" as I expect that this script will be run multiple times a day, and don't want the output file overwritten for multiple runs)
All the examples I found were very basic, like
with open('file.txt', 'w') as f: print('some text', file=f) but nothing explaining what to do when you want something printed to the file, then a bunch of calculations are done, then some more stuff is to be written, etc., so later I have lines like
if r_2 < 0.9: with open(outName, 'a') as f: print('WARNING: low r-squared value for linear regression. Try re-analyzing without using log10.', file=f) but I'm not sure if that's the right way to handle it.
My desired output would look like
qty log_qty probability probit 0 1.0 0.000000 0.07 3.524209 1 2.0 0.301030 0.15 3.963567 2 3.0 0.477121 0.24 4.293697 3 4.0 0.602060 0.32 4.532301 4 5.0 0.698970 0.40 4.746653 5 10.0 1.000000 0.60 5.253347 6 20.0 1.301030 0.87 6.126391 7 30.0 1.477121 0.90 6.281552 8 40.0 1.602060 0.95 6.644854 9 50.0 1.698970 0.98 7.053749 Linear regression using least-squares method yields: y = 2.062x + 3.353 with a corresponding R-squared value of 0.99066 The lower limit of detection (LLoD) at 95% CI is 39.4563. possibly adding a header at the top along the lines of "Probit analysis YYYY-MM-DD HH:MM:SS", but that's not important atm.
Code follows, sample data available here
import os import sys import math import pandas as pd import numpy as np import matplotlib.pyplot as plt import scipy import datetime from scipy.stats import norm from numpy.polynomial import Polynomial from tkinter import filedialog from tkinter import * # Initialize tkinter root = Tk() root.withdraw() # Function to check whether user-input column exists in data def checkInput(prompt,colCheck): while True: qVar = input(f"Enter the column header for {prompt}: ") if qVar in colCheck.columns: break print(f"Invalid input. Column '{qVar}' does not exist. Try again.") return qVar # Function to get the hitrate/probability of detection for a given quantity # Note: any values in df.result that cannot be parsed as a number will be converted to NaN def hitrate(qty, df): t_s = df[df.qty == qty].result t_s = t_s.apply(pd.to_numeric, args=('coerce',)).isna() return (len(t_s)-t_s.sum())/len(t_s) # Main function for generating plot and calculating LLoD-95 def regPlot(x, y, log): ts = time.time() dt = datetime.datetime.fromtimestamp(ts).strftime('%Y%m%d') tm = datetime.datetime.fromtimestamp(ts).strftime('%H:%M') params = {'mathtext.default': 'regular'} plt.rcParams.update(params) y95 = 6.6448536269514722 # probit score corresponding to 95% CI regFun = lambda m, x, b : (m*x) + b # equation for line regression = scipy.stats.linregress(x,y) r_2 = regression.rvalue*regression.rvalue #calculate R^2 # Solve linear equation for x at 95% CI log_llod = (y95 - regression.intercept) / regression.slope xmax = log_llod * 1.2 # Start plotting all the things! fig = plt.figure() ax = fig.add_subplot(111) ax.set_title('Limit of Detection') ax.set_ylabel('Probit score\n$(\sigma + 5)$') while True: if log == 'y': ax.set_xlabel('$log_{10}$(input quantity)') break elif log == 'n': ax.set_xlabel('input quantity') break raise ValueError('Error when calling regPlot(x,y,log) - User input invalid.') x_r = [0, xmax] y_r = [regression.intercept, regFun(regression.slope,x_r[1],regression.intercept)] ax.plot(x_r, y_r, '--k') # linear regression ax.plot(log_llod, y95, color='red', marker='o', markersize=8) # LLOD point ax.plot([0,xmax], [y95,y95], color='red', linestyle=':') # horiz. red line ax.plot([log_llod,log_llod], [regFun(regression.slope,x_r[0],regression.intercept),7.1], color='red', linestyle=':') # vert. red line ax.plot(x, y, 'bx') # actual (qty, probit) data points ax.grid() # grid ax.set_xlim(left=0) plt.ion() plt.show() plt.savefig(f"{dt} Probit Analysis at {tm}.png") plt.pause(0.001) return regression.slope, regression.intercept, r_2, regression.stderr, regression.intercept_stderr, log_llod def printResult(m,b,r2): print('\n Linear regression using least-squares method yields:\n') print('\t\ty = '+str("%.3f"%regression.slope)+'x + '+str("%.3f"%regression.intercept)+'\n') print('\twith a corresponding R-squared value of', str("%.5f"%r_2)+"\n") # Prompt user for data file and column headers, ask whether to use log(qty) print("In the directory prompt, select the .csv file containing data for analysis") path = filedialog.askopenfilename() #data = pd.read_csv(path, usecols=[conc, detect]) data = pd.read_csv(path) conc = checkInput('concentration/copy number', data) detect = checkInput('target detection', data) while True: logans = input("Analyze using log10(concentration/number of copies)? (y/n): ") if logans == 'y' or logans == 'n': break print('Invalid input. Please enter either y or n.') # Get timestamps for file names ts = time.time() dt = datetime.datetime.fromtimestamp(ts).strftime('%Y%m%d') tm = datetime.datetime.fromtimestamp(ts).strftime('%H:%M') outName = dt + " Probit Analysis " + tm +".txt" # Read the columns of data specified by the user and rename them for consistency data = data.rename(columns={conc:"qty", detect:"result"}) # Create list of unique qty values and initialize corresponding # log_10, probability, and probit vectors qtys = data['qty'].unique() log_qtys = [0] * len(qtys) prop = [0] * len(qtys) probit = [0] * len(qtys) # Calculate log10(qty), detection probability, and associated probit score for each unique qty for idx, val in enumerate(qtys): log_qtys[idx] = math.log10(val) prop[idx] = hitrate(val, data) probit[idx] = 5 + norm.ppf(prop[idx]) # Create a dataframe of quantaties, log(qtys), probabilities, and probit scores. hitTable = pd.DataFrame(np.vstack([qtys,log_qtys,prop,probit]).T, columns=['qty','log_qty','probability','probit']) # Convert infinite probit values to NaN and drop those rows. hitTable.probit.replace([np.inf,-np.inf],np.nan, inplace=True) hitTable.dropna(inplace=True) with open(outName, 'w') as f: print(hitTable, '\n', file=f) while True: if logans == 'y': m, b, r_2, stderr, int_stderr, log_llod = regPlot(hitTable.log_qty, hitTable.probit, logans) printResult(m,b,r_2) llod_95 = 10**log_llod if r_2 < 0.9: with open(outName, 'a') as f: print('WARNING: low r-squared value for linear regression. Try re-analyzing without using log10.', file=f) break elif logans == 'n': m, b, r_2, stderr, int_stderr, log_llod = regPlot(hitTable.qty, hitTable.probit, logans) printResult(m,b,r_2) llod_95 = log_llod if r_2 < 0.9: with open(outName, 'a') as f: print('WARNING: low r-squared value for linear regression. Try re-analyzing using log10.', file=f) break raise ValueError('Error when attempting to evaluate llod_95 - User input invalid.') with open(outName, 'a') as f: print("\nThe lower limit of detection (LLoD) at 95% CI is " + str("%.4f"%llod_95) + ".\n", file=f) #sys.stdout.close() EDIT: If possible, it'd be helpful to have the output written to a txt as well as printed in the CLI, so I don't have to double check the output every time during testing to make sure things are showing correctly. Also edited for formatting.
