{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using `bash` in bioinformatics\n",
    "\n",
    "The Unix commands shown in this notebook are quite intimidating at first. There's a lot of material becaue I try to cover all the stuff you willl need to construct a typicla bioinformatics pipeline. \n",
    "\n",
    "Try to understand the expressions and play aound to see what varying them will do. You do not need to memorize the commands now - you can always come back to this notebook for reference. You will also get more practice when you do the Bioiinformatics sessions to construct Unix pipelines for processing your RNA-seq data sets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hack to handle broken pipes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleanup () { \n",
    "    :\n",
    "}\n",
    "\n",
    "trap \"cleanup\" SIGPIPE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Treat unset variables as an error, and immediately exit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "set -u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.1. Files, directories and working with text\n",
    "\n",
    "Quickly reveiw what the following navigation and file/directory manipulation commands do:\n",
    "\n",
    "```bash\n",
    ">, >>, <, |, ls, cd, pwd, cp, mv, rm, mkdir, rmdir\n",
    "```\n",
    "\n",
    "This session dives quite deep into Unix shell scripting and invovles quite a bit of text processing. We review some basic commands for working with text in Unix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This is the file we will warm up with\n",
    "\n",
    "We create it using a HEREDOC that allows us to construct multi-line files and redirrect it to the file `junk.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat > junk.txt <<EOF \n",
    "1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n",
    "2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
    "3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n",
    "4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n",
    "EOF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We can see the contents of a file with `cat`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n",
      "2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n",
      "4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n"
     ]
    }
   ],
   "source": [
    "cat junk.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using `head` and `tail`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n",
      "2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n"
     ]
    }
   ],
   "source": [
    "head -3 junk.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n",
      "4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n"
     ]
    }
   ],
   "source": [
    "tail -3 junk.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n",
      "4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n"
     ]
    }
   ],
   "source": [
    "tail +3 junk.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pipling to chain commnads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n"
     ]
    }
   ],
   "source": [
    "head -3 junk.txt | tail -2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Searching within a file with `grep`\n",
    "\n",
    "- `grep` takes a REGEX (regular expression) as the query and a file (or files) to search\n",
    "- Normally grep shows lines that match the REGEX argument\n",
    "- The `-v` flag shows the lines that DO NOT match the REGEX argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2. GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "4. TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n"
     ]
    }
   ],
   "source": [
    "grep \"[2|4]\" junk.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n",
      "3. GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n"
     ]
    }
   ],
   "source": [
    "grep -v \"[2|4]\" junk.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Use a regular expression with `grep`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n"
     ]
    }
   ],
   "source": [
    "grep \"CCC.*CCC\" junk.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get rid of numbers with `cut`\n",
    "\n",
    "`cut` extracts columns separated by delimiters\n",
    "- here the delimiter is a space `-d ' '`\n",
    "- and we cut the second column `-f 2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n",
      "GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n",
      "TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n"
     ]
    }
   ],
   "source": [
    "cut -d ' ' -f 2 junk.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get rid of numbers with `grep`\n",
    "\n",
    "The `-o` flag to `grep` shows only the strings mathcing the regular expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n",
      "GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n",
      "TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n"
     ]
    }
   ],
   "source": [
    "grep -o '[A-Z]*' junk.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Convert DNA to RNA with `tr`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. CCCUAACCCCCUAACCCCCUAACCCUCAGUCGGGGAGGCGACAAUAGCUG\n",
      "2. GUCAUAUGUUCUGUACGUUAUUGGCCAACUGAUCAUACCUGAAUCGAGCC\n",
      "3. GAACCGGGAUUAUCAAAGACGAACAUGGUCGGGUCCUUGAACCAAACGAA\n",
      "4. UCUCCGUCCGCUGGCGUGUUUUUCUUUUCUCAAGUGGGCAAGUUACCCGG\n"
     ]
    }
   ],
   "source": [
    "cat junk.txt | tr T U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Reverse complement with `tr`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAGCTATTGTCGCCTCCCCGACTGAGGGTTAGGGGGTTAGGGGGTTAGGG\n",
      "GGCTCGATTCAGGTATGATCAGTTGGCCAATAACGTACAGAACATATGAC\n",
      "TTCGTTTGGTTCAAGGACCCGACCATGTTCGTCTTTGATAATCCCGGTTC\n",
      "CCGGGTAACTTGCCCACTTGAGAAAAGAAAAACACGCCAGCGGACGGAGA\n"
     ]
    }
   ],
   "source": [
    "grep -o '[ACTG]*' junk.txt |\n",
    "tr ACGT TGCA |\n",
    "rev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Split command over multiple lines\n",
    "\n",
    "Use `\\` to continue command to next line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAGCTATTGTCGCCTCCCCGACTGAGGGTTAGGGGGTTAGGGGGTTAGGG\n",
      "GGCTCGATTCAGGTATGATCAGTTGGCCAATAACGTACAGAACATATGAC\n",
      "TTCGTTTGGTTCAAGGACCCGACCATGTTCGTCTTTGATAATCCCGGTTC\n",
      "CCGGGTAACTTGCCCACTTGAGAAAAGAAAAACACGCCAGCGGACGGAGA\n"
     ]
    }
   ],
   "source": [
    "grep -o '[ACTG]*' junk.txt | tr \\\n",
    "ACGT TGCA | rev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually the line continuation character is not necesssary if you break at a pipe `|`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAGCTATTGTCGCCTCCCCGACTGAGGGTTAGGGGGTTAGGGGGTTAGGG\n",
      "GGCTCGATTCAGGTATGATCAGTTGGCCAATAACGTACAGAACATATGAC\n",
      "TTCGTTTGGTTCAAGGACCCGACCATGTTCGTCTTTGATAATCCCGGTTC\n",
      "CCGGGTAACTTGCCCACTTGAGAAAAGAAAAACACGCCAGCGGACGGAGA\n"
     ]
    }
   ],
   "source": [
    "grep -o '[ACTG]*' junk.txt | \n",
    "tr ACGT TGCA | \n",
    "rev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.2 Getting files\n",
    "\n",
    "We will download Chromosomes 7 and 8 of Crytococcus neoformans in FASTA format from NCBI Nucleotide. We can find their IDs through a [web query](https://www.ncbi.nlm.nih.gov/nuccore/?term=cryptococcus+neoformans+B-3501A).\n",
    "\n",
    "We find that the IDs are\n",
    "\n",
    "```\n",
    "CM000046.1 for chormosome 7\n",
    "CM000047.1 for chormosome 8\n",
    "```\n",
    "\n",
    "We will use the NCBI Entrez tool `efetch` to download the FASTA files to a folder `seqs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Trying to establish local installations of any missing Perl modules\n",
      "(as logged in /home/jovyan/edirect/setup-deps.log).\n",
      "Please be patient, as this step may take a little while.\n",
      "\n",
      "Entrez Direct has been successfully downloaded and installed.\n",
      "\n",
      "In order to complete the configuration process, please execute the following:\n",
      "\n",
      "  echo \"export PATH=\\${PATH}:/home/jovyan/edirect\" >> $HOME/.bashrc\n",
      "\n",
      "or manually edit the PATH variable assignment in your .bash_profile file.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "### setup efetch command\n",
    "cd ~\n",
    "#/bin/bash\n",
    "perl -MNet::FTP -e \\\n",
    "    '$ftp = new Net::FTP(\"ftp.ncbi.nlm.nih.gov\", Passive => 1);\n",
    "     $ftp->login; $ftp->binary;\n",
    "     $ftp->get(\"/entrez/entrezdirect/edirect.tar.gz\");'\n",
    "gunzip -c edirect.tar.gz | tar xf -\n",
    "rm edirect.tar.gz\n",
    "#builtin exit\n",
    "export PATH=${PATH}:$HOME/edirect >& /dev/null || setenv PATH \"${PATH}:$HOME/edirect\"\n",
    "./edirect/setup.sh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p seqs\n",
    "efetch -db=nuccore -format=fasta -id=CM000046.1 > seqs/CM000046.1.fa\n",
    "efetch -db=nuccore -format=fasta -id=CM000047.1 > seqs/CM000047.1.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa      \u001b[0m\u001b[01;31mcneo040623.scaffs.fa.gz\u001b[0m\n",
      "CM000047.1.fa      \u001b[01;31mCryptococcus_neoformans.ASM9104v1.43.gtf.gz\u001b[0m\n",
      "CM000047.1.fa.bak\n"
     ]
    }
   ],
   "source": [
    "ls seqs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">CM000046.1 Cryptococcus neoformans var. neoformans B-3501A chromosome 7, whole genome shotgun sequence\n",
      "TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAAGTTGTTTATTTAAGTAGTTT\n",
      "CTTTAATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAAATGCAACCTCTATAACCCC\n",
      "CTAACCCTTTTATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATA\n",
      "ACCCCCGAAAAGGATTGGTGTCGAGTTAGTAAAGGCGAGGGCTAGGACGCTGGTCTTGTGAGGAAGAGGA\n",
      "GGTTGGATTATTTTTTCTTTTCTTTAATAAGTTGTTTATTTAAGTAGTTTCTTTTATTCGGGTAACACAC\n",
      "ACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATAACCCCCATCCAAGGTGCCACGGTTGA\n",
      "GATTGAGCCCCACACACTGTCGGGACAAAGGAGAAACGACCTTCGGGTCAGAGGTTCCAGCGCTCTGGCC\n",
      "TTCACTGACTACGACCTGAAGGTTTACTCCCTCGGGGACCGAGACGCGAGGAGCACAGCCACCCCCAGCA\n",
      "CCCCCAACAGCAAACTGGCCGACTTCTGCTTGGACCGGTGCGTGAACTGGCTCGACAAGGTGGGTCAGGT\n"
     ]
    }
   ],
   "source": [
    "head seqs/CM000046.1.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">CM000047.1 Cryptococcus neoformans var. neoformans B-3501A chromosome 8, whole genome shotgun sequence\n",
      "ATGGCCTTCCTGACTACGACCTGAGGTTTACTCCCTCGGCGACCGAGACGCGAGGAGCACAGCCACCCCA\n",
      "AGCACCCCCAACAGCAAGCTGGCCGAATTCTGCTTGGACCGGTGCGTGAACTGGCTCGACAAGGTGGGTC\n",
      "AGGTCGTCTCGAAGAACGCTCCGAAAGTCACTGGTGGGGTCTTTAAACCGATCATCCTTTCCACTGGTGG\n",
      "CCTGATGAGCAGGAGCACAGCAGACGAATGGAAGGAGTGGAGGGAGGCGATGCCGGTGGGGGGGTTCGAG\n",
      "AAAATGGAGAAACGGATTGGTGTCGAGTTAGTAAAGGCGAGGGCTAGGACGCTGGTCTTGTGAGGAAGAG\n",
      "GAGGTTGGATTATTTTTTCTTTTCTTTAAAAAGTTGTTTATTTAAGTAGTTTCTTTAATTCGGGCAACCC\n",
      "ACACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATAACCCTCAATAGCTCTGGCTGGTAC\n",
      "AGATGGGCAGAAACTCACCTTGCTCCCACGGCATACATCCCTCCCAATAAGTATCTGCAGATGGCGAATA\n",
      "AACCCCATGCTCTGCGCAAATATATTGAAGGCTGACGACTTTGATGTTATGATCCACCCCATAAACGCTC\n"
     ]
    }
   ],
   "source": [
    "head seqs/CM000047.1.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.3 Compressing ard archiving files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `zip`, `gzip`, `bzip2`, `bgzip`\n",
    "\n",
    "In Bioinformatics the most common formats are `gzip` and `bgzip`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa      \u001b[0m\u001b[01;31mcneo040623.scaffs.fa.gz\u001b[0m\n",
      "CM000047.1.fa      \u001b[01;31mCryptococcus_neoformans.ASM9104v1.43.gtf.gz\u001b[0m\n",
      "CM000047.1.fa.bak\n"
     ]
    }
   ],
   "source": [
    "cd ~/seqs\n",
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compressing with `gzip`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;31mCM000046.1.fa.gz\u001b[0m   \u001b[01;31mcneo040623.scaffs.fa.gz\u001b[0m\n",
      "CM000047.1.fa.bak  \u001b[01;31mCryptococcus_neoformans.ASM9104v1.43.gtf.gz\u001b[0m\n",
      "\u001b[01;31mCM000047.1.fa.gz\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "gzip *fa\n",
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Uncompressing with `gzip`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa  cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "gunzip *gz\n",
    "ls "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compressing with `bgzip`\n",
    "\n",
    "`bgzip` (Blocked GNU Zip Format) is commonly used in Bioinformatics, and allows random access of the content. BAM, BCF and VCF file formats are typically `bgzip` compressed.\n",
    "\n",
    "Unfortunately `bgzip` does not take multiple files, and we have to use the `xargs` command to pipe in one file at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  -n, --max-args=MAX-ARGS      use at most MAX-ARGS arguments per command line\n"
     ]
    }
   ],
   "source": [
    "xargs --help | grep '\\-n,'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;31mCM000046.1.fa.gz\u001b[0m   \u001b[01;31mcneo040623.scaffs.fa.gz\u001b[0m\n",
      "CM000047.1.fa.bak  Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "\u001b[01;31mCM000047.1.fa.gz\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "ls *fa | xargs -n1 bgzip\n",
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note that `bgzip` files can be uncompressed with `gzip`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa  cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "gunzip *gz\n",
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Viewing compressed files\n",
    "\n",
    "Some utilities can work directly with gzipped files. One that we will be using extensively is `zcat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;31mCM000046.1.fa.gz\u001b[0m   \u001b[01;31mcneo040623.scaffs.fa.gz\u001b[0m\n",
      "CM000047.1.fa.bak  Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "\u001b[01;31mCM000047.1.fa.gz\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "ls *fa | xargs -n1 bgzip\n",
    "ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">CM000046.1 Cryptococcus neoformans var. neoformans B-3501A chromosome 7, whole genome shotgun sequence\n",
      "TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAAGTTGTTTATTTAAGTAGTTT\n",
      "CTTTAATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAAATGCAACCTCTATAACCCC\n",
      "CTAACCCTTTTATTCGGGTAACACACACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATA\n",
      "ACCCCCGAAAAGGATTGGTGTCGAGTTAGTAAAGGCGAGGGCTAGGACGCTGGTCTTGTGAGGAAGAGGA\n",
      "GGTTGGATTATTTTTTCTTTTCTTTAATAAGTTGTTTATTTAAGTAGTTTCTTTTATTCGGGTAACACAC\n",
      "ACGACAACCCAATAAATTAAACAACGAAAAATGCAACCTCTATAACCCCCATCCAAGGTGCCACGGTTGA\n",
      "GATTGAGCCCCACACACTGTCGGGACAAAGGAGAAACGACCTTCGGGTCAGAGGTTCCAGCGCTCTGGCC\n",
      "TTCACTGACTACGACCTGAAGGTTTACTCCCTCGGGGACCGAGACGCGAGGAGCACAGCCACCCCCAGCA\n",
      "CCCCCAACAGCAAACTGGCCGACTTCTGCTTGGACCGGTGCGTGAACTGGCTCGACAAGGTGGGTCAGGT\n"
     ]
    }
   ],
   "source": [
    "zcat CM000046.1.fa.gz | head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa  cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "gunzip *gz\n",
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Archival with `tar`\n",
    "\n",
    "Often we want to compress an entire direcotry (possibly with subdirectories). In Unix, this is most commonly done with `tar`. Flags to note\n",
    "\n",
    "- `c` compress\n",
    "- `x` extract \n",
    "- `z` gzip or gunzip\n",
    "- `f` files\n",
    "- `v` verbose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Usage: tar [OPTION...] [FILE]...\n",
      "GNU 'tar' saves many files together into a single tape or disk archive, and can\n",
      "restore individual files from the archive.\n",
      "\n",
      "Examples:\n",
      "  tar -cf archive.tar foo bar  # Create archive.tar from files foo and bar.\n",
      "  tar -tvf archive.tar         # List all files in archive.tar verbosely.\n",
      "  tar -xf archive.tar          # Extract all files from archive.tar.\n"
     ]
    }
   ],
   "source": [
    "tar --help | head -8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa  cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "cd \n",
    "ls seqs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seqs/\n",
      "seqs/CM000047.1.fa.bak\n",
      "seqs/CM000046.1.fa\n",
      "seqs/Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "seqs/CM000047.1.fa\n",
      "seqs/cneo040623.scaffs.fa\n",
      "\u001b[0m\u001b[01;31mseqs.tar.gz\u001b[0m\n",
      "\n",
      "seqs:\n",
      "CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa  cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "tar czvf seqs.tar.gz seqs\n",
    "ls seqs*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### View contents of archive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "drwxr-xr-x jovyan/users      0 2019-06-26 09:39 seqs/\n",
      "-rw-r--r-- jovyan/users 1194398 2019-06-26 09:30 seqs/CM000047.1.fa.bak\n",
      "-rw-r--r-- jovyan/users 1434513 2019-06-26 09:38 seqs/CM000046.1.fa\n",
      "-rw-r--r-- jovyan/users 27773163 2019-06-26 09:30 seqs/Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "-rw-r--r-- jovyan/users  1194378 2019-06-26 09:38 seqs/CM000047.1.fa\n",
      "-rw-r--r-- jovyan/users 18907235 2019-06-26 09:39 seqs/cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "tar tzvf seqs.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remove direcotry and regenerate from archive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;31mseqs.tar.gz\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "rm -rf seqs\n",
    "ls seqs*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seqs/\n",
      "seqs/CM000047.1.fa.bak\n",
      "seqs/CM000046.1.fa\n",
      "seqs/Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "seqs/CM000047.1.fa\n",
      "seqs/cneo040623.scaffs.fa\n",
      "\u001b[0m\u001b[01;31mseqs.tar.gz\u001b[0m\n",
      "\n",
      "seqs:\n",
      "CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa  cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "tar xzvf seqs.tar.gz\n",
    "ls seqs*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm seqs.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0.4 Checksums"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When working with genomic data, we deal with very large files. There is a small risk that these files will be corrupted over time or during data transfer. To ensure that files are not changed, we use a “checksum” function. This is a function that generates an long, essentially random number called a checksum that represents the contents of the file. When the file contents change, so will the checksum. In theory, there is a very small probability that two different files generate the same checksum, but in practice the probability is too small to worry about.\n",
    "\n",
    "There are several different algorithms for generating the checksums, and at least 3 Unix commands to do so, but they all work very similarly for our purposes.\n",
    "\n",
    "The strategy is:\n",
    "\n",
    "- Generate and store a checksum together with a data file whose integrity you care about\n",
    "- When you use or download the data, re-generate the checksum (using the same algorithm e.g. MD5) and compare with the checksum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa  CM000047.1.fa.bak     Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa  cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "ls ~/seqs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "md5sum ~/seqs/*fa > MD5_CHECKSUM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We check that the files are unchanged wth the `-c` flag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jovyan/seqs/CM000046.1.fa: OK\n",
      "/home/jovyan/seqs/CM000047.1.fa: OK\n",
      "/home/jovyan/seqs/cneo040623.scaffs.fa: OK\n"
     ]
    }
   ],
   "source": [
    "md5sum -c MD5_CHECKSUM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Suppsoe one of the files is corrupted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "echo 'Oops I did it again' >> seqs/CM000047.1.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jovyan/seqs/CM000046.1.fa: OK\n",
      "/home/jovyan/seqs/CM000047.1.fa: FAILED\n",
      "/home/jovyan/seqs/cneo040623.scaffs.fa: OK\n",
      "md5sum: WARNING: 1 computed checksum did NOT match\n"
     ]
    }
   ],
   "source": [
    "md5sum -c MD5_CHECKSUM  \\\n",
    "| cat # this is jusst so we can Run All"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using `diff` to find changes\n",
    "\n",
    "Let's rename the altered file, re-download the orginal and see what was changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "mv seqs/CM000047.1.fa seqs/CM000047.1.fa.bak"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "efetch -db=nuccore -format=fasta -id=CM000047.1 > seqs/CM000047.1.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- seqs/CM000047.1.fa\t2019-06-26 09:39:10.615210391 +0000\n",
      "+++ seqs/CM000047.1.fa.bak\t2019-06-26 09:39:09.187204139 +0000\n",
      "@@ -16820,3 +16820,4 @@\n",
      " CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG\n",
      " ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT\n",
      " AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG\n",
      "+Oops I did it again\n"
     ]
    }
   ],
   "source": [
    "diff -u seqs/CM000047.1.fa seqs/CM000047.1.fa.bak \\\n",
    "| cat # this is jusst so we can Run All"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explanaiton from [official documentation]\n",
    "\n",
    "Header\n",
    "\n",
    "```\n",
    "--- seqs/CM000047.1.fa\t2018-06-21 09:58:18.246185758 -0400\n",
    "+++ seqs/CM000047.1.fa.bak\t2018-06-21 09:55:55.967707764 -0400\n",
    "```\n",
    "\n",
    "Explanation\n",
    "\n",
    "```\n",
    "The unified output format starts with a two-line header, which looks like this:\n",
    "\n",
    "--- from-file from-file-modification-time\n",
    "+++ to-file to-file-modification-time\n",
    "The timestamp looks like ‘2002-02-21 23:30:39.942229878 -0800’ to indicate the date, time with fractional seconds, and time zone. The fractional seconds are omitted on hosts that do not support fractional timestamps.\n",
    "\n",
    "You can change the header’s content with the --label=label option. See Alternate Names.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hunk \n",
    "\n",
    "```\n",
    "@@ -16820,3 +16820,4 @@\n",
    " CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG\n",
    " ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT\n",
    " AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG\n",
    "+Oops I did it again\n",
    "```\n",
    "\n",
    "Explanation\n",
    "\n",
    "```\n",
    "@@ from-file-line-numbers to-file-line-numbers @@\n",
    " line-from-either-file\n",
    " line-from-either-file…\n",
    " \n",
    " If a hunk contains just one line, only its start line number appears. Otherwise its line numbers look like ‘start,count’. An empty hunk is considered to start at the line that follows the hunk.\n",
    "\n",
    "If a hunk and its context contain two or more lines, its line numbers look like ‘start,count’. Otherwise only its end line number appears. An empty hunk is considered to end at the line that precedes the hunk.\n",
    "\n",
    "The lines common to both files begin with a space character. The lines that actually differ between the two files have one of the following indicator characters in the left print column:\n",
    "\n",
    "‘+’\n",
    "A line was added here to the first file.\n",
    "\n",
    "‘-’\n",
    "A line was removed here from the first file.\n",
    " ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Verify\n",
    "\n",
    "We can confirm that a line was added at the position given by `diff`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG\n",
      "ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT\n",
      "AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG\n"
     ]
    }
   ],
   "source": [
    "cat seqs/CM000047.1.fa | tail +16820 | head -3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CCGTGTCGGCACGCATCTGATTTTTTTCTTATTATTCATAACAACGGCTACAGATCACGCAAGCGCATCG\n",
      "ACTCCGCTGCTGCGAAACACCACGAGTACTGCAACCATTGCTGCGGATGAGCTTAGGGGGTTAGGGGGTT\n",
      "AGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGGGTTAGGGG\n",
      "Oops I did it again\n"
     ]
    }
   ],
   "source": [
    "cat seqs/CM000047.1.fa.bak | tail +16820 | head -4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Working with FASTA files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Strucrure of a FASTA file\n",
    "\n",
    "From [NCBI FASTA specification](https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=BlastHelp)\n",
    "\n",
    "```\n",
    "A sequence in FASTA format begins with a single-line description, followed by lines of sequence data. The description line (defline) is distinguished from the sequence data by a greater-than (\">\") symbol at the beginning. It is recommended that all lines of text be shorter than 80 characters in length. An example sequence in FASTA format is:\n",
    "\n",
    "\t\t>P01013 GENE X PROTEIN (OVALBUMIN-RELATED)\n",
    "\t\tQIKDLLVSSSTDLDTTLVLVNAIYFKGMWKTAFNAEDTREMPFHVTKQESKPVQMMCMNNSFNVATLPAE\n",
    "\t\tKMKILELPFASGDLSMLVLLPDEVSDLERIEKTINFEKLTEWTNPNTMEKRRVKVYLPQMKIEEKYNLTS\n",
    "\t\tVLMALGMTDLFIPSANLTGISSAESLKISQAVHGAFMELSEDGIEMAGSTGVIEDIKHSPESEQFRADHP\n",
    "\t\tFLFLIKHNPTNTIVYFGRYWSP\n",
    "```\n",
    "\n",
    "Note that a FASTA file may contain mulitple sequences, each with its own description line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2019-06-26 09:39:11--  ftp://ftp.ensemblgenomes.org/pub/release-43/fungi/gtf/cryptococcus_neoformans/Cryptococcus_neoformans.ASM9104v1.43.gtf.gz\n",
      "           => ‘Cryptococcus_neoformans.ASM9104v1.43.gtf.gz’\n",
      "Resolving ftp.ensemblgenomes.org (ftp.ensemblgenomes.org)... 193.62.197.94\n",
      "Connecting to ftp.ensemblgenomes.org (ftp.ensemblgenomes.org)|193.62.197.94|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /pub/release-43/fungi/gtf/cryptococcus_neoformans ... done.\n",
      "==> SIZE Cryptococcus_neoformans.ASM9104v1.43.gtf.gz ... 1546239\n",
      "==> PASV ... done.    ==> RETR Cryptococcus_neoformans.ASM9104v1.43.gtf.gz ... done.\n",
      "Length: 1546239 (1.5M) (unauthoritative)\n",
      "\n",
      "Cryptococcus_neofor 100%[===================>]   1.47M  2.04MB/s    in 0.7s    \n",
      "\n",
      "2019-06-26 09:39:13 (2.04 MB/s) - ‘Cryptococcus_neoformans.ASM9104v1.43.gtf.gz’ saved [1546239]\n",
      "\n",
      "--2019-06-26 09:39:13--  http://www-sequence.stanford.edu/group/C.neoformans/files/cneo040623.scaffs.fa.gz\n",
      "Resolving www-sequence.stanford.edu (www-sequence.stanford.edu)... 171.66.102.5\n",
      "Connecting to www-sequence.stanford.edu (www-sequence.stanford.edu)|171.66.102.5|:80... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 5856251 (5.6M) [text/plain]\n",
      "Saving to: ‘cneo040623.scaffs.fa.gz’\n",
      "\n",
      "cneo040623.scaffs.f 100%[===================>]   5.58M  1.65MB/s    in 3.4s    \n",
      "\n",
      "2019-06-26 09:39:16 (1.65 MB/s) - ‘cneo040623.scaffs.fa.gz’ saved [5856251/5856251]\n",
      "\n",
      "/home/jovyan\n"
     ]
    }
   ],
   "source": [
    "cd seqs\n",
    "wget ftp://ftp.ensemblgenomes.org/pub/release-43/fungi/gtf/cryptococcus_neoformans/Cryptococcus_neoformans.ASM9104v1.43.gtf.gz\n",
    "wget http://www-sequence.stanford.edu/group/C.neoformans/files/cneo040623.scaffs.fa.gz\n",
    "cd -"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set up variables for files we will be using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DATA_DIR='/data/hts/2018/Granek_4837_180427A5'\n",
    "DATA_DIR='./seqs'\n",
    "FASTA='cneo040623.scaffs.fa.gz'\n",
    "FASTQ='10_MA_C_S3_L003_R1_001.fastq.gz'\n",
    "GTF='Cryptococcus_neoformans.ASM9104v1.43.gtf.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This is how you get the value of a shell variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./seqs\n"
     ]
    }
   ],
   "source": [
    "echo $DATA_DIR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### List files in data direcotry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa\n",
      "CM000047.1.fa\n",
      "CM000047.1.fa.bak\n"
     ]
    }
   ],
   "source": [
    "ls $DATA_DIR | head -3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inspect the contents of a zipped file.\n",
    "\n",
    "`zcat` is like `cat` but works on gzipped files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">chr01.b3501.040506 5 contigs (2263073 bp)  \n",
      "CCCTAACCCCCTAACCCCCTAACCCTCAGTCGGGGAGGCGACAATAGCTG\n",
      "GTCATATGTTCTGTACGTTATTGGCCAACTGATCATACCTGAATCGAGCC\n",
      "GAACCGGGATTATCAAAGACGAACATGGTCGGGTCCTTGAACCAAACGAA\n",
      "TCTCCGTCCGCTGGCGTGTTTTTCTTTTCTCAAGTGGGCAAGTTACCCGG\n",
      "GGATTTATACAACGATAAGAGGCTATCCGCTGACTATAAATTGTGTTAGC\n",
      "TGATCCAGATCCGCAGTTGGCACAAACCGATGTTTCCTTTTCCTCAGGTC\n",
      "CTGAATATTTCTCTTCCATCAAGGCTCCCAACCCCGAGGGGAGTATTTCC\n",
      "ACCCGATCTGATTCCAAGAGGTCCAGCGTCAACCAGGTACGTGAATGTTT\n",
      "TCTTCTTTAGTGGTGGGGACCCTATGCTTGTCCCACAGAACTGCACTGAA\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTA | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Show headers only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">chr01.b3501.040506 5 contigs (2263073 bp)  \n",
      ">chr02.b3501.040506 6 contigs (1607892 bp)  \n",
      ">chr03.b3501.040506 6 contigs (2025976 bp)  \n",
      ">chr04.b3501.040623 5 contigs (1762351 bp)  \n",
      ">chr05.b3501.040506 7 contigs (1432964 bp)  \n",
      ">chr06.b3501.040616 4 contigs (1411173 bp)  \n",
      ">chr07.b3501.040506 6 contigs (1356306 bp)  \n",
      ">chr08.b3501.040506 3 contigs (1153453 bp)  \n",
      ">chr09.b3501.040506 5 contigs (1072068 bp)  \n",
      ">chr10.b3501.040506 3 contigs (1017008 bp)  \n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTA   | grep '^>' | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Show headers with line numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1:>chr01.b3501.040506 5 contigs (2263073 bp)  \n",
      "45264:>chr02.b3501.040506 6 contigs (1607892 bp)  \n",
      "77423:>chr03.b3501.040506 6 contigs (2025976 bp)  \n",
      "117944:>chr04.b3501.040623 5 contigs (1762351 bp)  \n",
      "153193:>chr05.b3501.040506 7 contigs (1432964 bp)  \n",
      "181854:>chr06.b3501.040616 4 contigs (1411173 bp)  \n",
      "210079:>chr07.b3501.040506 6 contigs (1356306 bp)  \n",
      "237207:>chr08.b3501.040506 3 contigs (1153453 bp)  \n",
      "260278:>chr09.b3501.040506 5 contigs (1072068 bp)  \n",
      "281721:>chr10.b3501.040506 3 contigs (1017008 bp)  \n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTA| grep -n '^>' | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Count number of headers (one for each crhomosome)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTA  | grep -c '^>'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract first part of the sequence from chromosome 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">chr07.b3501.040506 6 contigs (1356306 bp)  \n",
      "TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAA\n",
      "GTTGTTTATTTAAGTAGTTTCTTTAATTCGGGTAACACACACGACAACCC\n",
      "AATAAATTAAACAACGAAAAAATGCAACCTCTATAACCCCCTAACCCTTT\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTA | sed -n '210079,210082p'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract first line of chromosome 7 and save to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcat $DATA_DIR/$FASTA | sed -n '210080p' > dna.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inspect the file we just saved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TGGTCTTATGAGGAAGAGGAGTTTGGATTATTTTTTCTTTTCTTTAAAAA\n"
     ]
    }
   ],
   "source": [
    "cat dna.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find complement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ACCAGAATACTCCTTCTCCTCAAACCTAATAAAAAAGAAAAGAAATTTTT\n"
     ]
    }
   ],
   "source": [
    "cat dna.txt | tr ACGT TGCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find reverse complement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TTTTTAAAGAAAAGAAAAAATAATCCAAACTCCTCTTCCTCATAAGACCA\n"
     ]
    }
   ],
   "source": [
    "cat dna.txt | tr ACGT TGCA | rev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Delete A and T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GGCGGGGGGGGGCC\n"
     ]
    }
   ],
   "source": [
    "cat dna.txt | tr -d 'AT' "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Only count non A and T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\n"
     ]
    }
   ],
   "source": [
    "cat dna.txt | tr -d 'AT' | tr -d '\\n' | wc -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bash: man: command not found\n"
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "man tr | grep '\\-d,' -A2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Count C and G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\n"
     ]
    }
   ],
   "source": [
    "grep '[CG]' -o dna.txt | wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bash: man: command not found\n"
     ]
    },
    {
     "ename": "",
     "evalue": "1",
     "output_type": "error",
     "traceback": []
    }
   ],
   "source": [
    "man grep | grep '\\-o,' -A2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Specify data directory as a shell variable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Working with FASTQ files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure of a FASTQ file\n",
    "\n",
    "From [Wikipedia](https://en.wikipedia.org/wiki/FASTQ_format)\n",
    "\n",
    "```\n",
    "A FASTQ file normally uses four lines per sequence.\n",
    "\n",
    "Line 1 begins with a '@' character and is followed by a sequence identifier and an optional description (like a FASTA title line).\n",
    "Line 2 is the raw sequence letters.\n",
    "Line 3 begins with a '+' character and is optionally followed by the same sequence identifier (and any description) again.\n",
    "Line 4 encodes the quality values for the sequence in Line 2, and must contain the same number of symbols as letters in the sequence.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Specify target file as a shell variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR='/data/hts2018_pilot/Granek_4837_180427A5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;31m10_MA_C_S3_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m26_RZ_C_S16_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m10_MA_C_S3_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m27_MA_P_S38_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m10_MA_C_S3_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m27_MA_P_S38_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m10_MA_C_S3_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m27_MA_P_S38_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m10_RZ_C_S11_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m27_MA_P_S38_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m10_RZ_C_S11_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m27_RZ_P_S46_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m10_RZ_C_S11_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m27_RZ_P_S46_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m10_RZ_C_S11_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m27_RZ_P_S46_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_MA_J_S20_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m27_RZ_P_S46_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_MA_J_S20_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_MA_C_S1_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_MA_J_S20_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_MA_C_S1_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_MA_J_S20_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_MA_C_S1_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_RZ_J_S28_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_MA_C_S1_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_RZ_J_S28_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_RZ_C_S9_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_RZ_J_S28_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_RZ_C_S9_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m11_RZ_J_S28_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_RZ_C_S9_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_MA_P_S36_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_RZ_C_S9_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_MA_P_S36_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_TOT_C_S17_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_MA_P_S36_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_TOT_C_S17_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_MA_P_S36_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_TOT_C_S17_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_RZ_P_S44_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m2_TOT_C_S17_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_RZ_P_S44_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_MA_P_S39_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_RZ_P_S44_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_MA_P_S39_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m12_RZ_P_S44_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_MA_P_S39_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_MA_J_S21_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_MA_P_S39_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_MA_J_S21_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_RZ_P_S47_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_MA_J_S21_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_RZ_P_S47_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_MA_J_S21_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_RZ_P_S47_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_RZ_J_S29_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m35_RZ_P_S47_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_RZ_J_S29_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m36_MA_J_S24_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_RZ_J_S29_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m36_MA_J_S24_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m13_RZ_J_S29_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m36_MA_J_S24_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_MA_C_S4_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m36_MA_J_S24_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_MA_C_S4_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m36_RZ_J_S32_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_MA_C_S4_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m36_RZ_J_S32_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_MA_C_S4_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m36_RZ_J_S32_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_RZ_C_S12_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m36_RZ_J_S32_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_RZ_C_S12_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m38_MA_P_S40_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_RZ_C_S12_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m38_MA_P_S40_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m14_RZ_C_S12_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m38_MA_P_S40_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_MA_C_S5_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m38_MA_P_S40_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_MA_C_S5_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m38_RZ_P_S48_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_MA_C_S5_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m38_RZ_P_S48_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_MA_C_S5_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m38_RZ_P_S48_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_RZ_C_S13_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m38_RZ_P_S48_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_RZ_C_S13_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_MA_J_S19_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_RZ_C_S13_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_MA_J_S19_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m15_RZ_C_S13_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_MA_J_S19_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_MA_P_S37_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_MA_J_S19_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_MA_P_S37_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_RZ_J_S27_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_MA_P_S37_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_RZ_J_S27_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_MA_P_S37_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_RZ_J_S27_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_RZ_P_S45_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_RZ_J_S27_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_RZ_P_S45_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_TOT_J_S34_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_RZ_P_S45_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_TOT_J_S34_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m16_RZ_P_S45_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m3_TOT_J_S34_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_MA_J_S18_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m3_TOT_J_S34_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_MA_J_S18_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_MA_J_S25_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_MA_J_S18_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_MA_J_S25_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_MA_J_S18_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_MA_J_S25_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_RZ_J_S26_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_MA_J_S25_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_RZ_J_S26_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_RZ_J_S33_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_RZ_J_S26_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_RZ_J_S33_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m1_RZ_J_S26_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_RZ_J_S33_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_MA_C_S6_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m40_RZ_J_S33_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_MA_C_S6_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m45_MA_P_S41_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_MA_C_S6_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m45_MA_P_S41_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_MA_C_S6_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m45_MA_P_S41_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_RZ_C_S14_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m45_MA_P_S41_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_RZ_C_S14_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m45_RZ_P_S49_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_RZ_C_S14_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m45_RZ_P_S49_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m21_RZ_C_S14_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m45_RZ_P_S49_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_MA_C_S7_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m45_RZ_P_S49_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_MA_C_S7_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m47_MA_P_S42_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_MA_C_S7_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m47_MA_P_S42_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_MA_C_S7_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m47_MA_P_S42_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_RZ_C_S15_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m47_MA_P_S42_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_RZ_C_S15_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m47_RZ_P_S50_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_RZ_C_S15_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m47_RZ_P_S50_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m22_RZ_C_S15_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m47_RZ_P_S50_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_MA_J_S22_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m47_RZ_P_S50_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_MA_J_S22_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_MA_P_S35_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_MA_J_S22_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_MA_P_S35_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_MA_J_S22_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_MA_P_S35_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_RZ_J_S30_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_MA_P_S35_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_RZ_J_S30_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_RZ_P_S43_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_RZ_J_S30_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_RZ_P_S43_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m23_RZ_J_S30_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_RZ_P_S43_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m24_MA_J_S23_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_RZ_P_S43_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m24_MA_J_S23_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_TOT_P_S51_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m24_MA_J_S23_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_TOT_P_S51_L002_R1_001.fastq.gz\u001b[0m\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[01;31m24_MA_J_S23_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_TOT_P_S51_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m24_RZ_J_S31_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m4_TOT_P_S51_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m24_RZ_J_S31_L002_R1_001.fastq.gz\u001b[0m  \u001b[01;31m9_MA_C_S2_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m24_RZ_J_S31_L003_R1_001.fastq.gz\u001b[0m  \u001b[01;31m9_MA_C_S2_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m24_RZ_J_S31_L004_R1_001.fastq.gz\u001b[0m  \u001b[01;31m9_MA_C_S2_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m26_MA_C_S8_L001_R1_001.fastq.gz\u001b[0m   \u001b[01;31m9_MA_C_S2_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m26_MA_C_S8_L002_R1_001.fastq.gz\u001b[0m   \u001b[01;31m9_RZ_C_S10_L001_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m26_MA_C_S8_L003_R1_001.fastq.gz\u001b[0m   \u001b[01;31m9_RZ_C_S10_L002_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m26_MA_C_S8_L004_R1_001.fastq.gz\u001b[0m   \u001b[01;31m9_RZ_C_S10_L003_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m26_RZ_C_S16_L001_R1_001.fastq.gz\u001b[0m  \u001b[01;31m9_RZ_C_S10_L004_R1_001.fastq.gz\u001b[0m\n",
      "\u001b[01;31m26_RZ_C_S16_L002_R1_001.fastq.gz\u001b[0m  Granek_4837_180427A5.checksum\n",
      "\u001b[01;31m26_RZ_C_S16_L003_R1_001.fastq.gz\u001b[0m  README.rtf\n"
     ]
    }
   ],
   "source": [
    "ls $DATA_DIR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Show the contents of a zipped file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "@NB501800:128:HTYKFBGX5:3:11401:19195:1019 1:N:0:ATTACTCG+AGGATAGG\n",
      "NGTAAGGGTCTGGGGGACGTTGTTGTGGATGGAAACAACCCTGGCCTCTCGGGGCTGAGCAGGGAGTCGCTTGAGC\n",
      "+\n",
      "#AAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEE/<<\n",
      "@NB501800:128:HTYKFBGX5:3:11401:24339:1020 1:N:0:ATTACTCG+AGGATAGG\n",
      "NCTCGCCTGGGGTGCTTGCGTGTCCCTGTCCTTGGCCCTGGTCACCGGGGAGATGCACATCATCTTCTTCTGGCTC\n",
      "+\n",
      "#AAAAEAAEEEEEEEEEEEEEEEE/EEE/EEEEEEEEEEEEEEEE<A/EEEE6EEE//EEEAEEEEEEE/EAEEEA\n",
      "@NB501800:128:HTYKFBGX5:3:11401:24140:1020 1:N:0:ATTACTCG+AGGATAGG\n",
      "NTGATCACGGGAGACAAGCTAGAGATGCAAAATGCCATAGTACGATCATCCCATACAAGATGGGTGAGAGAAAAAG\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTQ | head -n 10 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extracting just the sequence\n",
    "\n",
    "We need to extract every 4th line starting from the 2nd.  We cna use a tool called `sed` (stream editor) to do so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NGTAAGGGTCTGGGGGACGTTGTTGTGGATGGAAACAACCCTGGCCTCTCGGGGCTGAGCAGGGAGTCGCTTGAGC\n",
      "NCTCGCCTGGGGTGCTTGCGTGTCCCTGTCCTTGGCCCTGGTCACCGGGGAGATGCACATCATCTTCTTCTGGCTC\n",
      "NTGATCACGGGAGACAAGCTAGAGATGCAAAATGCCATAGTACGATCATCCCATACAAGATGGGTGAGAGAAAAAG\n",
      "NGGGCCTTGGTCCTGCCCTCAGCCTTGATAGGACCCTTGCCGTCGATGGGGTTACCCAAAGCGTCAACAACTCGGC\n",
      "NACCTTTCGGAGACCTCGGTGGGTCTTCCTGGGGAGACGGGTAGTACCCCATCGGGCGGTAACACCCTCGTAACCG\n",
      "NCCGAACTGTTGAAAAGAGCTTCCACGCTCCAAAGGCAGGGACAAGAGTAAAGAGCCACCACAGCCAATCTCCAAA\n",
      "NGCAAACTGAACGTCCTCAGTCTCCATTCGCTGGAAGTCCGATTCAGGACGGGCAGCGTCATCGATAGAGAATGGA\n",
      "NGCATTACCAGATGCACTGATAGCAATGATCGTCTCATCTTCCAAGTACTTGCTCACAGCGATACTATCACTCAAA\n",
      "GGCGATCTTCCTCTTGGGGAAGTGAACGTTCTTCTTGTGCTGCTTCTCGGTGATAGCAGAAGCCTCCTTGGTAGTC\n",
      "CGGCGATGATGAGAGTGAAGGGATCTTGCCACATGAGAGGGGATGATCGACGGGGACGAGGACGGATGGGGATGAT\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTQ | sed -n '2~4p;' | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Explanation**\n",
    "\n",
    "`sed -n '2~4p'`\n",
    "\n",
    "- `-n` means to suppress all output\n",
    "- `2~4p` means print (p) every 4th line starting at line 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Joining lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting sequence as a single string. We just need to delete newlline (\\n) characters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: We can break long commands uisng the line continuation character \\"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NGTAAGGGTCTGGGGGACGTTGTTGTGGATGGAAACAACCCTGGCCTCTCGGGGCTGAGCAGGGAGTCGCTTGAGCNCTCGCCTGGGGTGCTTGCGTGTCCCTGTCCTTGGCCCTGGTCACCGGGGAGATGCACATCATCTTCTTCTGGCTCNTGATCACGGGAGACAAGCTAGAGATGCAAAATGCCATAGTACGATCATCCCATACAAGATGGGTGAGAGAAAAAGNGGGCCTTGGTCCTGCCCTCAGCCTTGATAGGACCCTTGCCGTCGATGGGGTTACCCAAAGCGTCAACAACTCGGCNACCTTTCGGAGACCTCGGTGGGTCTTCCTGGGGAGACGGGTAGTACCCCATCGGGCGGTAACACCCTCGTAACCGNCCGAACTGTTGAAAAGAGCTTCCACGCTCCAAAGGCAGGGACAAGAGTAAAGAGCCACCACAGCCAATCTCCAAANGCAAACTGAACGTCCTCAGTCTCCATTCGCTGGAAGTCCGATTCAGGACGGGCAGCGTCATCGATAGAGAATGGANGCATTACCAGATGCACTGATAGCAATGATCGTCTCATCTTCCAAGTACTTGCTCACAGCGATACTATCACTCAAAGGCGATCTTCCTCTTGGGGAAGTGAACGTTCTTCTTGTGCTGCTTCTCGGTGATAGCAGAAGCCTCCTTGGTAGTCCGGCGATGATGAGAGTGAAGGGATCTTGCCACATGAGAGGGGATGATCGACGGGGACGAGGACGGATGGGGATGAT"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$FASTQ | \\\n",
    "sed -n '2~4p;' | \\\n",
    "head | \\\n",
    "tr -d '\\n' "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.Working wtih GTF files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure of a GTF file\n",
    "\n",
    "From [UCSD Genome Browser](http://genome.ucsc.edu/FAQ/FAQformat.html#format3)\n",
    "\n",
    "```\n",
    "The GTF (General Transfer Format) is an extension of GFF version 2 \n",
    "and used to represent transcription models. GFF (General Feature Format) \n",
    "consists of one line per feature, each containing 9 columns of data. \n",
    "\n",
    "Fields\n",
    "\n",
    "Fields are tab-separated. Also, all but the final field in each \n",
    "feature line must contain a value; \"empty\" columns are denoted \n",
    "with a '.'\n",
    "\n",
    "    seqname   - name of the chromosome or scaffold; chromosome names \n",
    "                without a 'chr' \n",
    "    source    - name of the program that generated this feature, or \n",
    "                the data source (database or project name)\n",
    "    feature   - feature type name. Current allowed features are\n",
    "                {gene, transcript, exon, CDS, Selenocysteine, start_codon,\n",
    "                stop_codon and UTR}\n",
    "    start     - start position of the feature, with sequence numbering \n",
    "                starting at 1.\n",
    "    end       - end position of the feature, with sequence numbering \n",
    "                starting at 1.\n",
    "    score     - a floating point value indiciating the score of a feature\n",
    "    strand    - defined as + (forward) or - (reverse).\n",
    "    frame     - one of '0', '1' or '2'. Frame indicates the number of base pairs\n",
    "                before you encounter a full codon. '0' indicates the feature \n",
    "                begins with a whole codon. '1' indicates there is an extra\n",
    "                base (the 3rd base of the prior codon) at the start of this feature.\n",
    "                '2' indicates there are two extra bases (2nd and 3rd base of the \n",
    "                prior exon) before the first codon. All values are given with\n",
    "                relation to the 5' end.\n",
    "    attribute - a semicolon-separated list of tag-value pairs (separated by a space), \n",
    "                providing additional information about each feature. A key can be\n",
    "                repeated multiple times.\n",
    "\n",
    "Attributes\n",
    "\n",
    "The following attributes are available. All attributes are semi-colon\n",
    "separated pairs of keys and values.\n",
    "\n",
    "- gene_id: The stable identifier for the gene\n",
    "- gene_version: The stable identifier version for the gene\n",
    "- gene_name: The official symbol of this gene\n",
    "- gene_source: The annotation source for this gene\n",
    "- gene_biotype: The biotype of this gene\n",
    "- transcript_id: The stable identifier for this transcript\n",
    "- transcript_version: The stable identifier version for this transcript\n",
    "- transcript_name: The symbold for this transcript derived from the gene name\n",
    "- transcript_source: The annotation source for this transcript\n",
    "- transcript_biotype: The biotype for this transcript\n",
    "- exon_id: The stable identifier for this exon\n",
    "- exon_version: The stable identifier version for this exon\n",
    "- exon_number: Position of this exon in the transcript\n",
    "- ccds_id: CCDS identifier linked to this transcript\n",
    "- protein_id: Stable identifier for this transcript's protein\n",
    "- protein_version: Stable identifier version for this transcript's protein\n",
    "- tag: A collection of additional key value tags\n",
    "- transcript_support_level: Ranking to assess how well a transcript is supported (from 1 to 5)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR=~/seqs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CM000046.1.fa         \u001b[0m\u001b[01;31mcneo040623.scaffs.fa.gz\u001b[0m\n",
      "CM000047.1.fa         Cryptococcus_neoformans.ASM9104v1.43.gtf\n",
      "CM000047.1.fa.bak     \u001b[01;31mCryptococcus_neoformans.ASM9104v1.43.gtf.gz\u001b[0m\n",
      "cneo040623.scaffs.fa\n"
     ]
    }
   ],
   "source": [
    "ls ${DATA_DIR}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inspect the GTF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#!genome-build ASM9104v1\n",
      "#!genome-version ASM9104v1\n",
      "#!genome-date 2016-06\n",
      "#!genome-build-accession GCA_000091045.1\n",
      "#!genebuild-last-updated 2016-06\n",
      "1\tena\tgene\t4839\t6236\t.\t+\t.\tgene_id \"CNA00010\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "1\tena\ttranscript\t4839\t6236\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\";\n",
      "1\tena\texon\t4839\t5178\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; exon_id \"AAW41410-1\"; exon_version \"1\";\n",
      "1\tena\tCDS\t4914\t5178\t.\t+\t0\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; protein_id \"AAW41410\";\n",
      "1\tena\tstart_codon\t4914\t4916\t.\t+\t0\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\";\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$GTF | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ignore comment lines version 1\n",
    "\n",
    "This uses `tail` to take all lines starting from line 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\tena\tgene\t4839\t6236\t.\t+\t.\tgene_id \"CNA00010\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "1\tena\ttranscript\t4839\t6236\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\";\n",
      "1\tena\texon\t4839\t5178\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; exon_id \"AAW41410-1\"; exon_version \"1\";\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$GTF | tail +6 | head -3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ingore comment lines 2\n",
    "\n",
    "This uses `grep` to ignore (`-v`) lines starting with `^#`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\tena\tgene\t4839\t6236\t.\t+\t.\tgene_id \"CNA00010\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "1\tena\ttranscript\t4839\t6236\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\";\n",
      "1\tena\texon\t4839\t5178\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; exon_id \"AAW41410-1\"; exon_version \"1\";\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$GTF | grep -v '^#' | head -3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find all unique feature types in GTF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CDS\n",
      "exon\n",
      "five_prime_utr\n",
      "gene\n",
      "start_codon\n",
      "stop_codon\n",
      "three_prime_utr\n",
      "transcript\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$GTF | grep -v '^#' | cut -f 3 | sort | uniq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get counts of all unique feature types in GTF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  44079 CDS\n",
      "  44271 exon\n",
      "   6738 five_prime_utr\n",
      "   6813 gene\n",
      "   6894 start_codon\n",
      "   2369 stop_codon\n",
      "   6718 three_prime_utr\n",
      "   7044 transcript\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$GTF | grep -v '^#' | cut -f 3 | sort | uniq -c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find the the attributes of the first 3 exons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\tena\texon\t4839\t5178\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; exon_id \"AAW41410-1\"; exon_version \"1\";\n",
      "1\tena\texon\t5241\t5377\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"2\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; exon_id \"AAW41410-2\"; exon_version \"1\";\n",
      "1\tena\texon\t5453\t5560\t.\t+\t.\tgene_id \"CNA00010\"; transcript_id \"AAW41410\"; exon_number \"3\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; exon_id \"AAW41410-3\"; exon_version \"1\";\n"
     ]
    }
   ],
   "source": [
    "zcat $DATA_DIR/$GTF | grep -v '^#' | awk '$3 == \"exon\"' | head -3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Store attributes of all genes in a file\n",
    "\n",
    "It is common to use Unix commands to do simple processing and save the results to a file, then use a scripting language such as R or Python to furhter process the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcat $DATA_DIR/$GTF | grep -v '^#' | awk '$3 == \"gene\"' | cut -f 9 > genes.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gene_id \"CNA00010\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00020\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00030\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00040\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00050\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00060\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n"
     ]
    }
   ],
   "source": [
    "head -6 genes.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Looping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop throgh first 10 lines\n",
    "\n",
    "- Creae a counter variable `LINE_NUM` and initialize it to 0\n",
    "- Redirect `genes.txt` to a while loop\n",
    "- Whlie loop reads each line into the variable `LINE`\n",
    "- The `[[ CONDITION ]]` tests if CONDITION is true\n",
    "- If LINE_NUM is less than `-lt` 10 print LINE and increment LINE_NUM\n",
    "- Note that `(( EXPR ))` is used to evaluate numeric expressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gene_id \"CNA00010\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00020\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00030\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00040\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00050\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00060\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00070\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00075\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00080\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n",
      "gene_id \"CNA00090\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n"
     ]
    }
   ],
   "source": [
    "LINE_NUM=0\n",
    "while read LINE\n",
    "do\n",
    "    if [[ $LINE_NUM -lt 10 ]]; then\n",
    "        echo $LINE\n",
    "        (( LINE_NUM++ ))\n",
    "    fi\n",
    "done < genes.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop through first 10 tRNAs\n",
    "\n",
    "Almost the same as previous program but we add a new condition to check tha the gene_biotype is tRNA using a REGULAr EXPRESSION via the match operator `=~`\n",
    "\n",
    "`.*gene_biotype[[:space:]]\\\"tRNA\\\".*` means\n",
    "\n",
    "- `.*`: Match any number of any character \n",
    "- `gene_biotype`: Match the string gene_biotype\n",
    "- `[[:space:]]`: Match a single space cahracter\n",
    "- `\\\"tRNA\\\"`: Math the string \"tRNA\" - we escape the double quotes to match linteral double quote characters\n",
    "- `.*`: Match any number of any character "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gene_id \"CNA02110\"; gene_name \"CNA02110\"; gene_source \"ena\"; gene_biotype \"tRNA\";\n",
      "gene_id \"EBG00005237496\"; gene_name \"tRNA\"; gene_source \"Rfam\"; gene_biotype \"tRNA\";\n",
      "gene_id \"CNA02190\"; gene_name \"CNA02190\"; gene_source \"ena\"; gene_biotype \"tRNA\";\n",
      "gene_id \"CNA02200\"; gene_name \"CNA02200\"; gene_source \"ena\"; gene_biotype \"tRNA\";\n",
      "gene_id \"CNA02280\"; gene_name \"CNA02280\"; gene_source \"ena\"; gene_biotype \"tRNA\";\n",
      "gene_id \"EBG00005237502\"; gene_name \"tRNA\"; gene_source \"Rfam\"; gene_biotype \"tRNA\";\n",
      "gene_id \"CNA02320\"; gene_name \"CNA02320\"; gene_source \"ena\"; gene_biotype \"tRNA\";\n",
      "gene_id \"CNA03230\"; gene_name \"CNA03230\"; gene_source \"ena\"; gene_biotype \"tRNA\";\n",
      "gene_id \"EBG00005237494\"; gene_name \"tRNA\"; gene_source \"Rfam\"; gene_biotype \"tRNA\";\n",
      "gene_id \"CNA03260\"; gene_name \"CNA03260\"; gene_source \"ena\"; gene_biotype \"tRNA\";\n"
     ]
    }
   ],
   "source": [
    "LINE_NUM=0\n",
    "while read LINE\n",
    "do\n",
    "    if [[ $LINE_NUM -lt 10 ]] && [[ $LINE =~ .*gene_biotype[[:space:]]\\\"tRNA\\\".* ]] ; then\n",
    "        echo $LINE\n",
    "        ((LINE_NUM++))\n",
    "    fi\n",
    "done < genes.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop through files in a directory and print the filename followed by the first line\n",
    "\n",
    "- Use a COMMAND EXPRESSION `$(EXPR)` to list files matching the glob pattern\n",
    "- We loop through each file using a `for` statement\n",
    "- For each file, extract the first line with `head -1` and save in variable `FIRST`\n",
    "- The `2>/dev/null` ignores output from the `stderr` stream\n",
    "- We use another COMMAND EXPRESSION to print the filename withou the directory (using `basename`)\n",
    "- We then print the first line storred in `FIRST`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR='/data/hts2018_pilot/Granek_4837_180427A5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11_MA_J_S20_L001_R1_001.fastq.gz\n",
      "@NB501800:128:HTYKFBGX5:1:11101:26760:1056 1:N:0:CGCTCATT+NGGATAGG\n",
      "\n",
      "11_MA_J_S20_L002_R1_001.fastq.gz\n",
      "@NB501800:128:HTYKFBGX5:2:11101:20327:1041 1:N:0:CGCTCATT+AGGATAGN\n",
      "\n",
      "11_MA_J_S20_L003_R1_001.fastq.gz\n",
      "@NB501800:128:HTYKFBGX5:3:11401:10896:1024 1:N:0:CGCTCATT+AGGATAGG\n",
      "\n",
      "11_MA_J_S20_L004_R1_001.fastq.gz\n",
      "@NB501800:128:HTYKFBGX5:4:11401:8259:1032 1:N:0:CGCTCATT+AGGATAGG\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for FILE in $(ls $DATA_DIR/*_MA_J_S20_*gz)\n",
    "do\n",
    "    FIRST=$(zcat $FILE 2>/dev/null | head -1)\n",
    "    echo $(basename $FILE)\n",
    "    echo $FIRST\n",
    "    echo\n",
    "done"
   ]
  }
 ],
 "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": 2
}