shell:
    """
    #module load gsnap
    gmap_build -D {params.gmap_db_dir} -d {params.gmap_db} -s none {input.mt_genome_fasta} &> {log}
    """

shell:
    """
    cat {input.mt_genome_fasta} {input.n_genome_fasta} > {output.mt_n_fasta}
    gmap_build -D {params.gmap_db_dir} -d {params.gmap_db} -s none {output.mt_n_fasta} &> {log}
    # rm {input.mt_genome_fasta}_{input.n_genome_fasta}.fasta
    """

shell:
    """

    mkdir -p {params.outDir}
    fastqc -t {threads} -o {params.outDir} {input} &> {log}

    """

run:
    if seq_type == "pe":
        print("PE mode")
        shell("gsnap -D {params.gmap_db_dir} -d {params.gmap_db} -o {params.uncompressed_output} -A sam --gunzip --nofails --pairmax-dna=500 --query-unk-mismatch=1 {params.RG_tag} -n 1 -Q -O -t {threads} {input[0]} {input[1]} &> {log} && gzip {params.uncompressed_output} &>> {log}")
    if seq_type == "se":
        print("SE mode")
        shell("gsnap -D {params.gmap_db_dir} -d {params.gmap_db} -o {params.uncompressed_output} -A sam --gunzip --nofails --pairmax-dna=500 --query-unk-mismatch=1 {params.RG_tag} -n 1 -Q -O -t {threads} {input[0]} &> {log} && gzip {params.uncompressed_output} &>> {log}")
    elif seq_type == "both":
        print("PE + SE mode")
        shell("gsnap -D {params.gmap_db_dir} -d {params.gmap_db} -o {params.uncompressed_output} -A sam --gunzip --nofails --pairmax-dna=500 --query-unk-mismatch=1 {params.RG_tag} -n 1 -Q -O -t {threads} {input[0]} {input[1]} {input[2]} &> {log} && gzip {params.uncompressed_output} &>> {log}")

run:
    sam_to_fastq(samfile=input.outmt_sam, outmt1=output.outmt1,
                 outmt2=output.outmt2, outmt=output.outmt)

run:
    if os.path.isfile(input.outmt):
        shell("gsnap -D {params.gmap_db_dir} -d {params.gmap_db} -o {params.uncompressed_output} --gunzip -A sam --nofails --query-unk-mismatch=1 -O -t {threads} {input.outmt} &> {log.logS} && gzip {params.uncompressed_output} &>> {log.logS}")
    else:
        open(output.outS, 'a').close()

run:
    if os.path.isfile(input.outmt1):
        shell("gsnap -D {params.gmap_db_dir} -d {params.gmap_db} -o {params.uncompressed_output} --gunzip -A sam --nofails --query-unk-mismatch=1 -O -t {threads} {input.outmt1} {input.outmt2} &> {log.logP} && gzip {params.uncompressed_output} &>> {log.logP}")
    else:
        open(output.outP, 'a').close()

run:
    filter_alignments(outmt=input.outmt,
                      outS=input.outS,
                      outP=input.outP,
                      OUT=output.sam,
                      ref_mt_fasta=params.ref_mt_fasta)

shell:
    """
    samtools sort -o {output.sorted_bam} -T {params.TMP} {input.bam} &> {log}
    # samtools sort -o {output.sorted_bam} -T ${{TMP}} {input.bam}
    """

shell:
    """
    samtools merge {output.merged_bam} {input} &> {log}
    samtools index {output.merged_bam} {output.merged_bam_index}
    """

shell:
    """
    samtools faidx {input.mt_n_fasta} &> {log}
    """

run:
    shell("picard CreateSequenceDictionary R={input.mt_n_fasta} O={output.genome_dict}")

shell:
    """
    samtools mpileup -B -f {params.genome_fasta} -o {output.pileup} {input.merged_bam} &> {log}
    """

run:
    mt_table_data = pileup2mt_table(pileup=input.pileup, ref_fasta=params.ref_mt_fasta)
    write_mt_table(mt_table_data=mt_table_data, mt_table_file=output.mt_table)

run:
    # function (and related ones) from mtVariantCaller
    # vcf_dict = mtvcf_main_analysis(sam_file = input.sam, mtable_file = input.mt_table, name2 = wildcards.sample)
    tmp_sam = os.path.split(input.merged_bam)[1].replace(".bam", ".sam")
    shell("samtools view {merged_bam} > {tmp_dir}/{tmp_sam}".format(merged_bam=input.merged_bam,
                                                                    tmp_dir=params.TMP,
                                                                    tmp_sam=tmp_sam))
    vcf_dict = mtvcf_main_analysis(sam_file="{tmp_dir}/{tmp_sam}".format(tmp_dir=params.TMP,
                                                                         tmp_sam=tmp_sam),
                                   mtable_file=input.mt_table, name2=wildcards.sample)
    # ref_genome_mt will be used in the VCF descriptive field
    # seq_name in the VCF data
    seq_name = get_seq_name(params.ref_mt_fasta)
    VCF_RECORDS = VCFoutput(vcf_dict, reference=wildcards.ref_genome_mt,
                            seq_name=seq_name, vcffile=output.single_vcf)
    bed_output(VCF_RECORDS, seq_name=seq_name, bedfile=output.single_bed)
    # fasta output
    #contigs = pileup2mt_table(pileup=input.pileup, fasta=params.ref_mt_fasta, mt_table=in.mt_table)
    mt_table_data = pileup2mt_table(pileup=input.pileup,
                                    ref_fasta=params.ref_mt_fasta)
    gapped_fasta = mt_table_handle2gapped_fasta(mt_table_data=mt_table_data)
    contigs = gapped_fasta2contigs(gapped_fasta=gapped_fasta)
    fasta_output(vcf_dict=vcf_dict, ref_mt=params.ref_mt_fasta,
                 fasta_out=output.single_fasta, contigs=contigs)

run:
    shell("bcftools merge {input.single_vcf_list} -O v -o {output.merged_vcf}")