workDir = '/home/nick/notebook/SIPSim/dev/bac_genome1210/'
SIPSimExe = '/home/nick/notebook/SIPSim/SIPSim'
genomeDir = '/var/seq_data/ncbi_db/genome/Feb2015/'

"""
Explanation: Description:

Simulating bacterial communities using all of the bacterial genomes available.

Workflow:

symlink all genomes
index all genomes

Setting variables
End of explanation
"""


import os,sys
import numpy as np
import pandas as pd
#from ggplot import *
import matplotlib.pyplot as plt

%load_ext rpy2.ipython
%matplotlib inline

%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)

if not os.path.isdir(workDir):
    os.mkdir(workDir)

genomeDir = os.path.join(workDir, 'genomes')
if not os.path.isdir(genomeDir):
    os.mkdir(genomeDir)

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


# using 'prokaryote.txt' file from the NCBI genome ftp site
!cd $genomeDir; \
    wget ftp://ftp.ncbi.nih.gov/genomes/GENOME_REPORTS/prokaryotes.txt

! printf "Number of lines in old (delta-GC dataset) version of prokaryotes.txt: "
! cd $genomeDir; \
    wc -l archive/deltaGC_dataset/prokaryotes.txt | perl -pe 's/ .+//'
! printf "Number of lines in new version of prokaryotes.txt: "
! cd $genomeDir; \
    wc -l prokaryotes.txt | perl -pe 's/ .+//'

os.chdir(genomeDir)

# loading table
proks = pd.read_csv('prokaryotes.txt', sep='\t')
print "Number of rows in full table: {}".format(proks.shape[0])

# filtering table
## just closed genomes (1 scaffold per chromosome)
### *Complete* in 'Status' field
proks = proks.loc[proks['Status'].isin(['Complete','Complete Genome'])]
print "Number of rows with a complete genome entry: {}".format(proks.shape[0])

"""
Explanation: Getting genomes from ncbi

Updating list of genomes from ncbi

Getting table of prokaryote genomes
End of explanation
"""


!cd $genomeDir; \
    seqDB_tools NCBIprokaryoteTableFilter \
    -t 30 prokaryotes.txt \
    > prokaryotes_filt.txt 2> NCBIprokaryoteTableFilter.log

!printf "Number of rows in filtered prokaryote table: "
!cd $genomeDir; wc -l prokaryotes_filt.txt | perl -pe 's/ .+//'

!wc -l NCBIprokaryoteTableFilter.log
!cd $genomeDir; \
    head -n 50 NCBIprokaryoteTableFilter.log

"""
Explanation: Filtering prokaryotes file

Only using 'complete' genomes.
Using the most recent version of each genome.
Adding taxonomy for each entry.
End of explanation
"""


!cd $genomeDir; \
    perl -ne 'print if /Rhizobium etli CFN |Lactobacillus gasseri ADH/' prokaryotes_filt.txt

!cd $genomeDir; \
    perl -ni -e 'print unless /Rhizobium etli CFN |Lactobacillus gasseri ADH/' prokaryotes_filt.txt

!cd $genomeDir; \
    perl -ne 'print if /Rhizobium etli CFN |Lactobacillus gasseri ADH/' prokaryotes_filt.txt

"""
Explanation: Manual filtering

Rhizobium_etli_CFN_42
no genome


Lactobacillus_gasseri_ADH
phage genome
End of explanation
"""


%%bash -s "$genomeDir"

cd $1

NY_misc_perl tblRandomByField -c 32 \
    -header <(perl -ne '@l=split /\t/; print unless \$l[9] eq "-"' prokaryotes_filt.txt) \
    > prokaryotes_filt_rand.txt

!cd $genomeDir; \
    wc -l prokaryotes_filt_rand.txt

!cd $genomeDir; \
    head -n 2 prokaryotes_filt_rand.txt

"""
Explanation: Random selection of 1 from each species

filtered out all genomes w/out accession numbers
End of explanation
"""


%%bash -s "$genomeDir"

perl -ne '@l=split /\t/; print if $l[25] =~/bacteria|superkingdom/i && $l[31] !~ /NA/' \
    prokaryotes_filt_rand.txt > prok-bac_filt_rand.txt
