shell:
    '(echo "$(( $(zcat {input.fastq1} | wc -l) / 4 ))" > {output.read_count}) 2> {log}'

shell:
    "(megahit -1 {input.fastq1} -2 {input.fastq2} {params.preset} --out-dir {params.outdir} -f && "
    " mv {params.outdir}/final.contigs.fa {output})"
    " > {log} 2>&1"

shell:
    "(megahit -r {input} {params.preset} --out-dir {params.outdir} -f && "
    " mv {params.outdir}/final.contigs.fa {output})"
    " > {log} 2>&1"

shell:
    "({wildcards.spadesflavor}.py -1 {input.fastq1} -2 {input.fastq2} -o {params.outdir} -t {threads} && "
    " if [ -f {params.outdir}/raw_contigs.fasta ]; then mv {params.outdir}/raw_contigs.fasta {output}; else mv {params.outdir}/contigs.fasta {output}; fi )"
    " > {log} 2>&1"

shell:
    "if [ -d {params.outdir} ]; then rm -Rf {params.outdir}; fi && "
    "(spades.py --corona -s {input} -o {params.outdir} -t {threads} && "
    " mv {params.outdir}/raw_contigs.fasta {output})"
    " > {log} 2>&1"

script:
    "../scripts/check_contigs.py"

shell:
    "(mkdir -p {params.outdir}/{wildcards.sample} && cd {params.outdir}/{wildcards.sample} &&"
    " ragoo.py ../../../../../{input.contigs} ../../../../../{input.reference} &&"
    " cd ../../../../../ && mv {params.outdir}/{wildcards.sample}/ragoo_output/ragoo.fasta {output})"
    " > {log} 2>&1"

script:
    "../scripts/ragoo-remove-chr0.py"

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

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

shell:
    "(zcat {input.left} > {output.left} &&"
    "zcat {input.right} > {output.right}) 2>{log}"

shell:
    "minimap2 --MD --eqx -ax asm5 {input} -o {output} 2> {log}"

script:
    "../scripts/assembly-benchmark-results.py"

script:
    "../scripts/summarize-non-cov2.py"

script:
    "../scripts/generate-mixtures.py"

script:
    "../scripts/evaluate-strain-call-error.py"

script:
    "../scripts/plot-caller-error.py"

shell:
    "(Trinity --left {input.fastq1} --max_memory 16G --right {input.fastq2} --CPU {threads} --seqType fq --output {params.outdir} && "
    "mv {params.outdir}.Trinity.fasta {output} ) > {log} 2>&1"

shell:
    """
    velveth {params.outdir} 21 -fastq.gz -shortPaired {input.fastq1} {input.fastq2} > {log} 2>&1
    velvetg {params.outdir} -ins_length 150 -exp_cov 10 >> {log} 2>&1
    mv {params.outdir}/contigs.fa {output} >> {log} 2>&1
    """

shell:
    "(cd {params.outdir} &&"
    " ragoo.py ../../../../../{input.contigs} ../../../../../{input.reference} &&"
    " cd ../../../../../ && mv {params.outdir}/ragoo_output/ragoo.fasta {output})"
    " > {log} 2>&1"

script:
    "../scripts/plot-assembly-comparison.py"

script:
    "../scripts/get-read-statistics.py"

script:
    "../scripts/plot-dependency-of-pangolin-call.py"

script:
    "../scripts/plot-pangolin-conflict.py"

script:
    "../scripts/collect_lineage_calls.py"

script:
    "../scripts/get_largest_contig.py"

script:
    "../scripts/select_random_lineages.py"

script:
    "../scripts/aggregate_read_calls.py"

script:
    "../scripts/benchmarking/compare-vcf.py"

script:
    "../scripts/benchmarking/aggregate-test-case-variants.py"

script:
    "../scripts/benchmarking/filter-test-case-variants.py"

script:
    "../scripts/benchmarking/get-test-case-variant-paths.py"

script:
    "../scripts/benchmarking/check-presence-of-test-case-variant-in-call.py"

script:
    "../scripts/masking.py"

script:
    "../scripts/plot-all-coverage.py"

script:
    "../scripts/plot-all-coverage.py"

script:
    "../scripts/quality-filter.py"

script:
    "../scripts/generate-high-quality-report.py"

script:
    "../scripts/generate-overview-table.py"

