shell:
    """
    zcat {input} > {output[0]}
    mkdir {output[1]}
    splitMfasta.pl {output[0]} \
    --outputpath={output[1]} --minsize={params.split_size}
    """

shell:
    "cat {input} | join_aug_pred.pl > {output}"

shell:
    """
    run_rcorrector.pl -1 {input[0]} -2 {input[1]} \
    -t {threads} \
    -od {params.outdir} \
    > {log} 2>&1
    """

shell:
    """
    samtools merge -@ {threads} {params.tmp_merge} {input}               
    samtools sort -@ {threads} -n -o {output} {params.tmp_merge}
    rm {params.tmp_merge}
    """

shell:
    """
    python /opt/agouti/agouti.py scaffold \
    -assembly {input.fa} \
    -bam {input.bam} \
    -gff {input.gff} \
    -outdir {params.outdir} \
    -minMQ {params.minMQ} -maxFracMM {params.maxFracMM} \
    > {log} 2>&1

    gzip -c {params.outdir}/agouti.agouti.fasta > {output}
    """

shell:
    """
    hisat2 -p {threads} \
    -x {params.index_prefix} \
    -1 {input.reads[0]} -2 {input.reads[1]} \
    2> {log} | \
    samtools sort -@ {threads} -O BAM -o {output[0]}
    samtools index {output}
    """

shell:
    """
    find {params.out_dir} -name '{params.to_rm}' -type d -delete
    cp /opt/gm_key_64 ~/.gm_key
    braker.pl \
    --genome {input.genome} \
    --prot_seq {input.prot_db} \
    --bam {params.list_bams} \
    --workingdir {params.out_dir} \
    --species {params.species} \
    --useexisting \
    --etpmode --softmasking --cores {threads} \
    --gff3 \
    --GENEMARK_PATH {params.genemark_path} \
    --PROTHINT_PATH {params.genemark_path}/ProtHint/bin \
    > {log} 2>&1
    """

shell:
    """
    gff3_to_fasta -g {input.gff} -f {input.fa} \
    -st pep -d simple -o {params.prefix} \
    > log 2>&1
    """

shell:
    """
    rm -r {params.outdir}
    mantis run \
    -i {input.pep} \
    -o {params.outdir} \
    -od '{params.org_detail}' \
    -km \
    -gff \
    -c {threads} \
    > {log} 2>&1
    """

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

shell:
    "wget -O {output} {params.link}"

shell:
    """
    java -jar {input.lrscaf} \
    -c {input.ref} \
    -a {input.paf} \
    -o {params.res}/res \
    -micl 1000 \
    -t mm \
    -i 0.3 \
    -p {threads} \
    > {log} 2>&1

    cp {params.res}/res/scaffolds.fasta {output}
    """

shell:
    """
    ragtag.py scaffold \
    {input.ref_coruscus} \
    {input.draft_assembly} \
    -o {params.out_dir} \
    --mm2-params '-x asm10 -t {threads}' -u \
    > {log} 2>&1

    cp {params.out_dir}/ragtag.scaffolds.fasta {output}
    """

shell:
    """
    minimap2 -t {threads} -ax map-ont \
    {input.ref} {input.ont_reads} 2> {log} | \
    samtools sort -@ {threads} | \
    samtools view -b -@ {threads} -o {output}
    samtools index {output}
    """

shell:
    """
    java -Xmx{params.java_mem}G -jar {input.pilon} \
    --genome {input.ref} \
    --frags {input.pe_bam} \
    --nanopore {input.ont_bam} \
    --targets {input.target} \
    --output {params.prefix} \
    --outdir {params.out_dir} \
    --changes --vcf --tracks --diploid --fix all \
    > {log} 2>&1
    """

shell:
    """
    ragtag.py scaffold \
    {input.ref_coruscus} \
    {input.draft_assembly} \
    -o {params.out_dir} \
    --mm2-params '-x asm10 -t {threads}' -u \
    > {log} 2>&1

    cp {params.out_dir}/ragtag.scaffolds.fasta {output}
    """

shell:
    "wget -O {output} {params.link}"

shell:
    """
    bwa-mem2 index {input}
    """

shell:
    """
    [ ! -e {input.ref}.0123 ] && bwa-mem2 index {input.ref}
    bwa-mem2 mem -t {threads} {input.ref} {input.R1} {input.R2} 2> {log} | \
    samtools sort -@ {threads} | \
    samtools view -b -@ {threads} -o {output}
    samtools index {output}
    """

