<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Selecting variants by number of unique barcodes
This notebook gets scores for the variants in an Experiment that are linked to multiple barcodes, and plots the relationship between each variant's score and number of unique barcodes.
Step1:  Modify the results_path variable in the next cell to match the output directory of your Enrich2-Example dataset.
Step2:  Open the Experiment HDF5 file.
Step3:  The pd.HDFStore.keys() method returns a list of all the tables in this HDF5 file.
Step4:  First we will work with the barcode-variant map for this analysis, stored in the "/main/barcodemap" table. The index is the barcode and it has a single column for the variant HGVS string.
Step5:  To find out how many unique barcodes are linked to each variant, we'll count the number of times each variant appears in the barcode-variant map using a Counter data structure. We'll then output the top ten variants by number of unique barcodes.
Step6:  Next we'll turn the Counter into a data frame.
Step7:  The data frame has the information we want, but it will be easier to use later if it's indexed by variant rather than row number.
Step8:  We'll use a cutoff to choose variants with a minimum number of unique barcodes, and store this subset in a new index. We'll also exclude the wild type by dropping the first entry of the index.
Step9:  We can use this index to get condition-level scores for these variants by querying the "/main/variants/scores" table. Since we are working with an Experiment HDF5 file, the data frame column names are a MultiIndex with two levels, one for experimental conditions and one for data values (see the pandas documentation for more information).
Step10:  There are fewer rows in multi_bc_scores than in multi_bc_variants because some of the variants were not scored in all replicate selections, and therefore do not have a condition-level score.
Now that we're finished getting data out of the HDF5 file, we'll close it.
Step11:  We'll add a column to the bc_counts data frame that contains scores from the multi_bc_scores data frame. To reference a column in a data frame with a MultiIndex, we need to specify all column levels.
Step12:  Many rows in bc_counts are missing scores (displayed as NaN) because those variants were not in multi_bc_scores. We'll drop them before continuing.
Step13:  Now that we have a data frame containing the subset of variants we're interested in, we can make a plot of score vs. number of unique barcodes. This example uses functions and colors from the Enrich2 plotting library.
<ASSISTANT_TASK:>
Python Code:
% matplotlib inline
from __future__ import print_function
import os.path
from collections import Counter
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
Explanation: Selecting variants by number of unique barcodes
This notebook gets scores for the variants in an Experiment that are linked to multiple barcodes, and plots the relationship between each variant's score and number of unique barcodes.
End of explanation
results_path = "/path/to/Enrich2-Example/Results/"
Explanation: Modify the results_path variable in the next cell to match the output directory of your Enrich2-Example dataset.
End of explanation
my_store = pd.HDFStore(os.path.join(results_path, "BRCA1_Example_exp.h5"))
Explanation: Open the Experiment HDF5 file.
End of explanation
my_store.keys()
Explanation: The pd.HDFStore.keys() method returns a list of all the tables in this HDF5 file.
End of explanation
bcm = my_store['/main/barcodemap']
bcm
Explanation: First we will work with the barcode-variant map for this analysis, stored in the "/main/barcodemap" table. The index is the barcode and it has a single column for the variant HGVS string.
End of explanation
variant_bcs = Counter(bcm['value'])
variant_bcs.most_common(10)
Explanation: To find out how many unique barcodes are linked to each variant, we'll count the number of times each variant appears in the barcode-variant map using a Counter data structure. We'll then output the top ten variants by number of unique barcodes.
End of explanation
bc_counts = pd.DataFrame(variant_bcs.most_common(), columns=['variant', 'barcodes'])
bc_counts
Explanation: Next we'll turn the Counter into a data frame.
End of explanation
bc_counts.index = bc_counts['variant']
bc_counts.index.name = None
del bc_counts['variant']
bc_counts
Explanation: The data frame has the information we want, but it will be easier to use later if it's indexed by variant rather than row number.
End of explanation
bc_cutoff = 10
multi_bc_variants = bc_counts.loc[bc_counts['barcodes']  >= bc_cutoff].index[1:]
multi_bc_variants
Explanation: We'll use a cutoff to choose variants with a minimum number of unique barcodes, and store this subset in a new index. We'll also exclude the wild type by dropping the first entry of the index.
End of explanation
multi_bc_scores = my_store.select('/main/variants/scores', where='index in multi_bc_variants')
multi_bc_scores
Explanation: We can use this index to get condition-level scores for these variants by querying the "/main/variants/scores" table. Since we are working with an Experiment HDF5 file, the data frame column names are a MultiIndex with two levels, one for experimental conditions and one for data values (see the pandas documentation for more information).
End of explanation
my_store.close()
Explanation: There are fewer rows in multi_bc_scores than in multi_bc_variants because some of the variants were not scored in all replicate selections, and therefore do not have a condition-level score.
Now that we're finished getting data out of the HDF5 file, we'll close it.
End of explanation
bc_counts['score'] = multi_bc_scores['E3', 'score']
bc_counts
Explanation: We'll add a column to the bc_counts data frame that contains scores from the multi_bc_scores data frame. To reference a column in a data frame with a MultiIndex, we need to specify all column levels.
End of explanation
bc_counts.dropna(inplace=True)
bc_counts
Explanation: Many rows in bc_counts are missing scores (displayed as NaN) because those variants were not in multi_bc_scores. We'll drop them before continuing.
End of explanation
fig, ax = plt.subplots()
enrich_plot.configure_axes(ax, xgrid=True)
ax.plot(bc_counts['barcodes'], 
        bc_counts['score'], 
        linestyle='none', marker='.', alpha=0.6,
        color=enrich_plot.plot_colors['bright5'])
ax.set_xlabel("Unique Barcodes")
ax.set_ylabel("Variant Score")
Explanation: Now that we have a data frame containing the subset of variants we're interested in, we can make a plot of score vs. number of unique barcodes. This example uses functions and colors from the Enrich2 plotting library.
End of explanation
<END_TASK>