shell:
    "rbt csv-report {input} --formatter {params.formatter} --pin-until {params.pin_until} {output} > {log} 2>&1"

script:
    "../scripts/generate-filter-overview.py"

shell:
    "rbt csv-report {input} --pin-until {params.pin_until} {output} > {log} 2>&1"

script:
    "../scripts/plot-lineages-over-time.py"

script:
    "../scripts/plot-variants-over-time.py"

script:
    "../scripts/aggregate-pangolin-calls-per-stage.py"

shell:
    "rbt csv-report {input} --pin-until {params.pin_until} {output} > {log} 2>&1"

script:
    "../scripts/collect-lineage-variants.py"

shell:
    "bcftools annotate -a {input.annotation} -c LINEAGES,SIGNATURES {input.calls} > {output} 2> {log}"

script:
    "../scripts/generate-lineage-variant-table.py"

shell:
    "nanoQC {input} -o {params.outdir} > {log} 2>&1"

shell:
    "echo $(( $(cat {input} | wc -l ) / 4)) > {output} 2> {log}"

shell:
    "notramp -a -r {input.reads} -p {input.primer} -g {input.ref_genome} -o {params.outdir}  2> {log}"

shell:
    "notramp -t --incl_prim -r {input.reads} -p {input.primer} -g {input.ref_genome} -o {params.outdir}  2> {log}"

shell:
    "bcftools consensus -f {input.fasta} {input.bcf} > {output} 2> {log}"

shell:
    "(cp {input} {output} && "
    ' sed -i "1s/.*/>{wildcards.sample}/" {output})'
    " 2> {log}"

script:
    "../scripts/update-sample-sheet.py"

script:
    "../scripts/vcf-to-fasta.py"

shell:
    "minimap2 --eqx -ax asm5 {input.assembly} {input.pseudoassembly} -o {output} 2> {log}"

script:
    "../scripts/aggregate-assembly-comparisons.py"

wrapper:
    "v1.15.1/bio/multiqc"

wrapper:
    "v1.15.1/bio/samtools/flagstat"

shell:
    "(samtools depth -aH -o {output} {input} && "
    " sed -i 's/{params.ref}.3/{wildcards.sample}/' {output})"
    " 2> {log}"

shell:
    "(kraken2 --db {input.db} --threads {threads} --unclassified-out {params.unclassified} "
    "--classified-out {params.classified} --report {output.report} --gzip-compressed "
    "--paired {input.reads} > {output.kraken_output}) 2> {log}"

shell:
    "(kraken2 --db {input.db} --threads {threads}"
    " --report {output.report} --gzip-compressed"
    " {input.reads} > {output.kraken_output})"
    "2> {log}"

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

script:
    "../scripts/extract-reads-of-interest.py"

shell:
    "(samtools sort  -@ {threads} -n {input} -o {output.bam_sorted}; "
    " samtools fastq -@ {threads} -0 {output.fq} {output.bam_sorted})"
    "> {log} 2>&1"

shell:
    "(kraken2 --db {input.db} --threads {threads} --report {output.report} --gzip-compressed "
    "--paired {input.reads} > {output.kraken_output}) 2> {log}"

shell:
    "(kraken2 --db {input.db} --threads {threads} --report {output.report} --gzip-compressed "
    "{input.reads} > {output.kraken_output}) 2> {log}"

shell:
    "ktImportTaxonomy -m 3 -t 5 -tax {input.taxonomy_database} -o {output} {input} 2> {log}"

wrapper:
    "v1.15.1/bio/samtools/sort"

script:
    "../scripts/bed-to-bedpe.py"

shell:
    "(cp {input.bam} {params.output_dir} &&"
    " cp {input.bai} {params.output_dir} &&"
    " cd {params.output_dir} &&"
    " bamclipper.sh -b {params.bam} -p {params.bed_path} -n {threads} -u 5 -d 5) "
    " > {params.cwd}/{log} 2>&1"

shell:
    "samtools fastq -@ {threads} {input.bam} -1 {output.fq1} -2 {output.fq2} 2> {log}"

shell:
    "samtools fastq -@ {threads} {input.bam} > {output} 2> {log}"

script:
    "../scripts/plot-primer-clipping.py"

wrapper:
    "v1.15.1/bio/bwa/index"

wrapper:
    "v1.15.1/bio/bwa/index"

