tnseq_process_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tn-seq Sequencing Data Processing Pipeline\n",
    "\n",
    "Jennifer Stiens\n",
    "2023\n",
    "j.j.stiens@gmail.com\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview of Process\n",
    "\n",
    "1. Download from genewiz and sanity check QC to make sure downloads intact\n",
    "2. FastQC (optional)\n",
    "3. Illumina adapter trimming and QC with Fastp\n",
    "4. Barcode from index read added to read heading of read 1 fastq files (read converted into 2-line reads, \".reads\")\n",
    "5. Trimming of the transposon tag\n",
    "6. Mapping with BWA-mem\n",
    "7. TA-site quantification, removal of PCR duplicates\n",
    "8. Generation of insertion files (.wig)\n",
    "\n",
    "This protocol uses primers based on those used in Mendum et al, 2019 (BMC Genomics) where the barcode is included in the i7 index read.\n",
    "The primers are included as a spreadsheet, "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Necessary modules and scripts\n",
    "\n",
    "tnseq_pro.py module contains all the python functions for processing the fastq and sam files.\n",
    "\n",
    "snakemake/tnseq/ directory contains the necessary snakemake scripts for fastp and mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download from genewiz and sanity checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "#check length of files to see if reads same in all files?\n",
    "cd fastq\n",
    "FILES=*.fastq.gz\n",
    "for file in $FILES; do wc -l $file; done >> sanity_check.txt\n",
    "\n",
    "# check files downloaded correctly\n",
    "for file in $FILES; do md5sum -c $file.md5; done >> md5_check.txt\n",
    "\n",
    "#head\n",
    "for file in $FILES; do echo $file; zcat $file | head -10 $file; done\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FastQC on all reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "cd fastq\n",
    "module load fastqc\n",
    "module load multiqc\n",
    "FILES=*.fastq.gz\n",
    "for f in $FILES; do fastqc ${f} -o fastqc; done\n",
    "cd fastqc\n",
    "multiqc .\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fastp for quality control and Illumina adapter trimming of read 1 only \n",
    "\n",
    "Don't allow automatic detection or will trim transposon tags which are needed for insertion site confirmation and position. This shouldn't actually trim that many files as most were probably trimmed by genewiz processing.\n",
    "\n",
    "config.yaml must be present in working directory"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "#config.yaml\n",
    "samples: [A1_R1_001, A2_R1_001, ...]\n",
    "PATH: <path to fastq file directory>\n",
    "genomeFile: ref_seqs/Mbovis_AF2122-97.fasta\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "cd fastp\n",
    "conda activate snakemake\n",
    "snakemake -np -s ~/snakemake/tnseq/fastp/snakefile.smk\n",
    "snakemake --cores 4 -s ~/snakemake/tnseq/fastp/snakefile.smk\n",
    "\n",
    "#on server\n",
    "#!nohup snakemake --cores 8 -s $my_path/snakemake/tnseq/fastp/snakefile.smk > nohup_fastp.out 2>&1 &\n",
    "\n",
    "#bash command for each file\n",
    "#fastp -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC -i fastq/<infile.fastq> -o trimmed/<outfile.fastq>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Move barcodes from P7 index read to the name portion of the header of R1 reads. \n",
    "\n",
    "P7 indexes were delivered as R2 fastq files from genewiz.\n",
    "\n",
    "Reads will be converted from trimmed reads in fastq format to fasta format (header and sequence only). This takes around 5 minutes per file on laptop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['tests/trimmed_reads/trimmed_test_5000_R1_001.fastq']\n",
      "tests/barcoded_reads/barcode_trimmed_test_5000_R1_001.reads\n"
     ]
    }
   ],
   "source": [
    "# iterate through files\n",
    "import scripts.tnseq_pro as tn\n",
    "\n",
    "#read all files in directory\n",
    "#trimmed_files = \"fastq/trimmed/\"\n",
    "trimmed_files = \"tests/trimmed_reads/\"\n",
    "barcoded_dir = \"tests/barcoded_reads\"\n",
    "tn.iterate_add_barcode(trimmed_files, barcoded_dir)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove transposon tag from reads with tag in first 20 bases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File being processed:  tests/barcoded_reads/barcode_trimmed_test_5000_R1_001.reads\n"
     ]
    }
   ],
   "source": [
    "import scripts.tnseq_pro as tn\n",
    "\n",
    "tn.iterate_tag_trim(\"tests/barcoded_reads\", \"tests/output/\")\n",
    "\n",
    "#tn.trim_tag_fastq(\"barcoded_output/<barcoded_fastq_file.reads>\", \"barcoded_output/\", tag=\"ACTTATCAGCCAACCTGTTA\", mismatch_max=2)\n",
    "#tn.trim_tag_fastq(\"tests/barcoded/barcode_trimmed_test_5000_R1_001.reads\", \"tests/barcoded\", tag=\"ACTTATCAGCCAACCTGTTA\", mismatch_max=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map tag-clipped, barcoded reads with BWA-mem\n",
    "\n",
    "Use snakemake pipeline to:\n",
    "\n",
    "1. map (bwa-mem)\n",
    "2. sort (samtools)\n",
    "3. index (samtools)\n",
    "4. filter for mapped reads only (samtools)\n",
    "5. compile read statistics (samtools)\n",
    "\n",
    "BWA generated index file (.fasta.fai) must be present in same file as reference fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "conda activate tnseq\n",
    "# have to index file first (generate .fasta.fai) for bwa-mem\n",
    "bwa-mem2 index ref_seqs/Mbovis_AF2122-97.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# with snakemake (maps, sorts, indexes and creates flagstats report)\n",
    "#check config.yaml file\n",
    "\n",
    "cd ~/tn_seq/menadione_tnseq/\n",
    "conda activate snakemake\n",
    "snakemake -np -s ~/snakemake/tnseq/mapping/snakefile.smk\n",
    "snakemake --cores 2 -s ~/snakemake/tnseq/mapping/snakefile.smk\n",
    "\n",
    "#on server\n",
    "snakemake -np -s $my_path/snakemake/tnseq/mapping/snakefile.smk\n",
    "nohup snakemake --cores 8 -s $my_path/snakemake/tnseq/mapping/snakefile.smk > nohup_map.out 2>&1 &"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sorted, indexed and filtered reads are in 'sorted_reads' directory along with \"flagstat\" files which can be used for mapping statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantification and filtering\n",
    "\n",
    "One function: 'sam_to_wig' is used for:\n",
    "\n",
    "1. Reducing reads to unique template counts (eliminating duplicates based on same barcode/insertion)\n",
    "2. Matching unique mapped reads to insertion site (TA-site)\n",
    "3. Quantifying number of reads per TA-site in genome and creating .wig file\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['tests/sorted_reads/mapped_test_5000_R1_001.sam']\n",
      "test_5000\n",
      "tests/sorted_reads/mapped_test_5000_R1_001.sam\n",
      "number of unique reads assigned to TA sites:  1145\n",
      "number of unique reads with no ta site match:  12\n"
     ]
    }
   ],
   "source": [
    "import scripts.tnseq_pro as tn\n",
    "tn.iterate_sam_to_wig(\"tests/sorted_reads\", \"tests/output\", \"ref_seqs/Mbovis_AF2122-97.fasta\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyse insertion files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File being processed:  tests/output/test_5000_insertions.wig\n"
     ]
    }
   ],
   "source": [
    "import scripts.tnseq_pro as tn  \n",
    "import glob\n",
    "wig_dir = \"tests/output\"\n",
    "wig_files = glob.glob(wig_dir + \"/*.wig\")\n",
    "for file in wig_files:\n",
    "        print(\"File being processed: \", file)\n",
    "        tn.analyze_dataset(file)\n",
    "        "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dna",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}