#!/usr/bin/env python
# encoding: utf-8

from pipe import *

# check if user ran make    
if not os.path.exists('bin'):
    ohno('bin/ doesn\'t exist. Did you run make?')

cutoff = 70
out    = 'out'
left_mates = 'data/s71.qseq'
right_mates = 'data/s73.qseq'

# Creates a simple function to prepend the output directory
# to the directory/filename you specify
d = get_outdir(out)

ohai('running pipeline!')

## MAKE DIRECTORIES
[ run('mkdir -p %s' % i) for i in [
    out,
    d('orfs'),
    d('anno'),
    d('refs'),
    d('tables') ] ]

trim_pairs(
    left_mates  = left_mates,
    right_mates = right_mates,
    out_left    = d('singletons_left.fastq'),
    out_right   = d('singletons_right.fastq'),
    out         = d('reads_trimmed.fastq'),
    cutoff      = cutoff
)

# ASSEMBLE W/ VELVET
kmers = {
     31: d('contigs_31'),
     51: d('contigs_51'),
     71: d('contigs_71')
}

for kmer in kmers:
    velvet(
        reads = [
            ('fastq', 'shortPaired', d('reads_trimmed.fastq')),
            ('fastq', 'short', d('singletons_left.fastq')),
            ('fastq', 'short', d('singletons_right.fastq'))
        ],
        outdir = kmers[kmer],
        kmer   = kmer
    )

# run final assembly
# type of assembly depends on whether or not contigs from first
# 3 assemblies are actually long.

velvet(
    reads    = [('fasta', 'long', d('contigs_%s/contigs.fa' % k)) for k in kmers],
    outdir   = d('final_contigs'),
    kmer     = 51
)

# PREDICT OPEN READING FRAMES
prodigal(
    input  = d('final_contigs/contigs.fa'),
    out    = d('orfs/predicted_orfs') # prefix*
)

# create table connecting seed and subsystems
# seed_sequence_number -> system;subsystem;subsubsystem;enzyme
# use this later to make functions tables
prepare_seed(
    seed = 'db/seed.fasta',
    peg  = 'db/subsystems2peg',
    role = 'db/subsystems2role',
    out  = 'db/seed_ss.txt'
)
 
## IDENTIFY ORFS WITH PHMMER
phmmer( 
    query = d('orfs/predicted_orfs.faa'),
    db    = 'db/seed.fasta',
    out   = d('anno/proteins.txt')
)

# flatten phmmer file (we only need top hit)
run('misc/flatten_phmmer.py out/anno/proteins.txt.table > out/anno/proteins_flat.txt',
    generates='out/anno/proteins_flat.txt')

## GET ORF COVERAGE

# reference assemble
reference_assemble( # clc specific
    reference = d('orfs/predicted_orfs.fna'),
    out       = d('refs/reads_versus_orfs.txt'),
    query     = [
        ('unpaired', d('reads_trimmed.fastq')),
        ('unpaired', d('singletons_left.fastq')),
        ('unpaired', d('singletons_right.fastq')) ],
)

# make coverage table
make_coverage_table( # clc specific
    reads     = 'reads.fasta',
    reference = d('orfs/predicted_orfs.fna'),
    clc_table = d('refs/reads_versus_orfs.txt'),
    phmmer    = d('anno/proteins_flat.txt'),
    out       = d('tables/orfs_coverage.txt'),
)

# make subsystems table from coverage table
make_subsystems_table(
    subsnames      = 'db/seed_ss.txt',
    coverage_table = d('tables/orfs_coverage.txt'),
    out            = d('tables/subsystems_coverage.txt'),
    reads_type     = 'fastq',
    reads          = [
        d('reads_trimmed.fastq'),
        d('singletons_left.fastq'),
        d('singletons_right.fastq') ]
)

# GET OTU COVERAGE
reference_assemble(
    reference = 'db/taxcollector.fa',
    out       = d('refs/reads_vs_taxcollector.txt'),
    query     = [
        ('unpaired', d('reads_trimmed.fastq')),
        ('unpaired', d('singletons_left.fastq')),
        ('unpaired', d('singletons_right.fastq')) ],
)

make_coverage_table(
    reference    = d('db/taxcollector.fa'),
    clc_table    = d('refs/16s_table.txt'),
    reads_format = 'fastq',
    reads        = [
        d('reads_trimmed.fastq'),
        d('singletons_left.fastq'),
        d('singletons_right.fastq') ]
)

# TODO make OTU abundancy matrices
# NOTE I should just make those scripts part of
# taxcollector :\

# estimate average genome size and use to normalize?

# make_subsystems_table()
# make_phylogeny_table()

@@ -2392,37 +2392,24 @@
  = %5B%0A       
- ('unpaired',
  d('reads_tr
@@ -2413,33 +2413,32 @@
 _trimmed.fastq')
-)
 ,%0A        ('unpa
@@ -2426,37 +2426,24 @@
 q'),%0A       
- ('unpaired',
  d('singleto
@@ -2449,33 +2449,32 @@
 ons_left.fastq')
-)
 ,%0A        ('unpa
@@ -2462,37 +2462,24 @@
 q'),%0A       
- ('unpaired',
  d('singleto
@@ -2490,25 +2490,24 @@
 ight.fastq')
-)
  %5D,%0A)%0A%0A# mak
@@ -2563,39 +2563,8 @@
 fic%0A
-    reads     = 'reads.fasta',%0A
     
@@ -2999,139 +2999,8 @@
 q',%0A
-    reads          = %5B%0A        d('reads_trimmed.fastq'),%0A        d('singletons_left.fastq'),%0A        d('singletons_right.fastq') %5D%0A
 )%0A%0A#
@@ -3149,37 +3149,24 @@
  = %5B%0A       
- ('unpaired',
  d('reads_tr
@@ -3170,33 +3170,32 @@
 _trimmed.fastq')
-)
 ,%0A        ('unpa
@@ -3183,37 +3183,24 @@
 q'),%0A       
- ('unpaired',
  d('singleto
@@ -3206,33 +3206,32 @@
 ons_left.fastq')
-)
 ,%0A        ('unpa
@@ -3227,21 +3227,8 @@
     
- ('unpaired',
  d('
@@ -3251,17 +3251,16 @@
 .fastq')
-)
  %5D,%0A)%0A%0Am
@@ -3259,24 +3259,28 @@
  %5D,%0A)%0A%0Amake_
+otu_
 coverage_tab
@@ -3403,343 +3403,6 @@
 q',%0A
-    reads        = %5B%0A        d('reads_trimmed.fastq'),%0A        d('singletons_left.fastq'),%0A        d('singletons_right.fastq') %5D%0A)%0A%0A# TODO make OTU abundancy matrices%0A# NOTE I should just make those scripts part of%0A# taxcollector :%5C%0A%0A# estimate average genome size and use to normalize?%0A%0A# make_subsystems_table()%0A# make_phylogeny_table(
 )%0A

make coverage table and make otu coverage table for illumina