wrapper:
    "v1.15.1/bio/bwa/mem"

wrapper:
    "v1.15.1/bio/picard/markduplicates"

wrapper:
    "v1.15.1/bio/samtools/calmd"

wrapper:
    "v1.15.1/bio/fastp"

shell:
    "curl -sSL https://www.ncbi.nlm.nih.gov/sars-cov-2/download-nuccore-ids > {output} 2> {log}"

shell:
    "((esearch -db nucleotide -query '{params.accession}' | "
    "efetch -format fasta > {output}) && [ -s {output} ]) 2> {log}"

shell:
    "(bigBedToBed http://hgdownload.soe.ucsc.edu/gbdb/wuhCor1/uniprot/unipChainCov2.bb"
    " -chrom=NC_045512v2 -start=0 -end=29903 {output})"
    "2>{log}"

shell:
    "(cut -f1-12 {input} | sed -e 's/ /-/g' | sed -e 's/NC_045512v2/NC_045512.2/g'"
    " | gt bed_to_gff3 -featuretype gene -thicktype transcript -blocktype CDS -o {output} -force /dev/stdin )"
    "2> {log}"

script:
    "../scripts/fix-protein-gff.py"

shell:
    "bgzip -c {input} > {output} 2> {log}"

shell:
    "curl -sSL https://raw.githubusercontent.com/W-L/ProblematicSites_SARS-CoV2/"
    "master/problematic_sites_sarsCov2.vcf | bgzip -c > {output} 2> {log}"

shell:
    "mkdir {output} && curl -SL https://genome-idx.s3.amazonaws.com/kraken/k2_standard_08gb_20220926.tar.gz | tar zxvf - -C {output} 2> {log}"

shell:
    "ktUpdateTaxonomy.sh {output} 2> {log}"

shell:
    "curl -SL -o {output} {params.human_genome} 2> {log}"

shell:
    "(mkdir -p {output} &&"
    " curl -L https://github.com/cov-lineages/pangoLEARN/archive/master.tar.gz |"
    " tar xvz --strip-components=1 -C {output})"
    " > {log} 2>&1"

shell:
    "(mkdir -p {output} &&"
    " curl -L https://github.com/cov-lineages/lineages/archive/master.tar.gz | "
    " tar xvz --strip-components=1 -C {output})"
    " > {log} 2>&1"

shell:
    "(curl -L -u $GISAID_API_TOKEN https://www.epicov.org/epi3/3p/resseq02/export/provision.json.xz |"
    " xz -d -T0 > {output})"
    " > {log} 2>&1"

shell:
    "sed -E 's/>(\S+)\\b/>{params.lineage}/;t' {input} > {output}"

script:
    "../scripts/get-strains-from-genbank.py"

script:
    "../scripts/extract-strains-from-gisaid-provision.py"

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

wrapper:
    "v1.15.1/bio/kallisto/index"

shell:
    "awk 'BEGIN {{ t=0.0;sq=0.0; n=0; }} ;NR%4==2 {{ n++;L=length($0);t+=L;sq+=L*L; }}END{{ m=t/n;printf(\"%f\\n\",m) ; }}' {input} && "
    "awk 'BEGIN {{ t=0.0;sq=0.0; n=0; }} ;NR%4==2 {{ n++;L=length($0);t+=L;sq+=L*L; }}END{{ m=t/n;printf(\"%f\\n\",sq/n-m*m) ; }}' {input} && "
    "awk 'BEGIN {{ t=0.0;sq=0.0; n=0; }} ;NR%4==2 {{ n++;L=length($0);t+=L;sq+=L*L; }}END{{ m=t/n;printf(\"%f\\n\",m) ; }}' {input} > {output.avg_read_length} && "
    "awk 'BEGIN {{ t=0.0;sq=0.0; n=0; }} ;NR%4==2 {{ n++;L=length($0);t+=L;sq+=L*L; }}END{{ m=t/n;printf(\"%f\\n\",sq/n-m*m) ; }}' {input} > {output.standard_deviation}"

wrapper:
    "v1.15.1/bio/tabix/index"

wrapper:
    "v1.15.1/bio/samtools/index"

shell:
    "bcftools index {input} 2> {log}"

shell:
    "bcftools sort -O b {input} -o {output} 2> {log}"

