Tool for fitting background noise of a flow cytometer
Fitting results:
Steps taken
- Determine approximate mean, height, and standard deviation
- Fit Gaussian function with approximate values as starting point
- Determine mean and standard deviation based on fit
- Calculate LoD: mean + 4*std
- pandas
- scipy
- numpy
- matplotlib
- flowio
import asyncio
from js import document, Uint8Array
from pyodide import create_proxy
import io
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import flowio
def gaussian(x, amplitude, mean, stddev):
return amplitude * np.exp(-(x - mean)**2 / (2 * stddev**2))
def fit(noise, factor=1):
bins_noise = noise.max() - noise.min()
bins, bin_edges = np.histogram(noise, bins=int(bins_noise / factor))
initial_amplitude = np.max(bins)
initial_std = np.std(noise)
initial_mean = np.mean(noise)
initial_guess = [initial_amplitude, initial_mean, initial_std]
params, covariance = curve_fit(gaussian, bin_edges[:-1], bins, p0=initial_guess)
amplitude_fit, mean_fit, stddev_fit = params
return bins, bin_edges, amplitude_fit, mean_fit, stddev_fit
async def process_file(event=None):
global events # Declare events as a global variable to be accessible in both functions
fileList = document.getElementById("myfile").files.to_py()
if not fileList:
return
bins_noise_factor = int(document.getElementById("bins_noise_input").value)
for f in fileList:
try:
array_buf = await f.arrayBuffer()
file_bytes = Uint8Array.new(array_buf).to_py()
file_content = io.BytesIO(file_bytes)
flow_data = flowio.FlowData(file_content, ignore_offset_error=True)
# Display FCS header
# document.getElementById("fcs_header").innerHTML = f"FCS Header: {flow_data.channels}"
# Get list of channels
channels = flow_data.channels
channel_select = document.getElementById("channel_select")
channel_select.innerHTML = "" # Clear existing options
for channel_number, channel_name in channels.items():
option = document.createElement("option")
option.value = str(int(channel_number) - 1)
option.text = channel_name['PnN']
channel_select.appendChild(option)
events = np.reshape(flow_data.events, (-1, flow_data.channel_count))
# Add event listener for channel selection
channel_select.addEventListener("change", create_proxy(lambda e: analyze_channel(int(e.target.value), bins_noise_factor)))
# Add event listener for bins noise factor
bins_input = document.getElementById("bins_noise_input")
bins_input.addEventListener("input", create_proxy(lambda e: analyze_channel(int(document.getElementById("channel_select").value), int(e.target.value))))
except Exception as e:
document.getElementById("fcs_header").innerHTML = f"Error: {e}"
def analyze_channel(channel_index, bins_noise_factor):
try:
if 'events' not in globals():
document.getElementById("print_output").innerHTML = "No FCS file loaded."
return
data = events[:, channel_index]
# Filter data
noise = data[(data > 1) & (data < 100000)]
# Fit Gaussian
bins, bin_edges, amplitude_fit, mean_fit, stddev_fit = fit(noise, factor=bins_noise_factor)
y_fit = gaussian(bin_edges[:-1], amplitude_fit, mean_fit, stddev_fit)
lod = mean_fit + (4 * stddev_fit)
document.getElementById("mean").innerHTML = f"mean (a.u.): {mean_fit}"
document.getElementById("std").innerHTML = f"std (a.u.): {stddev_fit}"
document.getElementById("lod").innerHTML = f"LOD (a.u.): {lod}"
plt.figure()
plt.xlabel(f"Channel (a.u.)")
plt.ylabel('Counts')
plt.xlim(np.min(noise), mean_fit + (10 * stddev_fit))
plt.plot(bin_edges[:-1], bins, label='Data')
plt.plot(bin_edges[:-1], y_fit, label='Gaussian fit')
plt.axvline(lod, c='red', linestyle='--', label='lod')
plt.legend()
from pyodide.http import pyfetch
import js
from io import BytesIO
import base64
buf = BytesIO()
plt.savefig(buf, format='png')
plt.close()
buf.seek(0)
img_str = base64.b64encode(buf.read()).decode("utf-8")
img_element = document.createElement("img")
img_element.src = "data:image/png;base64," + img_str
document.getElementById("graph-area").innerHTML = ""
document.getElementById("graph-area").appendChild(img_element)
except Exception as e:
document.getElementById("print_output").innerHTML = f"Error: {e}"
def main():
file_event = create_proxy(process_file)
e = document.getElementById("myfile")
e.addEventListener("change", file_event, False)
main()