<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Hail workshop
This notebook will introduce the following concepts
Step1:  Modes
Jupyter has two modes, a navigation mode and an editor mode.
Navigation mode
Step2:  This is a markdown cell, so even if something looks like code (as below), it won't get executed!
my_variable += 1
Step3:  Common gotcha
Step4:  Now we initialize Hail and set up Bokeh to display inline in the notebook.
Step5:  Download public 1000 Genomes data
The workshop materials are designed to work on a small (~20MB) downsampled chunk of the public 1000 Genomes dataset.
You can run these same functions on your computer or on the cloud!
Step6:  It is possible to call command-line utilities from Jupyter by prefixing a line with a !
Step7:  Part 1
Step8:  Read 1KG into Hail
We represent genetic data as a Hail MatrixTable, and name our variable mt to indicate this.
Step9:  What is a MatrixTable?
Let's describe it!
The describe method prints the schema, that is, the fields in the dataset and their types.
You can see
Step10:  count
MatrixTable.count returns a tuple with the number of rows (variants) and number of columns (samples).
Step11:  show
There is no mt.show() method, but you can show individual fields like the sample ID (s), or the locus (locus).
Step12:  <font color="brightred"><strong>Exercise
Step13:  Most of Hail's functionality is totally general-purpose!
Functions like summarize_variants are built out of Hail's general-purpose data manipulation functionality. We can use Hail to ask arbitrary questions about the data
Step14:  Or if we had travel data
Step15:  By sorting the result in Python, we can recover an interesting bit of biology...
Step16:  <font color="brightred"><strong>Question
Step17:  We can import it as a Hail Table with hl.import_table.
We call it "sa" for "sample annotations".
Step18:  While we can see the names and types of fields in the logging messages, we can also describe and show this table
Step19:  Add sample metadata into our 1KG MatrixTable
It's short and easy
Step20:  What's going on here?
Understanding what's going on here is a bit more difficult. To understand, we need to understand a few pieces
Step21:  <font color="brightred"><strong>Exercise
Step22:  <font color="brightred"><strong>Exercise
Step23:  Genotype QC
We explored GQ above, and analysts often set thresholds for GQ to filter entries (genotypes). Another useful metric is allele read balance.
This value is defined by
Step24:  Variant QC
Hail has the function hl.variant_qc to compute a list of useful statistics about variants from sequencing data.
Once again, Click the link above to see the documentation!
Step25:  Remove rare sites
Step26:  Remove sites far from Hardy-Weinberg equilbrium
Step27:  Part 3
Step28:  Two of the plots that analysts generally produce are a Manhattan plot and a Q-Q plot.
We'll start with the Manhattan plot
Step29:  Confounded!
The Q-Q plot indicates extreme inflation of p-values.
If you've done a GWAS before, you've probably included a few other covariates -- age, sex, and principal components.
Principal components are a measure of genetic ancestry, and can be used to control for population stratification.
We can compute principal components with Hail
Step30:  The eigenvalues reflect the amount of variance explained by each principal component
Step31:  The scores are the principal components themselves, computed per sample.
Step32:  The loadings are the contributions to each component for each variant.
Step33:  We can annotate the principal components back onto mt
Step34:  Principal components measure ancestry
Step35:  <font color="brightred"><strong>Question
Step36:  Part 4
Step37:  Create an interval key
Step38:  Annotate variants using these intervals
Step39:  Aggregate genotypes per gene
There is no hl.burden_test function -- instead, a burden test is the composition of two modular pieces of Hail functionality
Step40:  What is burden_mt?
It is a gene-by-sample matrix (compare to mt, a variant-by-sample matrix).
It has one row field, the gene.
It has one entry field, n_variants.
It has all the column fields from mt.
Run linear regression per gene
This should look familiar!
Step41:  Sorry, no hl.plot.manhattan for genes!
Instead, we can sort by p-value and print
Step42:  <font color="brightred"><strong>Exercise
Step43:  Infer sex from X-chromosome data
Step44:  Simulate genetic data
<ASSISTANT_TASK:>
Python Code:
print('Hello, world')
Explanation: Hail workshop
This notebook will introduce the following concepts:
Using Jupyter notebooks effectively
Loading genetic data into Hail
General-purpose data exploration functionality
Plotting functionality
Quality control of sequencing data
Running a Genome-Wide Association Study (GWAS)
Rare variant burden tests
Hail on Jupyter
From https://jupyter.org: 
"The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text. Uses include: data cleaning and transformation, numerical simulation, statistical modeling, data visualization, machine learning, and much more."
In the last year, the Jupyter development team released Jupyter Lab, an integrated environment for data, code, and visualizations. If you've used R Studio, this is the closest thing that works in Python (and many other languages!).
Why notebooks?
Part of what we think is so exciting about Hail is that it has coincided with a larger shift in the data science community.
Three years ago, most computational biologists at Broad analyzed genetic data using command-line tools, and took advantage of research compute clusters by explicitly using scheduling frameworks like LSF or Sun Grid Engine.
Now, they have the option to use Hail in interactive Python notebooks backed by thousands of cores on public compute clouds like Google Cloud, Amazon Web Services, or Microsoft Azure.
Using Jupyter
Running cells
Evaluate cells using SHIFT + ENTER. Select the next cell and run it.
End of explanation
# This is a code cell
my_variable = 5
Explanation: Modes
Jupyter has two modes, a navigation mode and an editor mode.
Navigation mode:
<font color="blue"><strong>BLUE</strong></font> cell borders
UP / DOWN move between cells
ENTER while a cell is selected will move to editing mode.
Many letters are keyboard shortcuts! This is a common trap.
Editor mode:
<font color="green"><strong>GREEN</strong></font> cell borders
UP / DOWN/ move within cells before moving between cells.
ESC will return to navigation mode.
SHIFT + ENTER will evaluate a cell and return to navigation mode.
Cell types
There are several types of cells in Jupyter notebooks. The two you will see here are Markdown (text) and Code.
End of explanation
print(my_variable)
Explanation: This is a markdown cell, so even if something looks like code (as below), it won't get executed!
my_variable += 1
End of explanation
import hail as hl
from bokeh.io import output_notebook, show
Explanation: Common gotcha: a code cell turns into markdown
This can happen if you are in navigation mode and hit the keyboard shortcut m while selecting a code cell.
You can either navigate to Cell &gt; Cell Type &gt; Code through the top menu, or use the keyboard shortcut y to turn it back to code.
Tips and tricks
Keyboard shortcuts:
SHIFT + ENTER to evaluate a cell
ESC to return to navigation mode
y to turn a markdown cell into code
m to turn a code cell into markdown
a to add a new cell above the currently selected cell
b to add a new cell below the currently selected cell
d, d (repeated) to delete the currently selected cell
TAB to activate code completion
To try this out, create a new cell below this one using b, and print my_variable by starting with print(my and pressing TAB!
Common gotcha: the state of your code seems wrong
Jupyter makes it easy to get yourself into trouble by executing cells out-of-order, or multiple times.
For example, if I declare x:
x = 5
Then have a cell that reads:
x += 1
And finally:
print(x)
If you execute these cells in order and once, I'll see the notebook print 6. However, there is nothing stopping you from executing the middle cell ten times, printing 16!
Solution
If you get yourself into trouble into this way, the solution is to clear the kernel (Python process) and start again from the top.
First, Kernel &gt; Restart &amp; Clear Output &gt; (accept dialog).
Second, Cell &gt; Run all above.
Set up our Python environment
In addition to Hail, we import a few methods from the bokeh plotting library. We'll see examples soon!
End of explanation
hl.init()
output_notebook()
Explanation: Now we initialize Hail and set up Bokeh to display inline in the notebook.
End of explanation
hl.utils.get_1kg('data/')
Explanation: Download public 1000 Genomes data
The workshop materials are designed to work on a small (~20MB) downsampled chunk of the public 1000 Genomes dataset.
You can run these same functions on your computer or on the cloud!
End of explanation
! ls -1 data/
Explanation: It is possible to call command-line utilities from Jupyter by prefixing a line with a !:
End of explanation
hl.import_vcf('data/1kg.vcf.bgz').write('data/1kg.mt', overwrite=True)
Explanation: Part 1: Explore genetic data with Hail
Import data from VCF
The Variant Call Format (VCF) is a common file format for representing genetic data collected on multiple individuals (samples).
Hail's import_vcf function can read this format.
However, VCF is a text format that is easy for humans but very bad for computers. The first thing we do is write to a Hail native file format, which is much faster!
End of explanation
mt = hl.read_matrix_table('data/1kg.mt')
Explanation: Read 1KG into Hail
We represent genetic data as a Hail MatrixTable, and name our variable mt to indicate this.
End of explanation
mt.describe()
Explanation: What is a MatrixTable?
Let's describe it!
The describe method prints the schema, that is, the fields in the dataset and their types.
You can see:
 - numeric types:
     - integers (int32, int64), e.g. 5
     - floating point numbers (float32, float64), e.g. 5.5 or 3e-8
 - strings (str), e.g. "Foo"
 - boolean values  (bool) e.g. True
 - collections:
     - arrays (array), e.g. [1,1,2,3]
     - sets (set), e.g. {1,3}
     - dictionaries (dict), e.g. {'Foo': 5, 'Bar': 10}
 - genetic data types:
     - loci (locus), e.g. [GRCh37] 1:10000 or [GRCh38] chr1:10024
     - genotype calls (call), e.g. 0/2 or 1|0
End of explanation
mt.count()
Explanation: count
MatrixTable.count returns a tuple with the number of rows (variants) and number of columns (samples).
End of explanation
mt.s.show(5)
mt.locus.show(5)
Explanation: show
There is no mt.show() method, but you can show individual fields like the sample ID (s), or the locus (locus).
End of explanation
hl.summarize_variants(mt)
Explanation: <font color="brightred"><strong>Exercise: </strong></font> show other fields
You can see the names of fields above. show() the first few values for a few of them, making sure to include at least one row field and at least one entry field. Capitalization is important.
To print fields inside the info structure, you must add another dot, e.g. mt.info.AN.
What do you notice being printed alongside some of the fields?
Hail has functions built for genetics
For example, hl.summarize_variants prints useful statistics about the genetic variants in the dataset.
End of explanation
mt.aggregate_rows(hl.agg.count_where(mt.alleles == ['A', 'T']))
Explanation: Most of Hail's functionality is totally general-purpose!
Functions like summarize_variants are built out of Hail's general-purpose data manipulation functionality. We can use Hail to ask arbitrary questions about the data:
End of explanation
snp_counts = mt.aggregate_rows(
    hl.array(hl.agg.counter(mt.alleles)))
snp_counts
Explanation: Or if we had travel data:
data.aggregate(
    hl.agg.count_where(data.departure_city == 'Boston')
)
The counter aggregator makes it possible to see distributions of categorical data, like alleles:
End of explanation
sorted(snp_counts,
       key=lambda x: x[1])
Explanation: By sorting the result in Python, we can recover an interesting bit of biology...
End of explanation
! head data/1kg_annotations.txt
Explanation: <font color="brightred"><strong>Question: </strong></font> What is interesting about this distribution?
<font color="brightred"><strong>Question: </strong></font> Why do the counts come in pairs?
Part 2: Annotation and quality control
Integrate sample information
We're building toward a genome-wide association test in part 3, but we don't just need genetic data to do a GWAS -- we also need phenotype data! Luckily, our hl.utils.get_1kg function also downloaded some simulated phenotype data.
This is a text file:
End of explanation
sa = hl.import_table('data/1kg_annotations.txt', 
                      impute=True, 
                      key='Sample')
Explanation: We can import it as a Hail Table with hl.import_table.
We call it "sa" for "sample annotations".
End of explanation
sa.describe()
sa.show()
Explanation: While we can see the names and types of fields in the logging messages, we can also describe and show this table:
End of explanation
mt = mt.annotate_cols(pheno = sa[mt.s])
Explanation: Add sample metadata into our 1KG MatrixTable
It's short and easy:
End of explanation
mt.describe()
Explanation: What's going on here?
Understanding what's going on here is a bit more difficult. To understand, we need to understand a few pieces:
1. annotate methods
In Hail, annotate methods refer to adding new fields. 
MatrixTable's annotate_cols adds new column fields.
MatrixTable's annotate_rows adds new row fields.
MatrixTable's annotate_entries adds new entry fields.
Table's annotate adds new row fields.
In the above cell, we are adding a new column field called "pheno". This field should be the values in our table sa associated with the sample ID s in our MatrixTable - that is, this is performing a join.
Python uses square brackets to look up values in dictionaries:
d = {'foo': 5, 'bar': 10}
d['foo']
You should think of this in much the same way - for each column of mt, we are looking up the fields in sa using the sample ID s.
End of explanation
mt = hl.sample_qc(mt)
mt.sample_qc.describe()
p = hl.plot.scatter(x=mt.sample_qc.r_het_hom_var,
                    y=mt.sample_qc.call_rate)
show(p)
Explanation: <font color="brightred"><strong>Exercise: </strong></font> Query some of these column fields using mt.aggregate_cols.
Some of the aggregators we used earlier:
 - hl.agg.counter
 - hl.agg.stats
 - hl.agg.count_where
Sample QC
We'll start with examples of sample QC.
Hail has the function hl.sample_qc to compute a list of useful statistics about samples from sequencing data.
Click the link above to see the documentation, which lists the fields and their descriptions.
End of explanation
mt = mt.filter_cols(mt.sample_qc.dp_stats.mean >= 4)
mt = mt.filter_cols(mt.sample_qc.call_rate >= 0.97)
Explanation: <font color="brightred"><strong>Exercise: </strong></font> Plot some other fields!
Modify the cell above. Remember hl.plot.histogram as well!
If you want to start getting fancy, you can plot more complicated expressions -- the ratio between two fields, for instance.
Filter columns using generated QC statistics
End of explanation
# call rate before filtering
mt.aggregate_entries(hl.agg.fraction(hl.is_defined(mt.GT)))
ab = mt.AD[1] / hl.sum(mt.AD)
filter_condition_ab = (
    hl.case()
    .when(mt.GT.is_hom_ref(), ab <= 0.1)
    .when(mt.GT.is_het(), (ab >= 0.25) & (ab <= 0.75))
    .default(ab >= 0.9) # hom-var
)
mt = mt.filter_entries(filter_condition_ab)
# call rate after filtering
mt.aggregate_entries(hl.agg.fraction(hl.is_defined(mt.GT)))
Explanation: Genotype QC
We explored GQ above, and analysts often set thresholds for GQ to filter entries (genotypes). Another useful metric is allele read balance.
This value is defined by:
$\quad AB = \dfrac{N_{alt}}{{N_{ref} + N_{alt}}}$
Where $N_{ref}$ is the number of reference reads and $N_{alt}$ is the number of alternate reads.
We want
End of explanation
mt = hl.variant_qc(mt)
mt.variant_qc.describe()
mt.variant_qc.AF.show()
Explanation: Variant QC
Hail has the function hl.variant_qc to compute a list of useful statistics about variants from sequencing data.
Once again, Click the link above to see the documentation!
End of explanation
mt = mt.filter_rows(hl.min(mt.variant_qc.AF) > 1e-6)
Explanation: Remove rare sites:
End of explanation
mt = mt.filter_rows(mt.variant_qc.p_value_hwe > 0.005)
# final variant and sample count
mt.count()
Explanation: Remove sites far from Hardy-Weinberg equilbrium:
End of explanation
gwas = hl.linear_regression_rows(y=mt.pheno.CaffeineConsumption, 
                                 x=mt.GT.n_alt_alleles(), 
                                 covariates=[1.0])
gwas.describe()
Explanation: Part 3: GWAS!
A GWAS is an independent association test performed per variant of a genetic dataset. We use the same phenotype and covariates, but test the genotypes for each variant separately. 
In Hail, the method we use is hl.linear_regression_rows.
We use the phenotype CaffeineConsumption as our dependent variable, the number of alternate alleles as our independent variable, and no covariates besides an intercept term (that's the 1.0).
End of explanation
p = hl.plot.manhattan(gwas.p_value)
show(p)
p = hl.plot.qq(gwas.p_value)
show(p)
Explanation: Two of the plots that analysts generally produce are a Manhattan plot and a Q-Q plot.
We'll start with the Manhattan plot:
End of explanation
pca_eigenvalues, pca_scores, pca_loadings = hl.hwe_normalized_pca(mt.GT, compute_loadings=True)
Explanation: Confounded!
The Q-Q plot indicates extreme inflation of p-values.
If you've done a GWAS before, you've probably included a few other covariates -- age, sex, and principal components.
Principal components are a measure of genetic ancestry, and can be used to control for population stratification.
We can compute principal components with Hail:
End of explanation
pca_eigenvalues
Explanation: The eigenvalues reflect the amount of variance explained by each principal component:
End of explanation
pca_scores.describe()
pca_scores.scores[0].show()
Explanation: The scores are the principal components themselves, computed per sample.
End of explanation
pca_loadings.describe()
Explanation: The loadings are the contributions to each component for each variant.
End of explanation
mt = mt.annotate_cols(pca = pca_scores[mt.s])
Explanation: We can annotate the principal components back onto mt:
End of explanation
p = hl.plot.scatter(mt.pca.scores[0], 
                    mt.pca.scores[1],
                    label=mt.pheno.SuperPopulation)
show(p)
Explanation: Principal components measure ancestry
End of explanation
gwas = hl.linear_regression_rows(
    y=mt.pheno.CaffeineConsumption, 
    x=mt.GT.n_alt_alleles(),
    covariates=[1.0, mt.pheno.isFemale, mt.pca.scores[0], mt.pca.scores[1], mt.pca.scores[2]])
p = hl.plot.qq(gwas.p_value)
show(p)
p = hl.plot.manhattan(gwas.p_value)
show(p)
Explanation: <font color="brightred"><strong>Question: </strong></font> Does your plot match your neighbors'?
If not, how is it different?
Control confounders and run another GWAS
End of explanation
gene_ht = hl.import_table('data/ensembl_gene_annotations.txt', impute=True)
gene_ht.show()
gene_ht.count()
Explanation: Part 4: Burden tests
GWAS is a great tool for finding associations between common variants and disease, but a GWAS can't hope to find associations between rare variants and disease. Even if we have sequencing data for 1,000,000 people, we won't have the statistical power to link a mutation found in only a few people to any disease.
But rare variation has lots of information - especially because statistical genetic theory dictates that rarer variants have, on average, stronger effects on disease per allele.
One possible strategy is to group together rare variants with similar predicted consequence. For example, we can group all variants that are predicted to knock out the function of each gene and test the variants for each gene as a group.
We will be running a burden test on our common variant dataset to demonstrate the technical side, but we shouldn't hope to find anything here -- especially because we've only got 10,000 variants!
Import gene data
We start by importing gene names and coordinates.
End of explanation
gene_ht = gene_ht.transmute(interval = hl.locus_interval(gene_ht['Chromosome'],
                                                         gene_ht['Gene start'],
                                                         gene_ht['Gene end'], 
                                                         reference_genome='GRCh37'))
gene_ht = gene_ht.key_by('interval')
Explanation: Create an interval key
End of explanation
mt = mt.annotate_rows(gene_info = gene_ht[mt.locus])
mt.gene_info.show()
Explanation: Annotate variants using these intervals
End of explanation
burden_mt = (
    mt
    .group_rows_by(gene = mt.gene_info['Gene name'])
    .aggregate(n_variants = hl.agg.count_where(mt.GT.n_alt_alleles() > 0))
)
burden_mt.describe()
Explanation: Aggregate genotypes per gene
There is no hl.burden_test function -- instead, a burden test is the composition of two modular pieces of Hail functionality:
group_rows_by / aggregate
hl.linear_regression_rows.
While this might be a few more lines of code to write than hl.burden_test, it means that you can flexibly specify the genotype aggregation however you like. Using other tools , you may have a few ways to aggregate, but if you want to do something different you are out of luck!
End of explanation
burden_results = hl.linear_regression_rows(
    y=burden_mt.pheno.CaffeineConsumption, 
    x=burden_mt.n_variants,
    covariates=[1.0, 
                burden_mt.pheno.isFemale, 
                burden_mt.pca.scores[0], 
                burden_mt.pca.scores[1], 
                burden_mt.pca.scores[2]])
Explanation: What is burden_mt?
It is a gene-by-sample matrix (compare to mt, a variant-by-sample matrix).
It has one row field, the gene.
It has one entry field, n_variants.
It has all the column fields from mt.
Run linear regression per gene
This should look familiar!
End of explanation
burden_results.order_by(burden_results.p_value).show()
Explanation: Sorry, no hl.plot.manhattan for genes!
Instead, we can sort by p-value and print:
End of explanation
ht = hl.identity_by_descent(mt).cache()
ht.describe()
ht.filter(ht.ibd.PI_HAT > 0.20).show()
Explanation: <font color="brightred"><strong>Exercise: </strong></font> Where along the genome can we find the top gene?
Part 5: Whirlwind tour
You've seen just a very small fraction of the functionality in Hail. Here are a few examples of things that a large library makes easy.
Find related individuals using IBD (identity by descent)
End of explanation
ht = hl.impute_sex(mt.GT).cache()
ht.show()
Explanation: Infer sex from X-chromosome data
End of explanation
sim_mt = hl.balding_nichols_model(n_populations=3,
                                  n_samples=1000, 
                                  n_variants=1000)
# simulate variant effects using spike-and-slab model
spike_prob = 0.2
sim_mt = sim_mt.annotate_rows(beta = hl.rand_bool(spike_prob) * hl.rand_norm(0, 1))
# compute risk scores from betas
sim_mt = sim_mt.annotate_cols(risk = hl.agg.sum(sim_mt.beta * sim_mt.GT.n_alt_alleles()) / sim_mt.count_rows())
show(hl.plot.histogram(sim_mt.risk))
Explanation: Simulate genetic data
End of explanation
<END_TASK>