cd ~/Desktop/SSUsearch/

mkdir -p ./workdir

#check seqfile files to process in data directory (make sure you still remember the data directory)
!ls ./data/test/data

"""
Explanation: Set up working directory
End of explanation
"""


Seqfile='./data/test/data/2d.fa'

"""
Explanation: README
This part of pipeline search for the SSU rRNA gene fragments, classify them, and extract reads aligned specific region. It is also heavy lifting part of the whole pipeline (more cpu will help).
This part works with one seqfile a time. You just need to change the "Seqfile" and maybe other parameters in the two cells bellow.
To run commands, click "Cell" then "Run All". After it finishes, you will see "*** pipeline runs successsfully :)" at bottom of this pape.
If your computer has many processors, there are two ways to make use of the resource:


Set "Cpu" higher number.


make more copies of this notebook (click "File" then "Make a copy" in menu bar), so you can run the step on multiple files at the same time.


(Again we assume the "Seqfile" is quality trimmed.)
Here we will process one file at a time; set the "Seqfile" variable to the seqfile name to be be processed
First part of seqfile basename (separated by ".") will be the label of this sample, so named it properly.
e.g. for "/usr/local/notebooks/data/test/data/1c.fa", "1c" will the label of this sample.
End of explanation
"""


Cpu='2'   # number of maxixum threads for search and alignment
Hmm='./data/SSUsearch_db/Hmm.ssu.hmm'   # hmm model for ssu
Gene='ssu'
Script_dir='./scripts'
Gene_model_org='./data/SSUsearch_db/Gene_model_org.16s_ecoli_J01695.fasta'
Ali_template='./data/SSUsearch_db/Ali_template.silva_ssu.fasta'

Start='577'  #pick regions for de novo clustering
End='727'
Len_cutoff='100' # min length for reads picked for the region

Gene_tax='./data/SSUsearch_db/Gene_tax.silva_taxa_family.tax' # silva 108 ref
Gene_db='./data/SSUsearch_db/Gene_db.silva_108_rep_set.fasta'

Gene_tax_cc='./data/SSUsearch_db/Gene_tax_cc.greengene_97_otus.tax' # greengene 2012.10 ref for copy correction
Gene_db_cc='./data/SSUsearch_db/Gene_db_cc.greengene_97_otus.fasta'

# first part of file basename will the label of this sample
import os
Filename=os.path.basename(Seqfile)
Tag=Filename.split('.')[0]

import os

New_path = '{}:{}'.format('~/Desktop/SSUsearch/external_tools/bin/', os.environ['PATH'])
Hmm=os.path.abspath(Hmm)
Seqfile=os.path.abspath(Seqfile)
Script_dir=os.path.abspath(Script_dir)
Gene_model_org=os.path.abspath(Gene_model_org)
Ali_template=os.path.abspath(Ali_template)
Gene_tax=os.path.abspath(Gene_tax)
Gene_db=os.path.abspath(Gene_db)
Gene_tax_cc=os.path.abspath(Gene_tax_cc)
Gene_db_cc=os.path.abspath(Gene_db_cc)

os.environ.update(
    {'PATH':New_path,
     'Cpu':Cpu, 
     'Hmm':os.path.abspath(Hmm), 
     'Gene':Gene, 
     'Seqfile':os.path.abspath(Seqfile), 
     'Filename':Filename, 
     'Tag':Tag, 
     'Script_dir':os.path.abspath(Script_dir), 
     'Gene_model_org':os.path.abspath(Gene_model_org), 
     'Ali_template':os.path.abspath(Ali_template), 
     'Start':Start, 
     'End':End,
     'Len_cutoff':Len_cutoff,
     'Gene_tax':os.path.abspath(Gene_tax), 
     'Gene_db':os.path.abspath(Gene_db), 
     'Gene_tax_cc':os.path.abspath(Gene_tax_cc), 
     'Gene_db_cc':os.path.abspath(Gene_db_cc)})

!echo "*** make sure: parameters are right"
!echo "Seqfile: $Seqfile\nCpu: $Cpu\nFilename: $Filename\nTag: $Tag"

cd workdir

mkdir -p $Tag.ssu.out

### start hmmsearch

%%bash
echo "*** hmmsearch starting"
time hmmsearch --incE 10 --incdomE 10 --cpu $Cpu \
  --domtblout $Tag.ssu.out/$Tag.qc.$Gene.hmmdomtblout \
  -o /dev/null -A $Tag.ssu.out/$Tag.qc.$Gene.sto \
  $Hmm $Seqfile
echo "*** hmmsearch finished"

!python $Script_dir/get-seq-from-hmmout.py \
    $Tag.ssu.out/$Tag.qc.$Gene.hmmdomtblout \
    $Tag.ssu.out/$Tag.qc.$Gene.sto \
    $Tag.ssu.out/$Tag.qc.$Gene

