# Change what's inside the quotes into your MBL STAMPS server username, or something else unique
# reminder! hit shift-ENTER to execute cell & move onto next cell
username="CHANGEME"

# configure paths etc.
import sys
import os
os.environ['PYTHONPATH'] = '/usr/local/share/khmer/python'
sys.path.insert(0, '/usr/local/share/khmer/python')

import khmer
import screed

cd /mnt

mkdir -p $username

cd $username

import random
random.seed(1)

x = ["A"] + ["G"] + ["C"] + ["T"]
x = x*1000
random.shuffle(x)
x = "".join(x)

y = ["A"] + ["G"] + ["C"] + ["T"]
y = y*1000
random.shuffle(y)
y = "".join(y)

print 'x is', x[:100]
print 'y is', y[:100]

outfp = open('metagenome.fa', 'w')
print >>outfp, ">x 1"
print >>outfp, x
print >>outfp, ">y 2"
print >>outfp, y
outfp.close()

!python /usr/local/share/2012-paper-diginorm/pipeline/make-biased-reads.py metagenome.fa | head -100000 > reads.fa

!/usr/local/share/khmer/scripts/normalize-by-median.py -k 20 -C 20 -x 1e8 reads.fa --savehash normC20k20.kh

!/usr/local/share/khmer/scripts/filter-abund.py normC20k20.kh reads.fa.keep

!/usr/local/share/khmer/scripts/normalize-by-median.py -k 20 -C 5 -x 1e8 reads.fa.keep.abundfilt

ls

!/usr/local/share/khmer/scripts/abundance-dist-single.py -x 1e8 -k 20 reads.fa.keep reads.fa.keep.hist
!/usr/local/share/khmer/scripts/abundance-dist-single.py -x 1e8 -k 20 reads.fa.keep.abundfilt reads.fa.ka.hist
!/usr/local/share/khmer/scripts/abundance-dist-single.py -x 1e8 -k 20 reads.fa.keep.abundfilt.keep reads.fa.kak.hist

dn1 = numpy.loadtxt('reads.fa.keep.hist')
abund = numpy.loadtxt('reads.fa.ka.hist')
dn2 = numpy.loadtxt('reads.fa.kak.hist')

plot(dn1[:,0], dn1[:,1], label='first round')
plot(abund[:,0], abund[:,1], label='trim errors')
plot(dn2[:,0], dn2[:,1], label='second round')
axis(ymax=2500)
legend()