search_fastq_for_sequences.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3d5086b6",
   "metadata": {},
   "source": [
    "# Seach fastq for sequences\n",
    "This notebook will help you determine if certain sequences exist in your fastq file."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afe8990e",
   "metadata": {},
   "source": [
    "### Usage\n",
    "Please provide file paths for a fastq file and a list of sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "204647be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run if dnaio is not installed\n",
    "!pip install dnaio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "d7433107",
   "metadata": {},
   "outputs": [],
   "source": [
    "### INPUTS ###\n",
    "fastq_file_path = \"data/FX4-230125-02_S4_L001_R2_001.fastq.gz\"\n",
    "sequence_list_path = \"data/list_of_sequences.csv\"\n",
    "output_file_path = \"data/FX4-230125-02_S4_L001_R2_001.output.not_parallel.csv\"\n",
    "njobs = None\n",
    "reverse = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "db56e44f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dnaio\n",
    "import csv\n",
    "from joblib import Parallel, delayed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "fe533154",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''Read and store list of sequences'''\n",
    "file = open(sequence_list_path, \"r\")\n",
    "sequence_list = file.read().split(\"\\n\")\n",
    "file.close()\n",
    "\n",
    "# Will only run if you would like to search for the reverse of your sequences\n",
    "if reverse:\n",
    "    reverse_sequences = []\n",
    "    for sequence in sequence_list:\n",
    "        reverse_sequences.append(sequence[::-1])\n",
    "    sequence_list = reverse_sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "4417c363",
   "metadata": {},
   "outputs": [],
   "source": [
    "def search_reads(read_number: int, read_sequence: str, sequence_list: list, output_file_path: str):\n",
    "    ''' Search a read sequence for sequences in list\n",
    "    Meant to be run in parallel\n",
    "    '''\n",
    "    for sequence in sequence_list:\n",
    "        if sequence in read_sequence:\n",
    "            # Open output file to append\n",
    "            output_file = open(output_file_path, \"a\")\n",
    "            output_file.write(f\"{read_number},{read_sequence},{sequence},{read_sequence.find(sequence)}\\n\")\n",
    "            output_file.close()\n",
    "            break #In theory there should only be one sequence per read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "b3569889",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''Search fastq file for sequences'''\n",
    "\n",
    "# Open output file and clear\n",
    "output_file = open(output_file_path, \"w\")\n",
    "output_file.write(\"read_number,read_sequence,pattern_found,nucleotide_index\\n\")\n",
    "output_file.close()\n",
    "\n",
    "# Set njobs if not set\n",
    "njobs = -1 if not njobs else njobs\n",
    "\n",
    "# Open fastq and run in parallel\n",
    "with dnaio.open(fastq_file_path) as f:\n",
    "    Parallel(n_jobs=njobs, backend=\"multiprocessing\")(delayed(search_reads)(read_number=i, read_sequence=read.sequence, sequence_list=sequence_list, output_file_path=output_file_path) for i, read in enumerate(f))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42107114",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "conda_python3",
   "language": "python",
   "name": "conda_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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}