shell:
    """
    for i in {{1..14}}; do
        grep '>' {input.ref} | sed 's/>//' | awk -v record=$i 'NR==record {{print $0}}' \
        > {params.out_dir}/target_$(printf '%02d' $i).txt
    done
    grep '>' {input.ref} | tail -n +15 | sed 's/>//' > {output[14]}
    """

shell:
    """
    java -Xmx{params.java_mem}G -jar {input.pilon} \
    --genome {input.ref} \
    --frags {input.pe_bam} \
    --targets {input.target} \
    --output {params.prefix} \
    --outdir {params.out_dir} \
    --changes --vcf --tracks --diploid --fix all \
    > {log} 2>&1
    """

shell:
    "cat {input} | bgzip -c > {output}"

shell:
    """
    cp {input.conf} {output.conf}
    SC=$( grep '>' {input.fa} | wc -l )
    SP=$( grep -v '>' {input.fa} | wc -m )
    sed -i "s/scaffold-count.*/scaffold-count:\ $SC/" {output.conf}
    sed -i "s/span.*/span:\ $SP/" {output.conf}
    sed -i "s/prefix.*/prefix:\ {wildcards.sample}_{wildcards.version}/" {output.conf}
    """

shell:
    """
    {params.blobtools_bin} add \
    --busco {input[2]} \
    --busco {input[1]} \
    {params.blobdir}
    """

shell:
    """
    cp {input[1]} {output[1]}
    cp -r {params.indir} results/blobtoolkit/blobdirs/ && \
    tar -czf {output[2]} {params.tardir} && \
    rm -r {params.tardir}
    """

shell:
    """
    cd {params.outdir}
    wget -c {params.metazoa} -O - | tar -xz
    wget -c {params.mollusca} -O - | tar -xz
    """

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

shell:
    """
    busco -f -m genome -i {params.fa} -o {params.outdir} \
    -q -c {threads} \
    -l {params.db} > {log} 2>&1
    cp -r {params.outdir} results/busco/ && rm -r {params.outdir}
    """

shell:
    """
    python workflow/scripts/summarize_buscos.py \
    -o {output} \
    --files {input}
    """

shell:
    "seqkit replace -p '\s.*$' -r '' {input} > {output}"

    shell:
        """
        echo "((((mgal_02,GCA900618805),medu_02),mtro_02),GCA017311375);" > {params.seqfile}
        echo "*GCA017311375 {input[0]}" >> {params.seqfile}
        echo "mtro_02 {input[1]}" >> {params.seqfile}
        echo "medu_02 {input[2]}" >> {params.seqfile}
        echo "mgal_02 {input[3]}" >> {params.seqfile}
        echo "GCA900618805 {input[4]}" >> {params.seqfile}

        cactus \
        {params.restart} \
	    --maxCores {threads} \
        results/cactusJobStore \
        {params.seqfile} \
        {output} \
        > {log} 2>&1
        """

shell:
    """
    paste <(zcat {input[0]}) <(zcat {input[1]} | cut -f 4-) | \
    bgzip -c > {output}
    """

shell:
    """
    python /opt/pcangsd/pcangsd.py \
    -beagle {input} \
    -minMaf {params.min_maf} \
    -iter 200 \
    -e {params.eigenvectors} \
    -o {params.prefix} \
    -sites_save \
    -snp_weights \
    -admix \
    -threads {threads} \
    > {log} 2>&1
    """

shell:
    """
    bwa-mem2 index {input} > {log} 2>&1
    """

shell:
    """
    bwa-mem2 mem -t {threads} \
    -R \"{params.rg_string}\" \
    {input.fa} \
    {input.fastqs} 2> {log} | \
    samtools sort -@ {threads} -o {output[0]} -
    samtools index {output[0]}
    """

shell:
    """
    bcftools mpileup \
    -f {input.ref} \
    -R {input.targets} \
    --redo-BAQ -a "FORMAT/AD,FORMAT/DP" \
    -Ou {input.bams} 2> {log} | \
    bcftools call -mG - -Ob --threads {threads} -o {output} \
    >> {log} 2>&1
    """

shell:
    """
    bcftools index -f {input.bcf1}
    bcftools index -f {input.bcf2}
    bcftools merge -m snps -Ob -o {output} -R {input.targets} \
    --threads {threads} {input.bcf1} {input.bcf2} > {log} 2>&1
    """

