# directory where you want the spacer blasting to be done
## CHANGE THIS!
workDir = "/home/nyoungb2/t/CLdb_Methanosarcina/"

"""
Explanation: Description:

This notebook goes through the analysis of blasting spacers vs a database containing potential protospacers.
Any blast datbase can be used.
A script is provided to easily blast all of the genomes in CLdb.

Before running this notebook:

run the Setup notebook

User-defined variables
End of explanation
"""


import os
from IPython.display import FileLinks
%load_ext rpy2.ipython

blastDir = os.path.join(workDir, 'arrayBlast')
if not os.path.isdir(blastDir):
    os.makedirs(blastDir)

"""
Explanation: Init
End of explanation
"""


!cd $blastDir; \
    CLdb -- array2fasta -cluster -cut 1 > spacers_cut1.fna

## getting number of spacer sequences
!printf 'Number of unique spacers:\t'
!cd $blastDir; \
    grep -c ">" spacers_cut1.fna; \
    echo; \
    head -n 6 spacers_cut1.fna

"""
Explanation: Blast vs all E. coli genomes in the CLdb
Spacer blast
Selecting spacer sequences
End of explanation
"""


!CLdb -- arrayBlast -h

# listing all arrayBlast subcommands
!CLdb -- arrayBlast --list

"""
Explanation: Blast

The array blast commands are sub-sub commands.
SO, that's why when you type: 'CLdb -- arrayBlast -h', you get then following:
End of explanation
"""


# spacer blast
## Really just a thin wrapper around blastn that will blast the spacers against all genomes in CLdb
### If you want to blast some else, just run the blast yourself (use 'perldoc CLdb_arrayBlast.pl' for more info)

!cd $blastDir; \
    CLdb -- arrayBlast -- run \
    -query spacers_cut1.fna \
    -fork 10 \
    > spacers_cut1_blastn.xml
    
# NOTE: using 10 CPUs    

"""
Explanation: Sub-sub-command help:
CLdb arrayBlast -- run -- -h
End of explanation
"""


# converting .xml to .srl
## This retains all of the info in the xml, but produces a smaller file that is faster to load and write.
## For large blast runs (xml files > 100 Mb), this will take a little while

!cd $blastDir; \
    CLdb -- arrayBlast -- xml2srl \
    spacers_cut1_blastn.xml \
    > spacers_cut1_blastn.srl

# checking output file
!cd $blastDir; \
    ls -thlc spacers_cut1_blastn.*

"""
Explanation: Note:

Using spacer clustering cutoff of 1, so just 'unique' spacer sequences
End of explanation
"""


# selecting DRs, so we can blast them.
## This is needed blast filter out spacers that hit CRISPR arrays

!cd $blastDir; \
    CLdb -- array2fasta \
    -cluster -cut 1 -r \
    > DRs_cut1.fna

## getting number of spacer sequences
!printf 'Number of unique DRs:\t'
!cd $blastDir; \
    grep -c ">" DRs_cut1.fna; \
    echo; \
    head -n 6 DRs_cut1.fna

# DR blast, just like spacer blast

!cd $blastDir; \
    CLdb -- arrayBlast -- run \
    -query DRs_cut1.fna \
    -fork 10 \
    > DRs_cut1_blastn.xml

# converting .xml to .srl
## This retains all of the info in the xml, but produces a smaller file that is faster to load and write.
## For large blast runs (xml files > 100 Mb), this will take a little while

!cd $blastDir; \
    CLdb -- arrayBlast -- xml2srl \
    DRs_cut1_blastn.xml \
    > DRs_cut1_blastn.srl


# checking output file
!cd $blastDir; \
    ls -thlc DRs_cut1_blastn.*

"""
Explanation: Note:

You could pipe the output from CLdb -- arrayBlast -- run directly into CLdb -- arrayBlast -- xml2srl
Example: CLdb -- arrayBlast -- run -query spacers_cut1.fna | CLdb -- arrayBlast -- xml2srl

DR Blast

BLASTing DR sequences against the same subjects.
This is needed to filter out spacer blast hits to CRISPR arrays.
End of explanation
"""


