diff --git a/onboarding_documentation/technical_documentation/analyst_tutorial/Final_analyst_tutorial.ipynb b/onboarding_documentation/technical_documentation/analyst_tutorial/Final_analyst_tutorial.ipynb new file mode 100644 index 0000000..0612ca4 --- /dev/null +++ b/onboarding_documentation/technical_documentation/analyst_tutorial/Final_analyst_tutorial.ipynb @@ -0,0 +1,421 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e87b9d2a", + "metadata": {}, + "source": [ + "## Analyst Tutorial: Ancestry Comparisong using PCA" + ] + }, + { + "cell_type": "markdown", + "id": "0f95cb81", + "metadata": {}, + "source": [ + "This Jupyter Notebook documents the analysis of ancestry comparison using Principal Component Analysis (PCA). The analysis is performed on a dataset generated by the Large Cohort Pipeline, which processed raw FASTQ files and generated PCA values. Additionally, we will compare the ancestry of our dataset to a public dataset by calculating PCA values for the public dataset.\n", + "\n", + "### Background\n", + "The Large Cohort Pipeline is a powerful tool for processing genomic data and extracting valuable insights. By submitting raw FASTQ files to the pipeline, we obtained PCA values that provide information about the genetic ancestry of our dataset. To gain further insights and compare our dataset with a public dataset, we will calculate PCA values for the public dataset as well.\n", + "\n", + "### Objectives\n", + "The main objectives of this analysis are as follows:\n", + "1. Compare the genetic ancestry of our dataset with the public dataset using PCA.\n", + "2. Visualize the PCA results to identify patterns and similarities/differences between the datasets.\n", + "3. Gain insights into the genetic diversity and population structure of the datasets.\n", + "4. Gain familiarity with using Hail to conduct analyses of large datasets\n", + "\n", + "### Methodology\n", + "We will the Hail for genomic data analysis, including PCA calculation and visualisation. The notebook consists of multiple cells, each containing code snippets and markdown cells for documentation.\n", + "\n", + "### Conclusion\n", + "By comparing the PCA values of our dataset with the public dataset, we aim to gain insights into the genetic ancestry and population structure of both datasets. The analysis will provide valuable information for understanding the genetic diversity and relationships between different populations.\n", + "\n", + "Let's proceed with the analysis and explore the fascinating world of genetic ancestry comparison using PCA!" + ] + }, + { + "cell_type": "markdown", + "id": "e642d8ed", + "metadata": {}, + "source": [ + "### Hail Background\n", + "**Hail**: Hail is an open-source, scalable framework for exploring and analyzing genomic data. It's built to handle large-scale genomic data and provides tools for data manipulation, statistical analysis, and visualisation. Hail is particularly powerful for genome-wide association studies (GWAS) and other genetic analyses.\n", + "\n", + "**Hail Tables**: Hail Tables are one of the two primary data structures in Hail. They are similar to SQL tables, R data.frames, or pandas DataFrames. A Hail Table is a distributed dataset, that can be used to hold various types of data, including numerical, categorical, and genomic data. Hail Tables support a variety of operations like filtering, aggregation, and joins.\n", + "\n", + "**Hail MatrixTables**: MatrixTables are the second primary data structure in Hail and are designed to hold large-scale, distributed genomic data. A MatrixTable consists of a two-dimensional matrix of entries, indexed by row keys (representing variant data) and column keys (representing sample data). Each cell in the matrix represents the genetic data for a specific variant and sample. Along with the matrix, there are also row, column, and global fields for storing related data. MatrixTables support a variety of operations like filtering, aggregation, linear algebra operations, and joins." + ] + }, + { + "cell_type": "markdown", + "id": "20e64b58", + "metadata": {}, + "source": [ + "Hail provides tutorials specifically designed to introduce individuals to working with Tables and MatrixTables:\n", + "- [Table tutorial](https://hail.is/docs/0.2/tutorials/03-tables.html)\n", + "- [MatrixTable tutorial](https://hail.is/docs/0.2/tutorials/07-matrixtable.html)\n", + "\n", + "Additionally, [here](https://hail.is/docs/0.2/tutorials-landing.html) is Hail's complete list of tutorials that provide users with the opportunity to explore more manipulations using their data structures" + ] + }, + { + "cell_type": "markdown", + "id": "ee7ab28a", + "metadata": {}, + "source": [ + "![Hail Matrix Table Diagram](Hail_matrix_table_diagram.png)\n", + "Source: https://hail.is/docs/0.2/_static/cheatsheets/hail_matrix_tables_cheat_sheet.pdf" + ] + }, + { + "cell_type": "markdown", + "id": "4cc981f2", + "metadata": {}, + "source": [ + "## Imports and Reading in Data from GCP" + ] + }, + { + "cell_type": "markdown", + "id": "7f93b1b4", + "metadata": {}, + "source": [ + "Before starting, we need to import the necessary packages. If you have followed the notebook setup properly, Hail should already be installed along with Python and some basic packages, and can be imported directly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d6ccf14-c2b7-4d1c-8c1f-642d2ed2aba5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import hail as hl\n", + "import pprint\n", + "from hail.plot import show\n", + "hl.plot.output_notebook()" + ] + }, + { + "cell_type": "markdown", + "id": "da901394", + "metadata": {}, + "source": [ + "This notebook should be set up inside the bucket that contains the data we are wanting to analyse, this way we have permissions to access the data. In the case of this tutorial it is the `bioheart` bucket. While some buckets/projects are capable of reading the contents of other buckets, not all of them are, so it is best to work from a notebook that is within the bucket containing our data." + ] + }, + { + "cell_type": "markdown", + "id": "25ca7847", + "metadata": {}, + "source": [ + "#### Reading in Data\n", + "We will now be reading in our data. If you remember in our config file we set the parameter `output_version` to be `5.0`. This has created a directory within the `large_cohort` directory of the `cpg-bioheart-test` bucket.\n", + "While we have mentioned `Tables` (`.ht`) and `MatrixTables` (`.mt`) above, we will also be working with `Variant Dataset` (`.vds`) files. \n", + "- `vds` or `Variant Dataset` ([documentation](https://hail.is/docs/0.2/vds/index.html))\n", + " - Contains variant information for every sample in the cohort\n", + " - Briefly; because we are dealing with cohort-level data there will be regions of genome that, for every sample in the cohort, is homozygous reference, as well as information where sample(s) are non-reference. To maintain both reference and non-reference information in a scalable manner, a `.vds` file is made up of two `MatrixTables`, the `reference_data` and `variant_data` `MatrixTable`'s.\n", + "- `sample_qc_ht`: this Hail Table is generated by the `SampleQC` stage of the `large cohort` pipeline. It contains the sex of each sample in the cohort and PCA values for each sample to be used in ancestry comparison as well as quality metrics.\n", + "- `relateds_to_drop_ht`: this Hail Table is generated by the `Relatedness` stage of the `large cohort` pipeline. It contains information about related samples in the cohort and which samples, if related, should be dropped from the cohort. For samples that are related, they are ranked according to their sample QC metrics and the sample with the worst metrics is dropped from the cohort so that there is only one sample left per family. See code [here](https://github.com/populationgenomics/production-pipelines/blob/824019c6eb0387d10aff047145e92583cd3e701c/cpg_workflows/large_cohort/relatedness.py#L78)\n", + "- `pop_ht`: a Hail Table containing the assigned global population labels. Hail uses a Random Forest machine learning model based on PCA scores to assign global population labels to samples. See code [here](https://github.com/populationgenomics/production-pipelines/blob/824019c6eb0387d10aff047145e92583cd3e701c/cpg_workflows/large_cohort/ancestry_pca.py#L215). Population label is assigned based on previously known ancestry of the HGDP and 1KG datasets (calculated from gnomAD).\n", + "- `vqsr.ht`: Hail Table resulting from the Variant Quality Score Recalibration process. This table contains information about the variants that were filtered out during the VQSR process.\n", + "- `freq_ht`: Hail Table containing information about the allele frequency of each variant in the cohort. This table is generated by the `VariantQC` stage of the `large cohort` pipeline. Generates frequency annotations of AF, AC, AN, InbreedingCoeff\n", + " - AF (Allele Frequency): This is the proportion of all alleles at a given location in a population that are of a particular type\n", + " - AC (Allele Count): This is the total number of times a specific allele appears in a population.\n", + " - AN (Allele Number): This is the total number of alleles at a given location in a population.\n", + " - InbreedingCoeff: This is a measure of the likelihood of genetic traits in a population being inherited due to the parents being closely related.\n", + "- `dense_mt`: Hail Matrix Table containing the genotype information for each sample in the cohort. This table is generated by the `Genotype` and subsequent `GenotypeQC` stages of the `large cohort` pipeline. This table contains the genotype information for a set of pre-determined QC sites (according to `gs://cpg-common-main/references/gnomad/v0/sample_qc/pre_ld_pruning_qc_variants.ht`) for each sample in the cohort. The table is dense, meaning that it contains all variants for all samples in the cohort." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe84e1cc-0b8d-473b-a239-24e7ae14187b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "vds = hl.vds.read_vds('gs://cpg-bioheart-test/vds/5-0.vds/')\n", + "sample_qc_ht = hl.read_table('gs://cpg-bioheart-test/large_cohort/5-0/sample_qc.ht/')\n", + "relateds_to_drop_ht = hl.read_table('gs://cpg-bioheart-test/large_cohort/5-0/relateds_to_drop.ht/')\n", + "pop_ht = hl.read_table('gs://cpg-bioheart-test/large_cohort/5-0/ancestry/inferred_pop.ht/')\n", + "vqsr_ht = hl.read_table('gs://cpg-bioheart-test/large_cohort/5-0/vqsr.ht/')\n", + "freq_ht = hl.read_table('gs://cpg-bioheart-test/large_cohort/5-0/frequencies.ht/')\n", + "dense_mt = hl.read_matrix_table('gs://cpg-bioheart-test/large_cohort/5-0/dense_subset.mt/')" + ] + }, + { + "cell_type": "markdown", + "id": "3197c3d7", + "metadata": {}, + "source": [ + "### Pre-processing and Filtering" + ] + }, + { + "cell_type": "markdown", + "id": "d522cce8", + "metadata": {}, + "source": [ + "As stated previously, a `hl.vds` object is a combination of two `MatrixTables`; `vds.reference_data` and `vds.variant_data`, which are both sparse representations of the data. `vds.reference_data` contains reference data as 'blocks' of information, meaning that reference data is not explicitly stored at each variant site but is represented in intervals ('blocks') where it is necessary. `vds.variant_data` is a sparse matrix of non-reference calls. Our analysis requires a dense representation of the data, meaning each position in the genome is represented in the data i.e. reference data is explicitly stored at each position and not in 'blocks', as well as non-reference calls being present. To achieve this, we will use the `to_dense_mt()` function to convert our `vds` into a `MatrixTable` with a dense representation of the data.\n", + "\n", + "Additionally, we use the `split_multi()` to split multi-allelic variants into bi-allelic variants. After splitting, each row represents a single variant." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4289f74d-60a4-4dfb-b17b-e198df09acff", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Row-level tables require a split+dense matrix table:\n", + "vds = hl.vds.split_multi(vds, filter_changed_loci=True)\n", + "mt = hl.vds.to_dense_mt(vds)" + ] + }, + { + "cell_type": "markdown", + "id": "94860935", + "metadata": {}, + "source": [ + "#### Sample vs Row level filtering\n", + "We will now filter our `MatrixTable` according to the outputs of the pipeline. \n", + "Any filtering that works using the column key is sample-level filtering and thus will be filtering on sample-level metrics. \n", + "For example: any samples that did not pass the `SampleQC` stage of the pipeline will be removed from the `MatrixTable` if they are not in the `sample_qc_ht` table. This is sample-level filtering because we are filtering on the sample-level metrics in the `sample_qc_ht` table.\n", + "`vqsr.ht` contains information about the variants that were filtered out during the VQSR process. We will use this table to filter out variants that did not pass the VQSR process. This is row-level filtering because we are filtering on the variant-level metrics in the `vqsr.ht` table." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e74f5e6-1c91-445e-af95-fd5ff131e912", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Hard-filtering samples and variants:\n", + "mt = mt.filter_cols(hl.len(sample_qc_ht[mt.col_key].filters) > 0, keep=False) # sample-level filtering\n", + "mt = mt.filter_cols(hl.is_defined(relateds_to_drop_ht[mt.col_key]), keep=False) # sample-level filtering\n", + "mt = mt.filter_rows(hl.len(vqsr_ht[mt.row_key].filters) > 0, keep=False) # row-level filtering" + ] + }, + { + "cell_type": "markdown", + "id": "18736b06", + "metadata": {}, + "source": [ + "#### Annotation\n", + "\n", + "Annotating our `MatrixTable` means adding `column fields` to the `MatrixTable`'s `columns` `Hail Table` that contain information about the samples and, vice versa, adding new `row fields` to the `rows` `Hail Table` that contains additional information about the calls. We will annotate our `MatrixTable` with the following information:\n", + "- Sample QC metrics from the `sample_qc_ht` table\n", + "- Sample population label from the `pop_ht` table\n", + "- Call information from the `vqsr` table\n", + "- Call frequency information from the `freq_ht` table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "510813e0-e22b-4be2-9d7b-c6f4e715be31", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Annotating samples and variants:\n", + "mt = mt.annotate_cols(**sample_qc_ht[mt.col_key])\n", + "mt = mt.rename({'filters' : 'sample_qc_filters'}) # to avoid conflict with vqsr column labels\n", + "mt = mt.annotate_cols(**pop_ht[mt.col_key])\n", + "mt = mt.annotate_rows(**vqsr_ht[mt.row_key])\n", + "mt = mt.annotate_rows(**freq_ht[mt.row_key])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6ae5753-b90a-48d5-a435-9d76929cb758", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# import annotations table with ancestry information\n", + "# even though we have an 'ancestry' output from the pipeline, we won't use it\n", + "table = (hl.import_table('annotations.txt', impute=True).key_by('Sample'))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "049c9540-1e6d-470e-877c-f17db3a215df", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# annotate matrix table columns with ancestry information\n", + "mt = mt.annotate_cols(superpopulation=table[mt.s].superpopulation, \n", + " population=table[mt.s].population)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f286ecc-0784-478c-b84e-76254f36aac2", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# plot ancestry\n", + "p = hl.plot.scatter(mt.pca_scores[0],\n", + " mt.pca_scores[1],\n", + " label=mt.superpopulation,\n", + " title='PCA', xlabel='PC1', ylabel='PC2')\n", + "show(p)" + ] + }, + { + "cell_type": "markdown", + "id": "f7c23590", + "metadata": {}, + "source": [ + "## Conducting PCA on Public Dataset\n", + "\n", + "In this section, we will be conducting a Principal Component Analysis (PCA) using Hail on the publicly available Human Genome Diversity Project - 1000 Genomes Project (HGDP-1KG) dataset. As we've seen, PCA can be used to identify and visualise genetic variation between individuals. It can help us understand the genetic structure of our dataset and identify potential population stratification, which is crucial for many genetic analyses.\n", + "\n", + "The HGDP-1KG dataset is a comprehensive collection of human genetic variation from globally diverse populations. It combines data from two significant projects: the Human Genome Diversity Project (HGDP) and the 1000 Genomes Project (1KG). This dataset is readily accessible through Hail's interface as it is hosted by `gnomAD`.\n", + "\n", + "Our goal in this section is to compare the ancestry of our dataset with the ancestry of the HGDP-1KG dataset. By projecting our samples onto the genetic variation patterns identified in the HGDP-1KG dataset, we can see where our samples fit in the global population structure. This can provide valuable insights into the genetic background of our samples and inform subsequent analyses.\n", + "\n", + "Let's get started with the PCA." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f86783b0-340c-4e23-a7c7-fd1b47127ad5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "hl.utils.get_1kg('data/')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d22599dd-007c-4fc4-b246-948ba059c680", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mt_1kg = hl.read_matrix_table('data/1kg.mt')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ddd6d2ba-0609-4752-99dc-47c33b5b2d23", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "table_1kg = (hl.import_table('data/1kg_annotations.txt', impute=True)\n", + " .key_by('Sample'))\n", + "mt_1kg = mt_1kg.annotate_cols(pheno = table_1kg[mt_1kg.s])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72024335-3aa6-4039-8c00-0b3c212f8a27", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Quality control\n", + "mt_1kg = hl.sample_qc(mt_1kg)\n", + "mt_1kg = mt_1kg.filter_cols((mt_1kg.sample_qc.dp_stats.mean >= 4) & (mt_1kg.sample_qc.call_rate >= 0.97))\n", + "\n", + "ab = mt_1kg.AD[1] / hl.sum(mt_1kg.AD)\n", + "\n", + "filter_condition_ab = ((mt_1kg.GT.is_hom_ref() & (ab <= 0.1)) |\n", + " (mt_1kg.GT.is_het() & (ab >= 0.25) & (ab <= 0.75)) |\n", + " (mt_1kg.GT.is_hom_var() & (ab >= 0.9)))\n", + "\n", + "fraction_filtered = mt_1kg.aggregate_entries(hl.agg.fraction(~filter_condition_ab))\n", + "print(f'Filtering {fraction_filtered * 100:.2f}% entries out of downstream analysis.')\n", + "mt_1kg = mt_1kg.filter_entries(filter_condition_ab)\n", + "\n", + "mt_1kg = hl.variant_qc(mt_1kg)\n", + "\n", + "mt_1kg = mt_1kg.filter_rows(mt_1kg.variant_qc.AF[1] > 0.01)\n", + "mt_1kg = mt_1kg.filter_rows(mt_1kg.variant_qc.p_value_hwe > 1e-6)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c999853d-4819-4bd3-9608-2330cbebcd54", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Calculate PCA scores\n", + "eigenvalues, pcs, _ = hl.hwe_normalized_pca(mt_1kg.GT)\n", + "mt_1kg = mt_1kg.annotate_cols(scores = pcs[mt_1kg.s].scores)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61b2ad7e-517f-4a3b-afc5-5952c00d4952", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "p = hl.plot.scatter(mt_1kg.scores[0],\n", + " mt_1kg.scores[1],\n", + " label=mt_1kg.pheno.SuperPopulation,\n", + " title='PCA', xlabel='PC1', ylabel='PC2')\n", + "show(p)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:python310]", + "language": "python", + "name": "conda-env-python310-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/onboarding_documentation/technical_documentation/analyst_tutorial/analyst_tutorial.md b/onboarding_documentation/technical_documentation/analyst_tutorial/analyst_tutorial.md new file mode 100644 index 0000000..f9e14c4 --- /dev/null +++ b/onboarding_documentation/technical_documentation/analyst_tutorial/analyst_tutorial.md @@ -0,0 +1,157 @@ +As an analyst you will be running analysis jobs via the `analysis-runner`, producing output, and then analysing the output. This tutorial will walk you through the process of running an analysis job, and then analysing the output. + +In this tutorial, we will begin with raw fastq files and process them using the `large cohort` pipeline. This pipeline performs several steps: + +1. Aligns the reads to the reference genome. +2. Conducts quality control checks. +3. Calls variants and combines these variant calls. +4. Performs an ancestry analysis. +5. Plots the results of the ancestry analysis. + +While the `large cohort` pipeline includes additional stages, for the purpose of this tutorial, we will only execute up to and including the Ancestry stage, which carries out Principal Component Analysis (PCA). + +For a description of what each stage of the `large cohort` pipeline does, please see the `large cohort` pipeline documentation [here](https://github.com/populationgenomics/production-pipelines#large-cohort-workflow) + +## Table of Contents +- [1. Getting test data](#getting-test-data) +- [2. Writing a config file](#writing-a-config-file) +- [3. Submitting an analysis job](#submitting-an-analysis-job) +- [4. Analysing the output](#analysing-the-output) + +## 1. Getting test data +As an analyst, it's crucial to initially run our analysis jobs on test data before executing them on the actual data. The primary reason for this approach is cost efficiency. Large datasets can be incredibly expensive to process and debug. By using test data, we can estimate the time and cost of running the job on the full dataset, allowing us to optimise resources. Additionally, this process helps us ensure that the pipeline is functioning correctly and producing the expected output. While these checks can also be performed on the main dataset, using test data allows for quicker iterations and adjustments. It's also a good practice to run the analysis job on test data whenever the pipeline is updated, to verify that the updates haven't introduced any issues. + +In this tutorial, we will be working with the bioheart project dataset. The BiohEART project aims to identify individuals at risk of atherosclerosis by detecting genetic markers in blood before the disease manifests. + +Our first step is to verify your access to the bioheart dataset. To do this, execute the following command: + +```bash +gsutil ls gs://cpg-bioheart-main-upload/ +``` + +The output of the above command should list a number of files ending in .fastq.gz. These are read files, which serve as the input to the pipeline. If instead you encounter an AccessDeniedException: 403, it means you do not have the necessary access to the bioheart dataset. The error message will look something like this: +``` +AccessDeniedException: 403 your.username@populationgenomics.org.au does not have storage.objects.list access to the Google Cloud Storage bucket. Permission 'storage.objects.list' denied on resource (or it may not exist). +``` + +### Subsetting the bucket +Our goal is to validate the pipeline's functionality and output by running it on a subset of a dataset. For this tutorial, we'll use a subset of the `bioheart` dataset. We'll employ the [create_test_subset.py](https://github.com/populationgenomics/metamist/blob/f6c226d08a8ee9875014d8c99cfe119742221efb/scripts/create_test_subset.py) script from the `metamist` repository to achieve this. It's important to note that this script is not specific to the bioheart project and can be used with any dataset. The script selects a random set of samples from the specified dataset and copies their read files to a test bucket. + +Before executing this script, ensure that you have installed all necessary dependencies, which include `gsutil` and `analysis-runner`. The `analysis-runner` is a command-line tool that facilitates the submission of analysis jobs to the cloud. It is instrumental in tracking analysis runs, guaranteeing the reproducibility of analyses, and simplifying the process of conducting analyses. In our setup, direct execution of scripts on production data is not permitted. Instead, we mandate that code be committed, reviewed, and subsequently executed. The `analysis-runner` is the tool that enables this process, acting as a bridge between the user and the cloud resources. It operates by allowing an analyst to submit a request for an analysis via the CLI. The `analysis-runner` server then processes this request, verifies the user's authorisation, creates a Batch pipeline from a specific commit in a GitHub repository, and initiates the analysis on behalf of the user. + +For more information on the analysis-runner please check the repo [here](https://github.com/populationgenomics/analysis-runner) and an exercise [here](https://github.com/populationgenomics/team-docs/tree/main/exercise-analysis-runner). + +If you haven't done so already, clone the `metamist` repository to your local machine. Once you've done that, navigate to your `metamist` repository. Make sure you're on the `main` branch of `metamist`, then execute the following command from the CLI: + +**TODO: We will need to set up more genomes in the actual bucket we're using for downstream ancestry analysis to work** +```bash +analysis-runner \ +--dataset bioheart --description "populate bioheart test subset" --output-dir "bioheart-test" \ +--access-level full \ +scripts/create_test_subset.py --project bioheart --samples XPG280371 XPG280389 XPG280397 XPG280405 XPG280413 --skip-ped +``` +**FOR REFERENCE: The above was taken from [this](https://centrepopgen.slack.com/archives/C03FA2M1MR9/p1700020527448029?thread_ts=1699935103.776929&cid=C03FA2M1MR9) Slack thread** + +Note that `create_test_subset.py` script takes `samples` as an input and not `sequencingGroups`. `samples` start with `XPG` while `sequencingGroups` start with `CPG`. + + +## 2. Writing a config file + +Now that we have test data available in the `bioheart-test` bucket, we can begin writing a config file. + +Config files, written in TOML format, are used to specify the parameters of an analysis job for the `large cohort` pipeline. They contain all the necessary parameters to run the job, such as images, references, specific stages of the pipeline to run, and the input data. You can find a comprehensive description of config files and their usage in the `team-docs` repo [here](https://github.com/populationgenomics/team-docs/blob/13755bd51356b50ce11e6be78a76e53ed0a3ccb1/cpg_utils_config.md). + +Please note that while config files are essential for the `large cohort` pipeline, they are not required to start the `analysis-runner`, which is a wrapper that can be used independently. + +When running stages in the `large cohort` pipeline, there are set config files that are used, one of which is the defaults file you can see [here](https://github.com/populationgenomics/production-pipelines/blob/0bcf9775206f10ee91ac197c8c178f844ecad447/cpg_workflows/defaults.toml). It's important to understand that this defaults file does not cover all possible parameters. It provides a base set of parameters, but it may not include every parameter that could be relevant for your specific analysis. User-defined parameters in your specific config file can override these defaults as necessary. You can find the default config file for the `large cohort` pipeline [here](https://github.com/populationgenomics/production-pipelines/blob/main/configs/defaults/large_cohort.toml), where you can see a description of some of the parameters and their uses. + +#### Guide: Writing a config file +If we want to run the `large cohort `pipeline on the `bioheart-test` dataset, we need to create a config file that is capable of doing this. Here's a step-by-step guide on how to do it: + +Since the default config file serves as our base, we only need to define the parameters that deviate from these defaults. This includes parameters that are left undefined in the default file due to their dataset-specific nature. In this case, we need to specify the input dataset, the sequencing type, and the output version. We also have the flexibility to specify the sequencing groups we want to run on. + +We can do the above by creating a new config file specifying the following paramters in the `[workflow]` of the `toml` file: + +```TOML +[workflow] +input_datasets = ['bioheart-test'] +sequencing_type = 'genome' +output_version = '1.0' # TODO: Do we want them to specify output_version? +only_sgs = [] # to be used to demonstrate how to run on a subset of samples +``` +- `input_datasets`: This parameter specifies the datasets to be loaded as inputs for the analysis. If this parameter is not provided, the datasets will be determined automatically based on the input of `analysis-runner`. It's an array that can contain multiple dataset names e.g. `input_datasets = ['bioheart', 'tob-wgs']`. +- `only_sgs`: This parameter is used to limit the analysis to specific `sequencingGroups`. If you want to run the analysis on a subset of sequencing groups, you can specify their IDs in this array. +- `sequencing_type`: This parameter is used to limit the data to a specific sequencing type. The default value is 'genome', but it can be changed to other sequencing types like 'exome' depending on the analysis. The `bioheart` dataset is WGS, so we will use 'genome' as the value for this parameter. +- `output_version`: This parameter is used to suffix the location of the workflow outputs (`get_workflow().prefix`) with a specific string. By default, the hash of all input paths will be used, though we can adopt a 'version' approach that can be useful for versioning your outputs, especially when you run the same analysis multiple times with different parameters or input data. + +- Please note that you do not need to change anything in the default `large cohort` config file. The changes you make should be in the new config file you're creating for the `bioheart-test` analysis that will eventually overwrite parameters in default files. + +Once the config file has been written, save it with an appropriate name (e.g. `bioheart-test.toml`). Instead of saving it on your local machine, push it to the `production-pipelines-configuration` repository, which is a private repository. In the following step, we will use the `analysis-runner` to submit the job and point it to the config file we have just pushed. + +## 3. Submitting an analysis job +Now that we have a config file, we can submit an analysis job. To do this we will use the `analysis-runner` tool. The `analysis-runner` tool is a command line tool that is used to submit analysis jobs to the cloud. + +Remember, we are using the `large cohort` pipeline to conduct an analysis of several genomes in the `bioheart` dataset. As the `large cohort` pipeline is a component of the `Production-pipelines` repository, we need to be within this repository to run it. This means that our current working directory should be the `Production-pipelines` repository when executing the pipeline. So our next step is to clone the `Production-pipelines` repository to our local machine. To do this, execute the following command in a directory of your choosing: + +```bash +git clone --recurse-submodules git@github.com:populationgenomics/production-pipelines.git +cd production-pipelines +``` + +Ensure you are on the `main` branch. + +#### Guide: Submit an analysis job +When working in any pipeline with many moving parts, it's important to be able to track the progress of your analysis. As mentioned previously, the `analysis-runner` is our tool of choice for this purpose. Therefore, in order to start our analysis, we point the `analysis-runner` to our config file with parameters specific to our analysis. + +The `analysis-runner` tool is a powerful utility for submitting analysis jobs. Detailed instructions on its usage can be found in the [official documentation](https://github.com/populationgenomics/analysis-runner) (also linked earlier). Before you can use it, ensure it's installed via `pip`. + +You have the option to either encapsulate the `analysis-runner` command within a bash script or run it directly from the command line, based on your preference. + +Let's break down the structure of an `analysis-runner` command: +- `--dataset`: This parameter specifies the dataset to be loaded as input for the analysis. +- `--description`: This parameter is used to provide a description of the analysis job. +- `--access-level`: This parameter is used to specify the access level of the analysis job. The access level can be either `test` or `full`. The `test` access level is used for test data, while the `full` access level is used for production data. +- `--output-dir`: This parameter is used to specify the output directory for the analysis job. The output directory is where the output files will be stored in the bucket of `dataset` provided (in this case `bioheart-test`) + - Because we specified the `access-level` as `test`, the output directory will be in the `bioheart-test` bucket. Had we specified `full` as the `access-level`, the output directory would be in the `bioheart-main` bucket. +- `--config`: This parameter is used to specify the config files to be used for the analysis job. The config files are used to specify the parameters of the analysis job. You can specify multiple config files by repeating this parameter. + - **IMPORTANT:** *The order of the config files is important. The config files are read in the order they are specified. This means that if you specify the same parameter in multiple config files, the value of the parameter in the last config file will be used. This is why the config file we created is listed last as an input.* +- `--image`: This parameter is used to specify the docker image to be used for the analysis job. The docker image is used to specify the environment in which the analysis job will be run. + - TODO: what else should we say about this? +- `main.py`: This is the Python script that initiates the pipeline. In our scenario, the `large cohort` pipeline is started by this script. +- `large_cohort`: This argument identifies the specific pipeline to execute. Here, it's set to `large_cohort`, indicating that we're running the `large cohort` pipeline, which is initiated by the `main.py` script. + +Putting it all together, we get the following: + +```bash +analysis-runner \ +--dataset bioheart \ +--description "Running large cohort up to Combiner stage on bioheart" \ +--output-dir "bioheart-test-YOUR-NAME" \ +--access-level test \ +--config configs/defaults/large_cohort.toml \ +--config path/to/you/config/config_filename.toml \ +--image australia-southeast1-docker.pkg.dev/cpg-common/images/cpg_workflows:latest \ +main.py large_cohort # python script to be execusted plus arguments (i.e. pipeline name) +``` + +The output of the `analysis-runner` command is link to the Hail Batch driver job. Something like this: +``` +Submitting analysis-runner@ for dataset "bioheart" +Request submitted successfully: https://batch.hail.populationgenomics.org.au/batches/ +``` +This driver job is responsible for setting up all the subsequent jobs to run the analysis. When this driver job is complete, a link to another Hail Batch service will be provided at the bottom of the `Logs`. It will look something like this: +``` +Submitted batch , see https://batch.hail.populationgenomics.org.au/batches/ +``` + +This is where all the individual jobs of the pipeline will be listed as well as their status. Once this secondary Batch is completed, the output will be available in `bioheart-test`. + +## 4. Analysing the output + +Now that we have some output to work with, lets begin analysing. + +Analysis is best done within a Jupyter Notebook. However, in order to interact with our data in the cloud using Hail and jupyter notebooks we need to do some set up first. + +Instructions on setting up a Jupyter Notebook in the cloud can be found [here](https://github.com/populationgenomics/team-docs/blob/main/notebooks.md). + +Once you have your notebook setup on GCP follow along with the instructions in this tutorials accompanying notebook, copying over code as required, good luck!