iXnos_Fig3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We used the ribosome profiling of *E. coli* from Mohammad et al. (2019) to train the neural network model. This requires iXnos (https://github.com/lareaulab/iXnos). The scripts align.sh and trim_linker.pl were modified to process Mohammad's data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```console\n",
    "# align.sh\n",
    "# bowtie-build\n",
    "# rsem-prepare-reference\n",
    "\n",
    "GENOME_DIR=$1\n",
    "PROC_DIR=$2\n",
    "\n",
    "for i in `awk 'BEGIN{FS=\"/\"} NR==1 {print $NF}' sra.id`; do \\\n",
    "\tfasterq-dump ${i} -e 24\n",
    "done\n",
    "\n",
    "if [ ! -f $PROC_DIR/trimmed.fastq ]; then\n",
    "        cat $PROC_DIR/*.fastq | $PROC_DIR/trim_linker.pl > $PROC_DIR/trimmed.fastq\n",
    "fi\n",
    "\n",
    "if [ ! -f $GENOME_DIR/ecol.ncrna.1.ebwt ]; then\n",
    "\tbowtie-build $GENOME_DIR/ecol.ncrna.fa $GENOME_DIR/ecol.ncrna\n",
    "fi\n",
    "\n",
    "if [ ! -f $GENOME_DIR/ecol.transcripts.13cds10.1.ebwt ]; then\n",
    "        bowtie-build $GENOME_DIR/ecol.transcripts.13cds10.fa $GENOME_DIR/ecol.transcripts.13cds10\n",
    "fi\n",
    "\n",
    "if [ ! -f $GENOME_DIR/ecol.transcripts.13cds10.idx.fa ]; then\n",
    "\trsem-prepare-reference $GENOME_DIR/ecol.transcripts.13cds10.fa $GENOME_DIR/ecol.transcripts.13cds10\n",
    "fi\n",
    "\n",
    "if [ ! -f $PROC_DIR/not_ncrna.fastq ]; then\n",
    "\tbowtie -v 2 -p 36 -S --un $PROC_DIR/not_ncrna.fastq \\\n",
    "        $GENOME_DIR/ecol.ncrna \\\n",
    "        $PROC_DIR/trimmed.fastq > $PROC_DIR/ncrna.sam 2> $PROC_DIR/ncrna.bowtiestats\n",
    "fi\n",
    "\n",
    "if [ ! -f $PROC_DIR/footprints.sam ]; then\n",
    "bowtie -a --norc -v 2 -p 36 -S --un $PROC_DIR/unmapped.fastq \\\n",
    "       $GENOME_DIR/ecol.transcripts.13cds10 \\\n",
    "       $PROC_DIR/not_ncrna.fastq > $PROC_DIR/footprints.sam 2> $PROC_DIR/footprints.bowtiestats\n",
    "fi\n",
    "\n",
    "rsem-calculate-expression --sam $PROC_DIR/footprints.sam $GENOME_DIR/ecol.transcripts.13cds10 $PROC_DIR/rsem 2> $PROC_DIR/rsem.stderr\n",
    "\n",
    "samtools view -h $PROC_DIR/rsem.transcript.bam > $PROC_DIR/rsem.transcript.sam\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```console\n",
    "# trim_linker.pl \n",
    "#!/usr/bin/perl\n",
    "use strict;\n",
    "\n",
    "# Trim off end sequence matching the linker from Mohammad's experiments:\n",
    "# CTGTAGGCACCATCAAT... (it can be longer than this)\n",
    "\n",
    "my $printstring = \"\";\n",
    "\n",
    "while(<>)\n",
    "{\n",
    "    my $header;\n",
    "    my $seq;\n",
    "    my $qheader;\n",
    "    my $qual;\n",
    "    \n",
    "    if (/^@/)\n",
    "    { \n",
    "\tprint $printstring;\n",
    "\t$printstring = \"\";\n",
    "\n",
    "\t$header = $_; chomp $header;\n",
    "\t$seq = <>; chomp $seq;\n",
    "\t$qheader = <>; chomp $qheader;\n",
    "\t$qual = <>; chomp $qual;\n",
    "    }\n",
    "\n",
    "    my $trimmed;\n",
    "\n",
    "    # no random barcodes at beginning\n",
    "    if ( ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCACCATCAAT/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCACCATCAA$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCACCATCA$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCACCATC$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCACCAT$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCACCA$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCACC$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCAC$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGCA$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGGC$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAGG$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTAG$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGTA$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTGT$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CTG$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)CT$/ ) or\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)C$/ ) or ## ...eats up anything ending in a C... but that's ok\n",
    "\t ( ($trimmed) = $seq =~ /^([ACGTN]+)$/ )  ## ...no evidence of linker..\n",
    "\t)\n",
    "    {\n",
    "    \n",
    "\tmy $length = length($trimmed);\n",
    "\n",
    "\tif ($length > 10)\n",
    "\t{\n",
    "\t    $printstring = $header . \"\\n\" . $trimmed . \"\\n\" . $qheader . \"\\n\";\n",
    "\t    $printstring .= substr( $qual, 0, $length ); ## fix this bug in other versions of the script!!\n",
    "\t    $printstring .= \"\\n\";\n",
    "\t}\n",
    "    }\n",
    "    \n",
    "    print $printstring;\n",
    "    $printstring = \"\";\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first prepared *E. coli* transcript files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```console\n",
    "# Get all CDS\n",
    "zcat Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.42.gtf.gz \\\n",
    "| gtfToGenePred stdin stdout | awk '$6!=$7' | genePredToBed stdin ecol.bed\n",
    "bedtools getfasta -fi Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.dna.toplevel.fa -bed ecol.bed -fo ecol.transcripts.txt -s -split -name -tab\n",
    "sed -i 's/(+)//;s/(-)//' ecol.transcripts.txt\n",
    "# Get padding of 13 and 10 at the 5' and 3' UTRs\n",
    "# Remove frameshift ORF AAC75929 (prfB)\n",
    "zcat Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.42.gtf.gz \\\n",
    "| gtfToGenePred stdin stdout | awk '$6!=$7' | genePredToBed stdin stdout \\\n",
    "| awk 'BEGIN{OFS=\"\\t\"} $10==1 {if($6~/+/){print $1,$2-13,$3+10,$4,$5,$6} \n",
    "else{print $1,$2-10,$3+13,$4,$5,$6}}' > ecol.13cds10.bed\n",
    "bedtools getfasta -fi Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.dna.toplevel.fa -bed ecol.13cds10.bed -fo ecol.transcripts.13cds10.txt -s -name -tab\n",
    "sed -i 's/(+)//;s/(-)//' ecol.transcripts.13cds10.txt\n",
    "\n",
    "# Get padding of 50 and 10 at the 5' and 3' UTRs for opening energy calculation\n",
    "# Remove frameshift ORF AAC75929 (prfB)\n",
    "zcat Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.42.gtf.gz \\\n",
    "| gtfToGenePred stdin stdout | awk '$6!=$7' | genePredToBed stdin stdout \\\n",
    "| awk 'BEGIN{OFS=\"\\t\"} $10==1 {if($6~/+/){print $1,$2-50,$3+10,$4,$5,$6} \n",
    "else{print $1,$2-10,$3+50,$4,$5,$6}}' > ecol.50cds10.bed\n",
    "bedtools getfasta \\\n",
    "-fi Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.dna.toplevel.fa \\\n",
    "-bed ecol.50cds10.bed \\\n",
    "-s -name -tab \\\n",
    "| sed 's/(+)//;s/(-)//' > ecol.transcripts.50cds10.txt\n",
    "join -t$'\\t' \\\n",
    "<(sort -k1,1 tx2gene.txt) \\\n",
    "<(sort -k1,1 ecol.transcripts.50cds10.txt) \\\n",
    "| awk '{print \">\" $2 \"\\n\" $3}' > ecol.transcripts.50cds10.fa\n",
    "\n",
    "# Get transcripts with ATG GTG TTG start codons\n",
    "ecol.transcripts.txt\n",
    "join -t$'\\t' \\\n",
    "<(awk '{print $2,$1}' ecol.transcripts.txt | awk '/^ATG|^GTG|^TTG/ {print $2}' | sort) \\\n",
    "<(sort -k1,1 ecol.transcripts.13cds10.txt) \\\n",
    "| awk '{print \">\" $1 \"\\n\" $2}' > ecol.transcripts.13cds10.fa\n",
    "awk 'BEGIN{RS=\">\"}NR>1{sub(\"\\n\",\"\\t\"); gsub(\"\\n\",\"\"); print $0}' ecol.transcripts.13cds10.fa \\\n",
    "| awk 'BEGIN{OFS=\"\\t\"}{print $1,\"13\",length($2),\"10\"}' > ecol.transcripts.13cds10.lengths.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We downloaded and examine the individual sra files. SRR7759806 and SRR7759807 showed the best triplet periodicity. They were merged for subsequent analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```console\n",
    "# Download sra files\n",
    "mkdir ~/src/iXnos/expts/mohammad\n",
    "cd ~/src/iXnos/expts/mohammad\n",
    "wget -O SRP158945_metadata.csv 'http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term=SRP158945'\n",
    "sed 's/,/\\t/g' SRP158945_metadata.csv | cut -f1,11 \\\n",
    "| awk '$1~/SRR/ {print \"/sra/sra-instant/reads/ByRun/sra/SRR/SRR775/\" $1 \"/\" $1 \".sra\"}' > sra.id\n",
    "ascp -l640M -T -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --user=anonftp --host=ftp.ncbi.nlm.nih.gov --mode=recv --file-list=sra.id .\n",
    "\n",
    "# Align reads\n",
    "for i in *.fastq; do \\\n",
    "  cd ~/src/iXnos/expts\n",
    "  mkdir ${i%%.*}\n",
    "  mv $i ${i%%.*}\n",
    "  cp align.sh trim_linker.pl ${i%%.*}\n",
    "  cd ${i%%.*}\n",
    "  ./align.sh ~/src/iXnos/genome_data/ .\n",
    "done\n",
    "\n",
    "# Examine plots\n",
    "for i in SRR7759805 SRR7759808 SRR7759809 SRR7759810 SRR7759811 SRR7759812 SRR7759813 SRR7759814; do \\\n",
    "  cd ~/src/iXnos/expts/mohammad;\n",
    "  mv $i process;\n",
    "  python ~/anaconda2/pkgs/iXnos/plot_frames.py;\n",
    "  mv process $i;\n",
    "  mv plots plots_$i;\n",
    "done\n",
    "\n",
    "# Merge SRR7759806 and SRR7759807\n",
    "cd ~/src/iXnos/expts/mohammad\n",
    "mkdir process\n",
    "samtools merge process/footprints.sam \\\n",
    "SRR7759806/footprints.sam \\\n",
    "SRR7759807/footprints.sam\n",
    "cd process\n",
    "rsem-calculate-expression --sam footprints.sam ~/src/iXnos/genome_data/ecol.transcripts.13cds10 rsem 2> rsem.stderr\n",
    "samtools view -h rsem.transcript.bam > rsem.transcript.sam\n",
    "python ~/anaconda2/pkgs/iXnos/plot_frames.py\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then train a neural network model using the processed ribosome profiling data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import lasagne as ls\n",
    "import iXnos.interface as inter\n",
    "import iXnos.optimizecodons as opt\n",
    "import iXnos.process as proc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/chunshenlim/src/iXnos/expts/mohammad//process/rsem.transcript.mapped.wts.sam'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expt_dir = \"/home/chunshenlim/src/iXnos/expts/mohammad/\"\n",
    "sam_fname = expt_dir + \"/process/rsem.transcript.sam\"\n",
    "inter.edit_sam_file(\n",
    "    expt_dir, sam_fname, filter_unmapped=True, \n",
    "    sam_add_simple_map_wts=False, RSEM_add_map_wts=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "expt_dir = \"/home/chunshenlim/src/iXnos/expts/mohammad/\"\n",
    "sam_fname = expt_dir + \"/process/rsem.transcript.mapped.wts.sam\"\n",
    "gene_len_fname = \"/home/chunshenlim/src/iXnos/genome_data/ecol.transcripts.13cds10.lengths.txt\"\n",
    "inter.do_frame_and_size_analysis(\n",
    "    expt_dir, sam_fname, gene_len_fname, min_size=13, \n",
    "    max_size=36, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "expt_dir = \"/home/chunshenlim/src/iXnos/expts/mohammad/\"\n",
    "sam_fname = expt_dir + \"/process/rsem.transcript.mapped.wts.sam\"\n",
    "gene_len_fname = \"/home/chunshenlim/src/iXnos/genome_data/ecol.transcripts.13cds10.lengths.txt\"\n",
    "# A site offset rules\n",
    "shift_dict = {25:{0:12, 1:False, 2:False}}\n",
    "# Number of codons to ignore at either end of CDS in model training\n",
    "cod_trunc_5p = 20\n",
    "cod_trunc_3p = 20\n",
    "# Minimum and maximum allowed footprint sizes\n",
    "min_fp_size = 25\n",
    "max_fp_size = 25\n",
    "# Choose number of genes in training and test sets\n",
    "num_tr_genes = 333\n",
    "num_te_genes = 167\n",
    "# Set coverage cutoffs for genes in training/test sets\n",
    "# Minimum footprint counts per gene\n",
    "min_cts_per_gene = 200\n",
    "# Minimum A site codons with nonzero counts\n",
    "min_cod_w_data = 100\n",
    "# Transcriptome file, e.g. for yeast\n",
    "gene_seq_fname = \"/home/chunshenlim/src/iXnos/genome_data/ecol.transcripts.13cds10.fa\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file /home/chunshenlim/src/iXnos/expts/mohammad//process/cts_by_codon.size.25.25.txt already exists\n",
      "file /home/chunshenlim/src/iXnos/expts/mohammad//process/outputs.size.25.25.txt already exists\n",
      "making file /home/chunshenlim/src/iXnos/expts/mohammad//process/tr_set_bounds.size.25.25.trunc.20.20.min_cts.200.min_cod.100.top.500.txt\n",
      "making file /home/chunshenlim/src/iXnos/expts/mohammad//process/te_set_bounds.size.25.25.trunc.20.20.min_cts.200.min_cod.100.top.500.txt\n",
      "610 genes meeting data cutoff\n"
     ]
    }
   ],
   "source": [
    "inter.process_sam_file(\n",
    "    expt_dir, sam_fname, gene_seq_fname, gene_len_fname,\n",
    "    shift_dict, cod_trunc_5p, cod_trunc_3p, min_fp_size,\n",
    "    max_fp_size, num_tr_genes, num_te_genes, min_cts_per_gene,\n",
    "    min_cod_w_data, paralog_groups_fname=False,\n",
    "    overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "name = \"high_salt\"\n",
    "expt_dir = \"/home/chunshenlim/src/iXnos/expts/mohammad/\"\n",
    "outputs_fname = expt_dir + \"/process/outputs.size.25.25.txt\"\n",
    "gene_seq_fname = \"/home/chunshenlim/src/iXnos/genome_data/ecol.transcripts.13cds10.fa\"\n",
    "gene_len_fname = \"/home/chunshenlim/src/iXnos/genome_data/ecol.transcripts.13cds10.lengths.txt\"\n",
    "tr_codons_fname = expt_dir + \"/process/tr_set_bounds.size.25.25.trunc.20.20.min_cts.200.min_cod.100.top.500.txt\"\n",
    "te_codons_fname = expt_dir + \"/process/te_set_bounds.size.25.25.trunc.20.20.min_cts.200.min_cod.100.top.500.txt\"\n",
    "rel_cod_idxs = range(-3,3)\n",
    "rel_nt_idxs = [cod_idx * 3 + i for cod_idx in rel_cod_idxs for i in range(3)]\n",
    "widths = [200]\n",
    "nonlinearity = \"tanh\"\n",
    "update_method = \"nesterov\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Making dense mlp\n",
      "nonnegative model\n"
     ]
    }
   ],
   "source": [
    "# Create neural network\n",
    "neural_net = inter.make_lasagne_feedforward_nn(\n",
    "    name, expt_dir, gene_seq_fname, gene_len_fname, tr_codons_fname,\n",
    "    te_codons_fname, outputs_fname, rel_cod_idxs=rel_cod_idxs,\n",
    "    rel_nt_idxs=rel_nt_idxs, nonlinearity=nonlinearity, widths=widths, \n",
    "    update_method=update_method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "learning rate: 0.01\n",
      "Epoch 1 took 4.049s\n",
      "  training loss:\t\t6.721418\n",
      "  test loss:\t\t6.753441\n",
      "learning rate: 0.00941176470588\n",
      "Epoch 2 took 3.730s\n",
      "  training loss:\t\t6.584628\n",
      "  test loss:\t\t6.733057\n",
      "learning rate: 0.00888888888889\n",
      "Epoch 3 took 3.662s\n",
      "  training loss:\t\t6.558858\n",
      "  test loss:\t\t6.786018\n",
      "learning rate: 0.00842105263158\n",
      "Epoch 4 took 3.472s\n",
      "  training loss:\t\t6.562880\n",
      "  test loss:\t\t6.683421\n",
      "learning rate: 0.008\n",
      "Epoch 5 took 3.565s\n",
      "  training loss:\t\t6.534845\n",
      "  test loss:\t\t6.669208\n",
      "learning rate: 0.00761904761905\n",
      "Epoch 6 took 3.391s\n",
      "  training loss:\t\t6.526851\n",
      "  test loss:\t\t6.671323\n",
      "learning rate: 0.00727272727273\n",
      "Epoch 7 took 3.328s\n",
      "  training loss:\t\t6.513343\n",
      "  test loss:\t\t6.639606\n",
      "learning rate: 0.00695652173913\n",
      "Epoch 8 took 3.707s\n",
      "  training loss:\t\t6.481373\n",
      "  test loss:\t\t6.625272\n",
      "learning rate: 0.00666666666667\n",
      "Epoch 9 took 3.471s\n",
      "  training loss:\t\t6.458457\n",
      "  test loss:\t\t6.679230\n",
      "learning rate: 0.0064\n",
      "Epoch 10 took 3.591s\n",
      "  training loss:\t\t6.439775\n",
      "  test loss:\t\t6.598577\n",
      "learning rate: 0.00615384615385\n",
      "Epoch 11 took 3.328s\n",
      "  training loss:\t\t6.408657\n",
      "  test loss:\t\t6.599809\n",
      "learning rate: 0.00592592592593\n",
      "Epoch 12 took 3.541s\n",
      "  training loss:\t\t6.401818\n",
      "  test loss:\t\t6.596228\n",
      "learning rate: 0.00571428571429\n",
      "Epoch 13 took 3.614s\n",
      "  training loss:\t\t6.377088\n",
      "  test loss:\t\t6.609076\n",
      "learning rate: 0.00551724137931\n",
      "Epoch 14 took 3.361s\n",
      "  training loss:\t\t6.356090\n",
      "  test loss:\t\t6.590171\n",
      "learning rate: 0.00533333333333\n",
      "Epoch 15 took 3.542s\n",
      "  training loss:\t\t6.348844\n",
      "  test loss:\t\t6.564938\n",
      "learning rate: 0.00516129032258\n",
      "Epoch 16 took 3.995s\n",
      "  training loss:\t\t6.317804\n",
      "  test loss:\t\t6.563464\n",
      "learning rate: 0.005\n",
      "Epoch 17 took 3.472s\n",
      "  training loss:\t\t6.256989\n",
      "  test loss:\t\t6.578596\n",
      "learning rate: 0.00484848484848\n",
      "Epoch 18 took 3.676s\n",
      "  training loss:\t\t6.284878\n",
      "  test loss:\t\t6.597509\n",
      "learning rate: 0.00470588235294\n",
      "Epoch 19 took 3.487s\n",
      "  training loss:\t\t6.274387\n",
      "  test loss:\t\t6.622429\n",
      "learning rate: 0.00457142857143\n",
      "Epoch 20 took 3.814s\n",
      "  training loss:\t\t6.262929\n",
      "  test loss:\t\t6.611515\n",
      "learning rate: 0.00444444444444\n",
      "Epoch 21 took 3.764s\n",
      "  training loss:\t\t6.225772\n",
      "  test loss:\t\t6.594164\n",
      "learning rate: 0.00432432432432\n",
      "Epoch 22 took 3.819s\n",
      "  training loss:\t\t6.217219\n",
      "  test loss:\t\t6.601016\n",
      "learning rate: 0.00421052631579\n",
      "Epoch 23 took 3.825s\n",
      "  training loss:\t\t6.216986\n",
      "  test loss:\t\t6.630048\n",
      "learning rate: 0.00410256410256\n",
      "Epoch 24 took 4.049s\n",
      "  training loss:\t\t6.196393\n",
      "  test loss:\t\t6.604311\n",
      "learning rate: 0.004\n",
      "Epoch 25 took 3.458s\n",
      "  training loss:\t\t6.177392\n",
      "  test loss:\t\t6.608903\n",
      "learning rate: 0.00390243902439\n",
      "Epoch 26 took 3.992s\n",
      "  training loss:\t\t6.155717\n",
      "  test loss:\t\t6.616409\n",
      "learning rate: 0.00380952380952\n",
      "Epoch 27 took 4.220s\n",
      "  training loss:\t\t6.144956\n",
      "  test loss:\t\t6.609216\n",
      "learning rate: 0.00372093023256\n",
      "Epoch 28 took 3.687s\n",
      "  training loss:\t\t6.121198\n",
      "  test loss:\t\t6.659405\n",
      "learning rate: 0.00363636363636\n",
      "Epoch 29 took 3.595s\n",
      "  training loss:\t\t6.116891\n",
      "  test loss:\t\t6.648084\n",
      "learning rate: 0.00355555555556\n",
      "Epoch 30 took 3.459s\n",
      "  training loss:\t\t6.108185\n",
      "  test loss:\t\t6.656893\n",
      "learning rate: 0.00347826086957\n",
      "Epoch 31 took 3.956s\n",
      "  training loss:\t\t6.097938\n",
      "  test loss:\t\t6.693715\n",
      "learning rate: 0.00340425531915\n",
      "Epoch 32 took 3.579s\n",
      "  training loss:\t\t6.034131\n",
      "  test loss:\t\t6.637483\n",
      "learning rate: 0.00333333333333\n",
      "Epoch 33 took 3.571s\n",
      "  training loss:\t\t6.065403\n",
      "  test loss:\t\t6.683349\n",
      "learning rate: 0.00326530612245\n",
      "Epoch 34 took 3.664s\n",
      "  training loss:\t\t6.058962\n",
      "  test loss:\t\t6.709966\n",
      "learning rate: 0.0032\n",
      "Epoch 35 took 3.734s\n",
      "  training loss:\t\t6.033556\n",
      "  test loss:\t\t6.735769\n",
      "learning rate: 0.00313725490196\n",
      "Epoch 36 took 3.776s\n",
      "  training loss:\t\t6.019363\n",
      "  test loss:\t\t6.756196\n",
      "learning rate: 0.00307692307692\n",
      "Epoch 37 took 3.534s\n",
      "  training loss:\t\t6.013037\n",
      "  test loss:\t\t6.710579\n",
      "learning rate: 0.00301886792453\n",
      "Epoch 38 took 4.529s\n",
      "  training loss:\t\t5.983729\n",
      "  test loss:\t\t6.706085\n",
      "learning rate: 0.00296296296296\n",
      "Epoch 39 took 3.665s\n",
      "  training loss:\t\t5.964023\n",
      "  test loss:\t\t6.771792\n",
      "learning rate: 0.00290909090909\n",
      "Epoch 40 took 3.597s\n",
      "  training loss:\t\t5.970914\n",
      "  test loss:\t\t6.737681\n",
      "learning rate: 0.00285714285714\n",
      "Epoch 41 took 3.918s\n",
      "  training loss:\t\t5.937856\n",
      "  test loss:\t\t6.830126\n",
      "learning rate: 0.00280701754386\n",
      "Epoch 42 took 3.578s\n",
      "  training loss:\t\t5.929502\n",
      "  test loss:\t\t6.757503\n",
      "learning rate: 0.00275862068966\n",
      "Epoch 43 took 3.575s\n",
      "  training loss:\t\t5.929956\n",
      "  test loss:\t\t6.889574\n",
      "learning rate: 0.00271186440678\n",
      "Epoch 44 took 3.963s\n",
      "  training loss:\t\t5.928693\n",
      "  test loss:\t\t6.850208\n",
      "learning rate: 0.00266666666667\n",
      "Epoch 45 took 3.792s\n",
      "  training loss:\t\t5.907813\n",
      "  test loss:\t\t6.757326\n",
      "learning rate: 0.00262295081967\n",
      "Epoch 46 took 4.088s\n",
      "  training loss:\t\t5.892130\n",
      "  test loss:\t\t6.794175\n",
      "learning rate: 0.00258064516129\n",
      "Epoch 47 took 4.145s\n",
      "  training loss:\t\t5.874660\n",
      "  test loss:\t\t6.791802\n",
      "learning rate: 0.00253968253968\n",
      "Epoch 48 took 3.638s\n",
      "  training loss:\t\t5.867027\n",
      "  test loss:\t\t6.824779\n",
      "learning rate: 0.0025\n",
      "Epoch 49 took 3.727s\n",
      "  training loss:\t\t5.852034\n",
      "  test loss:\t\t6.825814\n",
      "learning rate: 0.00246153846154\n",
      "Epoch 50 took 3.842s\n",
      "  training loss:\t\t5.843767\n",
      "  test loss:\t\t6.870696\n",
      "learning rate: 0.00242424242424\n",
      "Epoch 51 took 3.575s\n",
      "  training loss:\t\t5.827601\n",
      "  test loss:\t\t6.816054\n",
      "learning rate: 0.00238805970149\n",
      "Epoch 52 took 3.738s\n",
      "  training loss:\t\t5.800953\n",
      "  test loss:\t\t6.930156\n",
      "learning rate: 0.00235294117647\n",
      "Epoch 53 took 3.878s\n",
      "  training loss:\t\t5.813124\n",
      "  test loss:\t\t6.846173\n",
      "learning rate: 0.00231884057971\n",
      "Epoch 54 took 3.648s\n",
      "  training loss:\t\t5.796495\n",
      "  test loss:\t\t6.851785\n",
      "learning rate: 0.00228571428571\n",
      "Epoch 55 took 3.680s\n",
      "  training loss:\t\t5.786667\n",
      "  test loss:\t\t6.866404\n",
      "learning rate: 0.00225352112676\n",
      "Epoch 56 took 3.965s\n",
      "  training loss:\t\t5.776662\n",
      "  test loss:\t\t6.954507\n",
      "learning rate: 0.00222222222222\n",
      "Epoch 57 took 3.722s\n",
      "  training loss:\t\t5.759980\n",
      "  test loss:\t\t6.866863\n",
      "learning rate: 0.00219178082192\n",
      "Epoch 58 took 3.636s\n",
      "  training loss:\t\t5.722340\n",
      "  test loss:\t\t6.892121\n",
      "learning rate: 0.00216216216216\n",
      "Epoch 59 took 3.823s\n",
      "  training loss:\t\t5.736228\n",
      "  test loss:\t\t6.984660\n",
      "learning rate: 0.00213333333333\n",
      "Epoch 60 took 3.819s\n",
      "  training loss:\t\t5.720447\n",
      "  test loss:\t\t6.891850\n",
      "learning rate: 0.00210526315789\n",
      "Epoch 61 took 3.683s\n",
      "  training loss:\t\t5.714925\n",
      "  test loss:\t\t6.907398\n",
      "learning rate: 0.00207792207792\n",
      "Epoch 62 took 3.969s\n",
      "  training loss:\t\t5.695549\n",
      "  test loss:\t\t6.949372\n",
      "learning rate: 0.00205128205128\n",
      "Epoch 63 took 3.964s\n",
      "  training loss:\t\t5.689983\n",
      "  test loss:\t\t7.064771\n",
      "learning rate: 0.0020253164557\n",
      "Epoch 64 took 3.717s\n",
      "  training loss:\t\t5.679642\n",
      "  test loss:\t\t6.979869\n",
      "learning rate: 0.002\n",
      "Epoch 65 took 3.908s\n",
      "  training loss:\t\t5.655892\n",
      "  test loss:\t\t7.010615\n",
      "learning rate: 0.00197530864198\n",
      "Epoch 66 took 3.552s\n",
      "  training loss:\t\t5.654490\n",
      "  test loss:\t\t6.961313\n",
      "learning rate: 0.0019512195122\n",
      "Epoch 67 took 3.645s\n",
      "  training loss:\t\t5.650347\n",
      "  test loss:\t\t6.964334\n",
      "learning rate: 0.00192771084337\n",
      "Epoch 68 took 3.947s\n",
      "  training loss:\t\t5.625448\n",
      "  test loss:\t\t7.120770\n",
      "learning rate: 0.00190476190476\n",
      "Epoch 69 took 3.991s\n",
      "  training loss:\t\t5.622757\n",
      "  test loss:\t\t7.011893\n",
      "learning rate: 0.00188235294118\n",
      "Epoch 70 took 4.131s\n",
      "  training loss:\t\t5.587313\n",
      "  test loss:\t\t7.088116\n",
      "learning rate: 0.00186046511628\n",
      "Epoch 71 took 3.908s\n",
      "  training loss:\t\t5.583459\n",
      "  test loss:\t\t7.023372\n",
      "learning rate: 0.00183908045977\n",
      "Epoch 72 took 3.735s\n",
      "  training loss:\t\t5.587992\n",
      "  test loss:\t\t7.099008\n",
      "learning rate: 0.00181818181818\n",
      "Epoch 73 took 3.990s\n",
      "  training loss:\t\t5.571424\n",
      "  test loss:\t\t7.155041\n",
      "learning rate: 0.00179775280899\n",
      "Epoch 74 took 3.753s\n",
      "  training loss:\t\t5.568525\n",
      "  test loss:\t\t7.101817\n",
      "learning rate: 0.00177777777778\n",
      "Epoch 75 took 3.712s\n",
      "  training loss:\t\t5.547688\n",
      "  test loss:\t\t7.097993\n",
      "learning rate: 0.00175824175824\n",
      "Epoch 76 took 3.941s\n",
      "  training loss:\t\t5.540742\n",
      "  test loss:\t\t7.120319\n",
      "learning rate: 0.00173913043478\n",
      "Epoch 77 took 3.978s\n",
      "  training loss:\t\t5.528085\n",
      "  test loss:\t\t7.130242\n",
      "learning rate: 0.00172043010753\n",
      "Epoch 78 took 3.831s\n",
      "  training loss:\t\t5.511949\n",
      "  test loss:\t\t7.111835\n",
      "learning rate: 0.00170212765957\n",
      "Epoch 79 took 3.817s\n",
      "  training loss:\t\t5.487768\n",
      "  test loss:\t\t7.073235\n",
      "learning rate: 0.00168421052632\n",
      "Epoch 80 took 3.653s\n",
      "  training loss:\t\t5.497879\n",
      "  test loss:\t\t7.150865\n",
      "learning rate: 0.00166666666667\n",
      "Epoch 81 took 3.899s\n",
      "  training loss:\t\t5.482321\n",
      "  test loss:\t\t7.218219\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "learning rate: 0.00164948453608\n",
      "Epoch 82 took 3.691s\n",
      "  training loss:\t\t5.466751\n",
      "  test loss:\t\t7.245365\n",
      "learning rate: 0.00163265306122\n",
      "Epoch 83 took 3.721s\n",
      "  training loss:\t\t5.460148\n",
      "  test loss:\t\t7.213522\n",
      "learning rate: 0.00161616161616\n",
      "Epoch 84 took 3.909s\n",
      "  training loss:\t\t5.453848\n",
      "  test loss:\t\t7.203413\n",
      "learning rate: 0.0016\n",
      "Epoch 85 took 3.633s\n",
      "  training loss:\t\t5.437642\n",
      "  test loss:\t\t7.160859\n",
      "learning rate: 0.00158415841584\n",
      "Epoch 86 took 3.981s\n",
      "  training loss:\t\t5.431930\n",
      "  test loss:\t\t7.318815\n",
      "learning rate: 0.00156862745098\n",
      "Epoch 87 took 3.712s\n",
      "  training loss:\t\t5.416323\n",
      "  test loss:\t\t7.207446\n",
      "learning rate: 0.00155339805825\n",
      "Epoch 88 took 4.012s\n",
      "  training loss:\t\t5.399868\n",
      "  test loss:\t\t7.195398\n",
      "learning rate: 0.00153846153846\n",
      "Epoch 89 took 3.886s\n",
      "  training loss:\t\t5.385329\n",
      "  test loss:\t\t7.284057\n",
      "learning rate: 0.00152380952381\n",
      "Epoch 90 took 3.878s\n",
      "  training loss:\t\t5.357571\n",
      "  test loss:\t\t7.320951\n",
      "learning rate: 0.00150943396226\n",
      "Epoch 91 took 3.979s\n",
      "  training loss:\t\t5.364010\n",
      "  test loss:\t\t7.312585\n",
      "learning rate: 0.0014953271028\n",
      "Epoch 92 took 3.796s\n",
      "  training loss:\t\t5.348741\n",
      "  test loss:\t\t7.268694\n",
      "learning rate: 0.00148148148148\n",
      "Epoch 93 took 4.724s\n",
      "  training loss:\t\t5.338409\n",
      "  test loss:\t\t7.271367\n",
      "learning rate: 0.00146788990826\n",
      "Epoch 94 took 3.659s\n",
      "  training loss:\t\t5.320792\n",
      "  test loss:\t\t7.328790\n",
      "learning rate: 0.00145454545455\n",
      "Epoch 95 took 3.756s\n",
      "  training loss:\t\t5.313194\n",
      "  test loss:\t\t7.327226\n",
      "learning rate: 0.00144144144144\n",
      "Epoch 96 took 4.057s\n",
      "  training loss:\t\t5.302305\n",
      "  test loss:\t\t7.300800\n",
      "learning rate: 0.00142857142857\n",
      "Epoch 97 took 3.735s\n",
      "  training loss:\t\t5.257737\n",
      "  test loss:\t\t7.360095\n",
      "learning rate: 0.00141592920354\n",
      "Epoch 98 took 3.882s\n",
      "  training loss:\t\t5.279853\n",
      "  test loss:\t\t7.381225\n",
      "learning rate: 0.00140350877193\n",
      "Epoch 99 took 3.845s\n",
      "  training loss:\t\t5.261516\n",
      "  test loss:\t\t7.346964\n",
      "learning rate: 0.00139130434783\n",
      "Epoch 100 took 4.009s\n",
      "  training loss:\t\t5.262769\n",
      "  test loss:\t\t7.322045\n"
     ]
    }
   ],
   "source": [
    "# Run neural network for a desired Number of epochs\n",
    "num_epochs = 100\n",
    "neural_net.run_epochs(num_epochs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Making dense mlp\n",
      "nonnegative model\n"
     ]
    }
   ],
   "source": [
    "RANDOM_SEED = 12345\n",
    "RANDOM = np.random.RandomState(RANDOM_SEED)\n",
    "ls.random.set_rng(RANDOM)\n",
    "nn_dir = \"/home/chunshenlim/src/iXnos/expts/mohammad/lasagne_nn/high_salt/\"\n",
    "epoch = 19 # lowest test loss\n",
    "nt_feats = True\n",
    "my_nn = inter.load_lasagne_feedforward_nn(nn_dir, epoch)\n",
    "\n",
    "# Test for scoring\n",
    "cod_seq = \"ATGGCCTGTGATGAATTTGGTCACATTAAACTCAACCCTCAACGCTCCACTGTCTGGTATTAG\"\n",
    "rel_cod_idxs = range(-3,3)\n",
    "rel_nt_idxs = [cod_idx * 3 + i for cod_idx in rel_cod_idxs for i in range(3)]\n",
    "cds_score = opt.score_cod_seq_full(cod_seq, my_nn, rel_cod_idxs, rel_nt_idxs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17.00791147385095"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cds_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```console\n",
    "# Score the PSI:Biolog targets\n",
    "python ~/avoidance/src/Avoidance-v2/Random_forest/iXnos.py \\\n",
    "-i pET21_NESG.fa \\\n",
    "-d ~/avoidance/src/iXnos/expts/mohammad/lasagne_nn/high_salt \\\n",
    "-e 19 -p 32\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}