# filtering spacer hits to arrays
## This is based on the premise that DR hits adjacent to spacer hits signify that the spacer is hitting a CRISPR array

!cd $blastDir; \
    CLdb -- arrayBlast -- filter_arrayHits \
    spacers_cut1_blastn.srl \
    DRs_cut1_blastn.srl \
    > spacers_cut1_blastn_filt.srl
    
    
# checking output file
!cd $blastDir; \
    echo; \
    ls -thlc spacers_cut1_blastn_filt.srl

"""
Explanation: Filtering out spacer hits to CRISPRs

Based on DR hits falling adjacent to spacer hits
End of explanation
"""


!cd $blastDir; \
    CLdb -- arrayBlast -- srl2csv \
    spacers_cut1_blastn_filt.srl \
    > spacers_cut1_blastn_filt.txt
    
# assessing table of blast hits    
!cd $blastDir; \
    head spacers_cut1_blastn_filt.txt

#!printf "Number of blast hits: "
nhits = !cd $blastDir; wc -l spacers_cut1_blastn_filt.txt
nhits = nhits[0].split(' ')[0]
print "Number of blast hits: {}".format(int(nhits) - 1)

"""
Explanation: Converting blast .srl to .csv
End of explanation
"""


!cd $blastDir; \
    CLdb -- arrayBlast -- add_crRNA \
    < spacers_cut1_blastn_filt.srl \
    > spacers_cut1_blastn_filt_crDNA.srl
   
## checking output
!cd $blastDir; echo; \
    ls -thlc spacers_cut1_blastn_filt_crDNA.srl

"""
Explanation: Conclusions

We have some BLAST hits!
Let's add some information to the blast results

adding crRNA info

Adding crRNA region to .srl (actually the DNA template, refering to it as 'crDNA')
How much of the adjacent DR sequences are included in the crDNA is determined by the user (default: 10bp on either side)
End of explanation
"""


!cd $blastDir; \
    CLdb -- arrayBlast -- add_proto -workers 24 \
    < spacers_cut1_blastn_filt_crDNA.srl \
    > spacers_cut1_blastn_filt_crDNA_proto.srl   
    
## checking output
!cd $blastDir; echo; \
    ls -thlc spacers_cut1_blastn_filt_crDNA_proto.srl   

"""
Explanation: Adding protospacer info

Similar to adding crRNA info, but getting information from the blast database(s)
End of explanation
"""


# aligning crRNA and protospacer

!cd $blastDir; \
    CLdb -- arrayBlast -- align_proto \
    < spacers_cut1_blastn_filt_crDNA_proto.srl \
    > spacers_cut1_blastn_filt_crDNA_proto_aln.srl
    
## checking output
!cd $blastDir; echo; \
    ls -thlc spacers_cut1_blastn_filt_crDNA_proto_aln.srl

"""
Explanation: Aligning crRNA & protospacer

Alignment added to the .srl file.
These are individual alignments between the spacer and protospacer (using clustalw).
End of explanation
"""


# getting alignment of protospacer and crDNA and writing as a fasta

!cd $blastDir; \
    CLdb -- arrayBlast -- get_align \
    -outfmt taxon_name,subtype,evalue \
    < spacers_cut1_blastn_filt_crDNA_proto_aln.srl \
    > blastn_crDNA-proto_aln.fna
    

## checking output
!cd $blastDir; echo; \
    ls -thlc blastn_crDNA-proto_aln.fna

# let's take a look at the sequences 
!cd $blastDir; \
    head -n 10 blastn_crDNA-proto_aln.fna

"""
Explanation: Getting alignments (fasta)

You can parse out the alignments you want (parse by subtype, taxon_name, etc.)
OR you can parse them after writing (add subtype, taxon_name, e-value, etc. to the sequence name)
The alignments will be oriented to the crDNA (5'-3' for the crRNA)
See CLdb -- arrayBlast --perldoc -- get_align for info on what -outfmt does
End of explanation
"""


# file/directory
inFile = os.path.join(blastDir, 'blastn_crDNA-proto_aln.fna')
outDir = os.path.join(blastDir, 'blastn_crDNA-proto_aln_split')
if not os.path.isdir(outDir):
    os.makedirs(outDir)