"""
Explanation: Other parameters to set
End of explanation
"""


%%bash
echo "*** Starting mothur align"
cat  $Gene_model_org $Tag.ssu.out/$Tag.qc.$Gene > $Tag.ssu.out/$Tag.qc.$Gene.RFadded

# mothur does not allow tab between its flags, thus no indents here
time mothur "#align.seqs(candidate=$Tag.ssu.out/$Tag.qc.$Gene.RFadded, template=$Ali_template, threshold=0.5, flip=t, processors=$Cpu)"

rm -f mothur.*.logfile

"""
Explanation: Pass hits to mothur aligner
End of explanation
"""


!python $Script_dir/mothur-align-report-parser-cutoff.py \
    $Tag.ssu.out/$Tag.qc.$Gene.align.report \
    $Tag.ssu.out/$Tag.qc.$Gene.align \
    $Tag.ssu.out/$Tag.qc.$Gene.align.filter \
    0.5
    

!python $Script_dir/remove-gap.py $Tag.ssu.out/$Tag.qc.$Gene.align.filter $Tag.ssu.out/$Tag.qc.$Gene.align.filter.fa

"""
Explanation: Get aligned seqs that have > 50% matched to references
End of explanation
"""


!python $Script_dir/region-cut.py $Tag.ssu.out/$Tag.qc.$Gene.align.filter $Start $End $Len_cutoff

!mv $Tag.ssu.out/$Tag.qc.$Gene.align.filter."$Start"to"$End".cut.lenscreen $Tag.ssu.out/$Tag.forclust

"""
Explanation: Search is done here (the computational intensive part). Hooray!


\$Tag.ssu.out/\$Tag.qc.\$Gene.align.filter:      
    aligned SSU rRNA gene fragments


\$Tag.ssu.out/\$Tag.qc.\$Gene.align.filter.fa:   
    unaligned SSU rRNA gene fragments


Extract the reads mapped 150bp region in V4 (577-727 in E.coli SSU rRNA gene position) for unsupervised clustering
End of explanation
"""


%%bash
rm -f $Tag.ssu.out/$Tag.qc.$Gene.align.filter.silva_taxa_family*.taxonomy
mothur "#classify.seqs(fasta=$Tag.ssu.out/$Tag.qc.$Gene.align.filter.fa, template=$Gene_db, taxonomy=$Gene_tax, cutoff=50, processors=$Cpu)"
mv $Tag.ssu.out/$Tag.qc.$Gene.align.filter.silva_taxa_family*.taxonomy \
    $Tag.ssu.out/$Tag.qc.$Gene.align.filter.wang.silva.taxonomy

!python $Script_dir/count-taxon.py \
    $Tag.ssu.out/$Tag.qc.$Gene.align.filter.wang.silva.taxonomy \
    $Tag.ssu.out/$Tag.qc.$Gene.align.filter.wang.silva.taxonomy.count
!rm -f mothur.*.logfile

"""
Explanation: Classify SSU rRNA gene seqs using SILVA
End of explanation
"""


%%bash
rm -f $Tag.ssu.out/$Tag.qc.$Gene.align.filter.greengene_97_otus*.taxonomy
mothur "#classify.seqs(fasta=$Tag.ssu.out/$Tag.qc.$Gene.align.filter.fa, template=$Gene_db_cc, taxonomy=$Gene_tax_cc, cutoff=50, processors=$Cpu)"
mv $Tag.ssu.out/$Tag.qc.$Gene.align.filter.greengene_97_otus*.taxonomy \
    $Tag.ssu.out/$Tag.qc.$Gene.align.filter.wang.gg.taxonomy

!python $Script_dir/count-taxon.py \
    $Tag.ssu.out/$Tag.qc.$Gene.align.filter.wang.gg.taxonomy \
    $Tag.ssu.out/$Tag.qc.$Gene.align.filter.wang.gg.taxonomy.count
!rm -f mothur.*.logfile

# check the output directory
!ls $Tag.ssu.out

"""
Explanation: Classify SSU rRNA gene seqs with Greengene for copy  correction later
End of explanation
"""


!echo "*** pipeline runs successsfully :)"

"""
Explanation: This part of pipeline (working with one sequence file) finishes here. Next we will combine samples for community analysis (see unsupervised analysis).
Following are files useful for community analysis:

1c.577to727: aligned fasta file of seqs mapped to target region for de novo clustering
1c.qc.ssu.align.filter: aligned fasta file of all SSU rRNA gene fragments
1c.qc.ssu.align.filter.wang.gg.taxonomy: Greengene taxonomy (for copy correction)
1c.qc.ssu.align.filter.wang.silva.taxonomy: SILVA taxonomy
End of explanation
"""


