run:
        shell('prokka --force --outdir {params.outdir} --prefix {params.prefix} --locustag {params.str} --addgenes --increment 5 --centre NIOO-KNAW --genus {params.genus} --species {params.species} --str {params.str} --gcode 11 --cpus 20 --evalue 1e-03 --rfam {input.file}')

run:
    #Extract gene lengths
    shell("bioawk -c fastx '{{ print $name, length($seq) }}' < intermediate_files/combined_proteins/combined_proteins.fasta >intermediate_files/combined_proteins/length_genes.txt")

    #Run mmseq2 clustering to 0.95
    shell('mmseqs createdb {input} {params.mmseq_db}')
    shell('mmseqs cluster {params.mmseq_db} {params.mmseq_db_clu} {params.mmseq_temp} --min-seq-id 0.95 --cov-mode 0 -c 0.8')
    shell('mmseqs createtsv {params.mmseq_db} {params.mmseq_db} {params.mmseq_db_clu} {params.mmseq_tsv}')

    #Extract mmseq2 representatives
    shell('mmseqs createsubdb {params.mmseq_db_clu} {params.mmseq_db} {params.mmseq_rep}')
    shell('mmseqs convert2fasta {params.mmseq_rep} {params.mmseq_fasta}')


    #Run diamond
    shell('diamond makedb --in {params.mmseq_fasta} -d {params.diamond_db}')
    shell('diamond blastp -b 10 -p {THREADS} -c 1 -e 0.00001 --un {params.diamond_unaligned} -k 5000 -d {params.diamond_db} -q {params.mmseq_fasta} -o {params.diamond_results} --outfmt 6 qseqid sseqid pident bitscore')
    shell("awk -F'\t' '$3>20' {params.diamond_results} > {params.diamond_results_filtered}")
    shell("cut -f1,2,4 {params.diamond_results_filtered} > {params.diamond_results_clean}")

    ##Get unaaligned fasta headers
    shell('grep "^>" {params.diamond_unaligned}  > {params.diamond_unaligned_headers}')

    ### Run MCL
    shell('mcxload -abc {params.diamond_results_clean} --stream-mirror -o {params.mcl_data} -write-tab {params.mcl_tab}')
    shell('mcl {params.mcl_data} -I 3.0 -use-tab {params.mcl_tab}')
    shell('mv out.mcl_data.mci.I30 intermediate_files/clustering/')

    ##Generate matrix

    shell('Rscript src/MCL_merge.R {params.mmseq_tsv} {params.mcl_clusters} {params.diamond_unaligned_headers} {output.binary_table} {input} {output.fasta}')

run:
    shell('wget -P ./intermediate_files/PFAM/ ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam33.1/Pfam-A.hmm.gz')
    shell('gunzip intermediate_files/PFAM/Pfam-A.hmm.gz' )
    shell('hmmpress intermediate_files/PFAM/Pfam-A.hmm')
    shell('hmmsearch --tblout {output.pfam} --cpu 30 -E 1e-5 ./intermediate_files/PFAM/Pfam-A.hmm {input}')

run:
    shell('download_eggnog_data.py -y --data_dir {params.data}')
    shell('emapper.py -i {input} --cpu 25 -o intermediate_files/eggnog_annotation/eggnog_annotation --data_dir {params.data} --pident 30 --query_cover 50 --subject_cover 50 --report_orthologs')

shell:
    'set +u; source /opt/miniconda3/etc/profile.d/conda.sh; conda activate antismash_microlife; set -u; antismash --cb-general --cb-knownclusters --cb-subclusters --output-dir {params.out_dir} -c {params.threads} --asf --pfam2go --genefinding-tool prodigal --smcog-trees {input}'

library(stringr)
library(dplyr)
args = commandArgs(trailingOnly=TRUE)
#args[1] is the id2tags input
#args[2] is the prokka_headers.txt input
#args[3] is the output




id2tags<- read.delim2(args[1], header = F)
id2tags<- id2tags[,1:2]
colnames(id2tags)<- c('geneid', 'gene')

prokka_header <- read.delim(args[2], header = F)
prokka_header$prokkadescription <- prokka_header$V1
prokka_header$remove <- sapply(strsplit(prokka_header$V1," "), `[`, 1)

prokka_header$prokkadescription <- str_remove(prokka_header$prokkadescription, prokka_header$remove)
prokka_header$prokkadescription <- sub('.', '', prokka_header$prokkadescription)

prokka_header$remove <- sub('.', '', prokka_header$remove)

prokka_header <- prokka_header[,c(3,2)]
colnames(prokka_header)[1] <- 'gene'

##Merge id2tags with prokka headers
prokka_annotation <- merge(id2tags, prokka_header, by= 'gene', all.x = T)
prokka_annotation <- prokka_annotation[,2:3]
prokka_annotation<- prokka_annotation %>% distinct(geneid, .keep_all = TRUE)

write.table(prokka_annotation, args[3], row.names = F)

import argparse, os
import Bio
from Bio import SeqIO

def genbank2faa(infile, outfile):
    seq = SeqIO.parse(infile,'genbank')
    fastaseq = []
    if not os.path.exists(outfile):
        print( "Converting :%s"%(infile) )
        #partial CDS_counter
        partial = 0

        id_2_tags = "%s.tags"%(outfile)

        tagfile = open(id_2_tags,'w')

        for s in seq:

            features = s.features
            for f in features:
                if f.type == "source":
                    sourcefeature = f

            if 'organism' in sourcefeature.qualifiers:
                organism = sourcefeature.qualifiers['organism'][0]
                strain = ''

                if 'strain' in sourcefeature.qualifiers:
                    strain = sourcefeature.qualifiers['strain'][0]
                    strain = strain.replace(' ','').replace(':','').replace('|','')
                organism = organism + '_' + strain

            #else:
            #    organism = s.annotations['organism'].replace(' ','_').replace(':','_').replace('|','_')

            organism = organism.replace(' ','_').replace(':','_').replace('|','_')

            ltprefix = estimate_locustag_prefix(features)
            ltprefixnew = "%s_"%(ltprefix)


            moltype="C"
            for f in features:
                if f.type == "source":
                    if "plasmid" in f.qualifiers:
                        moltype="P"

            cdscount=0
            for f in features:
                if f.type == "CDS":
                    cdscount = cdscount + 1
                    if 'locus_tag' in f.qualifiers:
                        locus_tags = f.qualifiers['locus_tag']
                        locus_tag = locus_tags[-1]
                        if not '_' in locus_tag:
                           locus_tag = locus_tag.replace(ltprefix, ltprefixnew) 

                        if not isinstance(f.location.start, Bio.SeqFeature.ExactPosition):
                            # location is not exact.
                            partial = partial +1
                            locus_tag = "%s_%02d"%(locus_tag, partial)

                        prot = f.extract(s)
                        prot.seq = prot.seq.translate(table=11, to_stop=True)
                        prot.id = locus_tag

                        tagfile.write(locus_tag)
                        tagfile.write("\t")
                        tagfile.write(moltype)
                        tagfile.write("\t")
                        tagfile.write(organism)
                        tagfile.write("\t")
                        tagfile.write(prot.description)
                        tagfile.write("\n")

                        if 'product' in f.qualifiers:
                            prot.description = f.qualifiers['product'][0]
                        else:
                            prot.description = 'unknown'

                        fastaseq.append(prot)
            print( "Written :%s sequences from %s"%(len(fastaseq), s.description) )

        if partial>0:
            print("Warning sequence %s contains %d partial CDS sequences. This might inflate the number of clusters."%(infile, partial))

        SeqIO.write(fastaseq, outfile,'fasta')
    else:
        print( "Not Converting :%s"%(infile) )
        print( "%s exists"%(outfile) )


def genbank2ffn(infile, outfile):
    seq = SeqIO.parse(infile,'genbank')
    fastaseq = []
    print( "Converting :%s"%(infile) )
    #partial CDS_counter
    partial = 0


    for s in seq:
        features = s.features
        for f in features:
            if f.type == "gene":
                if 'locus_tag' in f.qualifiers:
                    locus_tags = f.qualifiers['locus_tag']
                    locus_tag = locus_tags[-1]
                    if not isinstance(f.location.start, Bio.SeqFeature.ExactPosition):
                        # location is not exact.
                        partial = partial +1
                        locus_tag = "%s_%02d"%(locus_tag, partial)

                    gene = f.extract(s)
                    gene.id = locus_tag
                    if 'gene' in f.qualifiers:
                        gene.description = f.qualifiers['gene'][0]
                    else:
                        gene.description = gene.id

                    fastaseq.append(gene)
        print( "Written :%s sequences from %s"%(len(fastaseq), s.description) )

    if partial>0:
        print("Warning sequence %s contains %d partial CDS sequences. This might inflate the number of clusters."%(infile, partial))

    SeqIO.write(fastaseq, outfile,'fasta')

def estimate_locustag_prefix(features):

    strings = [f.qualifiers['locus_tag'][-1] for f in features if 'locus_tag' in f.qualifiers]

    """ Find the longest string that is a prefix of all the strings.
    """
    if not strings:
        return ''
    prefix = strings[0]
    for s in strings:
        if len(s) < len(prefix):
            prefix = prefix[:len(s)]
        if not prefix:
            return ''
        for i in range(len(prefix)):
            if prefix[i] != s[i]:
                prefix = prefix[:i]
                break
    return prefix

if __name__=='__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument('-g', '--genbank', dest = 'genbank', help = 'genbank input file', required = True)
    parser.add_argument('-f', '--fasta', dest = 'fasta', help = 'fasta CDS output file', required = True)

    args = parser.parse_args()
    genbank2faa(args.genbank, args.fasta)

'R Script for Data Mining and Manipulation of Output Tables from BiG-SCAPE'

'Libraries to be Imported'
library(plyr)
library(dplyr)
library(stats)
library(tidyverse)
library(ggplot2)

#Set Working Directory
#setwd('/Results_I')
#ex:setwd("E:/NACTAR")

args = commandArgs(trailingOnly=TRUE)
"Import Data Tables from BiG-SCAPE"
#Clustering Tables
clustering_df_mix <- read.delim(args[1], header = TRUE) #clustering dataframe

#Network Tables
network_df <- read.delim(args[2], header = TRUE) #general network dataframe
network_df <- network_df [-c(4:11)] #removal of unwanted columns

#Annotation Tables
annotations_df_full<- read.delim(args[3], header = TRUE) #annotations dataframe

'Merge Files for Absence Presence Table'
annotations_merged_df_mix <- merge(annotations_df_full , clustering_df_mix , by.x= "BGC", by.y = "X.BGC.Name") #merge clustering and annotations
annotations_merged_df_mix  <- annotations_merged_df_mix [, c(1,8,6,5)]
#annotations_merged_df_mix  <- annotations_merged_df_mix [, c(1,8,4,5,6)] #removal of unwanted columns

colnames(annotations_merged_df_mix) <- c('BGC Name', 'GCF No', 'Organism', 'BGC Class')

'Extract Rows based on Logical Criteria (User Based Criteria) (Optional by User)'
#sum(network_df_mix$Raw.distance > 0.1)
network_df__distance_filtered <- network_df %>% filter(Raw.distance > 0.10)
network_df_bgc_filtered <- network_df[!grepl("BGC", network_df$Clustername.1),]



#Export Tables to CSV Format
# Write data to txt file: tab separated values
# sep = "\t"
write.table(network_df__distance_filtered, file = args[4], sep = "\t",
            row.names = FALSE, quote = FALSE)
# Write data to csv files:
# decimal point = "." and value separators = comma (",")

write.table(annotations_merged_df_mix, file = args[5], sep = '\t',
            row.names = FALSE, quote = FALSE)




###Obtain GCF annotations list
merged_annotations <- annotations_merged_df_mix
merged_annotations$GCF.No <- paste0( 'GCF',merged_annotations$GCF.No)
merged_annotations<- merged_annotations[,c(2,4)]
merged_annotations <- distinct(merged_annotations)

write.csv(merged_annotations, args[6], row.names = F)

'Python Script used to create an absence presence table of GCFs per Strain from BiG-SCAPE file'

#Import Libraries
import pandas as pd
import sys
import os
from string import digits

#Importing DataFrame
df1 = pd.read_csv(sys.argv[1], sep="\t")
df1.columns = ['BGC Name', 'GCF No', 'Organism', 'BGC Class']
#print(df1.isnull())
#print(df1.head())

'Method for Creating Absene Presence Table'
#Dataframe Copies
df2 = df1.copy() #Copy Original DataFrame
df2.columns = ['BGC Name', 'GCF No', 'Organism', 'BGC Class']
#DataFrame Manipulation
df2['GCF No'] = df1['BGC Class'].astype(str) + 'GCF' + df2['GCF No'].astype(str)
df1[['BGC','BGC2']] = df1['BGC Name'].str.split('_',n=1,expand=True)
df2['Genome'] = df1['Organism'].astype(str) + df1['BGC Name'].astype(str)
df2.set_index('GCF No')
df2 = df2[['GCF No','Genome']]
df2['Genome'] = df2['Genome'].str.replace(' ', '_').str.replace('Unclassified','_').str.replace('__','_').str.replace('.','')
path = os.path.join(os.path.abspath(os.getcwd()), 'intermediate_files/BiG-SCAPE')
df2.to_csv(os.path.join(path,r'abs_pres_table.csv'))