outPre = os.path.join(outDir, 'split')    

# spliting alignment
!split -l 4 -d -a 4 $inFile $outPre

# checking out files
!printf 'Number of files produced: '
!find $outDir -name "split*" | wc -l

"""
Explanation: The alignment is really just for each crDNA-protospacer pair (every 4 lines)
First, the crRNA (sequence name starts with "crRNA")
Second, the protospacer (sequence name start with "proto")
Example: 

~~~

crDNA
ATAGACA
proto
ATAGACA
~~~



Note: the lower case letters are sequence adjacent to the actual protospacer


If you wanted each crRNA-protospacer pair individually, you could use the split bash command like this:

WARNING: this produces a lot of files, so we will make them in a new directory
End of explanation
"""


!cd $blastDir; \
    CLdb -- arrayBlast -- get_PAM -PAM 1 4 -f - \
    < blastn_crDNA-proto_aln.fna \
    > blastn_crDNA-proto_aln_PAM.fna


## checking output
!cd $blastDir; echo; \
    head blastn_crDNA-proto_aln_PAM.fna

"""
Explanation: Getting PAMs

Looking at the alignments can be very useful, but it can help to extract what we specifically want and summarize.
Let's pull out the potential PAM regions from the alignment 
We are going to assume that the PAM is:
4bp long 
1bp away from the 5' end of the protospacer ('right' side of the alignment)
End of explanation
"""


!cd $blastDir; \
    egrep -v "^>" blastn_crDNA-proto_aln_PAM.fna | \
    sort | uniq -c | sort -k 1 -n -r 

"""
Explanation: Summary table of each unique PAM sequence
End of explanation
"""


!cd $blastDir; \
    CLdb -- arrayBlast -- get_SEEDstats \
    --SEED -8 -1 \
    -prefix blastn_crDNA-proto_aln \
    < blastn_crDNA-proto_aln.fna 
    
## checking output files
!cd $blastDir; echo; \
    ls -thlc blastn_crDNA-proto_aln_SEEDstats*.txt

# checking files
!cd $blastDir; \
    head -n 6 blastn_crDNA-proto_aln_SEEDstats-byPos.txt
    
print '-' * 50


!cd $blastDir; \
    head -n 6 blastn_crDNA-proto_aln_SEEDstats-summary.txt

"""
Explanation: Spacer (crDNA) - protospacer mismatch distribution

Let's look at potential spacer targeting by assessing mismatches of the crDNAs and protospacers.
This script will make 2 summary files of mismatches for each crDNA-protospacer alignment
Given that spacers can vary in length, the alignment is relative to a user-defined 'SEED' region.
I'm going to use the default SEED region: '-SEED -8 -1'
This SEED region covers the last 8 bp of the alignment
End of explanation
"""


%%R
library(ggplot2)
library(gridExtra)

%%R -i blastDir -w 800

# loading table
infile = file.path(blastDir, 'blastn_crDNA-proto_aln_SEEDstats-byPos.txt')

tbl = read.delim(infile)
tbl = tbl[tbl$region != 'protospacer', ]

# mismatches and counts, & normalized mismatches
p.mis = ggplot(tbl, aes(pos_rel_SEED, mismatches, group=region, fill=region)) +
    geom_bar(stat='identity', position='dodge') +
    labs(y='Number\nof mismatches', x='Alignment position relative to SEED start')
p.count = ggplot(tbl, aes(pos_rel_SEED, count, group=region, fill=region)) +  
    geom_bar(stat='identity', position='dodge') +
    labs(y='Number of alignments\nwith nucleotide\nat position', x='Alignment position relative to SEED start')
p.norm = ggplot(tbl, aes(pos_rel_SEED, mismatch_norm, group=region, fill=region)) +
  geom_bar(stat='identity', position='dodge') +
    labs(y='Mismatchs /\nalignments', x='Alignment position relative to SEED start')

grid.arrange(p.count, p.mis, p.norm, ncol=1)

"""
Explanation: Plotting mismatches
End of explanation
"""


