puf_processing_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Original PUF processing with R1 and R2 not being mixed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "pd.set_option('display.max_colwidth', None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "folders = ['Step1','Step2','Step3','Step4','Step5','Step6_Barcode','Step6_Barcode_Norm','Step6_Indel','Step6_Indel_Norm','Step7','Step8','Step8.5','Excel','Step9','Normalized_Matrix']\n",
    "for folder in folders:\n",
    "    if not os.path.exists(folder):\n",
    "        os.mkdir(folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Move all raw fastq files to Step1 folder\n",
    "# If unziped, comment out the following line\n",
    "# !gunzip ./Step1/*.gz\n",
    "files = [f'./Step1/{f}' for f in listdir('./Step1/') if isfile(join('./Step1/', f)) and f.endswith('.fastq')]\n",
    "for file in files:\n",
    "    os.system(f'awk \\'NR%4 ==2\\' {file} > ./Step2/{file[file.rfind(\"/\")+1:]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = [f'./Step2/{f}' for f in listdir('./Step2/') if isfile(join('./Step2/', f))]\n",
    "for file in files:\n",
    "    os.system(f'sed -e \"s/^.//g\" {file} > ./Step3/{file[file.rfind(\"/\")+1:]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "files_f = [f'./Step3/{f}' for f in listdir('./Step3/') if isfile(join('./Step3/', f)) and f.find(\"R1\") != -1]\n",
    "for file in files_f:\n",
    "    file1 = file\n",
    "    file2 = file.replace(\"R1\", \"R2\")\n",
    "    os.system(f\"paste -d '\\\\0' {file1} {file2} | cat> ./Step4/{file[file.rfind('/')+1:]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0-Pe_R1_001.fastq\n",
      "G0-Pe_R1_001.fastq\n",
      "F0-Pe-repeat_R1_001.fastq\n",
      "F0-Pe_R1_001.fastq\n",
      "H0-FT-repeat_R1_001.fastq\n",
      "G0-FT-repeat_R1_001.fastq\n",
      "H0-FT_R1_001.fastq\n",
      "G0-FT_R1_001.fastq\n",
      "F0-FT-repeat_R1_001.fastq\n",
      "F0-FT_R1_001.fastq\n",
      "G0-Pe-repeat_R1_001.fastq\n",
      "H0-Pe-repeat_R1_001.fastq\n"
     ]
    }
   ],
   "source": [
    "files = [f'./Step4/{f}' for f in listdir('./Step4/') if isfile(join('./Step4/', f))]\n",
    "for file in files:\n",
    "    print(file[file.rfind('/')+1:])\n",
    "    os.system(f'grep \"^CTTATATTCCCAGGGCCGGTTCGCGATCGCCCTGCAGG[A-Z][A-Z][A-Z][A-Z][A-Z]TAGTTATTAATGACTCACGGGGATTTCCAAGTCTCCACCCCATTGACGTCAATGGG.*CCGCCCTCGACCGCCTTGATTCTCATGGTCTGGGTGC[A-Z]*GTGGTGGTTGTTCACGGTGCCCT\" {file} > ./Step5/{file[file.rfind(\"/\")+1:]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0-Pe_R1_001.fastq\n",
      "G0-Pe_R1_001.fastq\n",
      "F0-Pe-repeat_R1_001.fastq\n",
      "F0-Pe_R1_001.fastq\n",
      "H0-FT-repeat_R1_001.fastq\n",
      "G0-FT-repeat_R1_001.fastq\n",
      "H0-FT_R1_001.fastq\n",
      "G0-FT_R1_001.fastq\n",
      "F0-FT-repeat_R1_001.fastq\n",
      "F0-FT_R1_001.fastq\n",
      "G0-Pe-repeat_R1_001.fastq\n",
      "H0-Pe-repeat_R1_001.fastq\n"
     ]
    }
   ],
   "source": [
    "#raw barcode and indel\n",
    "files = [f'./Step5/{f}' for f in listdir('./Step5/') if isfile(join('./Step5/', f))]\n",
    "for file in files:\n",
    "    print(file[file.rfind('/')+1:])\n",
    "    os.system(f'sed -e \"s/CTTATATTCCCAGGGCCGGTTCGCGATCGCCCTGCAGG\\(.*\\)TAGTTATTAATGACTCACGGGGATTTCCAAGTCTCCACCCCATTGACGTCAATGGG.*CCGCCCTCGACCGCCTTGATTCTCATGGTCTGGGTGC[A-Z]*GTGGTGGTTGTTCACGGTGCCCT[A-Z]*/\\\\1/\"  {file} > ./Step6_Barcode/{file[file.rfind(\"/\")+1:]}')\n",
    "    os.system(f'sed -e \"s/CTTATATTCCCAGGGCCGGTTCGCGATCGCCCTGCAGG[A-Z][A-Z][A-Z][A-Z][A-Z]TAGTTATTAATGACTCACGGGGATTTCCAAGTCTCCACCCCATTGACGTCAATGGG.*CCGCCCTCGACCGCCTTGATTCTCATGGTCTGGGTGC\\(.*\\)GTGGTGGTTGTTCACGGTGCCCT[A-Z]*/\\\\1/\"  {file} > ./Step6_Indel/{file[file.rfind(\"/\")+1:]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F0-FT-repeat_R1_001.fastq\n",
      "F0-FT_R1_001.fastq\n",
      "F0-Pe-repeat_R1_001.fastq\n",
      "F0-Pe_R1_001.fastq\n",
      "G0-FT-repeat_R1_001.fastq\n",
      "G0-FT_R1_001.fastq\n",
      "G0-Pe-repeat_R1_001.fastq\n",
      "G0-Pe_R1_001.fastq\n",
      "H0-FT-repeat_R1_001.fastq\n",
      "H0-FT_R1_001.fastq\n",
      "H0-Pe-repeat_R1_001.fastq\n",
      "H0-Pe_R1_001.fastq\n"
     ]
    }
   ],
   "source": [
    "#barcode indel combo\n",
    "files_barcode = [f'./Step6_Barcode/{f}' for f in listdir('./Step6_Barcode/') if isfile(join('./Step6_Barcode/', f))]\n",
    "files_barcode.sort()\n",
    "files_indel = [f'./Step6_Indel/{f}' for f in listdir('./Step6_Indel/') if isfile(join('./Step6_Indel/', f))]\n",
    "files_indel.sort()\n",
    "for file_barcode, file_indel in zip(files_barcode, files_indel):\n",
    "    print(file_barcode[file_barcode.rfind('/')+1:])\n",
    "    os.system(f\"paste -d '\\\\0' {file_barcode} {file_indel} | cat > ./Step7/{file_barcode[file_barcode.rfind('/')+1:]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F0-FT-repeat_R1_001.fastq\n",
      "F0-FT_R1_001.fastq\n",
      "F0-Pe-repeat_R1_001.fastq\n",
      "F0-Pe_R1_001.fastq\n",
      "G0-FT-repeat_R1_001.fastq\n",
      "G0-FT_R1_001.fastq\n",
      "G0-Pe-repeat_R1_001.fastq\n",
      "G0-Pe_R1_001.fastq\n",
      "H0-FT-repeat_R1_001.fastq\n",
      "H0-FT_R1_001.fastq\n",
      "H0-Pe-repeat_R1_001.fastq\n",
      "H0-Pe_R1_001.fastq\n"
     ]
    }
   ],
   "source": [
    "files_barcode = [f'./Step6_Barcode/{f}' for f in listdir('./Step6_Barcode/') if isfile(join('./Step6_Barcode/', f))]\n",
    "normalizer = [f'./Step5/{f}' for f in listdir('./Step5/') if isfile(join('./Step5/', f)) and f.find(\"R1\") != -1]\n",
    "files_barcode.sort()\n",
    "normalizer.sort()\n",
    "for file_barcode,normal in zip(files_barcode,normalizer):\n",
    "    print(file_barcode[file_barcode.rfind('/')+1:])\n",
    "    # In Step6_Barcode_Normalized, create frequency table for each barcode. Have an extra column for normalized frequency\n",
    "    df = pd.read_csv(file_barcode, sep='\\t', header=None)\n",
    "    normal = pd.read_csv(normal, header=None)\n",
    "    df.columns = ['barcode']\n",
    "    df['count'] = 1\n",
    "    df = df.groupby('barcode').count().reset_index()\n",
    "    # Set index to barcode\n",
    "    df = df.set_index('barcode')\n",
    "    # Sort from highest to lowest\n",
    "    df = df.sort_values(by='count', ascending=False)\n",
    "    df['normalized'] = df['count']/len(normal.index)*1000000\n",
    "    df.to_csv(f'./Step6_Barcode_Norm/{file_barcode[file_barcode.rfind(\"/\")+1:]}', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F0-FT-repeat_R1_001.fastq\n",
      "F0-FT_R1_001.fastq\n",
      "F0-Pe-repeat_R1_001.fastq\n",
      "F0-Pe_R1_001.fastq\n",
      "G0-FT-repeat_R1_001.fastq\n",
      "G0-FT_R1_001.fastq\n",
      "G0-Pe-repeat_R1_001.fastq\n",
      "G0-Pe_R1_001.fastq\n",
      "H0-FT-repeat_R1_001.fastq\n",
      "H0-FT_R1_001.fastq\n",
      "H0-Pe-repeat_R1_001.fastq\n",
      "H0-Pe_R1_001.fastq\n"
     ]
    }
   ],
   "source": [
    "files_indel = [f'./Step6_Indel/{f}' for f in listdir('./Step6_Indel/') if isfile(join('./Step6_Indel/', f))]\n",
    "normalizer = [f'./Step5/{f}' for f in listdir('./Step5/') if isfile(join('./Step5/', f)) and f.find(\"R1\") != -1]\n",
    "files_indel.sort()\n",
    "normalizer.sort()\n",
    "for file_indel,normal in zip(files_indel,normalizer):\n",
    "    print(file_indel[file_indel.rfind('/')+1:])\n",
    "    # In Step6_Indel_Normalized, create frequency table for each indel. Have an extra column for normalized frequency\n",
    "    df = pd.read_csv(file_indel, sep='\\t', header=None)\n",
    "    normal = pd.read_csv(normal, header=None)\n",
    "    df.columns = ['indel']\n",
    "    df['count'] = 1\n",
    "    df = df.groupby('indel').count().reset_index()\n",
    "    # Set index toindel \n",
    "    df = df.set_index('indel')\n",
    "    # Sort from highest to lowest\n",
    "    df = df.sort_values(by='count', ascending=False)\n",
    "    df['normalized'] = df['count']/len(normal.index)*1000000\n",
    "    df.to_csv(f'./Step6_Indel_Norm/{file_indel[file_indel.rfind(\"/\")+1:]}', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0-Pe_R1_001.fastq\n",
      "G0-Pe_R1_001.fastq\n",
      "F0-Pe-repeat_R1_001.fastq\n",
      "F0-Pe_R1_001.fastq\n",
      "H0-FT-repeat_R1_001.fastq\n",
      "G0-FT-repeat_R1_001.fastq\n",
      "H0-FT_R1_001.fastq\n",
      "G0-FT_R1_001.fastq\n",
      "F0-FT-repeat_R1_001.fastq\n",
      "F0-FT_R1_001.fastq\n",
      "G0-Pe-repeat_R1_001.fastq\n",
      "H0-Pe-repeat_R1_001.fastq\n"
     ]
    }
   ],
   "source": [
    "#indel filter\n",
    "files = [f'./Step7/{f}' for f in listdir('./Step7/') if isfile(join('./Step7/', f))]\n",
    "for file in files:\n",
    "    print(file[file.rfind('/')+1:])\n",
    "    os.system(f\"grep -v -x '.\\{{45\\}}' {file}  | cat >  ./Step8/{file[file.rfind('/')+1:]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n",
      "/var/folders/9m/qcg58ft12w72z4h1j7s196w00000gn/T/ipykernel_88831/4136541351.py:19: FutureWarning: save is not part of the public API, usage can give unexpected results and will be removed in a future version\n",
      "  writer.save()\n"
     ]
    }
   ],
   "source": [
    "normalizer = [f'./Step5/{f}' for f in listdir('./Step5/') if isfile(join('./Step5/', f)) and f.find(\"R1\") != -1]\n",
    "normalizer.sort()\n",
    "barcode_indels = [f'./Step8/{f}' for f in listdir('./Step8/') if isfile(join('./Step8/', f)) and f.find(\"R1\") != -1]\n",
    "barcode_indels.sort()\n",
    "for file,normal in zip(barcode_indels,normalizer):\n",
    "    writer = pd.ExcelWriter(f'./Excel/{file[file.rfind(\"/\")+1:]}.xlsx', engine='openpyxl') \n",
    "    sheet_names = ['Barcode','Indel','Barcode_Indel']\n",
    "    df = pd.read_csv(file, header=None)\n",
    "    normal = pd.read_csv(normal, header=None)\n",
    "    barcode_indel_freq = df.value_counts()\n",
    "    barcode_indel_str = pd.DataFrame(list(barcode_indel_freq.index))\n",
    "    barcode_indel_freq = barcode_indel_freq.reset_index(drop=True)\n",
    "    df = pd.DataFrame()\n",
    "    df['Barcode'] = barcode_indel_str[0].str[:5] \n",
    "    df['Indel'] = barcode_indel_str[0].str[5:]\n",
    "    df['Freq'] = barcode_indel_freq/len(normal.index)*1000000\n",
    "    df.to_excel(writer, sheet_name=f'0',index=False)\n",
    "    df.to_csv(f'./Step9/{file[file.rfind(\"/\")+1:]}',index=False)\n",
    "    writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reverse_complement(seq):\n",
    "    seq_dict = {'A':'T','T':'A','G':'C','C':'G','N':'N'}\n",
    "    return \"\".join([seq_dict[base] for base in reversed(seq)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "normed_freqs = [f'./Step9/{f}' for f in listdir('./Step9/') if isfile(join('./Step9/', f))]\n",
    "normed_freqs.sort()\n",
    "for file in normed_freqs:\n",
    "    table = pd.read_csv(file)\n",
    "    table['Indel'] = table['Indel'].apply(reverse_complement)\n",
    "    table_pivot = table.pivot_table(index=['Barcode'],columns=['Indel'], values=['Freq'], aggfunc='sum')\n",
    "    barcode = table[['Barcode','Freq']]\n",
    "    barcode = barcode.groupby('Barcode').sum()\n",
    "    barcode.sort_values(by='Freq', ascending=False, inplace=True)\n",
    "    indel = table[['Indel','Freq']]\n",
    "    indel = indel.groupby('Indel').sum()\n",
    "    indel.sort_values(by='Freq', ascending=False, inplace=True)\n",
    "    table_pivot = table_pivot.fillna(0)\n",
    "    table_pivot = table_pivot['Freq']\n",
    "    table_pivot = table_pivot.reindex(barcode.index, axis=0)\n",
    "    table_pivot = table_pivot.reindex(indel.index,axis=1)\n",
    "    #Write table_pivot to excel\n",
    "    table_pivot.to_csv(f'./Normalized_Matrix/{file[file.rfind(\"/\")+1:]}.csv',index=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9.15"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "224d2b38dd5c4bf9bb346b4b7ffc570356b5a9eca7c29d4f05461d7abc1d8f91"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}