wrapper:
    "v1.15.1/bio/samtools/faidx"

shell:
    "gzip --keep {input}"

wrapper:
    "v1.15.1/bio/vep/plugins"

wrapper:
    "v1.15.1/bio/vep/annotate"

script:
    "../scripts/delly.py"

shell:
    "(medaka_haploid_variant -i {input.sample} -r {input.ref} -o medaka_tmp/medaka/{wildcards.date}/{wildcards.reference}/{wildcards.sample}"
    " -t {threads} -m {params.model} && mv medaka_tmp/medaka/{wildcards.date}/{wildcards.reference}/{wildcards.sample}/medaka.annotated.vcf {output} &&"
    " rm -r medaka_tmp/medaka/{wildcards.date}/{wildcards.reference}/{wildcards.sample}) > {log} 2>&1"

shell:
    "(longshot -P 0 -F -A --no_haps --bam {input.bam} --ref {input.ref} --out {output} &&"
    " sed -i '2 i\##contig=<ID={params.reference_name}>' {output})"
    " 2> {log}"

    shell:
        "bcftools view -Oz -o {output} {input} 2> {log}"  # TODO -Oz generates a .vcf.gz!!


rule render_scenario:
    input:
        local(get_resource("scenario.yaml")),
    output:
        "results/{date}/scenarios/{sample}.yaml",
    log:
        "logs/{date}/render-scenario/{sample}.log",
    conda:
        "../envs/unix.yaml"

shell:
    "sed 's/sample:/{wildcards.sample}:/' {input} > {output}"

shell:
    "varlociraptor estimate alignment-properties {input.ref} --bam {input.bam} > {output} 2> {log}"

shell:
    "varlociraptor preprocess variants {params.extra} --candidates {input.candidates} "
    "{input.ref} --bam {input.bam} --max-depth {params.depth} --output {output} 2> {log}"

shell:
    "varlociraptor "
    "call variants {params.biases} generic --obs {wildcards.sample}={input.obs} "
    "--scenario {input.scenario} > {output} 2> {log}"

wrapper:
    "v1.15.1/bio/bcftools/concat"

wrapper:
    "v1.15.1/bio/vembrane/filter"

shell:
    "varlociraptor filter-calls control-fdr --mode local-smart {input} --var {wildcards.vartype} "
    "--events {params.events} --fdr {params.fdr} > {output} 2> {log}"

wrapper:
    "v1.15.1/bio/bcftools/concat"

shell:
    "rbt vcf-report {input.ref} --bams {params.bams} --vcfs {params.bcfs} --formats {params.format_field} "
    "--infos PROB_* -d {params.max_read_depth} -l {params.js_files} -- {output} 2> {log}"

script:
    "../scripts/ucsc_vcf.py"

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

import sys

sys.stderr = open(snakemake.log[0], "w")

import pandas as pd

pangolin_calls = []

assert len(snakemake.input) == len(snakemake.params.samples)
assert len(snakemake.input) == len(snakemake.params.stages)

for path, sample, stage in zip(
    snakemake.input, snakemake.params.samples, snakemake.params.stages
):
    pangolin_call = pd.read_csv(path)
    pangolin_call["sample"] = sample
    pangolin_call["stage"] = stage

    pangolin_calls.append(pangolin_call)

pangolin_calls_by_stage = pd.concat(pangolin_calls, axis=0, ignore_index=True)

failed_called = (pangolin_calls_by_stage["qc_status"] == "fail") | (
    pangolin_calls_by_stage["lineage"] == "None"
)
pangolin_calls_by_stage.loc[failed_called, "lineage"] = (
    "Call failed. Reason: " + pangolin_calls_by_stage.loc[failed_called, "note"]
)

pangolin_calls_by_stage = pangolin_calls_by_stage.pivot(
    index="sample", columns="stage", values="lineage"
).reset_index()

only_pseudo_assembly = ["scaffold", "polished", "masked-polished", "pseudo"]
only_consensus_assembly = [
    "scaffold",
    "polished",
    "masked-polished",
    "consensus",
    "masked-consensus",
]
both_fallbacks = list(
    set(only_pseudo_assembly).symmetric_difference(set(only_consensus_assembly))
)

columns = pangolin_calls_by_stage.columns.to_list()

if all(stage in columns for stage in both_fallbacks):
    sorted_columns = [
        "sample",
        "scaffold",
        "polished",
        "masked-polished",
        "consensus",
        "masked-consensus",
        "pseudo",
    ]
elif all(stage in columns for stage in only_pseudo_assembly):
    sorted_columns = ["sample", "scaffold", "polished", "masked-polished", "pseudo"]
elif all(stage in columns for stage in only_consensus_assembly):
    sorted_columns = [
        "sample",
        "scaffold",
        "polished",
        "masked-polished",
        "consensus",
        "masked-consensus",
    ]

pangolin_calls_by_stage = pangolin_calls_by_stage[sorted_columns]

pangolin_calls_by_stage.rename(
    columns={
        "sample": "Sample",
        "scaffold": "Scaffolded Seq",
        "polished": "Polished Seq",
        "masked-polished": "Masked Polished Seq",
        "pseudo": "Pseudo Seq",
        "consensus": "Consensus Seq",
        "masked-consensus": "Masked Consensus Seq",
    },
    inplace=True,
)

pangolin_calls_by_stage.fillna("Not called on", inplace=True)
# pangolin_calls_by_stage.replace({"None": "Call failed"}, inplace=True)
pangolin_calls_by_stage.sort_values(by="Sample", inplace=True)

print(pangolin_calls_by_stage)

pangolin_calls_by_stage.to_csv(snakemake.output[0], index=False)

import sys

sys.stderr = open(snakemake.log[0], "w")

import pandas as pd


def aggregate_calls(sm_input, sm_output):
    all_calls = []
    for i, file in enumerate(sm_input):
        call = pd.read_csv(file, sep="\t")
        call["sample"] = i
        all_calls.append(call)

    all_calls = pd.concat(all_calls, axis=0, ignore_index=True)
    all_calls.to_csv(sm_output, sep="\t", index=False)


if __name__ == "__main__":
    aggregate_calls(snakemake.input, snakemake.output[0])

sys.stderr = open(snakemake.log[0], "w")
import pysam

# see https://pysam.readthedocs.io/en/latest/api.html#pysam.AlignedSegment.cigartuples
BAM_CEQUAL = 7
BAM_CDIFF = 8
BAM_CSOFT_CLIP = 4
BAM_CHARD_CLIP = 5
BAM_CINS = 1
BAM_CDEL = 2

N_WINDOW = 100


def is_uncertain_region(record, rpos, rend, ref_fasta):
    refseq = ref_fasta.fetch(record.reference_name, rpos - N_WINDOW, rend + N_WINDOW)
    return ("N" * 10) in refseq


def get_edit_dist(record, ref_fasta):
    edit_dist = 0
    qpos = 0
    rpos = record.reference_start
    for op, length in record.cigartuples:
        if op == BAM_CEQUAL:
            # no edit
            qpos += length
            rpos += length
        else:
            if op == BAM_CDIFF:
                refseq = ref_fasta.fetch(record.reference_name, rpos, rpos + length)
                if not is_uncertain_region(record, rpos, rpos + length, ref_fasta):
                    # Only consider edits where the reference has a true nucleotide
                    # because IUPAC codes lead to mason simulating an N in the reads.
                    edit_dist += sum(base in "ACGT" for base in refseq)
                qpos += length
                rpos += length

            elif op == BAM_CSOFT_CLIP:
                edit_dist += length
                qpos += length
            elif op == BAM_CHARD_CLIP:
                edit_dist += length
            elif op == BAM_CINS:
                refseq = ref_fasta.fetch(record.reference_name, rpos, rpos + length)
                if not is_uncertain_region(record, rpos, rpos + length, ref_fasta):
                    # only count edit distance if the region behind the insert does not
                    # contain N stretches in the reference. Rationale: those regions are apparently
                    # hard to assemble, and we cannot properly simulate reads for them, so
                    # we should not count them in a benchmark.
                    edit_dist += length
                qpos += length
            elif op == BAM_CDEL:
                refseq = ref_fasta.fetch(record.reference_name, rpos, rpos + length)
                if not is_uncertain_region(record, rpos, rpos + length, ref_fasta):
                    # only count edit distance if the region behind the insert does not
                    # contain N stretches in the reference. Rationale: those regions are apparently
                    # hard to assemble, and we cannot properly simulate reads for them, so
                    # we should not count them in a benchmark.
                    edit_dist += length
                rpos += length
            else:
                raise ValueError(f"Unsupported CIGAR operation: {op}")
    return edit_dist


with open(snakemake.output[0], "w") as out:
    print(
        "Accession",
        "Contigs",
        "Total contigs",
        "Contig length",
        "Reference length",
        "Contig frac",
        "Cum frac",
        "Relevant edit dist",
        "Cigar string",
        "Edit frac",
        sep="\t",
        file=out,
    )
    sum_of_edit_dist = 0
    largest_no_contig = 0
    small_larg_cover = 1.0
    small_larg_cover_name = "asd"
    smallest_cum_frac = 1.0

    # for bam_file in snakemake_input:
    for bam_file, ref_fasta in zip(snakemake.input.bams, snakemake.input.refs):
        current_contig = 1

        ref_fasta = pysam.FastaFile(ref_fasta)

        with pysam.AlignmentFile(bam_file, "rb") as samfile:
            total_contigs = samfile.count()
            accession = samfile.get_reference_name(0)

            largest_no_contig = (
                total_contigs
                if total_contigs > largest_no_contig
                else largest_no_contig
            )

        with pysam.AlignmentFile(bam_file, "rb") as samfile:
            ref_lengths = samfile.lengths[0]
            largest_coverage = 0
            cum_frac = 0
            for read in samfile.fetch():
                query_alignment_length = read.query_alignment_length
                frac = round(query_alignment_length / ref_lengths, 2)

                # Calculate edit distance from CIGAR string, because NM tag counts matching Ns as edits.
                edit = get_edit_dist(read, ref_fasta)

                cum_frac += frac
                cum_frac = round(cum_frac, 2)
                sum_of_edit_dist = sum_of_edit_dist + edit
                edit_frac = round(edit / query_alignment_length, 5)
                largest_coverage = frac if frac > largest_coverage else largest_coverage

                print(
                    accession,
                    current_contig,
                    total_contigs,
                    query_alignment_length,
                    ref_lengths,
                    frac,
                    cum_frac,
                    edit,
                    read.cigarstring,
                    edit_frac,
                    sep="\t",
                    file=out,
                )

                current_contig += 1
            smallest_cum_frac = (
                cum_frac if smallest_cum_frac > cum_frac else smallest_cum_frac
            )
            small_larg_cover = (
                largest_coverage
                if small_larg_cover > largest_coverage
                else small_larg_cover
            )
            small_larg_cover_name = (
                accession
                if small_larg_cover >= largest_coverage
                else small_larg_cover_name
            )

    print("Largest number of contigs")
    print(largest_no_contig)
    print("Smallest largest coverage in " + str(small_larg_cover_name))
    print(small_larg_cover)
    print("Sum of edit distances")
    print(sum_of_edit_dist)
    print("Smallest cum. fraction of contigs")
    print(smallest_cum_frac)
    print("Largest number of contigs", file=out)
    print(largest_no_contig, file=out)
    print("Smallest largest coverage in " + str(small_larg_cover_name), file=out)
    print(small_larg_cover, file=out)
    print("Sum of edit distances", file=out)
    print(sum_of_edit_dist, file=out)
    print("Smallest cum. fraction of contigs", file=out)
    print(smallest_cum_frac, file=out)

import pandas as pd

# Function to create a bedpe file from a bed file
bed_list = []
with open(snakemake.input[0]) as f:
    line = f.readlines()
    for name in line:
        bed_list.append(name.split())
df_bed = pd.DataFrame(
    bed_list, columns=["chrom", "start", "end", "name", "score", "strand"]
)
df_sense = df_bed.loc[df_bed["strand"] == "+"]
df_antisense = df_bed.loc[df_bed["strand"] == "-"]
# The dataframes for the sense and antisense strands need to be set to the same index so they can be merged again for the bedpe file
df_sense.reset_index(inplace=True)
df_antisense.reset_index(inplace=True)
data = [
    df_sense["chrom"],
    df_sense["start"],
    df_sense["end"],
    df_antisense["chrom"],
    df_antisense["start"],
    df_antisense["end"],
]
headers = ["chrom1", "start1", "end1", "chrom2", "start2", "end2"]
df_bedpe = pd.concat(data, axis=1, keys=headers)
df_bedpe.to_csv(snakemake.output[0], header=None, index=None, sep="\t", mode="a")

