tutorial_HLAQCImputation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7080090d",
   "metadata": {},
   "source": [
    "# A Statistical Genetics Guide to Identifying HLA Alleles Driving Complex Disease"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fcaff7e",
   "metadata": {},
   "source": [
    "## Saisriram Gurajala, Saori Sakaue, Michelle Curtis\n",
    "## Updated: 02/10/2023"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a06102b8",
   "metadata": {},
   "source": [
    "# Dataset: Mock Genome Sequencing Array Data from the Human Genome Diversity Project"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a333bdb6",
   "metadata": {},
   "source": [
    "## Quality Control of Target Genotype Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c231493-2920-42e6-86b6-402ec69fb700",
   "metadata": {},
   "source": [
    "**The following explains steps and provides scripts for the 'Quality control of the target genotype data' section of Sakaue et. al.**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f62db4d",
   "metadata": {},
   "source": [
    "> **Important Note: Sample QC and Variant QC should be done on all chromosomes. However, we only used chromosome 6 in this vignette due to the size restriction of data uploads to GitHub. Please make sure to perform QC by using your PLINK genotype data on all chromosomes.**\n",
    "\n",
    "### 1. Deduplication of SNP data\n",
    "\n",
    "**Corresponds to 'Per-variant QC' subsection**\n",
    "\n",
    "> An initial step of variant qc comprises identification and removal of duplicated snps. This filtering is implemented below via a python script that parses the output of the plink --missing command for duplicates by position and lists the rsID of the duplicate with the highest genotype failure rate. The listed duplicates are then removed by the plink --exclude function. \n",
    "\n",
    "``` bash\n",
    "cd /data/\n",
    "\n",
    "#Initial missigness report of Genotype Data\n",
    "plink --bfile hgdp_chr6 --missing --out hgdp_chr6\n",
    "```\n",
    "``` bash\n",
    "#head commands to show missigness output \n",
    "head hgdp_chr6.lmiss\n",
    "```\n",
    "![Image](./images/hgdp_chr6.lmiss.png)\n",
    "\n",
    "```bash\n",
    "head hgdp_chr6.imiss\n",
    "```\n",
    "\n",
    "![Image](./images/hgdp_chr6.imiss.png)\n",
    "\n",
    "\n",
    "``` bash \n",
    "#Python script that parses --missing output and lists duplicated variants\n",
    "python ../scripts/get_duprem_var.py hgdp_chr6\n",
    "#Removal of identified duplicated snps\n",
    "plink --bfile hgdp_chr6  --exclude hgdp_chr6.remdup.snp --make-bed --out hgdp_chr6.dedup\n",
    "```\n",
    "\n",
    "### 2. Reverse/Forward strand flips \n",
    "\n",
    "**Corresponds to 'Per-variant QC' subsection**\n",
    "\n",
    "> An important step of variant quality control is to remove strands inconsistent with the reference panel, which is aligned to the forward strand. snpflip is a program identifying reverse and ambiguous strand SNPs. We recommend removal of both reverse and ambiguous strand SNPs, but here only remove ambiguous SNPs as no reverse SNPs were identified in this dataset. snpflip is available on [github](https://github.com/biocore-ntnu/snpflip), and as a python package that can be installed with the following command: `pip install snpflip`. The below commands remove the ambiguous SNPs as no reverse strand SNPs were found. \n",
    "\n",
    "``` bash\n",
    "\n",
    "#SNPFLIP to identifed snps with flipped strands\n",
    "snpflip --fasta-genome hg19.fa --bim-file hgdp_chr6.dedup.bim -o hgdp_chr6.dedup\n",
    "#Removal of amibigious and flipped snps\n",
    "plink --bfile hgdp_chr6.dedup --exclude hgdp_chr6.dedup.ambiguous --make-bed --out hgdp_chr6.dedup.ambstrandrem\n",
    "```\n",
    "\n",
    "#### **Optional: Removal of Palindromic SNPs**\n",
    "\n",
    "**Corresponds to 'Per-variant QC' subsection**\n",
    "\n",
    "> Target genotype alleles are usually aligned to the forward strand. This alignment is carried out by matching to the reference human genome sequence forward strand by position. If the alleles between the target and the reference genome are different (e.g., A/C in the reference but T/G in the target), alleles are flipped in the target dataset. However, when handling palindromic SNPs (i.e., SNPs with A/T or G/C alleles), population-derived AF and target dataset AF is compared to eliminate allele ambiguity. If AFs are largely different, we can flip alleles to be consistent with the population-derived AF. However, when the target and reference populations are different this strategy may be ineffective within the MHC. It may also be ineffective when there is large AF differences between cases and controls in case-control studies, or when the study sample size is too small (e.g., N < 100), to estimate AF accurately. Therefore, when strand information of palindromic SNPs is ambiguous in the target genotyping array or the genotyped data, it may be best to exclude all palindromic SNPs. Palindromic SNPs are included in the list of ambiguous SNPs outlined by snpflip, and were therefore removed using the above commands. \n",
    "The following commands will enable removal of AT/GC snps directly:\n",
    "\n",
    "```bash\n",
    "# Remove A/T C/G SNPs\n",
    "awk '(($5==\"C\"&&$6==\"G\")||($5==\"G\"&&$6==\"C\")||($5==\"A\"&&$6==\"T\")||($5==\"T\"&&$6==\"A\")) {print $2} ' hgdp_chr6.dedup.bim > AT_CG_SNPS.txt\n",
    "```\n",
    "```bash\n",
    "#head command to show hgdp_chr6.dedup.bim fields \n",
    "head hgdp_chr6.dedup.bim\n",
    "```\n",
    "![Image](./images/hgdp_dedup_bim_palindromic.png)\n",
    "\n",
    "```bash\n",
    "plink --bfile hgdp_chr6.dedup --exclude AT_CG_SNPS.txt --make-bed --out hgdp_chr6.dedup_AT_CG_removed\n",
    "```\n",
    "\n",
    "### 3. Initial SNP QC: Removal of very poor variants\n",
    "\n",
    "**Corresponds to 'Per-variant QC' subsection**\n",
    "\n",
    "> Prior to initial sample qc, very poor quality variants must be removed or individual genotype quality will be poorly characterized. Here we propose the use of 10% genotype call failure rate and a hardy-weinberg equilibrium threshold of 1e-10. \n",
    "\n",
    "``` bash\n",
    "#Variant missningness of 10% and hwe 1e-10 based filtering of extremely poor quality variants\n",
    "plink --bfile hgdp_chr6.ambstrandrem --geno 0.1 --hwe 1e-10 --make-bed --out hgdp_chr6.ambstrandrem.1stSNPQC\n",
    "```\n",
    "\n",
    "### 4. Sample QC\n",
    "\n",
    "**Corresponds to 'Per-individual QC' subsection**\n",
    "\n",
    "#### 4.1 Genotype Missingness Rate\n",
    "\n",
    "> Poor quality DNA samples often have high rates of genotyping failure, thereby leading to poor genotyping, imputation, and flawed association results. The below selection outlines how to filter individuals on per individual missing genotype rates using plink and R. \n",
    "\n",
    "##### Extract Missingness Report per Variant and Sample:\n",
    "\n",
    "\n",
    "```bash \n",
    "\n",
    "#Extract missingness report per variant and sample\n",
    "plink --bfile hgdp_chr6.ambstrandrem.1stSNPQC --missing --out hgdp_chr6.dedup.ambstrandrem.1stSNPQC\n",
    "\n",
    "```\n",
    "\n",
    "\n",
    "##### Identifying and Implementing a Per-Sample Missigness Threshold: \n",
    "\n",
    "    In R: \n",
    "    \n",
    "``` R\n",
    "library(ggplot2)\n",
    "name = 'hgdp_chr6.dedup.ambstrandrem.1stSNPQC.imiss'\n",
    "dat = read.table(name, header = TRUE)\n",
    "png('hgdp_chr6.dedup.ambstrandrem.1stSNPQC.imiss')\n",
    "ggplot(dat) +\n",
    "    geom_histogram(aes(x = F_MISS, y = ..count..)) +\n",
    "    theme_classic() +\n",
    "    geom_vline(xintercept = 0.015, color = \"red\", linetype = \"dashed\") +\n",
    "    labs(x = \"Genotype Missingness Rate\", y = \"Frequency\", title = \"Post 1ST SNP QC\")\n",
    "dev.off()\n",
    "```\n",
    "\n",
    "> An example figure outlining the distribution of missingness rate per sample post initial variant quality control is shown below. Since selection of thresholds is dataset dependent, the best way to select a threshold is to analyze the distribution. Selected thresholds should satisfy two conditions: retention of the bulk of the distribution and the exclusion of outliers. Here, a per-sample genotype missingness rate of 0.015 was chosen as it restricts the right tail of the visualized distribution.\n",
    "\n",
    "![Image](./images/hgdp_all_chr6.hg19.ba.only.GSA.dedup.ambstrandrem.1stSNPQC.imiss.png)\n",
    "\n",
    "\n",
    "\n",
    "``` bash\n",
    "#Remove samples with a missingness rate of > 0.015\n",
    "plink --bfile hgdp_chr6.dedup.ambstrandrem.1stSNPQC --mind 0.015 --make-bed --out hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind\n",
    "```\n",
    "\n",
    "#### 4.2 Removal of Samples with Outlier Heterozygosity \n",
    "\n",
    "> The distribution of mean heterozygosity (excluding sex chromosomes) across all individuals should be inspected to identify individuals with an excessive or reduced proportion of heterozygote genotypes, which may be indicative of DNA sample contamination or inbreeding, respectively. Below we outline a procedure to identify and remove these individuals using plink and R. \n",
    "\n",
    "##### Extracting Heterozygosity Statistics per Sample:\n",
    "\n",
    "``` bash\n",
    "#Using Plink to write per sample heterozygosity statistics\n",
    "plink --bfile hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind --het --out hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind\n",
    "```\n",
    "\n",
    "\n",
    "##### Implementing Outlier Heterozygosity Thresholds:\n",
    "\n",
    "> Here we use plink's --het output to find the heterozygosity rate per individual. The mean heterozygosity rate per individual is defined as follows: (N(NM) - O(HM)) / N(NM), where N(NM) is the rate of nonmissing genotypes per individual and O(HM) is the observed rate of homozygous genotypes per individual. The distribution of heterozygosity per individual is dataset dependent, and heterozygosity thresholds should therefore be selected in a cohort specific manner. We suggest visualizing the distribution and selecting thresholds based on the corresponding mean and standard deviation. \n",
    "    \n",
    "    In R: \n",
    "``` R\n",
    "library(ggplot2)\n",
    "#read in file \n",
    "name = 'hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het'\n",
    "dat = read.table(name, header = TRUE)\n",
    "#calculate heterozygosity rate per indvidual\n",
    "dat[, 'het'] = (dat[, 'N.NM.'] - dat[, 'O.HOM.'])/dat[, 'O.HOM.']\n",
    "m = mean(dat[, 'het'])\n",
    "sd = sd(dat[, 'het'])\n",
    "#Create sample ID column that matches fam naming convention\n",
    "dat[, 'sample'] = paste(dat[, 'FID'], dat[, 'IID'], sep = '_')\n",
    "#Plotting heterozygosity distribution as a histogram\n",
    "png(paste0(name, '.png'))\n",
    "ggplot(dat) +\n",
    "    geom_histogram(aes(x = het, y = ..count..), binwidth = 0.002) +\n",
    "    geom_vline(xintercept = c(m - (3 * sd), m, m + (3 * sd)), color = \"red\", linetype = \"dashed\") +\n",
    "    theme_classic() +\n",
    "    labs(x = \"Heterozygosity per Sample\", y = \"Frequency\", title = \"Outlier Heterozygosity Identification\")\n",
    "dev.off()\n",
    "high = dat[dat[, 'het'] > m + (3 * sd) | dat[, 'het'] < m - (3 * sd), ]\n",
    "write.table(high[, c(1:2)], paste0(name, \".het_outside_3sd.sample.txt\"), quote=F, col.names=T, row.names=F)\n",
    "```\n",
    "\n",
    "> An example figure where heterozygosity filtering may be required is shown below. The distribution is quite wide in this case, with three standard deviations in either direction spanning ~ 20-50% heterozygosity. Based on this, heterozygosity quality control is appropriate. \n",
    "\n",
    "![Image](./images/hgdp_all_chr6.hg19.ba.only.GSA.dedup.ambstrandrem.1stSNPQC.mind.het.png)\n",
    "\n",
    "> **Note: While the commands above lead to removal of both high and low heterozygosity rate individuals, we recommend only filtering out high heterozygosity individuals. Low heterozygosity is suggestive of inbreeding, which is not necessarily indicative of low sample quality. High heterozygosity can indicate sample contamination, and thereby requires removal. The outlined change will implement identification of only high heterozygosity individuals:**\n",
    "\n",
    "```R\n",
    "#original line\n",
    "high = dat[dat[, 'het'] > m + (3 * sd) | dat[, 'het'] < m - (3 * sd), ]\n",
    "#changed line\n",
    "high = dat[dat[, 'het'] > m + (3 * sd), ]\n",
    "```\n",
    "\n",
    "\n",
    "``` bash\n",
    "#Removal of outlier heterozygosity samples with plink\n",
    "plink --bfile hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind --remove hgdp_chr6.hg19.dedup.ambstrandrem.1stSNPQC.mind.het.het_outside_3sd.sample.txt --make-bed --out hgdp_chr6.hg19.dedup.ambstrandrem.1stSNPQC.mind.het\n",
    "```\n",
    "\n",
    "\n",
    "#### 4.3 Removal of Samples with Discordant Sex Information\n",
    "\n",
    "> Sample sex information is relevant in sex-specific analyses or imputation of the sex chromosome. Inaccurate assignment of sex (i.e. male recorded as female) may be a sign of sample swapping or other sample quality issue. Homozygosity rate across the non pseudo-autosomal region (non-PAR) of the X chromosome is used to impute sample sex assignment. \n",
    "\n",
    "##### Calculating X Chromosome Homozygosity per Sample:\n",
    "\n",
    ">**While our genotype data in this repository only includes mock genotype data on chromosome 6 due to size restrictions in GitHub, please assume that the data {bim,bed,fam} files are the plink genotype files with all chromosomes, including chromosome X.**\n",
    "\n",
    "``` bash\n",
    "#split off X chromosome PAR\n",
    "plink --bfile data --split-x hg19 --make-bed --out data.par_split\n",
    "#output sex check\n",
    "plink --bfile data.par_split --check-sex\n",
    "```\n",
    "\n",
    "> By default, plink's check-sex will assign samples as female if the F coefficient is < 0.2 or male if the F coefficient is > 0.8. To be most accurate, these cutoffs can be changed manually based on the distribution of the F coefficient, as visualized below.\n",
    "\n",
    "\n",
    "\n",
    "##### Implementing X Chromosome Homozygosity Thresholds:\n",
    "\n",
    "in R: \n",
    "\n",
    "```R\n",
    "library(ggplot2)\n",
    "name = 'plink.sexcheck'\n",
    "dat = read.table(name, header = TRUE)\n",
    "\n",
    "png(paste0(name, '.png'))\n",
    "ggplot(dat)+\n",
    "    geom_histogram(aes(x = F), bins = 30)+\n",
    "    theme_classic(20)+\n",
    "    labs(x = 'F coefficient', y = 'Frequency', title = 'X Chromosome Homozygosity')+\n",
    "    geom_vline(xintercept = c(0.4, 0.9), color = \"red\", linetype = \"dashed\")\n",
    "dev.off()\n",
    "```\n",
    "\n",
    "![Image](./images/plink.sexcheck.png)\n",
    "\n",
    "> We recommend removal of any samples where pedigree-defined sex is discordant with genotype-inferred sex assignment, as this suggests the sample may have been labeled incorrectly.\n",
    "\n",
    "``` bash\n",
    "#update sex in fam file and re-ouptut sex check\n",
    "plink --bfile data.par_split --impute-sex 0.4 0.85 --make-bed --out data.par_split.impute_sex\n",
    "#removal of samples with discordant sex information\n",
    "awk ' ($5==\"PROBLEM\") { print $1,$2 } ' data.par_split.impute_sex.sexcheck >> discordantsex.sample.txt\n",
    "plink --bfile data.par_split.impute_sex --remove discordantsex.sample.txt --make-bed --out data.par_split.impute_sex.concordant\n",
    "```\n",
    "\n",
    "\n",
    "#### 4.4 Removal of Highly Genetically Related Samples\n",
    "\n",
    "> A required feature of association studies with cohorts sampled from a general population is the removal of genetically identical individuals in the cohort. Genetically identical samples are detrimental to the overall power of further association tests, as the individual is overrepresented in the study cohort. Additionally, genetically identical samples can be the product of poor sample handling and experimental error. Below we outline steps to remove highly related individuals with plink. \n",
    "\n",
    "##### Preparing Data for IBD Computation:\n",
    "\n",
    "``` bash\n",
    "#per sample and per variant missigness rate \n",
    "plink --bfile hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het --missing --out hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het\n",
    "```\n",
    "\n",
    "> Identity-by-state (IBS) is a computation used to dissect the degree of genetic relatedness between individuals by considering the average proportion of shared alleles at common polymorphic sites. When implementing identity-by-state based computations of genetic relatedness, it is essential to first filter out regions of high linkage disequilibrium, like the MHC, in order to characterize individuals with independent SNPs. Identity-by-descent (IBD) is a metric derived from genome wide IBS computations and quantifies the degree of recent shared ancestry. The following commands computes IBD per pair of individuals while filtering out the MHC and other high LD SNPs, and identifies as well as removes highly related individuals. \n",
    "    \n",
    "``` bash\n",
    "HLA_START=24000000\n",
    "HLA_STOP=36000000\n",
    "\n",
    "#identify MHC and flanking snps \n",
    "gawk '$1==6 && $4 > '${HLA_START}' && $4 < '${HLA_STOP}'{print $2}' hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.bim > hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.hla.snps\n",
    "\n",
    "#removal of MHC snps, application of MAF filter and LD pruning \n",
    "plink --bfile hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het --exclude hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.hla.snps --maf 0.05 --indep 50 5 2 --out hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.noHLA.maf\n",
    "```\n",
    "\n",
    "\n",
    "##### Filtering of Samples with High Sample Relatedness:\n",
    "\n",
    "``` bash\n",
    "#exclusion of high ld snps, ibd computation using --genome \n",
    "plink --bfile hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het --exclude hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.noHLA.maf.prune.out --genome --out hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het\n",
    "\n",
    "#Identification of samples with IBD statistic > 0.9 \n",
    "awk '{if($10>0.9)print}' hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.genome > hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.ibd.rem.txt\n",
    "\n",
    "#head command to show hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.genome fields\n",
    "head hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.genome\n",
    "```\n",
    "![Image](./images/hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.genome.png)\n",
    "\n",
    "```bash\n",
    "#python script that extracts id of samples to be removed\n",
    "python ../scripts/get_remID.py hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.imiss hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.ibd.rem.txt hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.ibd.rem.fam\n",
    "\n",
    "#plink removal of samples\n",
    "plink --bfile hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het --remove hgdp_chr6.dedup.ambstrandrem.1stSNPQC.mind.het.ibd.rem.fam --make-bed --out hgdp_chr6.1stSNPQC.1stSampQC\n",
    "```\n",
    "\n",
    "\n",
    "### 5. Variant QC\n",
    "\n",
    "**Corresponds to 'Per-variant QC' subsection**\n",
    "\n",
    "#### 5.1 Per Variant Missingness Rate\n",
    "\n",
    "> Genotype call rate thresholds often range between 95-99% percent, but are selected on a cohort by cohort basis. The initial SNP filtering of 10% represented the minimum quality required to profile individuals, but variants used as input for association studies must be of higher quality. Removed variants have the potential to be disease-relevant, so we suggest a distribution based approach for threshold selection that seeks to reasonably differentiate SNPs of poor quality from the rest of genotyped variants with good quality. The below selection outlines steps to conduct this filtering using plink and R: \n",
    "\n",
    "\n",
    "``` bash\n",
    "#Generate per variant missingness rate\n",
    "plink --bfile  hgdp_chr6.1stSNPQC.1stSampQC --missing --out hgdp_chr6.1stSNPQC.1stSampQC\n",
    "```\n",
    "\n",
    "in R: \n",
    "\n",
    "``` R\n",
    "library(ggplot2)\n",
    "name = 'hgdp_chr6.1stSNPQC.1stSampQC.lmiss'\n",
    "#read in missingness report\n",
    "dat = read.table(name, header = TRUE)\n",
    "#Plot histogram of variant missingness rate \n",
    "png(paste0(name, '.png'))\n",
    "ggplot(dat) +\n",
    "    geom_histogram(aes(x = F_MISS, y = ..count..), binwidth = 0.01) +\n",
    "    theme_classic() +\n",
    "    geom_vline(xintercept = 0.02, color = \"red\", linetype = \"dashed\") +\n",
    "    labs(x = \"Variant Missigness Rate\", y = \"Frequency\", title = \"1st SNP and SAMPLE QC\")\n",
    "dev.off()\n",
    "```\n",
    "\n",
    "![Image](./images/hgdp_all_chr6.hg19.ba.only.GSA.1stSNPQC.1stSampQC.lmiss.png)\n",
    "\n",
    "> This threshold is identified by selecting a missingness rate that cuts off outlier variants while preserving the bulk of the distribution. Here we selected a per variant missingness rate of 0.02. \n",
    "    \n",
    "``` bash\n",
    "#plink filtering based on chosen --geno threshold\n",
    "plink --bfile hgdp_chr6.1stSNPQC.1stSampQC --geno 0.02 --make-bed --out hgdp_chr6.1stSNPQC.1stSampQC.miss\n",
    "```\n",
    "    \n",
    "#### 5.2 Allele Frequency Comparison with 1KG\n",
    "\n",
    "> A useful metric for divergence of SNP behavior from an expected population average is to match allele frequencies by variant between the target dataset and 1KG, and compute the difference in allele frequency for matched SNPs between the two datasets. Since the HGDP project represents genotypes from many ancestries, here we matched with 1KG's multipopulation average. For less admixed or diverse cohorts, we suggest matching to a corresponding 1KG population (i.e. if the target set is mostly African, we reccommend using data from the 1KG African cohorts for comparison). The 1KG frequency file is derived from the 1000 genomes data downloaded from the [ftp site](http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/), and commands to create it are listed below. \n",
    "    \n",
    "``` bash \n",
    "#reset language for join steps\n",
    "LANG=C\n",
    "\n",
    "#write hwe summary statistics\n",
    "plink --bfile hgdp_chr6.1stSNPQC.1stSampQC.miss --hardy --allow-no-sex --out hgdp_chr6.1stSNPQC.1stSampQC.miss.frq.diff\n",
    "\n",
    "#head command to show fields of hgdp_chr6.1stSNPQC.1stSampQC.miss.frq.diff.hwe\n",
    "head hgdp_chr6.1stSNPQC.1stSampQC.miss.frq.diff.hwe \n",
    "```\n",
    "![Image](./images/tmpfile.hwe.png)\n",
    "\n",
    "```bash\n",
    "tmpfile='hgdp_chr6.1stSNPQC.1stSampQC.miss.frq.diff'\n",
    "BED='hgdp_chr6.1stSNPQC.1stSampQC.miss'\n",
    "#computes allele frequency of A1 allele per marker\n",
    "cat ${tmpfile}.hwe | awk '/ALL/{split($6,V,\"/\"); if(V[1]+V[2]+V[3]==0){Freq=0}else{Freq=(V[2]/2+V[1])/(V[1]+V[2]+V[3])}print $2, $4, $5, Freq}' > ${tmpfile}.hwe.freq\n",
    "\n",
    "#head command to show fields of ${tmpfile}.hwe.freq\n",
    "head ${tmpfile}.hwe.freq\n",
    "```\n",
    "![Image](./images/tmpfile.hwe.freq.png)\n",
    "\n",
    "\n",
    "```bash\n",
    "#joins frequency with bim file for chr and pos, renames snps to match 1KG convention, lists chr_pos, A1, A2, and defined allele frequency\n",
    "cat ${tmpfile}.hwe.freq | sort -k1,1 | join - <(cat ${BED}.bim | awk '{print $2, $1 \"_\" $4}' | sort -k1,1) | awk '{print $5, $2, $3, $4}' | sort -k1,1 > ${tmpfile}.chr_pos_allele_freq\n",
    "\n",
    "#head command to show fields of ${tmpfile}.chr_pos_allele_freq\n",
    "head ${tmpfile}.chr_pos_allele_freq\n",
    "```\n",
    "![Image](./images/tmpfile.chr_pos_allele_freq.png)\n",
    "\n",
    "```bash\n",
    "#joins previously generated allele freq file with 1KG allele freq file, matches A1 and A2 alleles\n",
    "join ${tmpfile}.chr_pos_allele_freq 1KGp3v5.Ref.Frq.chr_pos_allele | awk '{if($2==$5 && $3==$6){print $1, $4, $7}else{if($2==$6 && $3==$5){print $1, 1-$4, $7}}}' > ${tmpfile}.plot\n",
    "\n",
    "#head command to show fields of ${tmpfile}.plot\n",
    "head ${tmpfile}.plot\n",
    "```\n",
    "\n",
    "![Images](./images/tmpfile.plot.png)\n",
    "\n",
    "\n",
    "```bash\n",
    "#counts number of alleles with > 20% allele frequency between target and 1KG\n",
    "cat ${tmpfile}.plot | awk '{gos=$2-$3;if(gos>0.20||gos<-0.20){COUNT++;}}END{print COUNT \" / \" NR}'\n",
    "\n",
    "rm ${tmpfile}.chr_pos_allele_freq\n",
    "rm ${tmpfile}.hwe.freq\n",
    "rm ${tmpfile}.hwe\n",
    "\n",
    "#Make 1KG reference file using similar commands to above\n",
    "plink --vcf ALL.chr6.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz --make-bed --out 1KGp3v5\n",
    "plink --bfile 1KGp3v5 --hardy --allow-no-sex --out 1KGp3v5.miss.frq.diff\n",
    "cat 1KGp3v5.miss.frq.diff.hwe | awk '/ALL/{split($6,V,\"/\"); if(V[1]+V[2]+V[3]==0){Freq=0}else{Freq=(V[2]/2+V[1])/(V[1]+V[2]+V[3])}print $2, $4, $5, Freq}' > 1KGp3v5.miss.frq.diff.hwe.freq \n",
    "cat 1KGp3v5.miss.frq.diff.hwe.freq | sort -k1,1 | join - <(cat 1KGp3v5.bim | awk '{print $2, $1 \"_\" $4}' | sort -k1,1) | awk '{print $5, $2, $3, $4}' | sort -k1,1 > 1KGp3v5.Ref.Frq.chr_pos_allele\n",
    "```\n",
    "\n",
    "> Choosing the allele frequency difference threshold is done via plotting, demonstrated below: \n",
    "\n",
    "in R: \n",
    "\n",
    "```R\n",
    "library(ggplot2)\n",
    "#read in data\n",
    "name = 'hgdp_chr6.1stSNPQC.1stSampQC.miss.frq.diff.plot'\n",
    "dat = read.table(name)\n",
    "#make an AF by AF plot comparing target and 1KG datasets\n",
    "png(paste0(name, '.png'))\n",
    "ggplot(dat, aes(x = V2, y = V3)) +\n",
    "    geom_point() +\n",
    "    theme_classic() +\n",
    "    labs(x = \"Target Dataset Freq\", y = \"1KG Freq\") +\n",
    "    geom_abline(intercept = c(0.20, -0.20), size = 1.5, linetype = \"dashed\", color = \"red\")\n",
    "dev.off()\n",
    "```\n",
    "\n",
    "![Image](./images/hgdp_all_chr6.hg19.ba.only.GSA.1stSNPQC.1stSampQC.miss.frq.diff.plot.png)\n",
    "\n",
    "> **Important Note: When the population does not exactly match between the target and the database or the reference, this strategy might be ineffective within the MHC. MHC is known to have most highly variable AF across populations. Thus, we could suggest using a liberal threshold in removing variants based on the AF differences.**\n",
    "    \n",
    "### 6. Extracting the MHC: Final Input for Imputation\n",
    "\n",
    "> Here we isolate the MHC using the command below.  \n",
    "    \n",
    "``` bash\n",
    "#isolate MHC\n",
    "plink --bfile hgdp_chr6.1stSNPQC.1stSampQC.miss --chr 6 --from-mb 28 --to-mb 34 --make-bed --out hgdp_chr6.final\n",
    "```\n",
    "\n",
    "### 7. Variant Naming Conventions: Target and Reference Genotype Data\n",
    "> **Note**: While the Michigan Imputation Server and Minimac3 match input genotype data to reference haplotypes by genomic position, SNP2HLA utilizes variant names defined in the second column in the input .bim plink file. We provide a python script to update the input .bim file to match variant naming conventions, and recommend this step to simplify downstream analysis. \n",
    "\n",
    "```bash\n",
    "#rename variants \n",
    "mv hgdp_chr6.final.bim hgdp_chr6.final.bim.old\n",
    "python ../scripts/rename_bim.py Tutorial_1KGonly.bim hgdp_chr6.final.bim\n",
    "\n",
    "#head command showing final bim file fields\n",
    "head hgdp_chr6.final.bim\n",
    "```\n",
    "\n",
    "![Images](./images/hgdp_chr6.final.bim.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5d7a243",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Tools for Genotype Phasing and HLA Imputation\n",
    "\n",
    "**The following explains steps and provides examples for the 'Tools for genotype phasing and HLA imputation' section of Sakaue et. al 2022**\n",
    "\n",
    "\n",
    "> Here we demonstrate the use of two well known tools for imputation of genetic variants used in several studies, SNP2HLA and Minimac3. Minimac3 requires pre-phasing of the input genotype data, whereas SNP2HLA will phase input data with BEAGLE. While there are many tools for phasing genotype data, we reccommend the use of either EAGLE or SHAPEIT. EAGLE is optimized for use with very large sample sizes, while SHAPEIT performs better with a smaller cohort size. While we provide example commands for both EAGLE and SHAPEIT, we recommend use of EAGLE for cohorts larger than 10000 individuals. \n",
    "    \n",
    "    \n",
    "### 1. Genotype Phasing\n",
    "\n",
    "#### 1.1 EAGLE \n",
    "\n",
    "##### Inputs and Implementation:\n",
    "    \n",
    " > [EAGLE](https://alkesgroup.broadinstitute.org/Eagle/#x1-30002) is available as a compressed archive (.tar.gz) from the following [link](https://alkesgroup.broadinstitute.org/Eagle/downloads/). Included with the EAGLE package is the `eagle` executable, and `/tables/genetic_map_hg19_withX.txt.gz` - the genetic maps required as an input. EAGLE can phase with or without a reference panel, but requires the following for imputation: the target genotype data in plink, .vcf, or .bcf format, and the previously mentioned genetic maps which are derived from HapMap's publically available genetic maps. The following command acts on the previously qcd genotype data.\n",
    "    \n",
    "```bash\n",
    "Eagle_v2.4.1/eagle --bfile hgdp_all_chr6.hg19.ba.only.GSA.final --geneticMapFile Eagle_v2.4.1/tables/genetic_map_hg19_withX.txt.gz --outPrefix hgdp_chr6.final.EAGLE.phased --numThreads 8\n",
    "```\n",
    "\n",
    "##### Output Formatting: \n",
    "\n",
    "> If plink inputs are used, EAGLE will output haplotypes in the Oxford HAPS/SAMPLE format which is compatible with SHAPEIT2. If a vcf/bcf output is desired from EAGLE directly, addition of the following option to the above command will result in compressed vcf (*.vcf.gz) output: `--vcfOutFormat=z`. If an uncompressed vcf is needed, addition of `--vcfOutFormat` to the above command will result in an .vcf output.  Conversion of the output to vcf/bcf file and optional compression can be achieved via the following command using SHAPEIT and samtools: \n",
    "\n",
    "```bash\n",
    "shapeit -convert --input-haps hgdp_chr6.final.EAGLE.phased --output-vcf hgdp_chr6.final.EAGLE.phased.vcf\n",
    "bgzip -c hgdp_chr6.final.EAGLE.phased.vcf > hgdp_chr6.final.EAGLE.phased.vcf.gz\n",
    "tabix hgdp_chr6.final.EAGLE.phased.vcf.gz\n",
    "```\n",
    "\n",
    "\n",
    "#### 1.2 SHAPEIT\n",
    "\n",
    "##### Inputs and Implementation: \n",
    "\n",
    "> [SHAPEIT](https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html#home) is our recommended tool to phase genotypes for this example dataset, and is downloadable as a precompiled binary in a compressed tar archive format. SHAPEIT requires a genetic map, derived from recombination rates found by the HapMap consortium, which can be downloaded [here](https://ftp.ncbi.nlm.nih.gov/hapmap/recombination/2011-01_phaseII_B37/). We have included the genetic map for chromosome 6, and the `shapeit` executable is available when downloading from the SHAPEIT website linked previously. SHAPEIT will accept input files in the following formats: PLINK PED/MAP, PLINK BED/BIM/FAM, OXFORD GEN/SAMPLE, and VCF. SHAPEIT additionally requires input data to be split by chromosome, as well as a corresponding genetic map for the chromosome being phased. Here we phased chromsome 6 of the qcd genotype data with the following command: \n",
    "    \n",
    "```bash\n",
    "shapeit -B hgdp_chr6.final -M genetic_map_chr6_combined_b37.txt -O hgdp_chr6.final.shapeit.phased --thread 8 --seed 0 --output-log hgdp_chr6.final.shapeit.phased.log\n",
    "```\n",
    "\n",
    "##### Output Formatting: \n",
    "\n",
    "> SHAPEIT outputs phased haplotypes in HAPS/SAMPLE format. The package includes functionality for conversion of HAPS/SAMPLE output to vcf/bcf for phased genotype parsing by other programs, such as PLINK. An example command for conversion, and compression of the converted vcf using samtools, is shown below:\n",
    "    \n",
    "```bash \n",
    "shapeit -convert --input-haps hgdp_chr6.final.shapeit.phased --output-vcf hgdp_chr6.final.shapeit.phased.vcf\n",
    "bgzip -c hgdp_chr6.hg19.final.shapeit.phased.vcf > hgdp_chr6.final.shapeit.phased.vcf.gz\n",
    "tabix hgdp_chr6.final.shapeit.phased.vcf.gz\n",
    "```\n",
    "\n",
    "\n",
    "### 2. HLA Imputation\n",
    "    \n",
    "> Two commonly available tools for imputation are SNP2HLA and Minimac3, as well as instructions on how to format input genotype data for imputation with the Michigan Imputation Server, are discussed in the following section. SNP2HLA is available as both a C shell script that can be implemented in linux or C shell based distributions, or as a module in the HLA-TAPAS package. Here we describe use of SNP2HLA in both contexts. \n",
    "\n",
    "    \n",
    "#### 2.1 SNP2HLA.csh\n",
    "\n",
    "**Corresponds to 'SNP2HLA' subsection**\n",
    "\n",
    "\n",
    "##### Additional Dependencies:\n",
    "\n",
    "> Implementation of the C shell version of SNP2HLA requires several dependencies post installation and unpacking of the SNP2HLA tarball from the [SNP2HLA website](https://software.broadinstitute.org/mpg/snp2hla/). Several additional dependencies are required to be present in the SNP2HLA directory in the package: the [plink/1.0.7](https://zzz.bwh.harvard.edu/plink/download.shtml) run file, and the beagle.jar, linkage2beagle.jar, beagle2linkage.jar files from [Beagle 3.0.4](http://faculty.washington.edu/browning/beagle/beagle.html#download). \n",
    "\n",
    "##### Inputs and Implementation: \n",
    "\n",
    "> The C shell distribution of SNP2HLA requires the target genotype data in plink bim/bed/fam format, and the reference panel in Beagle format (*.bgl.phased and .markers). Generating the correct reference panel inputs from the provided reference vcf file can be done with the vcf2beagle.jar command from beagle, an example of which is shown below:\n",
    "``` bash\n",
    "zcat Tutorial_1KGonly.bgl.phased.vcf.gz | java -jar vcf2beagle.jar ? Tutorial_1KGonly\n",
    "gunzip Tutorial_1KGonly.bgl.gz \n",
    "```\n",
    "\n",
    "> Imputing with SNP2HLA.csh is done with a command in the following format - `./SNP2HLAcsh input_plink_prefix reference_panel_prefix output_prefix ./plink`. The following command serves as an example of imputation with SNP2HLA.csh if housed in linux distribution: \n",
    "\n",
    "```bash\n",
    "cd /SNP2HLA_package_v1.0.3/SNP2HLA/\n",
    "\n",
    "./SNP2HLA.csh hgdp_chr6.final Tutorial_1KGonly hgdp_chr6.final.SNP2HLA.bash.imputed ./plink 2000 1000\n",
    "```\n",
    "##### Outputs:\n",
    "\n",
    "> SNP2HLA.csh outputs imputed genotype data in the following files: {output}.dosage (imputed allele dosage data, rows are markers and columns are individuals) which is recommended for downstream analysis, {output}.bgl.phased (imputed best guess genotypes in beagle phased format), {output}.[bim/bam/bed] (PLINK format best guess genotypes), {output}.bgl.gprobs (imputation posterior probabilities), and {output}.bgl.r2 (predicted r2 with true genotypes). \n",
    "\n",
    "\n",
    "#### 2.2 SNP2HLA.py\n",
    "\n",
    "**Corresponds to 'SNP2HLA' subsection**\n",
    "\n",
    "\n",
    "> SNP2HLA.py is a script derived through the [HLA-TAPAS package](https://github.com/immunogenomics/HLA-TAPAS) which we modified for edge cases and additional functionality. SNP2HLA.py is our recommended method of imputation with SNP2HLA due to ease of use and imputation speed compared to the previous C shell implementation. \n",
    "\n",
    "##### Additional Dependencies:\n",
    "\n",
    "> The following dependencies are required for imputation with SNP2HLA.py, and must be prepared in the 'dependency/' folder of the HLA-TAPAS directory: [PLINK v1.9b](https://www.cog-genomics.org/plink2/), [BEAGLE v4.1](https://faculty.washington.edu/browning/beagle/b4_1.html#download) (must be renamed to beagle.jar), [beagle2vcf.jar, linkage2beagle.jar, and vcf2beagle.jar](https://faculty.washington.edu/browning/beagle_utilities/utilities.html). \n",
    "\n",
    "\n",
    "##### Inputs and Implementation:\n",
    "\n",
    "> SNP2HLA.py requires input genotype data in PLINK .bim/.bed/.fam format, as well as several input reference panel formats. Reference panels must be prepared in the following formats: .bgl.phased.vcf.gz, .markers, .FRQ.frq, .bim. These files can be obtained from the HLA-TAPAS MakeReference module. The following commands will also produce all required inputs from a beagle phased, compressed vcf format reference panel file: \n",
    "\n",
    "``` bash\n",
    "plink --vcf Tutorial_1KGonly.bgl.phased.vcf.gz --keep-allele-order --make-bed --out Tutorial_1KGonly\n",
    "plink --bfile Tutorial_1KGonly --freq --out Tutorial_1KGonly.FRQ\n",
    "zcat Tutorial_1KGonly.bgl.phased.vcf.gz | grep -v \"#\" | awk '{print $3,$2,$4,$5}' > Tutorial_1KGonly.markers\n",
    "```\n",
    "\n",
    "\n",
    "> Following preparation of input files, the following command can be used to impute with SNP2HLA.py. We have provided SNP2HLA.py and associated dependencies in the /scripts/ folder of this github repository\n",
    "\n",
    "```bash\n",
    "cd /scripts/\n",
    "\n",
    "python ./SNP2HLA.py \\\n",
    "  -i hgdp_chr6.final \\\n",
    "  -o hgdp_chr6.final.SNP2HLApy.imputed \\\n",
    "  -rf Tutorial_1KGonly \\\n",
    "  --nthreads 10 \\\n",
    "  --java-mem=80g --tolerated-diff=0.5\n",
    "```\n",
    "##### Outputs:\n",
    "\n",
    "> SNP2HLA.py outputs a single phased and imputed vcf file: {output}.bgl.phased.vcf.gz, that can be recoded to a tabular, raw dosage format using plink. This can be accomplished via the following command, where {output}_allele_ids is a file with variant ids and the corresponding desired allele to be counted: \n",
    "\n",
    "```bash\n",
    "#Creates allele export file that specifies ALT allele to be counted\n",
    "zless -S hgdp_chr6.final.SNP2HLApy.imputed.bgl.phased.vcf.gz | grep -v \"#\" | awk '{print $3, $5}' > hgdp_chr6.final.SNP2HLApy.imputed_allele_ids\n",
    "#Exports to tabular format \n",
    "plink2 --vcf hgdp_chr6.final.SNP2HLApy.imputed.bgl.phased.vcf.gz dosage=DS --export A --export-allele hgdp_chr6.final.SNP2HLApy.imputed_allele_ids --out hgdp_chr6.final.SNP2HLApy.imputed\n",
    "\n",
    "#head command for allele export file fields\n",
    "head hgdp_chr6.final.SNP2HLApy.imputed_allele_ids\n",
    "```\n",
    "\n",
    "![Image](./images/hgdp_chr6.final.SNP2HLApy.imputed_allele_ids.png)\n",
    "\n",
    "```bash\n",
    "#command to show exported tabular output fields\n",
    "less -S hgdp_chr6.final.SNP2HLApy.imputed.raw \n",
    "```\n",
    "![Image](./images/hgdp_chr6.final.SNP2HLApy.imputed.raw.png)\n",
    "\n",
    "#### 2.3 Minimac3\n",
    "\n",
    "> We note that the original SNP2HLA implementation using BEAGLE does not scale to a large number of samples in the target dataset, especially at a study size larger 10,000. To address this, we also provide a pipeline using another representative imputation software, Minimac3, which can scale to hundreds of thousands or millions of individuals. [Minimac3](https://genome.sph.umich.edu/wiki/Minimac3_Usage#Introduction) is an updated implementation of the widely used imputation tools Minimac2 and Minimac1. It is available as both a standard executable, and as an openMP programming enabled version. We recommend use of the latter for faster parameter estimation and imputation.  \n",
    "\n",
    "##### Inputs and Implementation:\n",
    "\n",
    "> Input files to Minimac3 required for imputation include pre-phased gwas data and a reference panel. After phasing, protocols for which described in above section 1. Genotype Phasing, phased files must be converted to vcf format. Reference panels must be supplied in vcf format or m3vcf format. Minimac3 commands are formatted as follows: `minimac3 --refHaps reference_panel_vcf --haps target_genotype_vcf --prefix output_file_prefix --chr chromsome_imputed`. Additional options include --cpus for the number of cpus to be used for the imputation task, and --doseOutput for output of imputed dosage probabilities rather than best guess genotype calls. An example command for Minimac3 is displayed below: \n",
    "\n",
    "```bash\n",
    "Minimac3-omp --refHaps Tutorial_1KGonly.bgl.phased.vcf.gz --haps hgdp_chr6.final.EAGLE.phased.vcf.gz --prefix hgdp_chr6.final.EAGLE.phased.imputed --cpus 6 --chr 6 --doseOutput\n",
    "```\n",
    "\n",
    "##### Outputs: \n",
    "\n",
    "> Minimac3 outputs in .vcf format by default. The addition of the `--doseOutput` flag in the example command above results in dosage probabilities written to output files. The above command will output 6 files: {output}.rec (switchrate per marker), {output}.erate (error rate per marker), {output}.m3vcf.gz (best guess genotypes in m3vcf format), {output}.dose.gz (tabular dosage output of samples by markers), {output}.info (imputation statistics per marker), and {output}.dose.vcf.gz (vcf of dosage per marker per sample). For downstream analysis, we recommend conversion of {output}.dose.vcf.gz to tabular dosage format via the plink command in section '2.2 SNP2HLA.py', subsection 'Outputs:'. \n",
    "\n",
    "#### 2.4 Michigan Imputation Server\n",
    "\n",
    "**Corresponds to 'Michigan Imputation Server' subsection**\n",
    "\n",
    "\n",
    "> The [Michigan Imputation Server](https://imputationserver.sph.umich.edu/index.html#!pages/home) is a web-based tool for imputation that utilizes Minimac4 and provides a variety of reference panels for imputation. HLA imputation at one-field(two-digit) and two-field(four-digit) resolution with the Michigan Imputation Server can be accomplished with the following available references: [Four-digit Multi-ethnic HLA v1, and Four-digit Multi-ethnic HLA v2](https://imputationserver.readthedocs.io/en/latest/reference-panels/). We reccommend use of Four-digit Multi-ethnic HLA v2 due to higher imputation accuracy and better representation of East Asian populations. The following section describes how to prepare data for upload to the Michigan Imputation Server. \n",
    "\n",
    "##### Data Preparation\n",
    "\n",
    "> Input data must be provided in a sorted vcf format per chromosome, and can either be phased or unphased. If unphased, the 'Phasing' option when running a job must be selected. Sections of this vignette that discuss the necessary quality control, phasing, and file conversion steps are QUALITY CONTROL OF TARGET GENOTYPE DATA, and 1. Genotype Phasing. An outline of imputation with the Michigan Imputation Server is shown below: \n",
    "\n",
    "![Image](./images/SuppleFig2_MIS_usage.png)\n",
    "\n",
    "\n",
    "### 3. Runtime benchmarking\n",
    "\n",
    "> Below, we show the runtime benchmarking of HLA imputation using different platforms. **a.** For SNP2HLA, we used BEAGLE4 for imputation algorithm (Luo et al. Nat Genet 2021) with using 10 CPUs. For benchmarking of imputation by minmac4 (MM4), we used SHAPEIT2 as phasing algorithm with samples < 10,000 and EAGLE2 as phasing algorithm with samples > 5,000 as we described in the manuscript with using 10 CPUs. **b.** For Michigan Imputation Server, we uploaded the unphased genotype data and standard imputation pipeline with default setting was performed (with 1CPU).\n",
    "\n",
    "![Image](./images/SuppleFig4_imputation_runtime.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e81d3af",
   "metadata": {},
   "source": [
    "## Post-imputation QC\n",
    "\n",
    "**The following explains steps and provides examples for the 'Post-imputation QC' section of Sakaue et. al 2022**\n",
    "\n",
    "> Before supplying imputed data for association testing, the imputation results must be subjected to another round of quality control to ensure only high quality variants are included for statistical tests. We reccommend filtering on imputation R2, and while the choice of R2 is dataset dependent, R2 thresholds of 0.7  are widely used to filter out poorly imputed variants in single cohort studies. The below command serves as an example to quality control imputation results and output high quality variants in compressed vcf format using bcftools:\n",
    "\n",
    "```bash\n",
    "zcat hgdp_chr6.final.SHAPEIT.imputed.dose.vcf.gz | sed 's/PASS;GENOTYPED/PASS/' | bcftools view  -i 'R2>.7' -Oz -o hgdp_chr6.final.SHAPEIT.imputed.R20.7.dose.vcf.gz \n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}