script:
    "../scripts/btk_conta_extraction.py"

shell:
    """
    tail -n +2 {input[1]} | cut -d ',' -f 1 > {params.tmp_kept_list}

    {params.blobtools_bin} filter \
    --list {params.tmp_kept_list} \
    --output {output[0]} {params.blobdir}

    {params.blobtools_bin} filter \
    --list {params.tmp_kept_list} \
    --invert \
    --output {output[1]} {params.blobdir}

    rm {params.tmp_kept_list}
    """

shell:
    """
    seqkit grep -f <(tail -n +2 {input.kept} | cut -d ',' -f 1) {input[0]} \
    | seqkit replace -is -p "^N+|N+$" -r "" > {input.fa}_tmp

    seqkit fx2tab -n -i --gc --length -B N {input.fa}_tmp \
    | awk '($2 >= 1000 && $4 < 90) {{print $1}}' > {input.fa}_filt_list
    seqkit grep -f {input.fa}_filt_list {input.fa}_tmp > {input.fa}_tmp2

    seqkit sort --by-length -2 --reverse {input.fa}_tmp2 \
    | seqkit replace -p .+ -r "{params.scaff_prefix}_{{nr}}" --nr-width {params.nr_width} \
    | gzip -c > {output}

    rm {input.fa}_tmp* {input.fa}_filt_list
    """

shell:
    """
    wget {params.ftp_0} -O {output[0]}
    wget {params.ftp_1} -O {output[1]}
    wget {params.ftp_2} -O {output[2]}
    wget {params.ftp_3} -O {output[3]}
    """

shell:
    """
    cp {input} {output}
    """

shell:
    """
    agat_sp_manage_IDs.pl --gff {input} \
    --ensembl --prefix {params.prefix} \
    -o {output} > {log} 2>&1
    """

shell:
    """
    agat_convert_sp_gxf2gxf.pl -g {input} -gvi 3 -gvo 3 \
    -o {output} > {log} 2>&1
    """

shell:
    "cat {input} | gzip -c > {output}"

shell:
    "cp {input} {output}"

shell:
    """
    kmc_tools -t{threads} transform \
    {params.db} histogram \
    {output} -cx10000 \
    > {log} 2>&1
    """

shell:
    """
    /opt/genomescope2.0/genomescope.R -i {input} \
    -o {params.outdir} -n {params.name_prefix} \
    -k {params.k} -p 2 \
    > {log} 2>&1
    """

shell:
    """
    kat plot spectra-cn \
    -o {output} \
    -t {params.title} \
    {input} \
    > {log} 2>&1
    """

shell:
    """
    bwa-mem2 index {input} > {log} 2>&1
    """

shell:
    """
    bwa-mem2 mem -t {threads} {input.fa} \
    {input.fastqs} 2> {log} | \
    samtools sort -@ {threads} -o {output[0]} -
    samtools index {output[0]}
    """

shell:
    """
    samtools stats -@ {threads} {input} > {output[0]}
    """

shell:
    """
    seqkit rmdup -s -D {params} -i \
    -o {output} {input}
    """

shell:
    """
    seqkit replace -p .+ \
    -r "{params.scaff_prefix}_s{{nr}}" --nr-width {params.nr_width} \
    -o {output} {input}
    """

shell:
    """
    seqkit head -n 14 {input} | \
    seqkit replace -p "(.+)" -r "\$1 [location=chromosome] [chromosome=LG{{nr}}]" \
    --nr-width 2 \
    -o {output}

    seqkit fx2tab {input} | tail -n +15 | seqkit tab2fx | 
    gzip -c >> {output}
    """

shell:
    """
    paste \
    <(seqkit seq --name --only-id {input[0]}) \
    <(seqkit seq --name --only-id {input[1]}) \
    > {output}
    """

shell:
    """
    awk -F'\\t' 'NR==FNR{{a[$1]=$2; next}}{{id=$1; if(id in a) gsub(id,a[id])}} 1' \
    {input[0]} <(zcat {input[1]}) | \
    gzip -c > {output}
    """

shell:
    """
    cp {input[0]} {output[0]}
    cp {input[1]} {output[1]}
    """

shell:
    """
    /opt/proc10xG/process_10xReads.py \
    -o {params.out_prefix} \
    -1 {input[0]} -2 {input[1]} \
    > {log} 2>&1
    """

script:
    "../scripts/filter_barcodes.R"

shell:
    """
    /opt/proc10xG/filter_10xReads.py \
    -L {input.barcodes} \
    -1 {input[0]} -2 {input[1]} \
    -o {params.out_prefix} \
    > {log} 2>&1
    """

shell:
    """
    /opt/proc10xG/regen_10xReads.py \
    -1 {input[0]} -2 {input[1]} \
    -o {params.out_prefix} \
    > {log} 2>&1
    """

shell:
    """
    zcat {input.fa} > {params.tmp_fa}
    longranger mkref {params.tmp_fa} > {log} 2>&1
    mv refdata-{wildcards.sample}_v1.pseudohap {output[0]}
    """

shell:
    """
    longranger align \
    --id={params.run_id} \
    --fastqs={params.input_dir} \
    --sample={params.sample} \
    --reference={input[2]} \
    --localcores={threads} \
    --localmem={params.mem} \
    > {log} 2>&1

    rm -r {output[0]} \
    && cp -al {params.run_id} {output[0]} \
    && rm -r {params.run_id}
    """

shell:
    "samtools sort -@ {threads} -n -O BAM -o {output} {input} > {log} 2>&1"

shell:
    """
    cd {params.workdir}
    ngscstat {params.input} 2> ../../../{log}
    """

shell:
    """
    calcuts {input[0]} > {output} 2> {log}
    """

shell:
    "/opt/purge_dups/scripts/hist_plot.py "
    "-c {input[0]} "
    "{input[1]} "
    "{output}"

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

shell:
    "(minimap2 -t {threads} "
    "-xasm5 -DP {input} {input} "
    "| gzip -c - > {output}) 2> {log}"

shell:
    "purge_dups -2 -M{params.M} -E{params.E} -T {input.cutoffs} "
    "-c {input.basecov} {input.selfmap} > {output} "
    "2> {log}"

shell:
    """
    get_seqs {input.bed} {input.fa} \
    -p {params.prefix}
    gzip {params.prefix}.*.fa
    """

shell:
    "cp {input} {output}"

shell:
    "cp {input} {output}"

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

shell:
    """
    export LC_ALL=C
    BuildDatabase -name {params.db_name} {input} \
    > {log} 2>&1
    """

shell:
    """
    zcat {input} | \
    assembly_stats /dev/stdin \
    > {output}
    """

script:
    "../scripts/merge_stats.py"

shell:
    """
    d4tools stat -s hist {input} > {output}
    """

script:
    "../scripts/plot_coverage_hist.R"

shell:
    """
    supernova mkoutput \
    --style={params.style} \
    --asmdir={params.asm_dir} \
    --outprefix={params.outprefix} \
    --headers=full \
    > {log} 2>&1
    """

shell:
    """
    tar -cf - -C {params.input_dir}/outs assembly | zstdmt -T{threads} > {params.tmp_archive} 2> {log} && \
    rm -r {params.input_dir}/outs/assembly && \
    cp -r {params.input_dir} {params.output_dir} && rm -r {params.input_dir}
    """

shell:
    "bash workflow/scripts/collect_assembly_stats.sh"

echo "species,version,$(head -n1 results/supernova_assemblies/gallo_v1/outs/summary.csv)" \
> results/supernova_assemblies/supernova_assemblies_stats.csv

for S in 'gallo' 'edu' 'tros'; do
	for V in 'v1' 'v2'; do
		echo "${S},${V},$(tail -n +2 results/supernova_assemblies/${S}_${V}/outs/summary.csv)" \
		>> results/supernova_assemblies/supernova_assemblies_stats.csv
	done
done

import pandas as pd
import json
import re

out_df = pd.DataFrame()
for stat_file in snakemake.input:
    m = re.search('results/stats/(.+?)_(.+?).stats.json', stat_file)
    sample = m.group(1)
    version = m.group(2)
    with open(stat_file, 'r') as f:
        data = json.loads(f.read())
    data['C'] = data.pop('Contig Stats')
    data['S'] = data.pop('Scaffold Stats')
    tmpdf = pd.json_normalize(data)
    tmpdf.insert(0, "asm", [f"{sample}_{version}"])
    out_df = out_df.append(tmpdf)

out_df.to_csv(snakemake.output[0], sep=',', index=False)