In [1]:
% matplotlib inline
In [2]:
from __future__ import print_function
import os.path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from enrich2.variant import WILD_TYPE_VARIANT
import enrich2.plots as enrich_plot
pd.set_option("display.max_rows", 10) # rows shown when pretty-printing
Modify the results_path variable in the next cell to match the output directory of your Enrich2-Example dataset.
In [3]:
results_path = "/path/to/Enrich2-Example/Results/"
Open the Selection HDF5 file with the variants we are interested in.
In [5]:
my_store = pd.HDFStore(os.path.join(results_path, "Rep1_sel.h5"))
The pd.HDFStore.keys() method returns a list of all the tables in this HDF5 file.
In [6]:
my_store.keys()
Out[6]:
We will work with the "/main/variants/counts" table first. Enrich2 names the columns for counts c_n where n is the time point, beginning with 0 for the input library.
We can use a query to extract the subset of variants in the table that exceed the specified cutoff. Since we're only interested in variants, we'll explicitly exclude the wild type. We will store the data we extract in the variant_count data frame.
In [7]:
read_cutoff = 10
In [8]:
variant_counts = my_store.select('/main/variants/counts', where='c_0 > read_cutoff and index != WILD_TYPE_VARIANT')
variant_counts
Out[8]:
The index of the data frame is the list of variants that exceeded the cutoff.
In [9]:
variant_counts.index
Out[9]:
We can use this index to get the scores for these variants by querying the "/main/variants/scores" table. We'll store the result of the query in a new data frame named variant_scores, and keep only the score and standard error (SE) columns.
In [10]:
variant_scores = my_store.select('/main/variants/scores', where='index in variant_counts.index')
variant_scores = variant_scores[['score', 'SE']]
variant_scores
Out[10]:
Now that we're finished getting data out of the HDF5 file, we'll close it.
In [11]:
my_store.close()
To more easily explore the relationship between input count and score, we'll add a column to the variant_scores data frame that contains input counts from the variant_counts data frame.
In [12]:
variant_scores['input_count'] = variant_counts['c_0']
variant_scores
Out[12]:
Now that all the information is in a single data frame, we can make a plot of score vs. input count. This example uses functions and colors from the Enrich2 plotting library. Taking the log10 of the counts makes the data easier to visualize.
In [13]:
fig, ax = plt.subplots()
enrich_plot.configure_axes(ax, xgrid=True)
ax.plot(np.log10(variant_scores['input_count']),
variant_scores['score'],
linestyle='none', marker='.', alpha=0.6,
color=enrich_plot.plot_colors['bright4'])
ax.set_xlabel("log10(Input Count)")
ax.set_ylabel("Variant Score")
Out[13]:
In [ ]: