Author: Jonathan Friedman & Eugene Yurtsev Last Updated Date: 2013-September-11
Please read the instructions on the webpage (https://bitbucket.org/gorelab) on how to install these packages. Have fun!
In [80]:
%load_ext autoreload
%autoreload 2
import FlowCytometryTools
from FlowCytometryTools import FCPlate, FCMeasurement
import os
import glob
from GoreUtilities import graph
(Note: The lines above with the % mark instruct python to reload the dependencies if any of them changed.)
In [81]:
datadir = os.path.join(FlowCytometryTools.__path__[0], 'tests', 'data', 'Plate01')
datafile = os.path.join(datadir, 'CFP_Well_A4.fcs')
print 'Loading file from path: {0}'.format(datafile)
sample = FCMeasurement(ID='Test Plate', datafile=datafile)
In [82]:
# You can insert your own directory instead of datadir here.
datadir = os.path.join(FlowCytometryTools.__path__[0], 'tests', 'data', 'Plate01')
print 'Loading files from directory path: {0}'.format(datadir)
# Load the files
plate = FCPlate.from_dir(ID='Test Plate', path=datadir)
In [83]:
plate = FCPlate.from_dir(ID='Test Plate', path=datadir).transform('hlog', channels=('Y2-A', 'B1-A'))
# The line above is equivalent to the two steps below.
# plate = FCPlate.from_dir(ID='Test Plate', path=datadir)
# plate = plate.transform('hlog', channels=('Y2-A', 'B1-A'))
In [84]:
print plate
In [85]:
well = plate['A3']
print well
Before we can do any plotting, we need to figure out what information was collected by the flow cytometer.
In [86]:
print well.channel_names
In [87]:
print well.channels
If you want to understand what the different fields mean, you need to read the specification for the FCS format. To check which FCS format was used do the following:
In [88]:
well.meta['__header__']['FCS format']
Out[88]:
There is more meta data buried inside of each FCS file. This meta data is stored inside of a dictionary called meta.
In [89]:
print type(well.meta)
print well.meta.keys()
In [90]:
print well.meta['$TOT']
In [91]:
print well.meta['$DATE']
You'll find that some basic plotting functionality has already been implemented for you. Let's check it out by first plotting 1d and 2d histograms on single wells, and then we'll move on to plots on an entire 96-well plate.
In [92]:
figure(figsize=(10, 10))
well.plot('B1-A', bins=100);
In [93]:
figure(figsize=(10, 10))
well.plot('B1-A', bins=100, alpha=0.7, color='green', normed=1);
In [18]:
figure(figsize=(10, 10))
well.plot(['B1-A', 'Y2-A']);
This plot function accepts all the same arguments as does the matplotlib pcolormesh function. For example, we can change the default colormap used for plotting.
In [94]:
well.plot(['B1-A', 'Y2-A'], cmap=cm.Oranges, colorbar=False);
Also, instead of plotting histograms, we can choose to plot a scatter plot by setting the argument kind='scatter'. In this case, the plot function will behave like the matplotlib scatter function and pass all additional arguments to scatter.
In [95]:
well.plot(['B1-A', 'Y2-A'], kind='scatter', color='red', s=1, alpha=0.3);
You can check which wells have data loaded in them simply by using the print command
In [96]:
print plate
Many of the wells in the plate above are empty. Use the function dropna() to compactify the plate.
In [97]:
plate = plate.dropna()
print plate
plate.shape
Out[97]:
In [100]:
plate.plot('Y2-A', bins=100, autolabel=False, ylim=(0, 1000), xlim=(-1000, 10000));
In [101]:
figure(figsize=(10, 10))
plate.plot(['B1-A', 'Y2-A']);
In [102]:
figure(figsize=(10, 10))
plate.plot(['B1-A', 'Y2-A'], ids=['A3', 'C7']);
In [103]:
from FlowCytometryTools import ThresholdGate, PolyGate
In [104]:
y2_gate = ThresholdGate(1000.0, 'Y2-A', region='above')
Important: You must install wx for python in order to use this GUI.
In [9]:
plate['A3'].view()
In [29]:
poly_gate = PolyGate([(7307.6923076923067, 5809.523809523811), (3483.5164835164837, 1455.7823129251719), (5901.0989010989015, -3714.2857142857138),
(8230.7692307692305, 367.34693877551035), (8230.7692307692305, 367.34693877551035)], ['HDR-T', 'FSC-A'], region='in')
In [105]:
figure(figsize=(10, 4))
ax1 = subplot(1, 2, 1)
title('Sample with gates drawn but not applied')
plate['A3'].plot(['B1-A', 'Y2-A'], gates=[y2_gate], apply_gates=False);
ax2 = subplot(1, 2, 2)
title('Sample with gates drawn and applied')
plate['A3'].plot(['B1-A', 'Y2-A'], gates=[y2_gate], apply_gates=True);
graph.autoscale_subplots([ax1, ax2])
In [106]:
b1_gate = ThresholdGate(2000.0, 'B1-A', region='above')
plate['A3'].plot(('B1-A','Y2-A'), gates=[y2_gate, b1_gate], apply_gates=False);
In [108]:
figure(figsize=(8, 8))
#suptitle('All gates plotted but not applied')
plate.plot(('B1-A','Y2-A'), gates=[y2_gate, b1_gate], apply_gates=False);
figure(figsize=(8, 8))
#suptitle('B1-A gate applied')
plate.plot(('B1-A','Y2-A'), gates=[b1_gate], apply_gates=True);
figure(figsize=(8, 8))
#suptitle('Y2-A gate applied')
plate.plot(('B1-A','Y2-A'), gates=[y2_gate], apply_gates=True);
In [109]:
original_well = plate['A3']
gated_well = original_well.gate(y2_gate, )
figure(figsize=(4, 10))
ax = subplot(2, 1, 1)
title('Original well')
original_well.plot('Y2-A', bins=300, color='gray', ax=ax)
xlim(-1000, 10000)
ax = subplot(2, 1, 2)
title('Gated well')
gated_well.plot('Y2-A', bins=300, color='y', ax=ax)
ylim(0, 400)
xlim(-1000, 10000)
Out[109]:
To make calculations, you need to write a function that accepts as input a well and returns the desired calculation. As an example, we will create a function that counts the total number of events.
In [110]:
well = plate['A3']
To access the underlying data in the well, use the get_data function. The returned data is a pandas data frame. That's right. You have the power of pandas at your hands!
In [111]:
data = well.get_data()
print data
In [112]:
data['Y2-A']
Out[112]:
In [113]:
data['Y2-A'].describe()
Out[113]:
In [114]:
count = data.shape[0]
print 'The count of events is: {0}'.format(count)
In [115]:
def count_events(well):
data = well.get_data()
count = data.shape[0]
return count
In [116]:
print well.apply(count_events)
In [117]:
plate.apply(count_events)
Out[117]:
In [124]:
def count_in_gate(well):
gate = ThresholdGate(1000.0, 'Y2-A', region='above')
new_gate = well.gate(gate)
return new_gate.get_data().shape[0]
print plate.apply(count_in_gate)
gate = [ThresholdGate(1000.0, 'Y2-A', region='above')]
plate.plot('Y2-A', gates=gate, color='g');