# path to raw files
## CHANGE THIS!
rawFileDir = "~/perl/projects/CLdb/data/Methanosarcina/"
# directory where the CLdb database will be created
## CHANGE THIS!
workDir = "~/t/CLdb_Methanosarcina/"

# viewing file links
import os
import zipfile
import csv
from IPython.display import FileLinks
# pretty viewing of tables
## get from: http://epmoyer.github.io/ipy_table/
from ipy_table import *

rawFileDir = os.path.expanduser(rawFileDir)
workDir = os.path.expanduser(workDir)

"""
Explanation: This notebook describes the setup of CLdb with a set of Methanosarcina genomes.

Dataset Notes

The CRISPR systems were classified according to Vestergaard G, Garrett RA, Shah SA. (2014). CRISPR adaptive immune systems of Archaea. RNA Biol 11: 156–167


General Notes

This notebook assumed that you have CLdb in your PATH
End of explanation
"""


FileLinks(rawFileDir)

"""
Explanation: The required files are in '../ecoli_raw/':

a loci table
array files
genome nucleotide sequences
genbank (preferred) or fasta format

Let's look at the provided files for this example:
End of explanation
"""


!CLdb -h

"""
Explanation: Checking that CLdb is installed in PATH
End of explanation
"""


# this makes the working directory
if not os.path.isdir(workDir):
    os.makedirs(workDir)

# unarchiving files in the raw folder over to the newly made working folder
files = ['array.zip','loci.zip', 'accessions.txt.zip']
files = [os.path.join(rawFileDir, x) for x in files]
for f in files:
    if not os.path.isfile(f):
        raise IOError, 'Cannot find file: {}'.format(f)
    else:
        zip = zipfile.ZipFile(f)
        zip.extractall(path=workDir)         

print 'unzipped raw files:'        
FileLinks(workDir)      

"""
Explanation: Setting up the CLdb directory
End of explanation
"""


# making genbank directory
genbankDir = os.path.join(workDir, 'genbank')
if not os.path.isdir(genbankDir):
    os.makedirs(genbankDir)    

# downloading genomes
!cd $genbankDir; \
    CLdb -- accession-GI2fastaGenome -format genbank -fork 9 < ../accessions.txt
    
# checking files
!cd $genbankDir; \
    ls -thlc *.gbk

"""
Explanation: Downloading the genome genbank files. Using the 'GIs.txt' file

GIs.txt is just a list of GIs and taxon names.
End of explanation
"""


!CLdb -- makeDB -h

"""
Explanation: Creating/loading CLdb of E. coli CRISPR data
End of explanation
"""


!cd $workDir; \
    CLdb -- makeDB -r -drop
    
CLdbFile = os.path.join(workDir, 'CLdb.sqlite')
print 'CLdb file location: {}'.format(CLdbFile)

"""
Explanation: Making CLdb sqlite file
End of explanation
"""


s = 'DATABASE = ' + CLdbFile
configFile = os.path.join(os.path.expanduser('~'), '.CLdb')

with open(configFile, 'wb') as outFH:
    outFH.write(s)
    
print 'Config file written: {}'.format(configFile)

# checking that the config is set
!CLdb --config-params

"""
Explanation: Setting up CLdb config

This way, the CLdb script will know where the CLdb database is located.
Otherwise, you would have to keep telling the CLdb script where the database is.
End of explanation
"""


lociFile = os.path.join(workDir, 'loci', 'loci.txt')

# reading in file
tbl = []
with open(lociFile, 'rb') as f:
    reader = csv.reader(f, delimiter='\t')
    for row in reader:
        tbl.append(row)

# making table
make_table(tbl)
apply_theme('basic')

"""
Explanation: Loading loci

The next step is loading the loci table.
This table contains the user-provided info on each CRISPR-CAS system in the genomes.
Let's look at the table before loading it in CLdb

Checking out the CRISPR loci table
End of explanation
"""


!CLdb -- loadLoci -h

!CLdb -- loadLoci < $lociFile

"""
Explanation: Notes on the loci table:
* As you can see, not all of the fields have values. Some are not required (e.g., 'fasta_file').
* You will get an error if you try to load a table with missing values in required fields.
* For a list of required columns, see the documentation for CLdb -- loadLoci -h.
Loading loci info into database
End of explanation
"""


# This is just a quick summary of the database 
## It should show 10 loci for the 'loci' rows
!CLdb -- summary

"""
Explanation: Notes on loading

A lot is going on here:
Various checks on the input files
Extracting the genome fasta sequence from each genbank file 
The genome fasta is required


Loading of the loci information into the sqlite database

Notes on the command

Why didn't I use the 'required' -database flag for CLdb -- loadLoci???
I didn't have to use the -database flag because it is provided via the .CLdb config file that was previously created.
End of explanation
"""


# an example array file (obtained from CRISPRFinder)
arrayFile = os.path.join(workDir, 'array', 'Methanosarcina_acetivorans_C2A_1.txt')
!head $arrayFile

"""
Explanation: The summary doesn't show anything for spacers, DRs, genes or leaders!
That's because we haven't loaded that info yet...
Loading CRISPR arrays

The next step is to load the CRISPR array tables.
These are tables in 'CRISPRFinder format' that have CRISPR array info.
Let's take a look at one of the array files before loading them all.
End of explanation
"""


# loading CRISPR array info
!CLdb -- loadArrays 

# This is just a quick summary of the database 
!CLdb -- summary

"""
Explanation: Note: the array file consists of 4 columns:

spacer start
spacer sequence
direct-repeat sequence
direct-repeat stop

All extra columns ignored!
End of explanation
"""


geneDir = os.path.join(workDir, 'genes')
if not os.path.isdir(geneDir):
    os.makedirs(geneDir)

!cd $geneDir; \
    CLdb -- getGenesInLoci 2> CAS.log > CAS.txt
    
# checking output    
!cd $geneDir; \
    head -n 5 CAS.log; \
    echo -----------; \
    tail -n 5 CAS.log; \
    echo -----------; \
    head -n 5 CAS.txt

# loading gene table into the database
!cd $geneDir; \
    CLdb -- loadGenes < CAS.txt 

"""
Explanation: Note: The output should show 75 spacer & 85 DR entries in the database
Loading CAS genes

Technically, all coding seuqences in the region specified in the loci table (CAS_start, CAS_end) will be loaded.
This requires 2 subcommands:
The 1st gets the gene info
The 2nd loads the info into CLdb
End of explanation
"""


!CLdb -- setSenseStrand 

"""
Explanation: Setting array sense strand

The strand that is transcribed needs to be defined in order to have the correct sequence for downstream analyses (e.g., blasting spacers and getting PAM regions)
The sense (reading) strand is defined by (order of precedence):
The leader region (if defined; in this case, no).
Array_start,Array_end in the loci table
The genome negative strand will be used if array_start > array_end
End of explanation
"""


!CLdb -- clusterArrayElements -s -r

"""
Explanation: Spacer and DR clustering

Clustering of spacer and/or DR sequences accomplishes:
A method of comparing within and between CRISPRs
A reducing redundancy for spacer and DR blasting
End of explanation
"""


# summary
!cd $workDir; \
    CLdb -- summary -name -subtype > summary.txt

# checking output
!cd $workDir; \
     cat summary.txt

"""
Explanation: Database summary
End of explanation
"""