printf "Number of bacterial genomes: "
wc -l prok-bac_filt_rand.txt


perl -ne '@l=split /\t/; print if $l[25] =~/archaea|superkingdom/i && $l[31] !~ /NA/' \
    prokaryotes_filt_rand.txt > prok-arc_filt_rand.txt
printf "Number of archaeal genomes: "
wc -l prok-arc_filt_rand.txt

"""
Explanation: Parsing by domain

Also, must have a genus taxonomic classification
End of explanation
"""


%%bash -s "$genomeDir"

cd $1

seqDB_tools NCBIprokaryoteTableFilter prok-bac_filt_rand.txt -t 30 -w -d prok-bac-genomes

!cd $genomeDir; \
    seqDB_tools NCBIprokaryoteTableFilter prok-arc_filt_rand.txt -t 30 -w -d prok-arc-genomes

"""
Explanation: Downloading genomes
End of explanation
"""


bacGenomeDir = os.path.join(genomeDir, 'prok-bac-genomes')

genomeLocalDir = os.path.join(workDir, 'genomes')

if not os.path.isdir(genomeLocalDir):
    os.mkdir(genomeLocalDir)

%%bash -s "$genomeLocalDir" "$bacGenomeDir"
# renaming sequences
cd $1

# making sure each sequence is unique
find $2 -name "*fasta" |\
    perl -pe 's/.+\///' | \
    xargs -P 24 -I % bash -c \
    "/var/seq_data/ncbi_db/genome/genome_rename.pl < $2/% > %"

!printf "Number of bacterial genomes: "
!find $genomeLocalDir -name "*fasta" | wc -l

# list of all genomes files and their associated names
!cd $genomeLocalDir; \
    find . -name "*fasta" | \
    perl -pe 's/.+\///' | \
    perl -pe 's/(.+)(\.[^.]+)/\$1\t\$1\$2/' > genomes_all_list.txt

"""
Explanation: Copying and renaming genomes

Need to have each scaffold/chromosome name as 1 string be MFEprimer is terrible
End of explanation
"""


# determining the largest genome
!cd $genomeLocalDir; \
    find . -name "*fasta" | \
    xargs -I % seq_tools fasta_info % --fn --tn --tl --tgc --sn --sl --sgc \
    > genome_info.txt

"""
Explanation: Checking genome info

Making sure there are no outliers that need to be filtered out
End of explanation
"""


%%R -i genomeLocalDir

setwd(genomeLocalDir)

tbl = read.delim('genome_info.txt', sep='\t', header=FALSE)
colnames(tbl) = c('file_name', 'total_seqs', 'total_length', 'total_GC', 'seq_name', 'seq_length', 'seq_GC')

%%R -w 500 -h 350
# genome lengths
ggplot(tbl, aes(total_length)) +
    geom_histogram(binwidth=10000) +
    labs(x='Genome length') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )

%%R
summary(tbl$total_length)

%%R -w 500 -h 350
# genome GC
ggplot(tbl, aes(total_GC)) +
    geom_histogram(binwidth=1) +
    labs(x='Genome G+C content') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )

%%R
summary(tbl$total_GC)

%%R -w 500 -h 350
# genome GC
ggplot(tbl, aes(total_GC)) +
    geom_density(fill='red', alpha=0.3) +
    scale_x_continuous(limits=c(10,90)) +
    labs(x='Genome G+C content') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )

%%R -w 500 -h 350
# sequences per genome
ggplot(tbl, aes(total_seqs)) +
    geom_histogram(binwidth=1) +
    labs(x='Chromosomes/contigs') +
    theme_bw() +
    theme(
        text = element_text(size=16)
    )

"""
Explanation: loading table and plotting
End of explanation
"""


%%R
tbl %>%
    filter(total_GC < 28) %>%
    arrange(total_GC) %>%
    mutate(taxon = gsub('.fasta', '', file_name)) %>%
    select(total_GC, taxon)

"""
Explanation: Genomes with low G+C
End of explanation
"""


!cd $workDir; \
    $SIPSimExe indexGenomes genomes/genomes_all_list.txt \
    --fp ./genomes/ --np 30 > genomes/index_log.txt

"""
Explanation: Indexing genomes

Needed for in-silico PCR
End of explanation
"""


