shell:
    "cat {input} > {output} 2> {log}"

shell:
    """
    awk -v RS='>[^\\n]+\\n' 'length() <= 2000 {{printf "%s", prt $0}} {{prt = RT}}' {input} > {output} 2> {log}
    """

shell:
    "seqtk seq -C {input} | seqtk rename - contig_ > {output} 2> {log}"

shell:
    "toAmos -s {input} -o {output} 2> {log}"

shell:
    """
    cd .snakemake/conda 
    find . > pathfiles.txt
    grep '\minimus2$' pathfiles.txt > minimus2_location.txt
    python ../../src/scripts/minimus2_replacement.py 
    touch ../../.minimus_replaced.txt
    """

shell:
    "minimus2 `file={input.file}; echo ${{file%.*}}` -D OVERLAP=100 MINID=95"

shell:
    "cat {input} > {output} 2> {log}"

shell: "cat {input.short} {input.merged} > {output} 2> {log}"

shell: "cut_up_fasta.py {input} -c 10000 -o 0 --merge_last -b contigs_10K.bed > {output}"

shell: "python src/scripts/assembly_rmdup.py"

shell: "python src/scripts/mmgenome_metabat.py"

    shell: "seqtk seq -L {params.length} {input}  > {output}"
"""

rule concatenate:
    input: expand("scratch/assembly/megahit/{treatment}/{kmers}/final.contigs.fa", treatment=config["treatment"], kmers=config["assembly-klist"])
    output: "scratch/binning/vamb/catalogue.fna.gz"
    conda: "../../../envs/vamb.yaml"
    log: "logs/vamb/concatenate.log"
    shell:     
        "concatenate.py {output} {input} 2> {log}"
    #Just cat with extras to make it more suitable to VAMB

rule vamb_map:
    input: 
        catalogue="scratch/binning/vamb/catalogue.fna.gz",
        forward="scratch/host_filtering/{sample}_R1.fastq" if config['host_removal'] \
             else "scratch/filter/{sample}_R1.fq",
        rev="scratch/host_filtering/{sample}_R2.fastq" if config['host_removal'] \
             else "scratch/filter/{sample}_R2.fq"
    output: temp("scratch/binning/vamb/{sample}.bam")
    conda: "../../../envs/minimap2.yaml"
    log: "logs/vamb/vamb_map_{sample}.log"
    shell:
        """

run:
    with open('{input.clusters}', 'w') as file:
        vamb.cluster.write_clusters(file, filtered_bins)
    keptcontigs = set.union(*filtered_bins.values())
    with open('{input.contigs}', 'rb') as file:
        fastadict = vamb.vambtools.loadfasta(file, keep=keptcontigs)
    bindir = '{params.outdir}'
    vamb.vambtools.write_bins(bindir, filtered_bins, fastadict, maxbins=500)

shell: 
   "seqtk seq -L {params.length} {input}  > {output} 2> {log}"

script:
    "../../../scripts/antismash_get_bgcs.py"

shell: "cp {input} {params.outdir} 2> {log}"

shell:
    """
    emapper.py --annotate_hits_table {input.diamond} --no_file_comments --cpu {threads} --data_dir {params.db} -o {input.seq}
    """

shell:
    """
    input=({input.forward})
    if [ "${{input##*.}}" == "bz2" ]; then
        pbzip2 -p{threads} -dc {input.forward}  > {output.forward} 2> {log.forward}
        pbzip2 -p{threads} -dc {input.rev}  > {output.rev} 2> {log.rev}
    elif [ "${{input##*.}}" == "gz" ]; then
        pigz -p {threads} -dc {input.forward}  > {output.forward} 2> {log.forward}
        pigz -p {threads} -dc {input.rev}  > {output.rev} 2> {log.rev}
    else
        cp {input.forward} {output.forward}
        cp {input.rev} {output.rev}
    fi
    """

shell: "bbmap.sh ref={params.ref} in1={input.forward} in2={input.rev} outu1={output.forward} outu2={output.rev} t={threads} 2> {log}"

shell: "bwa mem -t {threads} {params.refindex} {input.forward} {input.rev} -o {output} 2> {log}"

shell: "samtools view -b -f 4 {input} > {output} 2> {log}"

shell: "samtools sort {input} > {output}"

shell: "bamToFastq -i {input}  -fq {output.forward} -fq2 {output.rev} 2> {log}"

shell: "samtools view -b -F 4 {input} > {output} 2> {log}"

shell: "samtools sort {input} > {output}"

shell: "bamToFastq -i {input}  -fq {output.forward} -fq2 {output.rev} 2> {log}"

shell: "CAT add_names -i {input} -o {output} -t {params.tax} --only_official > {log} 2>&1"

shell: "CAT summarise -c {input.contigs} -i {input.classification} -o {output} > {log} 2>&1"

shell:
    "cat {input} > {output}"

shell: 
   "seqtk seq -L {params.length} {input}  > {output}"

shell:
    "cd-hit-est -i {input} -o {output} -T 90 -M 500000 -c 0.99 -n 10 > {log}"

shell:
    """awk '/^>/ {{print ">contig_" ++i; next}}{{print}}' < {input} > {output}"""

shell:
    "toAmos -s {input} -o {output}"

shell:
    "minimus2 `file={input}; echo ${{file%.*}}` -D OVERLAP=100 MINID=95"

shell:
    "cat {input} > {output}"

shell:
    """
    groopm parse {params.db} {input.contigs} scratch/coverm/bamfiles/*.bam -t {threads}
    groopm core {params.db} -t {threads}
    groopm refine {params.db} -t {threads}
    groopm recruit {params.db} -t {threads}
    groopm extract {params.db} {input.contigs} -o {params.outdir} -t {threads}
    """

shell: "bamm parse -c {output} -m pmean -b {input}"

shell: """echo "contig\tLength\t{params.sample}" > {output} && tail -n +2 {input} >> {output}"""

run:
    shell("grep -v '^#' {input.prediction} | cut -f1 -d ' ' > {output.posorfs}")
    shell("perl ~/install/mmgenome/scripts/extract.using.header.list.pl -l {output.posorfs} -s {input.orfs} -o {output.faa}")

shell: """
    blastp -query {input.faa} -db /data/db/blast/nr/20150311/nr -evalue 1e-5 -num_threads {threads} -max_target_seqs 5 -outfmt 5 -out {output.blast}
    # Here we need to run MEGAN first
    java -Xmx32G -Djava.awt.headless=true -Duser.language=en -Duser.region=US -cp '/data/tools/MEGAN/6.10.8/jars/MEGAN.jar:/data/tools/MEGAN/6.10.8/jars/data.jar' megan.tools.Blast2LCA -i {output.blast} -f BlastXML -ms 50 -me 0.01 -top 50 -a2t /data/db/megan/prot_acc2tax-Oct2017X1.abin
    """

run:
    import re
    import sys

    minscore = int(params.minscore)
    out = open(output.taxonomy, 'w')

    filteredtax = {}
    # Parse the taxonomy string and apply score filter
    for line in open(input.taxonomy):
        taxdict = {}
        s = line.split(';')
        for i in range(2,len(s),2):                
            try:
                level, value = s[i].strip().split('__')
                if int(s[i+1]) >= 80:
                    taxdict[level] = s[i].strip()
            except:
                pass
        taxonomy = ["root","cellular organisms"]
        for level in ["d","p","c","o","f","g","s"]:
            taxonomy.append(taxdict.setdefault(level, "unclassified"))

        taxstring = ";".join(taxonomy)
        genestring = ";".join(s[0:1])

        out.write("%s\t%s\n" % (genestring,taxstring))
    out.close()

shell: "perl ~/src/orochi/scripts/hmm.majority.vote.pl -i {input.taxonomy} -o {output.essential}"

script:
    "../../mmgenome_load_data.R"

with open('scratch/assembly/megahit/minimus2/all.merged.contigs.fasta') as f:
    lines = f.readlines()

rango = list(range(0, int(round(len(lines)/2)), 2))
headers = list()
seqs = list()
for i in rango:
    headers.append(lines[i])
    seqs.append(lines[i+1])

filtered_headers = list()
for i in range(len(headers)):
    prueba = headers[i].split(" ")
    filtered_headers.append(prueba[0])

setOfElems = set()
filtered_seqs = list()
filtered_complete_headers = list()
repetidos_descartados = list()
for elem in range(len(filtered_headers)):
    if filtered_headers[elem] in setOfElems:
        repetidos_descartados.append(filtered_headers[elem])
    else:
        setOfElems.add(filtered_headers[elem]) 
        filtered_seqs.append(seqs[elem])
        filtered_complete_headers.append(headers[elem])

with open('scratch/assembly/megahit/minimus2/filtered_assembly.fasta', "a") as new_file:
    for i in range(len(setOfElems)):
        new_file.write(filtered_complete_headers[i] + '\n')
        new_file.write(filtered_seqs[i] + '\n')

import os
with open('minimus2_location.txt', 'r') as f:
    lineas = f.readlines()
    path_to_read = lineas[0]

with open('minimus2', 'w') as new:
    pass

with open('minimus2', 'a') as new:
    with open(path_to_read[0:-1], 'r') as f:
        lineas = f.readlines()
        i = 0
        while i < 45:
            new.write(lineas[i])
            i = i+1
        correct_line = lineas[43]
    correct_path = correct_line.split("=")
    path = correct_path[1]
    path = path[0:-1]
    new.write('DELTAFILTER = ' + path + '/delta-filter' + '\n')
    new.write('SHOWCOORDS = ' + path + '/show-coords' + '\n')

with open('minimus2', 'a') as new:
    with open(path_to_read[0:-1], 'r') as f:
        lineas = f.readlines()
        i = 47
        while i > 46 and i < 57:
            new.write(lineas[i])
            i = i+1
    new.write('20: $(NUCMER) --maxmatch -c $(OVERLAP) $(REFSEQ) $(QRYSEQ) -p $(PREFIX)' + '\n')

with open('minimus2', 'a') as new:
    with open(path_to_read[0:-1], 'r') as f:
        lineas = f.readlines()
        i = 58
        while i > 57 and i < 85:
            new.write(lineas[i])
            i = i+1

os.chmod("minimus2", 0o777)
os.replace("minimus2", path_to_read[0:-1])