Unit2-RNAseq/3.2_sam_and_htseq.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RNAseq notebook 3.2: SAM files and read counting\n",
    "This notebook features examples on how to work with the sequence alignment map file, and how to derive gene counts from mapped reads.  \n",
    "**Notes**\n",
    "- full SAM files tend to be large so bash manipulation can take some time (typically minutes)\n",
    "- not all SAM attributes will be found in all SAM files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Find where a specific read aligned**  \n",
    "Example: 1000th read"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "%%bash  \n",
    "cd ../Downloads/hisat_out  \n",
    "grep SRR5454079.1000 ./SRR5454079.sam | head -n 1\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Print the first few aligned reads**\n",
    "```bash\n",
    "%%bash\n",
    "cd ../Downloads/hisat_out\n",
    "awk '/^SRR/' SRR5454079.sam | head\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Print the first few reads that aligned to the mitochondrial genome**  \n",
    "```bash\n",
    "%%bash  \n",
    "cd ../Downloads/hisat_out  \n",
    "awk '{if ($3 == \"MT\") {print}}' SRR5454079.sam | head\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises**   \n",
    "1) How many reads mapped to chromosome 20?  \n",
    "2) Find the 76th read in the SAM file. Where did it map in the human genome? Now use blastn to map the read. Do the results agree with each other?  \n",
    "3) Inspect the reference genome details in the SAM header. Beyond chromosomes, what else is included in the reference?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Check how many reads were uncounted due to multimapping (alignment not unique)**\n",
    "```bash\n",
    "%%bash\n",
    "cd Unit2-RNAseq/data\n",
    "tail SRR5454102_genecounts.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Check how many Ensembl genes have zero expression**  \n",
    "Spot and correct the mistake\n",
    "```bash\n",
    "%%bash\n",
    "cd Unit2-RNAseq/data\n",
    "awk '{if ($2 == 0) print }' SRR5454102_genecounts.txt  | wc -l\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**HOMEWORK**  \n",
    "1) Use awk to check the number of columns in the SAM file for all rows and print only the unique column counts. *HINT*: revisit Unit 1   \n",
    "2) Count how many reads from SRR5454079 mapped to chromosome 20 with 2 soft-clipped bases at the start of the read. *HINT*: Consult the SAM documentation on CIGAR strings.  \n",
    "3) Using the human transcriptome annotations by Ensembl, calculate counts per gene in the bam file for SRR5454079 and print the first 10 lines (use -s reverse)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}