s4U - Notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b6df2232",
   "metadata": {},
   "source": [
    "# Outline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e36bc3f6",
   "metadata": {},
   "source": [
    "# I. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b3a6a2a",
   "metadata": {},
   "source": [
    "# II. Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eed7144e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Following script saved as preprocessing.sh\n",
    "# Command: sbatch preprocessing.sh experiment keyfile\n",
    "\n",
    "###################################################################################################################\n",
    "#! /bin/bash\n",
    "\n",
    "#SBATCH --job-name wil_works\n",
    "#SBATCH --ntasks=4\n",
    "#SBATCH --cpus-per-task=4\n",
    "#SBATCH --partition=ccr\n",
    "#SBATCH --time=24:00:00\n",
    "#SBATCH --gres=lscratch:500 \n",
    "#SBATCH --mem=50g\n",
    "\n",
    "\n",
    "experiment=$1\n",
    "keyfile=$2\n",
    "reference=$3\n",
    "\n",
    "\n",
    "# Extract tarball\n",
    "tar -xvf ${experiment}.tar\n",
    "\n",
    "\n",
    "# Base Calling\n",
    "module load bcl2fastq/2.20.0 || exit 1\n",
    "bcl2fastq --runfolder-dir ${experiment} -o ${experiment}_fastq_dir -r 4 -w 4 -p 14 --barcode-mismatches 0\n",
    "\n",
    "# Concatenate all reads in one FASTQ file\n",
    "zcat ${experiment}_fastq_dir/*.fastq.gz > ${experiment}.fastq\n",
    "\n",
    "# Demultiplex FASTQ file\n",
    "perl /data/BatistaLab_NGS/Scripts/icSHAPE-master/scripts/splitFastq.pl -U ${experiment}.fastq -b 6:6 -l CGTGAT:1::ACATCG:2::GCCTAA:3::TGGTCA:4::CACTGT:5::ATTGGC:6::GATCTG:7::TCAAGT:8::CTGATC:9::AAGCTA:10::GTAGCC:11::TACAAG:12::TTGACT:13::GGAACT:14::TGACAT:15::GGACGG:16::GCGGAC:17::TTTCAC:18::GGCCAC:19::CGAAAC:20::CGTACG:21::CCACTC:22::ATCAGT:23::AGGAAT:24 -d ${experiment}.splitFastq\n",
    "\n",
    "# Compress FASTQ file\n",
    "gzip -f ${experiment}.fastq\n",
    "\n",
    "# Remove PCR duplicates\n",
    "cd ${experiment}.splitFastq/\n",
    "for file in *fastq; do perl /data/BatistaLab_NGS/Scripts/icSHAPE-master/scripts/readCollapse.pl -U ${file} ; done \n",
    "\n",
    "# Remove adapters\n",
    "for file in *collapsed.fastq; do perl /data/BatistaLab_NGS/Scripts/icSHAPE-master/scripts/trimming.pl -U ${file} -o ${file}.trimm.fastq -l 17 -t 0 -c phred33 -a /data/BatistaLab_NGS/Scripts/s4U_pipeline/TruSeq2-PE.fa -m 20; done \n",
    "cd ..\n",
    "\n",
    "# Rename FASTQ files according to KeyFile\n",
    "while IFS=$'\\t' read -r -a myArray\n",
    "do\n",
    " cp  ${experiment}.splitFastq/${myArray[0]}.collapsed.fastq.trimm.fastq ${experiment}.splitFastq/${myArray[1]}.fastq\n",
    "done < ${keyfile}\n",
    "\n",
    "rm ${experiment}.splitFastq/*collapsed.fastq.trimm.fastq\n",
    "\n",
    "# Prepare reference genome and  Align FASTQ files\n",
    "module load bowtie\n",
    "\n",
    "if ${reference} == \"Ecoli\"\n",
    "then\n",
    "  bowtie2-build /data/BatistaLab_NGS/Scripts/s4U_pipeline/transcriptome/transcriptome_uniquetRNA.fa Ecoli_transcriptome\n",
    "  for file in ${experiment}.splitFastq/*rep*.fastq ; do bowtie2 -p 16 -x Ecoli_transcriptome -U ${file} -S ${file}.bam ; done\n",
    "fi\n",
    "\n",
    "if ${reference} == \"Paeruginosa\"\n",
    "then\n",
    "  bowtie2-build /data/BatistaLab_NGS/Scripts/s4U_pipeline/transcriptome/Pseudomonas_aeruginosa_transcriptome_tRNAunique.fa PA_transcriptome\n",
    "  for file in ${experiment}.splitFastq/*rep*.fastq ; do bowtie2 -p 16 -x PA_transcriptome -U ${file} -S ${file}.bam ; done\n",
    "fi\n",
    "\n",
    "\n",
    "module load samtools\n",
    "\n",
    "# Keep mapped reads\n",
    "for file in ${experiment}.splitFastq/*fastq.bam ; do samtools view -F 4 -b -o ${file}.mapped.bam ${file} ; done\n",
    "\n",
    "# Sort reads \n",
    "for file in ${experiment}.splitFastq/*.mapped.bam ; do samtools sort -o ${file}.sorted.bam ${file} ; done\n",
    "\n",
    "# Create index\n",
    "for file in ${experiment}.splitFastq/*.sorted.bam ; do samtools index ${file} ; done\n",
    "\n",
    "# ??\n",
    "    for file in ${experiment}.splitFastq/*fastq.bam.mapped.bam; do samtools view ${file} |  grep -E \"@|NM:\" |  grep -v \"XS:\" | cut -f 3|  sort | uniq -c | awk '{printf(\"%s\\t%s\\n\", $2, $1)}' > ${file}.count ; done\n",
    "\n",
    "# Rename COUNT files according to KeyFile\n",
    "while IFS=$'\\t' read -r -a myArray\n",
    "do\n",
    " cp ${experiment}.splitFastq/${myArray[1]}.fastq.bam.mapped.bam.count ${experiment}.splitFastq/${myArray[1]}.count\n",
    "done < ${keyfile}\n",
    "\n",
    "####################################################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f32b98a",
   "metadata": {},
   "source": [
    "# III. Misincorporation HeatMaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accfc891",
   "metadata": {},
   "outputs": [],
   "source": [
    "#parse_pileups.py\n",
    "\n",
    "import sys\n",
    "import re\n",
    "\n",
    "infile = open(sys.argv[1], 'rt')\n",
    "infile.readline()\n",
    "outfile = open(sys.argv[2], 'w+')\n",
    "\n",
    "# List tRNAs ID from reference\n",
    "tRNAs = ['EBT00001514652','EBT00001514586','EBT00001514664','EBT00001514749','EBT00001514687','EBT00001514602',\\\n",
    "         'EBT00001514635','EBT00001514579','EBT00001514747','EBT00001514671','EBT00001514668','EBT00001514736',\\\n",
    "         'EBT00001514738','EBT00001514720','EBT00001514654','EBT00001514578','EBT00001514721','EBT00001514692',\\\n",
    "         'EBT00001514633','EBT00001514638','EBT00001514676','EBT00001514740','EBT00001514753','EBT00001514655',\\\n",
    "         'EBT00001514732','EBT00001514606','EBT00001514688','EBT00001514601','EBT00001514622','EBT00001514729',\\\n",
    "         'EBT00001514735','EBT00001514716','EBT00001514697','EBT00001514650','EBT00001514699','EBT00001514691',\\\n",
    "         'EBT00001514582','EBT00001514702','EBT00001514593','EBT00001514626','EBT00001514730','EBT00001514592',\\\n",
    "         'EBT00001514696','EBT00001514711','EBT00001514647','EBT00001514636','EBT00001514709','EBT00001514679',\\\n",
    "         'EBT00001514624','EBT00001514739','EBT00001514751','EBT00001514722','EBT00001514600','EBT00001514728',\\\n",
    "         'EBT00001514597','EBT00001514667','EBT00001514651','EBT00001514715','EBT00001514737','EBT00001514613',\\\n",
    "         'EBT00001514723','EBT00001514693','EBT00001514637','EBT00001514632','EBT00001514733','EBT00001514743',\\\n",
    "         'EBT00001514689','EBT00001514584','EBT00001514690','EBT00001514678','EBT00001514596','EBT00001514707',\\\n",
    "         'EBT00001514627','EBT00001514591','EBT00001514734','EBT00001514641','EBT00001514612','EBT00001514598',\\\n",
    "         'EBT00001514684','EBT00001514719','EBT00001514665','EBT00001514621','EBT00001514583','EBT00001514604',\\\n",
    "         'EBT00001514683','EBT00001514642', # E coli\n",
    "\n",
    "         'EBT00051365344','EBT00051365099','EBT00051370091','EBT00051369932','EBT00051369986','EBT00051367609',\\\n",
    "         'EBT00051369609','EBT00051366328','EBT00051367914','EBT00051366603','EBT00051366399','EBT00051365908',\\\n",
    "         'EBT00051369492','EBT00051366465','EBT00051366923','EBT00051368550','EBT00051368683','EBT00051365503',\\\n",
    "         'EBT00051369841','EBT00051365018','EBT00051367353','EBT00051365422','EBT00051365808','EBT00051368205',\\\n",
    "         'EBT00051369258','EBT00051365670','EBT00051364934','EBT00051365187','EBT00051369657','EBT00051366843',\\\n",
    "         'EBT00051366536','EBT00051368279','EBT00051367683','EBT00051368081','EBT00051367446','EBT00051370188',\\\n",
    "         'EBT00051369313','EBT00051369436','EBT00051370040','EBT00051370146','EBT00051367983','EBT00051367160',\\\n",
    "         'EBT00051366161'] # P a\n",
    "\n",
    "# Extract relevant information from pileup files for tRNA only\n",
    "for line in infile:\n",
    "    if line[0] != '#':\n",
    "        array = line.strip().split('\\t')\n",
    "        ID, position, ref, alt, INFO = array[0], array[1], array[3], array[4], array[7]\n",
    "\n",
    "        if ID in tRNAs:\n",
    "            alt = alt.split(',')\n",
    "            x = re.search(\"DP=([0-9]+);AD=([0-9,\\,]+)\", INFO)\n",
    "            DP, AD = x.group(1), x.group(2)\n",
    "            AD = AD.split(',')\n",
    "\n",
    "            outfile.write(ID+'\\t'+position+'\\t'+DP+'\\t'+ref+'\\t'+str(float(AD[0])/float(DP))+'\\t')\n",
    "            \n",
    "            for i in range(0, len(alt)):\n",
    "                outfile.write(alt[i]+'\\t'+str(float(AD[i+1])/float(DP))+'\\t')\n",
    "\n",
    "            outfile.write('\\n')\n",
    "            \n",
    "            \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cebac58",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Following script saved as pileup_analysis.sh\n",
    "# Command: sbatch pileup_analysis.sh experiment\n",
    "\n",
    "###################################################################################################################\n",
    "#! /bin/bash\n",
    "\n",
    "#SBATCH --job-name wil_works\n",
    "#SBATCH --ntasks=4\n",
    "#SBATCH --cpus-per-task=4\n",
    "#SBATCH --partition=ccr\n",
    "#SBATCH --time=24:00:00\n",
    "#SBATCH --gres=lscratch:500 \n",
    "#SBATCH --mem=50g\n",
    "\n",
    "experiment=$1\n",
    "reference=$2\n",
    "\n",
    "# Make pileups\n",
    "module load samtools\n",
    "#for file in ${experiment}.splitFastq/*bam.sorted.bam ; do bcftools mpileup -d 1000000000  ${file} -f transcriptome_uniquetRNA.fa -a INFO/AD > ${file}.pileup; done\n",
    "for file in ${experiment}.splitFastq/*bam.sorted.bam ; do bcftools mpileup -d 1000000000  ${file} -f Pseudomonas_aeruginosa_transcriptome_tRNAunique.fa -a INFO/AD > ${file}.pileup; done\n",
    "\n",
    "if ${reference} == \"Ecoli\"\n",
    "then\n",
    "  for file in ${experiment}.splitFastq/*bam.sorted.bam ; do bcftools mpileup -d 1000000000  ${file} -f transcriptome_uniquetRNA.fa -a INFO/AD > ${file}.pileup; done\n",
    "fi\n",
    "\n",
    "if ${reference} == \"Paeruginosa\"\n",
    "then\n",
    "  for file in ${experiment}.splitFastq/*bam.sorted.bam ; do bcftools mpileup -d 1000000000  ${file} -f Pseudomonas_aeruginosa_transcriptome_tRNAunique.fa -a INFO/AD > ${file}.pileup; done\n",
    "fi\n",
    "\n",
    "# Keep only tRNA\n",
    "for file in ${experiment}.splitFastq/*.pileup ; do python3 parse_pileups.py ${file} ${file}.out ; done\n",
    "\n",
    "# Rename files\n",
    "for file in ${experiment}.splitFastq/*.pileup.out ; do mv \"${file}\" \"`echo ${file} | sed 's/.fastq.bam.mapped.bam.sorted.bam//'`\"; done\n",
    "\n",
    "###################################################################################################################\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ba9959d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#summarize_pileups.py\n",
    "#https://www.geeksforgeeks.org/how-to-read-multiple-text-files-from-folder-in-python/\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import pandas as pd\n",
    "\n",
    "# Folder Path \n",
    "path = sys.argv[1]\n",
    "\n",
    "# Change the directory \n",
    "os.chdir(path) \n",
    "\n",
    "all_rows = {}\n",
    "\n",
    "# Read text files\n",
    "def read_text_file(file_path): \n",
    "    with open(file_path, 'r') as f: \n",
    "\n",
    "        \n",
    "        rows = {}\n",
    "        \n",
    "        for line in f:\n",
    "            \n",
    "            array = line.strip().split('\\t')\n",
    "            ID = array[0]+'|'+array[1]\n",
    "            \n",
    "            DP = array[2]\n",
    "\n",
    "            variants = ''\n",
    "            \n",
    "            for field in array[3:]:\n",
    "                if field != '<*>':\n",
    "                    variants = variants+field+','\n",
    "                else:\n",
    "                    break\n",
    "\n",
    "            if ID not in rows:\n",
    "                rows[ID] = ''\n",
    "\n",
    "                \n",
    "            rows[ID] = rows[ID] + 'DP='+DP+','+ variants\n",
    "\n",
    "            \n",
    "    return(rows)\n",
    "\n",
    "                \n",
    "\n",
    "            \n",
    "# iterate through all files\n",
    "for file in os.listdir(): \n",
    "    # Check whether file is in text format or not \n",
    "    if file.endswith(\"pileup.out\"):\n",
    "    \n",
    "        column = file.replace(\".pileup.out\", \"\")\n",
    "\n",
    "        file_path = f\"{path}{file}\"\n",
    "  \n",
    "        # call read text file function \n",
    "        all_rows[column] =  read_text_file(file_path)\n",
    "\n",
    "        \n",
    "final_rows = {}\n",
    "\n",
    "columns = []\n",
    "for column in all_rows:\n",
    "    columns.append(column)\n",
    "    \n",
    "    for row in all_rows[column]:\n",
    "        if row not in final_rows:\n",
    "            final_rows[row] = []\n",
    "            \n",
    "        final_rows[row].append(all_rows[column][row])\n",
    "        \n",
    "df = pd.DataFrame.from_dict(final_rows, orient='index', columns=columns)\n",
    "df.to_csv('pileup_summary.csv')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e478fca",
   "metadata": {},
   "source": [
    "```\n",
    "\n",
    "mkdir summarize_pileups/\n",
    "cp ${experiment}.splitFastq/*.pileup.out summarize_pileups/\n",
    "cd summarize_pileups/\n",
    "python3 summarize_pileups.py ./\n",
    "\n",
    "\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "273622a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# misincorporation_heatmap_Ecoli.py\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "translate = {'EBT00001514652' : 'alaT', 'EBT00001514586' : 'alaU', 'EBT00001514664' : 'alaV',\\\n",
    "             'EBT00001514749' : 'alaW', 'EBT00001514687' : 'alaX', 'EBT00001514602' : 'argQ',\\\n",
    "             'EBT00001514635' : 'argU', 'EBT00001514579' : 'argV', 'EBT00001514747' : 'argW',\\\n",
    "             'EBT00001514671' : 'argX', 'EBT00001514668' : 'argY', 'EBT00001514736' : 'argZ',\\\n",
    "             'EBT00001514738' : 'asnT', 'EBT00001514720' : 'asnU', 'EBT00001514654' : 'asnV',\\\n",
    "             'EBT00001514578' : 'asnW', 'EBT00001514721' : 'aspT', 'EBT00001514692' : 'aspU',\\\n",
    "             'EBT00001514633' : 'aspV', 'EBT00001514638' : 'cysT', 'EBT00001514676' : 'glnU',\\\n",
    "             'EBT00001514740' : 'glnV', 'EBT00001514753' : 'glnW', 'EBT00001514655' : 'glnX',\\\n",
    "             'EBT00001514732' : 'gltT', 'EBT00001514606' : 'gltU', 'EBT00001514688' : 'gltV',\\\n",
    "             'EBT00001514601' : 'gltW', 'EBT00001514622' : 'glyT', 'EBT00001514729' : 'glyU',\\\n",
    "             'EBT00001514735' : 'glyV', 'EBT00001514716' : 'glyW', 'EBT00001514697' : 'glyX',\\\n",
    "             'EBT00001514650' : 'glyY', 'EBT00001514699' : 'hisR', 'EBT00001514691' : 'ileT',\\\n",
    "             'EBT00001514582' : 'ileU', 'EBT00001514702' : 'ileV', 'EBT00001514593' : 'ileX',\\\n",
    "             'EBT00001514626' : 'ileY', 'EBT00001514730' : 'leuP', 'EBT00001514592' : 'leuQ',\\\n",
    "             'EBT00001514696' : 'leuT', 'EBT00001514711' : 'leuU', 'EBT00001514647' : 'leuV',\\\n",
    "             'EBT00001514636' : 'leuW', 'EBT00001514709' : 'leuX', 'EBT00001514679' : 'leuZ',\\\n",
    "             'EBT00001514624' : 'lysQ', 'EBT00001514739' : 'lysT', 'EBT00001514751' : 'lysV',\\\n",
    "             'EBT00001514722' : 'lysW', 'EBT00001514600' : 'lysY', 'EBT00001514728' : 'lysZ',\\\n",
    "             'EBT00001514597' : 'metT', 'EBT00001514667' : 'metU', 'EBT00001514651' : 'metV',\\\n",
    "             'EBT00001514715' : 'metW', 'EBT00001514737' : 'metY', 'EBT00001514613' : 'metZ',\\\n",
    "             'EBT00001514723' : 'pheU', 'EBT00001514693' : 'pheV', 'EBT00001514637' : 'proK',\\\n",
    "             'EBT00001514632' : 'proL', 'EBT00001514733' : 'proM', 'EBT00001514743' : 'selC',\\\n",
    "             'EBT00001514689' : 'serT', 'EBT00001514584' : 'serU', 'EBT00001514690' : 'serV',\\\n",
    "             'EBT00001514678' : 'serW', 'EBT00001514596' : 'serX', 'EBT00001514707' : 'thrT',\\\n",
    "             'EBT00001514627' : 'thrU', 'EBT00001514591' : 'thrV', 'EBT00001514734' : 'thrW',\\\n",
    "             'EBT00001514641' : 'trpT', 'EBT00001514612' : 'tyrT', 'EBT00001514598' : 'tyrU',\\\n",
    "             'EBT00001514684' : 'tyrV', 'EBT00001514719' : 'valT', 'EBT00001514665' : 'valU',\\\n",
    "             'EBT00001514621' : 'valV', 'EBT00001514583' : 'valW', 'EBT00001514604' : 'valX',\\\n",
    "             'EBT00001514683' : 'valY', 'EBT00001514642' : 'valZ'}\n",
    "\n",
    "\n",
    "infile = pd.read_csv('pileup_summary.csv')\n",
    "\n",
    "\n",
    "ID_POS = infile.iloc[:, 0]\n",
    "\n",
    "experiment = sys.argv[1]\n",
    "\n",
    "Rep1 = infile[experiment+'_rep1']\n",
    "Rep2 = infile[experiment+'_rep2']\n",
    "Rep3 = infile[experiment+'_rep3']\n",
    "\n",
    "dic_depths = {}\n",
    "dic = {}\n",
    "\n",
    "for i in range(0, len(ID_POS)):\n",
    "    \n",
    "    ID, POS = ID_POS[i].split('|')\n",
    "    \n",
    "    DP_depths = []\n",
    "    DP_refs = []\n",
    "    \n",
    "    for Rep in [Rep1, Rep2, Rep3]:\n",
    "        \n",
    "        if 'nan' in str(Rep[i]):\n",
    "            \n",
    "            DP_depths.append(np.nan)\n",
    "            DP_refs.append(np.nan)\n",
    "                \n",
    "        else:\n",
    "            DP = Rep[i].split(',')[0].split('=')[1]\n",
    "            \n",
    "            DP_ref = Rep[i].split(',')[2]\n",
    "\n",
    "            if int(DP) > 1:\n",
    "                DP_refs.append(float(DP_ref))\n",
    "                DP_depths.append(int(DP))\n",
    "            else:\n",
    "                DP_refs.append(np.nan)  \n",
    "                DP_depths.append(np.nan)\n",
    "    \n",
    "    DP_ref_mean = np.nanmean(DP_refs)\n",
    "    \n",
    "    ID = translate[ID]\n",
    "    \n",
    "    if ID not in dic:\n",
    "        dic_depths[ID] = []\n",
    "        dic[ID] = []\n",
    "        \n",
    "    dic_depths[ID].append(DP_depths)\n",
    "    dic[ID].append(DP_ref_mean)\n",
    "    \n",
    "\n",
    "df = pd.DataFrame.from_dict(dic, orient='index')\n",
    "df.to_csv(experiment+'_HeatMap.csv')\n",
    "\n",
    "df_depths = pd.DataFrame.from_dict(dic_depths, orient='index')\n",
    "df_depths.to_csv(experiment+'_HeatMap_depths.csv')\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48b8da13",
   "metadata": {},
   "source": [
    "```\n",
    "for experiment in M63p06_captured M63p06_input M63p2_captured M63p2_input M63m06_captured M63m06_input LB06_captured LB06_input LB2_captured LB2_input Spent0perc_captured Spent0perc_input Spent25perc_captured Spent25perc_input Spent50perc_captured Spent50perc_input dThiI06_input dThiI2_input ; do python3 misincorporation_heatmap_Ecoli.py ${experiment} ; done\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c08929b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# misincorporation_heatmap_PA.py\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "\n",
    "infile = pd.read_csv('PA_pileup_summary.csv')\n",
    "\n",
    "\n",
    "ID_POS = infile.iloc[:, 0]\n",
    "\n",
    "experiment = sys.argv[1]\n",
    "\n",
    "Rep1 = infile[experiment+'_rep1']\n",
    "Rep2 = infile[experiment+'_rep2']\n",
    "Rep3 = infile[experiment+'_rep3']\n",
    "\n",
    "dic_depths = {}\n",
    "dic = {}\n",
    "\n",
    "for i in range(0, len(ID_POS)):\n",
    "    \n",
    "    ID, POS = ID_POS[i].split('|')\n",
    "    \n",
    "    DP_depths = []\n",
    "    DP_refs = []\n",
    "    \n",
    "    for Rep in [Rep1, Rep2, Rep3]:\n",
    "        \n",
    "        if 'nan' in str(Rep[i]):\n",
    "            \n",
    "            DP_depths.append(np.nan)\n",
    "            DP_refs.append(np.nan)\n",
    "                \n",
    "        else:\n",
    "            DP = Rep[i].split(',')[0].split('=')[1]\n",
    "            \n",
    "            DP_ref = Rep[i].split(',')[2]\n",
    "\n",
    "            if int(DP) > 1:\n",
    "                DP_refs.append(float(DP_ref))\n",
    "                DP_depths.append(int(DP))\n",
    "            else:\n",
    "                DP_refs.append(np.nan)  \n",
    "                DP_depths.append(np.nan)\n",
    "    \n",
    "    DP_ref_mean = np.nanmean(DP_refs)\n",
    "        \n",
    "    if ID not in dic:\n",
    "        dic_depths[ID] = []\n",
    "        dic[ID] = []\n",
    "        \n",
    "    dic_depths[ID].append(DP_depths)\n",
    "    dic[ID].append(DP_ref_mean)\n",
    "    \n",
    "\n",
    "df = pd.DataFrame.from_dict(dic, orient='index')\n",
    "df.to_csv(experiment+'_HeatMap.csv')\n",
    "\n",
    "df_depths = pd.DataFrame.from_dict(dic_depths, orient='index')\n",
    "df_depths.to_csv(experiment+'_HeatMap_depths.csv')\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e69bb82",
   "metadata": {},
   "source": [
    "```\n",
    "\n",
    "for experiment in PA_capture PA_input ; do python3 misincorporation_heatmap_PA.py ${experiment} ; done\n",
    "\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "980a8843",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c88d6cb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "\n",
    "##########################################################################################\n",
    "# https://www.r-bloggers.com/2015/09/passing-arguments-to-an-r-script-from-command-lines/\n",
    "# heatmap.R\n",
    "\n",
    "library(ggplot2)\n",
    "library(heatmaply)\n",
    "library(plotly)\n",
    "library(orca)\n",
    "\n",
    "\n",
    "args = commandArgs(trailingOnly=TRUE)\n",
    "\n",
    "\n",
    "# test if there is at least one argument: if not, return an error\n",
    "if (length(args)==0) {\n",
    "  stop(\"At least one argument must be supplied (input file).n\", call.=FALSE)\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "data <- read.csv(args[1], header=TRUE)\n",
    "data <- cbind(data[,1], data[,2:76])\n",
    "colnames(data) <- c('ID', seq(1,length(data[1,])-1))\n",
    "\n",
    "toplot <- data[-1]\n",
    "rownames(toplot) <- data$ID\n",
    "toplot <- toplot[order(rownames(toplot)), ]\n",
    "\n",
    "\n",
    "p <- heatmaply(toplot, \n",
    "        dendrogram = \"none\",\n",
    "        file = paste0(args[1],'.pdf'),\n",
    "        #plot_method = \"plotly\",\n",
    "        xlab = \"\", ylab = \"\", \n",
    "        main = args[1],\n",
    "        #scale = \"row\",\n",
    "        margins = c(60,100,40,20),\n",
    "        grid_color = \"white\",\n",
    "        grid_width = 0.00001,\n",
    "        titleX = FALSE,\n",
    "        hide_colorbar = TRUE,\n",
    "        branches_lwd = 0.1,\n",
    "        fontsize_row = 5, fontsize_col = 5,\n",
    "        labRow = rownames(toplot),\n",
    "        column_text_angle = 90,\n",
    "        )\n",
    "#############################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d88f3d32",
   "metadata": {},
   "source": [
    "```\n",
    "conda create -n heatmaps python=3.6\n",
    "conda activate heatmaps\n",
    "conda config --env --add channels bioconda\n",
    "conda config --env --add channels conda-forge\n",
    "conda config --env --set channel_priority strict\n",
    "conda config --show-sources\n",
    "\n",
    "conda install -c conda-forge r-base\n",
    "conda install -c conda-forge r-ggplot2\n",
    "conda install -c conda-forge r-heatmaply\n",
    "conda install -c plotly plotly\n",
    "\n",
    "R\n",
    "install.packages(\"orca\")\n",
    "\n",
    "for experiment in M63p06_captured M63p06_input M63p2_captured M63p2_input M63m06_captured M63m06_input LB06_captured LB06_input LB2_captured LB2_input Spent0perc_captured Spent0perc_input Spent25perc_captured Spent25perc_input Spent50perc_captured Spent50perc_input dThiI06_input dThiI2_input ; do Rscript heatmap.R ${experiment}_HeatMap.csv ; done\n",
    "\n",
    "\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d77eaf7c",
   "metadata": {},
   "source": [
    "Example of heatmap <img src=\"M63p06_captured_HeatMap.csv.pdf\" width=1200 height=800 />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e41ea51c",
   "metadata": {},
   "source": [
    "```\n",
    "for experiment in PA_capture PA_input ; do Rscript heatmap.R ${experiment}_HeatMap.csv ; done\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5250c50",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56aedcae",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7e4b3649",
   "metadata": {},
   "source": [
    "# IV. Termination HeatMaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "593d36e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Following script saved as termination_analysis.sh\n",
    "# Command: sbatch termination_analysis.sh experiment\n",
    "\n",
    "###################################################################################################################\n",
    "#! /bin/bash\n",
    "\n",
    "#SBATCH --job-name wil_works\n",
    "#SBATCH --ntasks=4\n",
    "#SBATCH --cpus-per-task=4\n",
    "#SBATCH --partition=ccr\n",
    "#SBATCH --time=24:00:00\n",
    "#SBATCH --gres=lscratch:500 \n",
    "#SBATCH --mem=50g\n",
    "\n",
    "experiment=$1\n",
    "reference=$2\n",
    "\n",
    "module load bedtools\n",
    "\n",
    "if ${reference} == \"Ecoli\"\n",
    "then\n",
    "  for file in ${experiment}.splitFastq/*bam.sorted.bam ; do bedtools genomecov -d -5 -ibam  ${file} -g transcriptome_uniquetRNA.fa > ${file}_5end.txt ; bedtools genomecov -d -3 -ibam ${file} -g transcriptome_uniquetRNA.fa > ${file}_3end.txt ; done\n",
    "fi\n",
    "\n",
    "if ${reference} == \"Paeruginosa\"\n",
    "then\n",
    "  for file in ${experiment}.splitFastq/*bam.sorted.bam ; do bedtools genomecov -d -5 -ibam  ${file} -g Pseudomonas_aeruginosa_transcriptome_tRNAunique.fa > ${file}_5end.txt ; bedtools genomecov -d -3 -ibam ${file} -g transcriptome_uniquetRNA.fa > ${file}_3end.txt ; done\n",
    "fi\n",
    "\n",
    "\n",
    "###################################################################################################################\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "956d6244",
   "metadata": {},
   "outputs": [],
   "source": [
    "#summarize_terminations.py\n",
    "#https://www.geeksforgeeks.org/how-to-read-multiple-text-files-from-folder-in-python/\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "tRNAs = ['EBT00001514652','EBT00001514586','EBT00001514664','EBT00001514749','EBT00001514687','EBT00001514602',\\\n",
    "         'EBT00001514635','EBT00001514579','EBT00001514747','EBT00001514671','EBT00001514668','EBT00001514736',\\\n",
    "         'EBT00001514738','EBT00001514720','EBT00001514654','EBT00001514578','EBT00001514721','EBT00001514692',\\\n",
    "         'EBT00001514633','EBT00001514638','EBT00001514676','EBT00001514740','EBT00001514753','EBT00001514655',\\\n",
    "         'EBT00001514732','EBT00001514606','EBT00001514688','EBT00001514601','EBT00001514622','EBT00001514729',\\\n",
    "         'EBT00001514735','EBT00001514716','EBT00001514697','EBT00001514650','EBT00001514699','EBT00001514691',\\\n",
    "         'EBT00001514582','EBT00001514702','EBT00001514593','EBT00001514626','EBT00001514730','EBT00001514592',\\\n",
    "         'EBT00001514696','EBT00001514711','EBT00001514647','EBT00001514636','EBT00001514709','EBT00001514679',\\\n",
    "         'EBT00001514624','EBT00001514739','EBT00001514751','EBT00001514722','EBT00001514600','EBT00001514728',\\\n",
    "         'EBT00001514597','EBT00001514667','EBT00001514651','EBT00001514715','EBT00001514737','EBT00001514613',\\\n",
    "         'EBT00001514723','EBT00001514693','EBT00001514637','EBT00001514632','EBT00001514733','EBT00001514743',\\\n",
    "         'EBT00001514689','EBT00001514584','EBT00001514690','EBT00001514678','EBT00001514596','EBT00001514707',\\\n",
    "         'EBT00001514627','EBT00001514591','EBT00001514734','EBT00001514641','EBT00001514612','EBT00001514598',\\\n",
    "         'EBT00001514684','EBT00001514719','EBT00001514665','EBT00001514621','EBT00001514583','EBT00001514604',\\\n",
    "         'EBT00001514683','EBT00001514642', # E coli\n",
    "\n",
    "         'EBT00051365344','EBT00051365099','EBT00051370091','EBT00051369932','EBT00051369986','EBT00051367609',\\\n",
    "         'EBT00051369609','EBT00051366328','EBT00051367914','EBT00051366603','EBT00051366399','EBT00051365908',\\\n",
    "         'EBT00051369492','EBT00051366465','EBT00051366923','EBT00051368550','EBT00051368683','EBT00051365503',\\\n",
    "         'EBT00051369841','EBT00051365018','EBT00051367353','EBT00051365422','EBT00051365808','EBT00051368205',\\\n",
    "         'EBT00051369258','EBT00051365670','EBT00051364934','EBT00051365187','EBT00051369657','EBT00051366843',\\\n",
    "         'EBT00051366536','EBT00051368279','EBT00051367683','EBT00051368081','EBT00051367446','EBT00051370188',\\\n",
    "         'EBT00051369313','EBT00051369436','EBT00051370040','EBT00051370146','EBT00051367983','EBT00051367160',\\\n",
    "         'EBT00051366161'] # P a\n",
    "\n",
    "\n",
    "\n",
    "# Folder Path \n",
    "path = sys.argv[1]\n",
    "\n",
    "# Change the directory \n",
    "os.chdir(path) \n",
    "\n",
    "\n",
    "\n",
    "all_rows = {}\n",
    "\n",
    "# Read text File\n",
    "def read_text_file(file_path): \n",
    "    with open(file_path, 'r') as f: \n",
    "        \n",
    "        rows = {}\n",
    "        \n",
    "        for line in f:\n",
    "            \n",
    "            array = line.strip().split('\\t')\n",
    "            ID = array[0]+'|'+array[1]\n",
    "            \n",
    "            if array[0] in tRNAs:\n",
    "                Term_count = array[2]\n",
    "                \n",
    "                if ID not in rows:\n",
    "                    rows[ID] = ''\n",
    "\n",
    "                rows[ID] = Term_count\n",
    "\n",
    "    return(rows)\n",
    "\n",
    "                \n",
    "\n",
    "            \n",
    "# iterate through all file \n",
    "for file in os.listdir(): \n",
    "    # Check whether file is in text format or not \n",
    "    if file.endswith(\"bam_5end.txt\"):\n",
    "    \n",
    "        column = file.replace(\".fastq.bam.mapped.bam.sorted.bam_5end.txt\", \"\")\n",
    "\n",
    "        file_path = f\"{path}{file}\"\n",
    "  \n",
    "        # call read text file function \n",
    "        all_rows[column] =  read_text_file(file_path)\n",
    "\n",
    "        \n",
    "final_rows = {}\n",
    "\n",
    "columns = []\n",
    "for column in all_rows:\n",
    "    columns.append(column)\n",
    "    \n",
    "    for row in all_rows[column]:\n",
    "        if row not in final_rows:\n",
    "            final_rows[row] = []\n",
    "            \n",
    "        final_rows[row].append(all_rows[column][row])\n",
    "        \n",
    "df = pd.DataFrame.from_dict(final_rows, orient='index', columns=columns)\n",
    "df.to_csv('termination_summary.csv')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a45c61f3",
   "metadata": {},
   "source": [
    "```\n",
    "mkdir summarize_terminations/\n",
    "cp ${experiment}.splitFastq/*bam_5end.txt summarize_terminations/\n",
    "cd summarize_terminations/\n",
    "python3 summarize_terminations.py ./\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "815f333d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# termination_heatmap_Ecoli.py\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "translate = {'EBT00001514652' : 'alaT', 'EBT00001514586' : 'alaU', 'EBT00001514664' : 'alaV',\\\n",
    "             'EBT00001514749' : 'alaW', 'EBT00001514687' : 'alaX', 'EBT00001514602' : 'argQ',\\\n",
    "             'EBT00001514635' : 'argU', 'EBT00001514579' : 'argV', 'EBT00001514747' : 'argW',\\\n",
    "             'EBT00001514671' : 'argX', 'EBT00001514668' : 'argY', 'EBT00001514736' : 'argZ',\\\n",
    "             'EBT00001514738' : 'asnT', 'EBT00001514720' : 'asnU', 'EBT00001514654' : 'asnV',\\\n",
    "             'EBT00001514578' : 'asnW', 'EBT00001514721' : 'aspT', 'EBT00001514692' : 'aspU',\\\n",
    "             'EBT00001514633' : 'aspV', 'EBT00001514638' : 'cysT', 'EBT00001514676' : 'glnU',\\\n",
    "             'EBT00001514740' : 'glnV', 'EBT00001514753' : 'glnW', 'EBT00001514655' : 'glnX',\\\n",
    "             'EBT00001514732' : 'gltT', 'EBT00001514606' : 'gltU', 'EBT00001514688' : 'gltV',\\\n",
    "             'EBT00001514601' : 'gltW', 'EBT00001514622' : 'glyT', 'EBT00001514729' : 'glyU',\\\n",
    "             'EBT00001514735' : 'glyV', 'EBT00001514716' : 'glyW', 'EBT00001514697' : 'glyX',\\\n",
    "             'EBT00001514650' : 'glyY', 'EBT00001514699' : 'hisR', 'EBT00001514691' : 'ileT',\\\n",
    "             'EBT00001514582' : 'ileU', 'EBT00001514702' : 'ileV', 'EBT00001514593' : 'ileX',\\\n",
    "             'EBT00001514626' : 'ileY', 'EBT00001514730' : 'leuP', 'EBT00001514592' : 'leuQ',\\\n",
    "             'EBT00001514696' : 'leuT', 'EBT00001514711' : 'leuU', 'EBT00001514647' : 'leuV',\\\n",
    "             'EBT00001514636' : 'leuW', 'EBT00001514709' : 'leuX', 'EBT00001514679' : 'leuZ',\\\n",
    "             'EBT00001514624' : 'lysQ', 'EBT00001514739' : 'lysT', 'EBT00001514751' : 'lysV',\\\n",
    "             'EBT00001514722' : 'lysW', 'EBT00001514600' : 'lysY', 'EBT00001514728' : 'lysZ',\\\n",
    "             'EBT00001514597' : 'metT', 'EBT00001514667' : 'metU', 'EBT00001514651' : 'metV',\\\n",
    "             'EBT00001514715' : 'metW', 'EBT00001514737' : 'metY', 'EBT00001514613' : 'metZ',\\\n",
    "             'EBT00001514723' : 'pheU', 'EBT00001514693' : 'pheV', 'EBT00001514637' : 'proK',\\\n",
    "             'EBT00001514632' : 'proL', 'EBT00001514733' : 'proM', 'EBT00001514743' : 'selC',\\\n",
    "             'EBT00001514689' : 'serT', 'EBT00001514584' : 'serU', 'EBT00001514690' : 'serV',\\\n",
    "             'EBT00001514678' : 'serW', 'EBT00001514596' : 'serX', 'EBT00001514707' : 'thrT',\\\n",
    "             'EBT00001514627' : 'thrU', 'EBT00001514591' : 'thrV', 'EBT00001514734' : 'thrW',\\\n",
    "             'EBT00001514641' : 'trpT', 'EBT00001514612' : 'tyrT', 'EBT00001514598' : 'tyrU',\\\n",
    "             'EBT00001514684' : 'tyrV', 'EBT00001514719' : 'valT', 'EBT00001514665' : 'valU',\\\n",
    "             'EBT00001514621' : 'valV', 'EBT00001514583' : 'valW', 'EBT00001514604' : 'valX',\\\n",
    "             'EBT00001514683' : 'valY', 'EBT00001514642' : 'valZ'}\n",
    "\n",
    "\n",
    "infile = pd.read_csv('termination_summary.csv')\n",
    "\n",
    "\n",
    "ID_POS = infile.iloc[:, 0]\n",
    "\n",
    "experiment = sys.argv[1]\n",
    "\n",
    "Rep1 = infile[experiment+'_rep1']\n",
    "Rep2 = infile[experiment+'_rep2']\n",
    "Rep3 = infile[experiment+'_rep3']\n",
    "\n",
    "\n",
    "dic = {}\n",
    "\n",
    "for i in range(0, len(ID_POS)):\n",
    "    \n",
    "    ID, POS = ID_POS[i].split('|')\n",
    "    \n",
    "    Terminations = []\n",
    "\n",
    "    \n",
    "    for Rep in [Rep1, Rep2, Rep3]:\n",
    "        \n",
    "        if 'nan' in str(Rep[i]):\n",
    "            \n",
    "            Terminations.append(np.nan)\n",
    "                \n",
    "        else:\n",
    "            Terminations.append(Rep[i])\n",
    "\n",
    "    Terminations_mean = np.nanmean(Terminations)\n",
    "    \n",
    "    ID = translate[ID]\n",
    "    \n",
    "    if ID not in dic:\n",
    "        dic[ID] = []\n",
    "        \n",
    "    dic[ID].append(Terminations_mean)\n",
    "    \n",
    "    \n",
    "    \n",
    "dic_normalized = {}\n",
    "\n",
    "for ID in dic:\n",
    "    dic_normalized[ID] = []\n",
    "    \n",
    "    for i in range(0, len(dic[ID])):\n",
    "\n",
    "        dic_normalized[ID].append(1- float(dic[ID][i])/ sum(dic[ID][:i+1]))\n",
    "    \n",
    "    \n",
    "\n",
    "df = pd.DataFrame.from_dict(dic_normalized, orient='index')\n",
    "\n",
    "\n",
    "df.to_csv(experiment+'_termination_HeatMap.csv')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3bc5882",
   "metadata": {},
   "source": [
    "```\n",
    "for experiment in M63p06_captured M63p06_input M63p2_captured M63p2_input M63m06_captured M63m06_input LB06_captured LB06_input LB2_captured LB2_input Spent0perc_captured Spent0perc_input Spent25perc_captured Spent25perc_input Spent50perc_captured Spent50perc_input dThiI06_input dThiI2_input ; do python3 termination_heatmap_Ecoli.py ${experiment} ; done\n",
    "\n",
    "for experiment in M63p06_captured M63p06_input M63p2_captured M63p2_input M63m06_captured M63m06_input LB06_captured LB06_input LB2_captured LB2_input Spent0perc_captured Spent0perc_input Spent25perc_captured Spent25perc_input Spent50perc_captured Spent50perc_input dThiI06_input dThiI2_input ; do Rscript heatmap.R ${experiment}_termination_HeatMap.csv ; done\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f932661d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# termination_heatmap_PA.py\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "\n",
    "\n",
    "infile = pd.read_csv('PA_termination_summary.csv')\n",
    "\n",
    "\n",
    "ID_POS = infile.iloc[:, 0]\n",
    "\n",
    "experiment = sys.argv[1]\n",
    "\n",
    "Rep1 = infile[experiment+'_rep1']\n",
    "Rep2 = infile[experiment+'_rep2']\n",
    "Rep3 = infile[experiment+'_rep3']\n",
    "\n",
    "\n",
    "dic = {}\n",
    "\n",
    "for i in range(0, len(ID_POS)):\n",
    "    \n",
    "    ID, POS = ID_POS[i].split('|')\n",
    "    \n",
    "    Terminations = []\n",
    "\n",
    "    \n",
    "    for Rep in [Rep1, Rep2, Rep3]:\n",
    "        \n",
    "        if 'nan' in str(Rep[i]):\n",
    "            \n",
    "            Terminations.append(np.nan)\n",
    "                \n",
    "        else:\n",
    "            Terminations.append(Rep[i])\n",
    "\n",
    "    Terminations_mean = np.nanmean(Terminations)\n",
    "        \n",
    "    if ID not in dic:\n",
    "        dic[ID] = []\n",
    "        \n",
    "    dic[ID].append(Terminations_mean)\n",
    "    \n",
    "    \n",
    "    \n",
    "dic_normalized = {}\n",
    "\n",
    "for ID in dic:\n",
    "    dic_normalized[ID] = []\n",
    "    \n",
    "    for i in range(0, len(dic[ID])):\n",
    "\n",
    "        dic_normalized[ID].append(1- float(dic[ID][i])/ sum(dic[ID][:i+1]))\n",
    "    \n",
    "    \n",
    "\n",
    "df = pd.DataFrame.from_dict(dic_normalized, orient='index')\n",
    "\n",
    "\n",
    "df.to_csv(experiment+'_termination_HeatMap.csv')\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48c9c1ba",
   "metadata": {},
   "source": [
    "```\n",
    "\n",
    "for experiment in PA_capture PA_input ; do python3 termination_heatmap_PA.py ${experiment} ; done\n",
    "\n",
    "for experiment in PA_capture PA_input ; do Rscript heatmap.R ${experiment}_termination_HeatMap.csv ; done\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "968d5b06",
   "metadata": {},
   "source": [
    "# DeSeq2 Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c44afffc",
   "metadata": {},
   "source": [
    "## Keep only tRNAs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "443c162d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#parse_tRNA.py\n",
    "\n",
    "import sys\n",
    "import re\n",
    "\n",
    "infile = open(sys.argv[1], 'rt')\n",
    "outfile = open(sys.argv[2], 'w+')\n",
    "\n",
    "\n",
    "tRNAs = ['EBT00001514652','EBT00001514586','EBT00001514664','EBT00001514749','EBT00001514687','EBT00001514602',\\\n",
    "         'EBT00001514635','EBT00001514579','EBT00001514747','EBT00001514671','EBT00001514668','EBT00001514736',\\\n",
    "         'EBT00001514738','EBT00001514720','EBT00001514654','EBT00001514578','EBT00001514721','EBT00001514692',\\\n",
    "         'EBT00001514633','EBT00001514638','EBT00001514676','EBT00001514740','EBT00001514753','EBT00001514655',\\\n",
    "         'EBT00001514732','EBT00001514606','EBT00001514688','EBT00001514601','EBT00001514622','EBT00001514729',\\\n",
    "         'EBT00001514735','EBT00001514716','EBT00001514697','EBT00001514650','EBT00001514699','EBT00001514691',\\\n",
    "         'EBT00001514582','EBT00001514702','EBT00001514593','EBT00001514626','EBT00001514730','EBT00001514592',\\\n",
    "         'EBT00001514696','EBT00001514711','EBT00001514647','EBT00001514636','EBT00001514709','EBT00001514679',\\\n",
    "         'EBT00001514624','EBT00001514739','EBT00001514751','EBT00001514722','EBT00001514600','EBT00001514728',\\\n",
    "         'EBT00001514597','EBT00001514667','EBT00001514651','EBT00001514715','EBT00001514737','EBT00001514613',\\\n",
    "         'EBT00001514723','EBT00001514693','EBT00001514637','EBT00001514632','EBT00001514733','EBT00001514743',\\\n",
    "         'EBT00001514689','EBT00001514584','EBT00001514690','EBT00001514678','EBT00001514596','EBT00001514707',\\\n",
    "         'EBT00001514627','EBT00001514591','EBT00001514734','EBT00001514641','EBT00001514612','EBT00001514598',\\\n",
    "         'EBT00001514684','EBT00001514719','EBT00001514665','EBT00001514621','EBT00001514583','EBT00001514604',\\\n",
    "         'EBT00001514683','EBT00001514642', # E coli\n",
    "\n",
    "         'EBT00051365344','EBT00051365099','EBT00051370091','EBT00051369932','EBT00051369986','EBT00051367609',\\\n",
    "         'EBT00051369609','EBT00051366328','EBT00051367914','EBT00051366603','EBT00051366399','EBT00051365908',\\\n",
    "         'EBT00051369492','EBT00051366465','EBT00051366923','EBT00051368550','EBT00051368683','EBT00051365503',\\\n",
    "         'EBT00051369841','EBT00051365018','EBT00051367353','EBT00051365422','EBT00051365808','EBT00051368205',\\\n",
    "         'EBT00051369258','EBT00051365670','EBT00051364934','EBT00051365187','EBT00051369657','EBT00051366843',\\\n",
    "         'EBT00051366536','EBT00051368279','EBT00051367683','EBT00051368081','EBT00051367446','EBT00051370188',\\\n",
    "         'EBT00051369313','EBT00051369436','EBT00051370040','EBT00051370146','EBT00051367983','EBT00051367160',\\\n",
    "         'EBT00051366161'] # P a\n",
    "\n",
    "\n",
    "\n",
    "for line in infile:\n",
    "    if line[0] != '#':\n",
    "        array = line.strip().split('\\t')\n",
    "        ID, count = array\n",
    "        \n",
    "        if ID in tRNAs:\n",
    "\n",
    "            outfile.write(ID+'\\t'+count+'\\n')\n",
    "\n",
    "            \n",
    "            \n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09abc978",
   "metadata": {},
   "source": [
    "```\n",
    "for file in *count; do python parse_tRNA.py ${file} ${file}.tRNA ; done\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dde6c9e6",
   "metadata": {},
   "source": [
    "## Single-Factor analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42f902b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "#SingleFactorDeSeq2_tRNA.r\n",
    "\n",
    "library(dplyr)\n",
    "\n",
    "args = commandArgs(trailingOnly=TRUE)\n",
    "path <- args[1]\n",
    "current <- args[2]\n",
    "other <- args[3]\n",
    "\n",
    "################################################\n",
    "# PART 1: GATHER ALL INPUTS INTO ONE DATAFRAME #\n",
    "###############################################\n",
    "\n",
    "# LIST ALL INPUT FILES\n",
    "temp <- list.files(path = path, pattern=\"*.count.tRNA\")\n",
    "\n",
    "# START DATAFRAME WITH FIRST INPUT\n",
    "df <- as.data.frame(read.csv(paste0(path,temp[1]), header=FALSE, sep='\\t'))\n",
    "cols <- c('ID',substr(temp[1], 1, nchar(temp[1])-11))\n",
    "condition <- substr(temp[1], 1, nchar(temp[1])-16)\n",
    "\n",
    "# FILL DATAFRAME WITH THE REMAINING INPUTS\n",
    "for (file in temp[2:length(temp)]) {\n",
    "    cols <- c(cols,substr(file, 1, nchar(file)-11))\n",
    "    condition <- c(condition,substr(file, 1, nchar(file)-16))\n",
    "    df2 <- as.data.frame(read.csv(paste0(path,file), header=FALSE, sep='\\t'))\n",
    "    df <- full_join(df, df2, by='V1')\n",
    "}\n",
    "\n",
    "condition <- factor(condition)\n",
    "warnings()\n",
    "df[is.na(df)] <- 0\n",
    "colnames(df) <- cols\n",
    "write.table(df,\"s4U_counts_table.tRNA.txt\", quote = FALSE, row.names = FALSE)\n",
    "\n",
    "\n",
    "##################################\n",
    "# PART 2: SINGLE FACTOR ANALYSIS #\n",
    "##################################\n",
    "\n",
    "# READ DATA INTO MATRIX\n",
    "data <- read.table(\"s4U_counts_table.tRNA.txt\", header=T, row.names=1)\n",
    "countdata <- as.matrix(data)\n",
    "\n",
    "# MAKE MATRIX A DESEQ2 OBJECT\n",
    "library(DESeq2)\n",
    "coldata <- data.frame(row.names=colnames(countdata), condition)\n",
    "print(coldata)\n",
    "dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)\n",
    "print(summary(dds))\n",
    "dds\n",
    "\n",
    "#filter the dataset - remove all zeros, influences padj\n",
    "nrow(dds)\n",
    "dds <- dds[ rowSums(counts(dds)) > 1, ]\n",
    "nrow(dds)\n",
    "\n",
    "# Run the DESeq pipeline\n",
    "dds <- DESeq(dds)\n",
    "res <- results(dds)\n",
    "print(res)\n",
    "\n",
    "\n",
    "conditions <- c(paste0(current,'_captured'), paste0(current,'_input'), paste0(other,'_captured'), paste0(other,'_input'))\n",
    "\n",
    "for (current in conditions) { conditions <- conditions[conditions != current]\n",
    "                             \n",
    "        for (other in conditions) {\n",
    "            \n",
    "            print(current)\n",
    "            print(other)\n",
    "    \n",
    "            current_res = results(dds, contrast=c(\"condition\",current,other))\n",
    "            current_res\n",
    "            print(current_res)\n",
    "            summary(current_res)\n",
    "            current_ress_Orderedlfc <- current_res[order(current_res$log2FoldChange),]\n",
    "            head(current_res_Orderedlfc)\n",
    "\n",
    "            # Exploring and exporting results\n",
    "\n",
    "            pdf(file = paste0(current,'_VS_',other,'.pdf'),   # The directory you want to save the file in\n",
    "            width = 10, # The width of the plot in inches\n",
    "            height = 10) # The height of the plot in inches\n",
    "\n",
    "            plotMA(current_res, ylim=c(-10,10))\n",
    "            dev.off()\n",
    "\n",
    "            identify(current_res$baseMean, current_res$log2FoldChange, labels = row.names(current_res))\n",
    "            write.csv(as.data.frame(current_res), file= paste0(current,'_VS_',other,'.tRNA.csv'))\n",
    "\n",
    "            # Volcano plot\n",
    "            pdf(file = paste0(current,'_VS_',other,'.volcano.tRNA.pdf'),   # The directory you want to save the file in\n",
    "                width = 10, # The width of the plot in inches\n",
    "                height = 10) # The height of the plot in inches\n",
    "            par(mar=c(5,5,5,5), cex=1.0, cex.main=1.4, cex.axis=1.4, cex.lab=1.4)\n",
    "            topT <- as.data.frame(current_res)\n",
    "            #Adjusted P values (FDR Q values)\n",
    "            with(topT, plot(log2FoldChange, -log10(padj), pch=20, main=paste0(current,' VS ',other,'.volcano'), cex=1.0, xlab=bquote(~Log[2]~fold~change), ylab=bquote(~-log[10]~Q~value)))\n",
    "            with(subset(topT, padj<0.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(padj), pch=20, col=\"red\", cex=0.5))\n",
    "            #Add lines for absolute FC>2 and P-value cut-off at FDR Q<0.05\n",
    "            abline(v=0, col=\"black\", lty=3, lwd=1.0)\n",
    "            abline(v=-1, col=\"black\", lty=4, lwd=2.0)\n",
    "            abline(v=1, col=\"black\", lty=4, lwd=2.0)\n",
    "            abline(h=-log10(max(topT$pvalue[topT$padj<0.05], na.rm=TRUE)), col=\"black\", lty=4, lwd=2.0)\n",
    "            dev.off()\n",
    "        }\n",
    "    }\n",
    "            \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d3c00a4",
   "metadata": {},
   "source": [
    "```\n",
    "\n",
    "Rscript SingleFactorDeSeq2_tRNA.r LBvsM63analysis_tRNA/ LB06 LB2\n",
    "Rscript SingleFactorDeSeq2_tRNA.r LBvsM63analysis_tRNA/ M63m06 M63p06\n",
    "Rscript SingleFactorDeSeq2_tRNA.r LBvsM63analysis_tRNA/ M63p06 M63p2\n",
    "Rscript SingleFactorDeSeq2_tRNA.r LBvsM63analysis_tRNA/ LB06 M63p06\n",
    "Rscript SingleFactorDeSeq2_tRNA.r LBvsM63analysis_tRNA/ LB2 M63p2\n",
    "\n",
    "\n",
    "\n",
    "Rscript SingleFactorDeSeq2_tRNA.r Spent_analysis_tRNA/ Spent0perc Spent25perc\n",
    "Rscript SingleFactorDeSeq2_tRNA.r Spent_analysis_tRNA/ Spent0perc Spent50perc\n",
    "Rscript SingleFactorDeSeq2_tRNA.r Spent_analysis_tRNA/ Spent25perc Spent50perc\n",
    "\n",
    "\n",
    "\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11a86c09",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "#SingleFactorDeSeq2_PA_tRNA.r\n",
    "\n",
    "library(dplyr)\n",
    "\n",
    "args = commandArgs(trailingOnly=TRUE)\n",
    "path <- args[1]\n",
    "current <- args[2]\n",
    "other <- args[3]\n",
    "\n",
    "################################################\n",
    "# PART 1: GATHER ALL INPUTS INTO ONE DATAFRAME #\n",
    "###############################################\n",
    "\n",
    "# LIST ALL INPUT FILES\n",
    "temp <- list.files(path = path, pattern=\"*.count.tRNA\")\n",
    "\n",
    "# START DATAFRAME WITH FIRST INPUT\n",
    "df <- as.data.frame(read.csv(paste0(path,temp[1]), header=FALSE, sep='\\t'))\n",
    "cols <- c('ID',substr(temp[1], 1, nchar(temp[1])-11))\n",
    "condition <- substr(temp[1], 1, nchar(temp[1])-16)\n",
    "\n",
    "# FILL DATAFRAME WITH THE REMAINING INPUTS\n",
    "for (file in temp[2:length(temp)]) {\n",
    "    cols <- c(cols,substr(file, 1, nchar(file)-11))\n",
    "    condition <- c(condition,substr(file, 1, nchar(file)-16))\n",
    "    df2 <- as.data.frame(read.csv(paste0(path,file), header=FALSE, sep='\\t'))\n",
    "    df <- full_join(df, df2, by='V1')\n",
    "}\n",
    "\n",
    "condition <- c(rep(\"PA_input\",3), rep(\"PA_capture\",3))\n",
    "condition <- factor(condition)\n",
    "warnings()\n",
    "df[is.na(df)] <- 0\n",
    "colnames(df) <- cols\n",
    "write.table(df,\"s4U_counts_table.tRNA.txt\", quote = FALSE, row.names = FALSE)\n",
    "\n",
    "\n",
    "##################################\n",
    "# PART 2: SINGLE FACTOR ANALYSIS #\n",
    "##################################\n",
    "\n",
    "# READ DATA INTO MATRIX\n",
    "data <- read.table(\"s4U_counts_table.tRNA.txt\", header=T, row.names=1)\n",
    "countdata <- as.matrix(data)\n",
    "\n",
    "# MAKE MATRIX A DESEQ2 OBJECT\n",
    "library(DESeq2)\n",
    "coldata <- data.frame(row.names=colnames(countdata), condition)\n",
    "print(coldata)\n",
    "dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)\n",
    "print(summary(dds))\n",
    "dds\n",
    "\n",
    "#filter the dataset - remove all zeros, influences padj\n",
    "nrow(dds)\n",
    "dds <- dds[ rowSums(counts(dds)) > 1, ]\n",
    "nrow(dds)\n",
    "\n",
    "# Run the DESeq pipeline\n",
    "dds <- DESeq(dds)\n",
    "res <- results(dds)\n",
    "print(res)\n",
    "\n",
    "\n",
    "conditions <- c(paste0(current,'_input'), paste0(other,'_captured'))\n",
    "\n",
    "for (current in conditions) { conditions <- conditions[conditions != current]\n",
    "                             \n",
    "        for (other in conditions) {\n",
    "            \n",
    "            print(current)\n",
    "            print(other)\n",
    "    \n",
    "            current_res = results(dds)\n",
    "            current_res\n",
    "            print(current_res)\n",
    "            summary(current_res)\n",
    "            current_ress_Orderedlfc <- current_res[order(current_res$log2FoldChange),]\n",
    "            head(current_res_Orderedlfc)\n",
    "\n",
    "            # Exploring and exporting results\n",
    "\n",
    "            pdf(file = paste0(current,'_VS_',other,'.pdf'),   # The directory you want to save the file in\n",
    "            width = 10, # The width of the plot in inches\n",
    "            height = 10) # The height of the plot in inches\n",
    "\n",
    "            plotMA(current_res, ylim=c(-10,10))\n",
    "            dev.off()\n",
    "\n",
    "            identify(current_res$baseMean, current_res$log2FoldChange, labels = row.names(current_res))\n",
    "            write.csv(as.data.frame(current_res), file= paste0(current,'_VS_',other,'.tRNA.csv'))\n",
    "\n",
    "            # Volcano plot\n",
    "            pdf(file = paste0(current,'_VS_',other,'.volcano.tRNA.pdf'),   # The directory you want to save the file in\n",
    "                width = 10, # The width of the plot in inches\n",
    "                height = 10) # The height of the plot in inches\n",
    "            par(mar=c(5,5,5,5), cex=1.0, cex.main=1.4, cex.axis=1.4, cex.lab=1.4)\n",
    "            topT <- as.data.frame(current_res)\n",
    "            #Adjusted P values (FDR Q values)\n",
    "            with(topT, plot(log2FoldChange, -log10(padj), pch=20, main=paste0(current,' VS ',other,'.volcano'), cex=1.0, xlab=bquote(~Log[2]~fold~change), ylab=bquote(~-log[10]~Q~value)))\n",
    "            with(subset(topT, padj<0.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(padj), pch=20, col=\"red\", cex=0.5))\n",
    "            #Add lines for absolute FC>2 and P-value cut-off at FDR Q<0.05\n",
    "            abline(v=0, col=\"black\", lty=3, lwd=1.0)\n",
    "            abline(v=-1, col=\"black\", lty=4, lwd=2.0)\n",
    "            abline(v=1, col=\"black\", lty=4, lwd=2.0)\n",
    "            abline(h=-log10(max(topT$pvalue[topT$padj<0.05], na.rm=TRUE)), col=\"black\", lty=4, lwd=2.0)\n",
    "            dev.off()\n",
    "        }\n",
    "    }\n",
    "            \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54871961",
   "metadata": {},
   "source": [
    "```\n",
    "Rscript SingleFactorDeSeq2_tRNA.r PA_analysis/ PA PA\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8574f910",
   "metadata": {},
   "source": [
    "## Multi-Factor Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3193b86",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "#MultiFactorDeSeq2_tRNA.r\n",
    "\n",
    "library(dplyr)\n",
    "\n",
    "args = commandArgs(trailingOnly=TRUE)\n",
    "path <- args[1]\n",
    "current <- args[2]\n",
    "other <- args[3]\n",
    "\n",
    "################################################\n",
    "# PART 1: GATHER ALL INPUTS INTO ONE DATAFRAME #\n",
    "###############################################\n",
    "\n",
    "# LIST ALL INPUT FILES\n",
    "pattern1 = paste0(current,\"*.count.tRNA\", collapse = \"|\")\n",
    "pattern2 = paste0(other,\"*.count.tRNA\", collapse = \"|\")\n",
    "temp <- list.files(path = path, pattern=glob2rx(pattern1))\n",
    "temp <- c(temp, list.files(path = path, pattern=glob2rx(pattern2)))\n",
    "\n",
    "\n",
    "# START DATAFRAME WITH FIRST INPUT\n",
    "df <- as.data.frame(read.csv(paste0(path,temp[1]), header=FALSE, sep='\\t'))\n",
    "cols <- c('ID',substr(temp[1], 1, nchar(temp[1])-6))\n",
    "condition <- substr(temp[1], 1, nchar(temp[1])-11)\n",
    "\n",
    "# FILL DATAFRAME WITH THE REMAINING INPUTS\n",
    "for (file in temp[2:length(temp)]) {\n",
    "    cols <- c(cols,substr(file, 1, nchar(file)-6))\n",
    "    condition <- c(condition,substr(file, 1, nchar(file)-11))\n",
    "    df2 <- as.data.frame(read.csv(paste0(path,file), header=FALSE, sep='\\t'))\n",
    "    df <- full_join(df, df2, by='V1')\n",
    "}\n",
    "\n",
    "condition <- factor(condition)\n",
    "warnings()\n",
    "df[is.na(df)] <- 0\n",
    "colnames(df) <- cols\n",
    "write.table(df,\"s4U_counts_table.tRNA.txt\", quote = FALSE, row.names = FALSE)\n",
    "\n",
    "\n",
    "##################################\n",
    "# PART 2: MULTI FACTOR ANALYSIS #\n",
    "##################################\n",
    "\n",
    "# READ DATA INTO MATRIX\n",
    "data <- read.table(\"s4U_counts_table.tRNA.txt\", header=T, row.names=1)\n",
    "countdata <- as.matrix(data)\n",
    "\n",
    "# MAKE MATRIX A DESEQ2 OBJECT\n",
    "library(DESeq2)\n",
    "#coldata <- data.frame(row.names=colnames(countdata), condition)\n",
    "#dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)\n",
    "#dds\n",
    "\n",
    "#filter the dataset - remove all zeros, influences padj\n",
    "#nrow(dds)\n",
    "#dds <- dds[ rowSums(counts(dds)) > 1, ]\n",
    "#nrow(dds)\n",
    "\n",
    "# Run the DESeq pipeline\n",
    "#dds <- DESeq(dds)\n",
    "#res <- results (dds)\n",
    "\n",
    "\n",
    "condition <- colnames(countdata)\n",
    "\n",
    "MF_condition <- c()\n",
    "MF_assay <- c()\n",
    "MF_replicate <- c()\n",
    "\n",
    "for (row in condition) {\n",
    "\n",
    "        MF_condition <- c(MF_condition, unlist(strsplit(as.character(row), '_'))[1])\n",
    "        MF_assay <- c(MF_assay, unlist(strsplit(as.character(row), '_'))[2])\n",
    "        MF_replicate <- c(MF_replicate, unlist(strsplit(as.character(row), '_'))[3])\n",
    "        }\n",
    "\n",
    "MF_coldata <- data.frame(row.names=colnames(countdata), MF_condition, MF_assay, MF_replicate)\n",
    "#MF_coldata <- MF_coldata[15:21,]\n",
    "MF_coldata <- MF_coldata[,1:2]\n",
    "\n",
    "print(MF_coldata)\n",
    "\n",
    "ddsMF <- DESeqDataSetFromMatrix(countData=countdata, colData=MF_coldata, design=~ MF_assay + MF_condition + MF_assay:MF_condition)\n",
    "\n",
    "\n",
    "print(ddsMF)\n",
    "\n",
    "nrow(ddsMF)\n",
    "ddsMF <- ddsMF[ rowSums(counts(ddsMF)) > 1, ]\n",
    "nrow(ddsMF)\n",
    "\n",
    "\n",
    "# Diferential expression analysis\n",
    "\n",
    "current <- args[2]\n",
    "other <- args[3]\n",
    "\n",
    "ddsMF <- DESeq(ddsMF, test=\"LRT\", reduced= ~ MF_assay + MF_condition) #creates the analysis\n",
    "resultsNames(ddsMF)\n",
    "\n",
    "\n",
    "print('MILESTONE')\n",
    "\n",
    "#normalized_counts <- counts(ddsMF, normalized=TRUE)\n",
    "#write.table(normalized_counts, file=\"normalized_counts.ddMF.txt\", sep=\",\", quote=F, col.names=NA)\n",
    "\n",
    "resMF = results(ddsMF, contrast=c(\"MF_condition\",current,other))\n",
    "resMF\n",
    "summary(resMF)\n",
    "resMF_Orderedlfc <- resMF[order(resMF$log2FoldChange),]\n",
    "head(resMF_Orderedlfc)\n",
    "\n",
    "plot.new()\n",
    "#identify(resMF$baseMean, resMF$log2FoldChange, labels = row.names(resMF))\n",
    "#write.csv(as.data.frame(resMF), file= paste0(current,'_VS_',other,'.MF.csv'))\n",
    "\n",
    "write.csv(as.data.frame(resMF_Orderedlfc), file= paste0(current,'_VS_',other,'.MF.tRNA.csv'))\n",
    "\n",
    "# Volcano plot\n",
    "\n",
    "\n",
    "pdf(file = paste0(current,'_VS_',other,'.volcano.tRNA.pdf'),   # The directory you want to save the file in\n",
    "    width = 10, # The width of the plot in inches\n",
    "    height = 10) # The height of the plot in inches\n",
    "par(mar=c(5,5,5,5), cex=1.0, cex.main=1.4, cex.axis=1.4, cex.lab=1.4)\n",
    "topT <- as.data.frame(resMF)\n",
    "#Adjusted P values (FDR Q values)\n",
    "with(topT, plot(log2FoldChange, -log10(padj), pch=20, main=paste0(current,' VS ',other,'volcano'), cex=1.0, xlab=bquote(~Log[2]~fold~change), ylab=bquote(~-log[10]~Q~value)))\n",
    "with(subset(topT, padj<0.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(padj), pch=20, col=\"red\", cex=0.5))\n",
    "#Add lines for absolute FC>2 and P-value cut-off at FDR Q<0.05\n",
    "abline(v=0, col=\"black\", lty=3, lwd=1.0)\n",
    "abline(v=-1, col=\"black\", lty=4, lwd=2.0)\n",
    "abline(v=1, col=\"black\", lty=4, lwd=2.0)\n",
    "abline(h=-log10(max(topT$pvalue[topT$padj<0.05], na.rm=TRUE)), col=\"black\", lty=4, lwd=2.0)\n",
    "dev.off()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}