docs/example/bi-parental-pedigree.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8bbf4feb-901d-4c5d-bcde-00d6d3be5f8e",
   "metadata": {},
   "source": [
    "# Pedigree informed micro-haplotype calling in an autopolyploid bi-parental cross\n",
    "\n",
    "*(Last updated for MCHap version 0.10.0)*\n",
    "\n",
    "This notebook demonstrates pedigree informed micro-haplotype calling in a small bi-parental cross.\n",
    "\n",
    "**WARNING: The `call-pedigree` tool is considered highly experimental**!!!\n",
    "\n",
    "The following topics are introduced:\n",
    "\n",
    "- **Pooled de novo haplotype assembly with `mchap assemble`**\n",
    "  - Input files\n",
    "  - Sample pooling\n",
    "- **Re-calling genotypes with `mchap call`**\n",
    "  - Specifying a population prior on allele frequencies\n",
    "- **Re-calling genotypes with `mchap call-pedigree`**\n",
    "  - Specifying pedigree to inform haplotype calling\n",
    "\n",
    "**Software requirements:**\n",
    "\n",
    "This notebook uses the [bash-kernel](https://github.com/takluyver/bash_kernel) for Jupyter which can be installed with `pip install bash_kernel`.\n",
    "Alternatively, the code from this notebook may be run in a unix-like bash environment.\n",
    "\n",
    "In addition to using the MCHap software, these examples also require the `bgzip` and `tabix` tools which are part of [`htslib`](https://github.com/samtools/htslib).\n",
    "\n",
    "**Data sources:**\n",
    "\n",
    "The bi-parental population used within this notebook is a small subset of the population published by:\n",
    "J Tahir, C Brendolise, S Hoyte, M Lucas, S Thomson, K Hoeata, C McKenzie, A Wotton, K Funnell, E Morgan, D Hedderley, D Chagné, P M Bourke, J McCallum, S E Gardiner,* and L Gea \n",
    "“QTL Mapping for Resistance to Cankers Induced by Pseudomonas syringae pv. actinidiae (Psa) in a Tetraploid Actinidia chinensis Kiwifruit Population” \n",
    "Pathogens 2020, 9, 967; doi:10.3390/pathogens9110967\n",
    "\n",
    "- Raw sequences: [10.5281/zenodo.4285665](https://zenodo.org/record/4285666) and [10.5281/zenodo.4287636](https://zenodo.org/record/4287637)\n",
    "- Reference genome: DOI [10.5281/zenodo.5717386](https://zenodo.org/record/5717387) (chromosome 1 only)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6917e307-5ed1-4fd5-8fbe-3df393e24c5c",
   "metadata": {},
   "source": [
    "## De novo haplotype assembly with `mchap assemble`\n",
    "\n",
    "*For more background information on* `mchap assemble` *see the online [documentation](https://github.com/PlantandFoodResearch/MCHap/blob/master/docs/assemble.rst).*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef434b09-0305-49b2-9ec1-6190a604fa76",
   "metadata": {},
   "source": [
    "### Input files\n",
    "\n",
    "The required input files have been organised by file type:\n",
    "\n",
    "- `input/bam` BAM alignment files for the parents and progeny samples\n",
    "- `input/bed` Target loci for haplotype assembly\n",
    "- `input/vcf` VCF file of \"basis\" SNVs\n",
    "- `input/fasta`: Reference genome"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a3a629d-6067-4375-bea0-d81a9f80de3d",
   "metadata": {},
   "source": [
    "In this example we use one BAM file per sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3dda15df-9fa3-4b38-9b0c-c16df023dd99",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "convert_to_cram.sh       \u001b[0m\u001b[01;32mprogeny006.loci.bam\u001b[0m      \u001b[01;32mprogeny013.loci.bam.bai\u001b[0m\n",
      "\u001b[01;32mparent1.loci.bam\u001b[0m         \u001b[01;32mprogeny006.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny014.loci.bam\u001b[0m\n",
      "\u001b[01;32mparent1.loci.bam.bai\u001b[0m     \u001b[01;32mprogeny007.loci.bam\u001b[0m      \u001b[01;32mprogeny014.loci.bam.bai\u001b[0m\n",
      "\u001b[01;32mparent2.loci.bam\u001b[0m         \u001b[01;32mprogeny007.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny015.loci.bam\u001b[0m\n",
      "\u001b[01;32mparent2.loci.bam.bai\u001b[0m     \u001b[01;32mprogeny008.loci.bam\u001b[0m      \u001b[01;32mprogeny015.loci.bam.bai\u001b[0m\n",
      "\u001b[01;32mprogeny001.loci.bam\u001b[0m      \u001b[01;32mprogeny008.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny016.loci.bam\u001b[0m\n",
      "\u001b[01;32mprogeny001.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny009.loci.bam\u001b[0m      \u001b[01;32mprogeny016.loci.bam.bai\u001b[0m\n",
      "\u001b[01;32mprogeny002.loci.bam\u001b[0m      \u001b[01;32mprogeny009.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny017.loci.bam\u001b[0m\n",
      "\u001b[01;32mprogeny002.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny010.loci.bam\u001b[0m      \u001b[01;32mprogeny017.loci.bam.bai\u001b[0m\n",
      "\u001b[01;32mprogeny003.loci.bam\u001b[0m      \u001b[01;32mprogeny010.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny018.loci.bam\u001b[0m\n",
      "\u001b[01;32mprogeny003.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny011.loci.bam\u001b[0m      \u001b[01;32mprogeny018.loci.bam.bai\u001b[0m\n",
      "\u001b[01;32mprogeny004.loci.bam\u001b[0m      \u001b[01;32mprogeny011.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny019.loci.bam\u001b[0m\n",
      "\u001b[01;32mprogeny004.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny012.loci.bam\u001b[0m      \u001b[01;32mprogeny019.loci.bam.bai\u001b[0m\n",
      "\u001b[01;32mprogeny005.loci.bam\u001b[0m      \u001b[01;32mprogeny012.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny020.loci.bam\u001b[0m\n",
      "\u001b[01;32mprogeny005.loci.bam.bai\u001b[0m  \u001b[01;32mprogeny013.loci.bam\u001b[0m      \u001b[01;32mprogeny020.loci.bam.bai\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "ls input/bam"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "872b7113-9697-4b89-9d9b-ff2b17bdd6aa",
   "metadata": {},
   "source": [
    "The BED file specifies the genomic coordinates of assembly targets. The fourth column (loci ID) of the BED file is optional but, if present, the loci IDs will be included in the output VCF file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e0990cdb-39a8-4c2c-bd6f-ff5056f6243f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr1\t17590\t17709\tlocus001\n",
      "chr1\t568848\t568967\tlocus012\n",
      "chr1\t684808\t684927\tlocus016\n",
      "chr1\t809104\t809223\tlocus019\n"
     ]
    }
   ],
   "source": [
    "cat input/bed/targets4.bed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c86c9d66-930d-4d91-8fca-4bd38574312a",
   "metadata": {},
   "source": [
    "The input VCF file is used to specify basis SNVs for haplotype alleles. These may be multi-allelic SNVs. Any sample data within this VCF file will be ignored:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae59e2e6-5892-4eda-b571-4c0df85ca3e1",
   "metadata": {},
   "source": [
    "### Identify putative SNVs\n",
    "\n",
    "MCHap includes the `find-snvs` tool for identifying putative SNVs.\n",
    "This is a fast but simplistic approach to identifying potential basis variants for assembling into haplotypes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1442f65e-2749-4879-b819-aa9646dd12d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "mchap find-snvs \\\n",
    "    --bam input/bam/*.bam \\\n",
    "    --reference input/fasta/chr1.fa.gz \\\n",
    "    --targets input/bed/targets4.bed \\\n",
    "    --ind-maf 0.1 \\\n",
    "    --ind-mad 3 \\\n",
    "    --min-ind 2 \\\n",
    "    | bgzip > putative_snvs.vcf.gz\n",
    "tabix -p vcf putative_snvs.vcf.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15dc29b1-0e25-4408-a003-84f20b6d8e5b",
   "metadata": {},
   "source": [
    "Notes:\n",
    "- The `--targets` parameter should be a bed file defining genomic interval to search for putative SNVs\n",
    "- The `--ind-mad` parameter specifies a (per sample) minor allele depths \n",
    "- The `--ind-maf` parameter specifies a (per sample) minor allele frequency (calculated from allele depths)\n",
    "- The `--min-ind` parameter specifies the minimum number of sample required to meet the above conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9b6230e-42a3-4b20-b020-d101670ca7a0",
   "metadata": {},
   "source": [
    "The output of `mchap find-snvs` is a VCF without genotype calls.\n",
    "Instead, sample allele depths are reported.\n",
    "The total allele depths and mean of individual frequencies are also reported in an `INFO` field:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4ff6bbf6-5705-4a6f-ab6b-8728da8ea211",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##fileformat=VCFv4.3\n",
      "##fileDate=20240910\n",
      "##source=mchap v0.9.4.dev72+g719498f.d20240710\n",
      "##commandline=\"/home/cfltxm/mambaforge/envs/mchap/bin/mchap find-snvs --bam input/bam/parent1.loci.bam input/bam/parent2.loci.bam input/bam/progeny001.loci.bam input/bam/progeny002.loci.bam input/bam/progeny003.loci.bam input/bam/progeny004.loci.bam input/bam/progeny005.loci.bam input/bam/progeny006.loci.bam input/bam/progeny007.loci.bam input/bam/progeny008.loci.bam input/bam/progeny009.loci.bam input/bam/progeny010.loci.bam input/bam/progeny011.loci.bam input/bam/progeny012.loci.bam input/bam/progeny013.loci.bam input/bam/progeny014.loci.bam input/bam/progeny015.loci.bam input/bam/progeny016.loci.bam input/bam/progeny017.loci.bam input/bam/progeny018.loci.bam input/bam/progeny019.loci.bam input/bam/progeny020.loci.bam --reference input/fasta/chr1.fa.gz --targets input/bed/targets4.bed --ind-maf 0.1 --ind-mad 3 --min-ind 2\"\n",
      "##reference=file:input/fasta/chr1.fa.gz\n",
      "##contig=<ID=chr1,length=21898217>\n",
      "##INFO=<ID=REFMASKED,Number=0,Type=Flag,Description=\"Reference allele is masked\">\n",
      "##INFO=<ID=AD,Number=R,Type=Integer,Description=\"Total read depth for each allele\">\n",
      "##INFO=<ID=ADMF,Number=R,Type=Float,Description=\"Mean of sample allele frequencies calculated from read depth\">\n",
      "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n",
      "##FORMAT=<ID=AD,Number=R,Type=Integer,Description=\"Read depth for each allele\">\n",
      "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tparent1\tparent2\tprogeny001\tprogeny002\tprogeny003\tprogeny004\tprogeny005\tprogeny006\tprogeny007\tprogeny008\tprogeny009\tprogeny010\tprogeny011\tprogeny012\tprogeny013\tprogeny014\tprogeny015\tprogeny016\tprogeny017\tprogeny018\tprogeny019\tprogeny020\n",
      "chr1\t17591\t.\tG\tA\t.\t.\tAD=449,81;ADMF=0.842,0.158\tGT:AD\t.:16,7\t.:17,3\t.:25,4\t.:6,4\t.:18,3\t.:20,8\t.:17,1\t.:15,2\t.:26,4\t.:26,5\t.:18,6\t.:17,2\t.:18,1\t.:25,3\t.:32,5\t.:19,3\t.:16,6\t.:20,3\t.:23,6\t.:30,2\t.:21,1\t.:24,2\n",
      "chr1\t17602\t.\tG\tA\t.\t.\tAD=478,78;ADMF=0.849,0.151\tGT:AD\t.:18,9\t.:17,3\t.:25,4\t.:4,4\t.:19,3\t.:21,7\t.:16,1\t.:16,2\t.:28,4\t.:29,5\t.:22,6\t.:15,2\t.:19,1\t.:29,3\t.:32,2\t.:19,3\t.:18,6\t.:23,3\t.:24,6\t.:33,2\t.:25,1\t.:26,1\n",
      "chr1\t17612\t.\tA\tG\t.\t.\tAD=483,81;ADMF=0.846,0.154\tGT:AD\t.:18,8\t.:18,3\t.:28,3\t.:4,4\t.:19,3\t.:22,8\t.:18,1\t.:15,3\t.:29,4\t.:27,5\t.:22,7\t.:17,1\t.:18,1\t.:29,4\t.:35,2\t.:19,3\t.:18,7\t.:23,2\t.:21,6\t.:34,3\t.:23,1\t.:26,2\n",
      "chr1\t17616\t.\tA\tG\t.\t.\tAD=479,81;ADMF=0.843,0.157\tGT:AD\t.:17,8\t.:16,3\t.:27,3\t.:4,4\t.:18,3\t.:23,8\t.:17,2\t.:13,3\t.:31,4\t.:27,5\t.:21,5\t.:17,2\t.:20,1\t.:28,3\t.:34,2\t.:20,3\t.:18,7\t.:23,2\t.:24,6\t.:34,3\t.:22,1\t.:25,3\n",
      "chr1\t17624\t.\tC\tT\t.\t.\tAD=536,60;ADMF=0.895,0.103\tGT:AD\t.:24,4\t.:20,0\t.:24,7\t.:7,2\t.:20,5\t.:26,3\t.:19,0\t.:13,4\t.:30,7\t.:36,0\t.:28,0\t.:20,0\t.:22,4\t.:32,0\t.:37,0\t.:18,7\t.:19,6\t.:26,0\t.:30,1\t.:30,10\t.:27,0\t.:28,0\n",
      "chr1\t17625\t.\tG\tC\t.\t.\tAD=519,79;ADMF=0.856,0.144\tGT:AD\t.:18,8\t.:17,3\t.:29,3\t.:5,4\t.:22,3\t.:21,8\t.:17,2\t.:16,2\t.:37,4\t.:30,5\t.:21,4\t.:16,2\t.:23,1\t.:30,3\t.:36,2\t.:22,3\t.:19,7\t.:25,2\t.:26,6\t.:36,3\t.:27,1\t.:26,3\n",
      "chr1\t17629\t.\tG\tT\t.\t.\tAD=534,74;ADMF=0.865,0.135\tGT:AD\t.:20,8\t.:17,3\t.:30,3\t.:5,4\t.:22,3\t.:21,8\t.:19,1\t.:16,3\t.:36,2\t.:29,4\t.:23,5\t.:17,1\t.:26,1\t.:29,3\t.:41,2\t.:23,3\t.:18,7\t.:27,2\t.:26,5\t.:34,3\t.:27,1\t.:28,2\n",
      "chr1\t17655\t.\tG\tA\t.\t.\tAD=537,89;ADMF=0.849,0.15\tGT:AD\t.:22,10\t.:18,3\t.:28,3\t.:9,5\t.:25,3\t.:18,6\t.:25,3\t.:17,3\t.:35,3\t.:23,4\t.:20,6\t.:15,3\t.:28,3\t.:28,4\t.:36,4\t.:25,2\t.:22,5\t.:24,2\t.:27,8\t.:33,3\t.:30,1\t.:29,5\n"
     ]
    }
   ],
   "source": [
    "zcat putative_snvs.vcf.gz | head -n 20"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee86267b-e05f-47e1-878e-56ab416eb3c9",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Pooled assembly\n",
    "\n",
    "For this tutorial we will jump straight into a pooled assembly. A more beginner friendly example can be found in the standard MCHap [bi-parental example notebook](https://github.com/PlantandFoodResearch/MCHap/blob/master/docs/example/bi-parental.ipynb).\n",
    "\n",
    "Sample pools can be defined using a tabular file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "599ec80a-e0ea-4c96-9d29-300fa7666299",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parent1\tPOOL\n",
      "parent2\tPOOL\n",
      "progeny001\tPOOL\n",
      "progeny002\tPOOL\n",
      "progeny003\tPOOL\n",
      "progeny004\tPOOL\n",
      "progeny005\tPOOL\n",
      "progeny006\tPOOL\n",
      "progeny007\tPOOL\n",
      "progeny008\tPOOL\n",
      "progeny009\tPOOL\n",
      "progeny010\tPOOL\n",
      "progeny011\tPOOL\n",
      "progeny012\tPOOL\n",
      "progeny013\tPOOL\n",
      "progeny014\tPOOL\n",
      "progeny015\tPOOL\n",
      "progeny016\tPOOL\n",
      "progeny017\tPOOL\n",
      "progeny018\tPOOL\n",
      "progeny019\tPOOL\n",
      "progeny020\tPOOL\n"
     ]
    }
   ],
   "source": [
    "cat input/pools/sample-pools.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2198a916-264f-4215-8ebe-272c538df7b1",
   "metadata": {},
   "source": [
    "In the file shown above we assign all of the samples to a pool called \"POOL\".\n",
    "However, you can imagine a more complex scheme if out example data contained more than one bi-parental cross.\n",
    "\n",
    "We can then run the pooled assembly with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f9abbbf8-f143-482c-a26b-65fa92f9cc10",
   "metadata": {},
   "outputs": [],
   "source": [
    "mchap assemble \\\n",
    "        --bam input/bam/*.bam \\\n",
    "        --targets input/bed/targets4.bed \\\n",
    "        --variants putative_snvs.vcf.gz \\\n",
    "        --reference input/fasta/chr1.fa.gz \\\n",
    "        --sample-pool input/pools/sample-pools.txt \\\n",
    "        --ploidy 8 \\\n",
    "        --report AFP AOP \\\n",
    "        | bgzip > pooled_assembly.vcf.gz\n",
    "tabix -p vcf pooled_assembly.vcf.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05f4350d-9ca9-4949-a787-a9022aad6558",
   "metadata": {},
   "source": [
    "Notes:\n",
    "- The `--bam` argument may be file with a list of BAM paths ([documentation](https://github.com/PlantandFoodResearch/MCHap/blob/master/docs/assemble.rst#analyzing-many-samples))\n",
    "- The `--ploidy` and `--inbreeding` values can be specified per sample using a simple tabular file ([documentation](https://github.com/PlantandFoodResearch/MCHap/blob/master/docs/assemble.rst#sample-parameters))\n",
    "- The `--inbreeding` argument will default to `0` which may be unrealistic and it's often better to guess a \"sensible\" value ([documentation](https://github.com/PlantandFoodResearch/MCHap/blob/master/docs/assemble.rst#sample-parameters))\n",
    "- The `--sample-pool` argument is (optionally) used to define sample pools ([documentation](https://github.com/PlantandFoodResearch/MCHap/blob/master/docs/assemble.rst#sample-pooling))\n",
    "- The `--report AFP AOP` line tells MCHap to report posterior allele frequencies (`AFP`) and posterior allele occurance (`AOP`)\n",
    "- The output of `mchap assemble` is piped into `bgzip` and the resulting compressed VCF file is saved as `assemble.vcf.gz`\n",
    "- The `tabix` tool is then used to index the compressed VCF file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30af9ce0-fb91-44ed-a52f-437c66dcd98b",
   "metadata": {},
   "source": [
    "Look at the VCF header information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "364bad0c-ff38-491d-bd6d-86d8e723ada1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#fileformat=VCFv4.3\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#fileDate=20240910\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#source=mchap v0.9.4.dev72+g719498f.d20240710\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#phasing=None\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#commandline=\"/home/cfltxm/mambaforge/envs/mchap/bin/mchap assemble --bam input/bam/parent1.loci.bam input/bam/parent2.loci.bam input/bam/progeny001.loci.bam input/bam/progeny002.loci.bam input/bam/progeny003.loci.bam input/bam/progeny004.loci.bam input/bam/progeny005.loci.bam input/bam/progeny006.loci.bam input/bam/progeny007.loci.bam input/bam/progeny008.loci.bam input/bam/progeny009.loci.bam input/bam/progeny010.loci.bam input/bam/progeny011.loci.bam input/bam/progeny012.loci.bam input/bam/progeny013.loci.bam input/bam/progeny014.loci.bam input/bam/progeny015.loci.bam input/bam/progeny016.loci.bam input/bam/progeny017.loci.bam input/bam/progeny018.loci.bam input/bam/progeny019.loci.bam input/bam/progeny020.loci.bam --targets input/bed/targets4.bed --variants putative_snvs.vcf.gz --reference input/fasta/chr1.fa.gz --sample-pool input/pools/sample-pools.txt --ploidy 8 --report AFP AOP\"\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#randomseed=42\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#contig=<ID=chr1,length=21898217>\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FILTER=<ID=PASS,Description=\"All filters passed\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FILTER=<ID=NOA,Description=\"No observed alleles at locus\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FILTER=<ID=AF0,Description=\"All alleles have prior allele frequency of zero\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=AN,Number=1,Type=Integer,Description=\"Total number of alleles in called genotypes\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=UAN,Number=1,Type=Integer,Description=\"Total number of unique alleles in called genotypes\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=AC,Number=A,Type=Integer,Description=\"Allele count in genotypes, for each ALT allele, in the same order as listed\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=REFMASKED,Number=0,Type=Flag,Description=\"Reference allele is masked\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=NS,Number=1,Type=Integer,Description=\"Number of samples with data\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=MCI,Number=1,Type=Integer,Description=\"Number of samples with incongruent Markov chain replicates\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=DP,Number=1,Type=Integer,Description=\"Combined depth across samples\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=RCOUNT,Number=1,Type=Integer,Description=\"Total number of observed reads across all samples\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=END,Number=1,Type=Integer,Description=\"End position on CHROM\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=NVAR,Number=1,Type=Integer,Description=\"Number of input variants within assembly locus\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=SNVPOS,Number=.,Type=Integer,Description=\"Relative (1-based) positions of SNVs within haplotypes\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=AFP,Number=R,Type=Float,Description=\"Posterior mean allele frequencies\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#INFO=<ID=AOP,Number=R,Type=Float,Description=\"Posterior probability of allele occurring across all samples\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype quality\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=SQ,Number=1,Type=Integer,Description=\"Genotype support quality\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Read depth\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=RCOUNT,Number=1,Type=Integer,Description=\"Total count of read pairs within haplotype interval\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=RCALLS,Number=1,Type=Integer,Description=\"Total count of read base calls matching a known variant\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=MEC,Number=1,Type=Integer,Description=\"Minimum error correction\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=MECP,Number=1,Type=Float,Description=\"Minimum error correction proportion\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=GPM,Number=1,Type=Float,Description=\"Genotype posterior mode probability\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=SPM,Number=1,Type=Float,Description=\"Genotype support posterior mode probability\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=MCI,Number=1,Type=Integer,Description=\"Replicate Markov-chain incongruence, 0 = none, 1 = incongruence, 2 = putative CNV\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=AFP,Number=R,Type=Float,Description=\"Posterior mean allele frequencies\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[K#FORMAT=<ID=AOP,Number=R,Type=Float,Description=\"Posterior probability of allele occurring\">\n",
      "\u001b[01;31m\u001b[K#\u001b[m\u001b[KCHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tPOOL\n"
     ]
    }
   ],
   "source": [
    "zcat pooled_assembly.vcf.gz | grep \"^#\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "274d7199-45ca-44aa-b312-deca02e811b7",
   "metadata": {},
   "source": [
    "Note that the above VCF contains a single sample named \"POOL\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1a303b1-74d6-4565-b20e-7a91068414d2",
   "metadata": {},
   "source": [
    "Look at `locus001`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0bc2a7d8-e364-4571-86f4-5c39d345dd07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr1\t17591\t\u001b[01;31m\u001b[Klocus001\u001b[m\u001b[K\tGTTATTGGACAGTGACGATGGAGTGATTGCTGGCGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGTACAC\tATTATTGGACAATGACGATGGGGTGGTTGCTGGCCCAGTCCGCCAGCACCACCACCACCAAGTCAACATGTCGGACATTTATGGGGTGGTGCCACGAAACCTGATCACAAATGGCGCAC,GTTATTGGACAGTGACGATGGAGTGATTGCTGGTGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGCACAC\t.\tPASS\tAN=8;UAN=3;AC=1,1;NS=1;MCI=0;DP=619;RCOUNT=847;END=17709;NVAR=13;SNVPOS=1,12,22,26,34,35,39,65,73,96,104,115,116;AFP=0.75,0.125,0.125;AOP=1,1,1\tGT:GQ:SQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:SPM:MCI:AFP:AOP\t0/0/0/0/0/0/1/2:60:60:619:847:8019:17:0.002:1:1:0:0.75,0.125,0.125:1,1,1\n"
     ]
    }
   ],
   "source": [
    "zcat pooled_assembly.vcf.gz | grep \"locus001\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "559d087d-b500-42d3-8989-b083126c4d2c",
   "metadata": {},
   "source": [
    "### Individual genotype calling\n",
    "\n",
    "We will start by calling genotypes using a prior on allele frequencies derived from the population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "84e164f5-fa39-4881-801c-6d856ae24de9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mchap call \\\n",
    "        --bam input/bam/*.loci.bam \\\n",
    "        --haplotypes pooled_assembly.vcf.gz \\\n",
    "        --ploidy 4 \\\n",
    "        --inbreeding 0.01 \\\n",
    "        --prior-frequencies AFP \\\n",
    "        | bgzip > individual_calling.vcf.gz\n",
    "tabix -p vcf individual_calling.vcf.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c56b657-c879-4808-8b06-5b6527d2f388",
   "metadata": {},
   "source": [
    "Notes:\n",
    "- The `--prior-frequencies AFP` command tells MCHap to use the posterior allele frequencies from the assembly (`AFP`) as the priors for calling genotypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "08f9a746-cf8f-4aa2-a42e-553f99511ef5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr1\t17591\t\u001b[01;31m\u001b[Klocus001\u001b[m\u001b[K\tGTTATTGGACAGTGACGATGGAGTGATTGCTGGCGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGTACAC\tATTATTGGACAATGACGATGGGGTGGTTGCTGGCCCAGTCCGCCAGCACCACCACCACCAAGTCAACATGTCGGACATTTATGGGGTGGTGCCACGAAACCTGATCACAAATGGCGCAC,GTTATTGGACAGTGACGATGGAGTGATTGCTGGTGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGCACAC\t.\tPASS\tAN=88;UAN=3;AC=22,12;NS=22;MCI=0;DP=619;RCOUNT=847;END=17709;NVAR=13;SNVPOS=1,12,22,26,34,35,39,65,73,96,104,115,116\tGT:GQ:SQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:SPM:MCI\t0/0/1/2:60:60:30:43:383:2:0.005:1:1:0\t0/0/0/1:60:60:22:28:280:2:0.007:1:1:0\t0/0/1/2:60:60:31:41:401:1:0.002:1:1:0\t0/0/1/2:23:60:14:26:176:0:0:0.995:1:0\t0/0/1/2:60:60:25:31:320:0:0:1:1:0\t0/0/1/2:60:60:31:46:408:0:0:1:1:0\t0/0/0/1:60:60:22:33:291:0:0:1:1:0\t0/0/1/2:30:60:21:26:272:1:0.004:0.999:1:0\t0/0/1/2:60:60:38:50:497:1:0.002:1:1:0\t0/0/0/1:60:60:33:46:424:1:0.002:1:1:0\t0/0/0/1:24:60:28:41:358:4:0.011:0.996:1:0\t0/0/0/1:28:60:19:27:251:0:0:0.998:1:0\t0/0/1/2:60:60:27:35:347:0:0:1:1:0\t0/0/0/1:60:60:30:37:387:3:0.008:1:1:0\t0/0/0/1:60:60:39:58:511:0:0:1:1:0\t0/0/1/2:60:60:27:38:346:0:0:1:1:0\t0/0/1/2:60:60:25:34:326:1:0.003:1:1:0\t0/0/0/1:60:60:26:35:341:0:0:1:1:0\t0/0/1/2:60:60:34:45:441:0:0:1:1:0\t0/0/1/2:60:60:37:50:483:0:0:1:1:0\t0/0/0/1:7:7:30:39:386:3:0.008:0.786:0.786:0\t0/0/0/1:60:60:30:38:390:0:0:1:1:0\n"
     ]
    }
   ],
   "source": [
    "zcat individual_calling.vcf.gz | grep \"locus001\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c81e397a-87e2-4881-8224-d536156f9a22",
   "metadata": {},
   "source": [
    "## Pedigree informed genotype calling\n",
    "\n",
    "**WARNING: This the `call-pedigree` program is highly experimental!**\n",
    "\n",
    "Next we look at pedigree informed genotype calling.\n",
    "In MCHap, a pedigree is defined with a simple tabular format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "378412d6-e4e4-441b-b37b-3ddef7eb53e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parent1\t.\t.\n",
      "parent2\t.\t.\n",
      "progeny001\tparent1\tparent2\n",
      "progeny002\tparent1\tparent2\n",
      "progeny003\tparent1\tparent2\n",
      "progeny004\tparent1\tparent2\n",
      "progeny005\tparent1\tparent2\n",
      "progeny006\tparent1\tparent2\n",
      "progeny007\tparent1\tparent2\n",
      "progeny008\tparent1\tparent2\n",
      "progeny009\tparent1\tparent2\n",
      "progeny010\tparent1\tparent2\n",
      "progeny011\tparent1\tparent2\n",
      "progeny012\tparent1\tparent2\n",
      "progeny013\tparent1\tparent2\n",
      "progeny014\tparent1\tparent2\n",
      "progeny015\tparent1\tparent2\n",
      "progeny016\tparent1\tparent2\n",
      "progeny017\tparent1\tparent2\n",
      "progeny018\tparent1\tparent2\n",
      "progeny019\tparent1\tparent2\n",
      "progeny020\tparent1\tparent2\n"
     ]
    }
   ],
   "source": [
    "cat input/pedigree/pedigree.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eef38c1-ac25-4f77-9a35-0c7991fe7561",
   "metadata": {},
   "source": [
    "We can run `mchap call-pedigree` with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8c892051-df00-4389-a3e8-313115773b17",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/cfltxm/Repo/MCHap/mchap/application/call_pedigree.py:45: ExperimentalFeatureWarning: THIS PROGRAM IS HIGHLY EXPERIMENTAL!!!\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "mchap call-pedigree \\\n",
    "        --bam input/bam/*.loci.bam \\\n",
    "        --haplotypes pooled_assembly.vcf.gz \\\n",
    "        --ploidy 4 \\\n",
    "        --sample-parents input/pedigree/pedigree.txt \\\n",
    "        --gamete-error 0.1 \\\n",
    "        --prior-frequencies AFP \\\n",
    "        | bgzip > pedigree_calling.vcf.gz\n",
    "tabix -p vcf pedigree_calling.vcf.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "77763113-9a7e-4cca-adb0-d01c0c855d6b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[01;31m\u001b[K#CHROM\u001b[m\u001b[K\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tparent1\tparent2\tprogeny001\tprogeny002\tprogeny003\tprogeny004\tprogeny005\tprogeny006\tprogeny007\tprogeny008\tprogeny009\tprogeny010\tprogeny011\tprogeny012\tprogeny013\tprogeny014\tprogeny015\tprogeny016\tprogeny017\tprogeny018\tprogeny019\tprogeny020\n",
      "chr1\t17591\t\u001b[01;31m\u001b[Klocus001\u001b[m\u001b[K\tGTTATTGGACAGTGACGATGGAGTGATTGCTGGCGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGTACAC\tATTATTGGACAATGACGATGGGGTGGTTGCTGGCCCAGTCCGCCAGCACCACCACCACCAAGTCAACATGTCGGACATTTATGGGGTGGTGCCACGAAACCTGATCACAAATGGCGCAC,GTTATTGGACAGTGACGATGGAGTGATTGCTGGTGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGCACAC\t.\tPASS\tAN=88;UAN=3;AC=22,12;NS=22;MCI=0;DP=619;RCOUNT=847;END=17709;NVAR=13;SNVPOS=1,12,22,26,34,35,39,65,73,96,104,115,116\tGT:GQ:SQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:SPM:MCI:PEDERR\t0/0/1/2:60:60:30:43:383:2:0.005:1:1:0:0\t0/0/0/1:60:60:22:28:280:2:0.007:1:1:0:0\t0/0/1/2:60:60:31:41:401:1:0.002:1:1:0:0\t0/0/1/2:33:60:14:26:176:0:0:1:1:0:0\t0/0/1/2:60:60:25:31:320:0:0:1:1:0:0\t0/0/1/2:60:60:31:46:408:0:0:1:1:0:0\t0/0/0/1:60:60:22:33:291:0:0:1:1:0:0\t0/0/1/2:60:60:21:26:272:1:0.004:1:1:0:0\t0/0/1/2:60:60:38:50:497:1:0.002:1:1:0:0\t0/0/0/1:60:60:33:46:424:1:0.002:1:1:0:0\t0/0/0/1:24:60:28:41:358:4:0.011:0.996:1:0:0\t0/0/0/1:25:60:19:27:251:0:0:0.996:1:0:0\t0/0/1/2:60:60:27:35:347:0:0:1:1:0:0\t0/0/0/1:60:60:30:37:387:3:0.008:1:1:0:0\t0/0/0/1:60:60:39:58:511:0:0:1:1:0:0\t0/0/1/2:60:60:27:38:346:0:0:1:1:0:0\t0/0/1/2:60:60:25:34:326:1:0.003:1:1:0:0\t0/0/0/1:60:60:26:35:341:0:0:1:1:0:0\t0/0/1/2:60:60:34:45:441:0:0:1:1:0:0\t0/0/1/2:60:60:37:50:483:0:0:1:1:0:0\t0/0/0/1:5:5:30:39:386:3:0.008:0.683:0.683:0:0\t0/0/0/1:33:60:30:38:390:0:0:1:1:0:0\n"
     ]
    }
   ],
   "source": [
    "zcat pedigree_calling.vcf.gz | grep \"#CHROM\"\n",
    "zcat pedigree_calling.vcf.gz | grep \"locus001\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7cbb074-6857-4205-813f-eef9302d99d6",
   "metadata": {},
   "source": [
    "### Pedigree imputation\n",
    "\n",
    "If the pedigree names an individual that does not have an alignment file then MCHap will attempt to impute that individuals genotype.\n",
    "\n",
    "We can create a list of alingment files that does not include `parent1`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1539e88b-b5bb-4a18-a4a2-2127071a4218",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input/bam/parent2.loci.bam\n",
      "input/bam/progeny001.loci.bam\n",
      "input/bam/progeny002.loci.bam\n",
      "input/bam/progeny003.loci.bam\n",
      "input/bam/progeny004.loci.bam\n",
      "input/bam/progeny005.loci.bam\n",
      "input/bam/progeny006.loci.bam\n",
      "input/bam/progeny007.loci.bam\n",
      "input/bam/progeny008.loci.bam\n",
      "input/bam/progeny009.loci.bam\n",
      "input/bam/progeny010.loci.bam\n",
      "input/bam/progeny011.loci.bam\n",
      "input/bam/progeny012.loci.bam\n",
      "input/bam/progeny013.loci.bam\n",
      "input/bam/progeny014.loci.bam\n",
      "input/bam/progeny015.loci.bam\n",
      "input/bam/progeny016.loci.bam\n",
      "input/bam/progeny017.loci.bam\n",
      "input/bam/progeny018.loci.bam\n",
      "input/bam/progeny019.loci.bam\n",
      "input/bam/progeny020.loci.bam\n"
     ]
    }
   ],
   "source": [
    "ls input/bam/*.bam | tail -n 21 > bams.txt\n",
    "cat bams.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e18bbd4-cde3-4900-bf55-2eb0b46180c6",
   "metadata": {},
   "source": [
    "We then run `mchap call-pedigree` using that list of BAM files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c64875bc-a781-4f57-975e-0acc1b937b0d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/cfltxm/Repo/MCHap/mchap/application/call_pedigree.py:45: ExperimentalFeatureWarning: THIS PROGRAM IS HIGHLY EXPERIMENTAL!!!\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "mchap call-pedigree \\\n",
    "        --bam bams.txt \\\n",
    "        --haplotypes pooled_assembly.vcf.gz \\\n",
    "        --ploidy 4 \\\n",
    "        --sample-parents input/pedigree/pedigree.txt \\\n",
    "        --gamete-error 0.1 \\\n",
    "        --prior-frequencies AFP \\\n",
    "        | bgzip > pedigree_imputing.vcf.gz\n",
    "tabix -p vcf pedigree_imputing.vcf.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0df7d27-7070-40a2-9435-ad82fbd5f56a",
   "metadata": {},
   "source": [
    "We can see that `parent1` has been added as the last sample on the list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "722e11c6-3bd1-4b23-942c-8110e33f387d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[01;31m\u001b[K#CHROM\u001b[m\u001b[K\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tparent2\tprogeny001\tprogeny002\tprogeny003\tprogeny004\tprogeny005\tprogeny006\tprogeny007\tprogeny008\tprogeny009\tprogeny010\tprogeny011\tprogeny012\tprogeny013\tprogeny014\tprogeny015\tprogeny016\tprogeny017\tprogeny018\tprogeny019\tprogeny020\tparent1\n",
      "chr1\t17591\t\u001b[01;31m\u001b[Klocus001\u001b[m\u001b[K\tGTTATTGGACAGTGACGATGGAGTGATTGCTGGCGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGTACAC\tATTATTGGACAATGACGATGGGGTGGTTGCTGGCCCAGTCCGCCAGCACCACCACCACCAAGTCAACATGTCGGACATTTATGGGGTGGTGCCACGAAACCTGATCACAAATGGCGCAC,GTTATTGGACAGTGACGATGGAGTGATTGCTGGTGCAGGCCGCCAGCACCACCACCACCAAGTCGACATGTCCGACATTTATGGGGTGGTGCCACAAAACCTGCTCACAAATGGCACAC\t.\tPASS\tAN=88;UAN=3;AC=22,12;NS=22;MCI=0;DP=589;RCOUNT=804;END=17709;NVAR=13;SNVPOS=1,12,22,26,34,35,39,65,73,96,104,115,116\tGT:GQ:SQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:SPM:MCI:PEDERR\t0/0/0/1:16:60:22:28:280:2:0.007:0.976:1:0:0\t0/0/1/2:60:60:31:41:401:1:0.002:1:1:0:0\t0/0/1/2:28:60:14:26:176:0:0:0.998:1:0:0\t0/0/1/2:60:60:25:31:320:0:0:1:1:0:0\t0/0/1/2:60:60:31:46:408:0:0:1:1:0:0\t0/0/0/1:60:60:22:33:291:0:0:1:1:0:0\t0/0/1/2:60:60:21:26:272:1:0.004:1:1:0:0\t0/0/1/2:60:60:38:50:497:1:0.002:1:1:0:0\t0/0/0/1:60:60:33:46:424:1:0.002:1:1:0:0\t0/0/0/1:25:60:28:41:358:4:0.011:0.997:1:0:0\t0/0/0/1:28:60:19:27:251:0:0:0.998:1:0:0\t0/0/1/2:60:60:27:35:347:0:0:1:1:0:0\t0/0/0/1:60:60:30:37:387:3:0.008:1:1:0:0\t0/0/0/1:60:60:39:58:511:0:0:1:1:0:0\t0/0/1/2:60:60:27:38:346:0:0:1:1:0:0\t0/0/1/2:33:60:25:34:326:1:0.003:1:1:0:0\t0/0/0/1:60:60:26:35:341:0:0:1:1:0:0\t0/0/1/2:60:60:34:45:441:0:0:1:1:0:0\t0/0/1/2:60:60:37:50:483:0:0:1:1:0:0\t0/0/0/1:5:5:30:39:386:3:0.008:0.664:0.664:0:0\t0/0/0/1:60:60:30:38:390:0:0:1:1:0:0\t0/0/1/2:3:3:0:0:0:0:.:0.528:0.537:0:0\n"
     ]
    }
   ],
   "source": [
    "zcat pedigree_imputing.vcf.gz | grep \"#CHROM\"\n",
    "zcat pedigree_imputing.vcf.gz | grep \"locus001\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3864bb6-408d-452a-952e-44d6fc486305",
   "metadata": {},
   "source": [
    "Note that the genotype for `parent1` is the same as above, but the quality scores are much lower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4e3fcfc0-d33a-4acc-8b82-d0132782f1fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "individual_calling.vcf.gz\n",
      "parent1\n",
      "0/0/1/2:60:60:30:43:383:2:0.005:1:1:0\n",
      "pedigree_calling.vcf.gz\n",
      "parent1\n",
      "0/0/1/2:60:60:30:43:383:2:0.005:1:1:0:0\n",
      "pedigree_imputing.vcf.gz\n",
      "parent1\n",
      "0/0/1/2:3:3:0:0:0:0:.:0.528:0.537:0:0\n"
     ]
    }
   ],
   "source": [
    "echo individual_calling.vcf.gz\n",
    "zcat individual_calling.vcf.gz | grep \"#CHROM\" | cut -f10\n",
    "zcat individual_calling.vcf.gz | grep \"locus001\" | cut -f10\n",
    "echo pedigree_calling.vcf.gz\n",
    "zcat pedigree_calling.vcf.gz | grep \"#CHROM\" | cut -f10\n",
    "zcat pedigree_calling.vcf.gz | grep \"locus001\" | cut -f10\n",
    "echo pedigree_imputing.vcf.gz\n",
    "zcat pedigree_imputing.vcf.gz | grep \"#CHROM\" | cut -f31\n",
    "zcat pedigree_imputing.vcf.gz | grep \"locus001\" | cut -f31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22c1a82c-a190-4786-85b2-a5b556a555d1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}