library(reshape2)
library(tidyverse)
library(stringr)
library(dplyr)
library(readr)
library(Biostrings)
library(data.table)
library(parallel)

cl <- makeCluster(30)
args = commandArgs(trailingOnly=TRUE)

tsv = data.table(read_table(args[1], col_names = F))
colnames(tsv) <- c('V1', 'V2')

n_clusters_mmseq = length(unique(tsv$V1))


tsv2 <- tsv[, .( V2 = str_c(V2, collapse = ',') ), by = V1]

#write.table(tsv2,'diamond/mmseq2_clusters.txt', row.names = F, quote = F)


'Open MCL output, collapse it and clean it\n'

mcl <- data.table(read_csv( args[2], col_names = F))

#Load unaligned headers 
unalin <- read.delim2(args[3], header = F)
unalin <- data.frame(X1 = str_remove_all(unalin$V1, '>'))

mcl <- dplyr::bind_rows(mcl, unalin)


#Name clusters

mcl$clusters <- seq(nrow(mcl))
mcl$clusters <- sprintf("cluster_%06d", mcl$clusters)


mcl_melt <- mcl[, c(X1 = strsplit(X1, "\t")), by = clusters]

'Join mmseq2 clusters with MCL clusters\n'

M <- merge(mcl_melt, tsv2, by.x= 'X1', by.y = 'V1')
M <- M[,c(2,3)]
colnames(M)[2] <- 'descriptions'
M <- data.table(M)
M <- M[, .( descriptions = str_c(descriptions, collapse = ',') ), by = clusters]
M <- M[order(clusters)]
M$gene <- sapply(strsplit(M$descriptions,","), `[`, 1)

#write.table(M, 'clusters_info.txt', row.names = F, sep = '\t', quote = F)

'Create 1/0 matrix\n' 

bacteria_name <- function(string){
  b <- stringr::str_split(string, ',')[[1]]
  c <- sapply(stringr::str_split(b,"[|]"), `[`, 1)
  d <- stringr::str_c(c, collapse = ",")
  return(d)
}

MEGAMATRIX <- M
MEGAMATRIX$descriptions <- parSapply(cl, MEGAMATRIX$descriptions , FUN = bacteria_name)
MEGAMATRIX$gene <- NULL



tbl <- table(splitstackshape:::cSplit(MEGAMATRIX, "descriptions", sep = ",", direction = "long"))
tbl <- data.table(as.data.frame(tbl))
tbl <-dcast.data.table(tbl, clusters ~ descriptions, value.var = 'Freq')

#lista <- as.list(MEGAMATRIX$descriptions)
#binary_matrix <- as.data.frame((splitstackshape:::charMat(listOfValues = lista, fill = 0L)))
#binary_matrix <- cbind(clusters = MEGAMATRIX$clusters, binary_matrix)
binary_matrix <- merge(tbl, M, by = 'clusters')



####Extract sequence of representative proteins
#system("bioawk -c fastx '{ print $name, length($seq) }' < intermediate_files/combined_proteins/combined_proteins.fasta >intermediate_files/combined_proteins/length_genes.txt")
gene_lengths <- read.table('intermediate_files/combined_proteins/length_genes.txt', header = F)

##Function to extract longest gene

extract_long_genes <- function(character, gene_length){
  a <- as.character(stringr::str_split(character, ',')[[1]])
  sub <- data.table::data.table(gene_length[gene_length$V1 %in% a,])
  longest_gene <- sub[which.max(sub$V2)]$V1
  return(longest_gene)
}

'Obtain representative genes\n'

genes_description <- binary_matrix$descriptions

new_rep <- parSapply(cl, genes_description, extract_long_genes, gene_length = gene_lengths)

binary_matrix$gene <- new_rep

write.table(binary_matrix, args[4], row.names = F)

'Obtain representative sequences\n'

representative_genes <- binary_matrix$gene
fasta <-  Biostrings::readAAStringSet(args[5])
subset <- fasta[representative_genes]
writeXStringSet(subset, args[6])