import sys

sys.stderr = open(snakemake.log[0], "w")

import pandas as pd

all_variants = []
for path in snakemake.input:
    all_variants.append(pd.read_csv(path, sep="\t"))

pd.concat(all_variants).to_csv(snakemake.output[0], index=False, sep="\t")

import sys

from pysam import VariantFile

# sys.stderr = open(snakemake.log[0], "w")


tests_cases = []

for path, pos, variant, test_case in zip(
    snakemake.input,
    snakemake.params.poses,
    snakemake.params.variants,
    snakemake.params.test_cases,
):
    with VariantFile(path, "rb") as variant_file:
        for record in variant_file.fetch():
            if int(pos) == record.pos:
                tests_cases.append(test_case)

with open(snakemake.output[0], "w") as f:
    for path in tests_cases:
        print(f"{path}", file=f)

import sys

sys.stderr = open(snakemake.log[0], "w")

import pandas as pd
from pysam import VariantFile

AA_ALPHABET_TRANSLATION = {
    "Gly": "G",
    "Ala": "A",
    "Leu": "L",
    "Met": "M",
    "Phe": "F",
    "Trp": "W",
    "Lys": "K",
    "Gln": "Q",
    "Glu": "E",
    "Ser": "S",
    "Pro": "P",
    "Val": "V",
    "Ile": "I",
    "Cys": "C",
    "Tyr": "Y",
    "His": "H",
    "Arg": "R",
    "Asn": "N",
    "Asp": "D",
    "Thr": "T",
}


def phred_to_prob(phred):
    if phred is None:
        return 0
    return 10 ** (-phred / 10)


def get_AA_variant(record):
    for ann in record.info["ANN"]:
        ann = ann.split("|")
        hgvsp = ann[11]
        feature = ann[3]
        if hgvsp:
            _enssast_id, alteration = hgvsp.split(":", 1)
            _prefix, alteration = alteration.split(".", 1)
            for triplet, amino in AA_ALPHABET_TRANSLATION.items():
                alteration = alteration.replace(triplet, amino)
            hgvsp = f"{feature}:{alteration}"
    return hgvsp


def extract_data(variant_file: VariantFile):
    variants = []
    for record in variant_file.fetch():
        prob_not_present = phred_to_prob(record.info["PROB_ABSENT"][0]) + phred_to_prob(
            record.info["PROB_ARTIFACT"][0]
        )
        variant = get_AA_variant(record)
        if variant != "":
            variants.append(
                {
                    "chrom": record.chrom,
                    "pos": record.pos,
                    "variant": variant,
                    # "ANN" : record.info["ANN"],
                    "vaf": record.samples[0]["AF"][0],
                    "prob_not_present": prob_not_present,
                }
            )
    data = pd.DataFrame(variants)
    data["prob_present"] = 1 - data["prob_not_present"]
    data.drop(columns="prob_not_present", inplace=True)
    data.drop_duplicates(subset=["chrom", "pos", "variant"], inplace=True)
    return data


found_variants_list = []
for illumina_bcf, ont_bcf in zip(snakemake.input.illumina_bcf, snakemake.input.ont_bcf):
    with VariantFile(illumina_bcf, "rb") as illumina_file:
        illumina_variants = extract_data(illumina_file)

    with VariantFile(ont_bcf, "rb") as ont_file:
        ont_variants = extract_data(ont_file)

    found_variants = pd.merge(
        illumina_variants,
        ont_variants,
        how="outer",
        left_on=["chrom", "pos", "variant"],
        right_on=["chrom", "pos", "variant"],
        suffixes=("_illumina", "_ont"),
    )
    found_variants_list.append(found_variants)

found_variants = pd.concat(found_variants_list)

# found_variants.fillna({'vaf_illumina':0.0, 'vaf_ont':0.0, 'prob_present_illumina':0.0, 'prob_present_ont':0.0}, inplace=True)
found_variants["difference"] = found_variants["vaf_illumina"].astype(
    float
) - found_variants["vaf_ont"].astype(float)
found_variants.sort_values(by="difference", inplace=True, ascending=False)

found_variants.drop(columns=["difference"], inplace=True)
found_variants["test_case"] = snakemake.wildcards.test_case
found_variants.to_csv(snakemake.output[0], index=False, sep="\t")