text:
using DataFrames, CSV, Glob, Logging, StatsBase

"""
    Nucleotide_Essentials.potential_mismatches

Returns an ```Vector{Any}``` of potential barcodes with a single nucleotide change, including both deletions and substitutions

Supported keyword arguments include:

* ```Path::String```: The full or relative path to a .fastq file
* ```mismatch::Int64```: The number of altered nucleotides to include (1 is only supported at this time)

# Example:
```julia
potential_mismatches("GCGT", 1)
17-element Vector{Any}:
"GCGT"
"CCGT"
"ACGT"
"TCGT"
"GGGT"
"GAGT"
"GTGT"
"GCCT"
"GCAT"
"GCTT"
"GCGG"
"GCGC"
"GCGA"
"CGT"
"GGT"
"GCT"
"GCG"
```
"""
function potential_mismatches(og_barcode::Union{String15, String}, mismatch::Int64)
    potentials = []
    barcode = collect(og_barcode)
    if mismatch == 1
        for spot in 1:length(barcode)
            barcode = collect(og_barcode)
            barcode[spot] = 'G'
            push!(potentials, String(barcode))
            barcode = collect(og_barcode)
            barcode[spot] = 'C'
            push!(potentials, String(barcode))
            barcode = collect(og_barcode)
            barcode[spot] = 'A'
            push!(potentials, String(barcode))
            barcode = collect(og_barcode)
            barcode[spot] = 'T'
            push!(potentials, String(barcode))
        end
        for spot in 1:length(barcode)
            # Now handle potential deletions in the barcode
            barcode = collect(og_barcode)
            deleteat!(barcode, spot)
            del_barcode = join(barcode)
            push!(potentials, String(barcode))
        end
        potentials = unique(potentials)
    end
    return potentials
end

"""
    Nucleotide_Essentials.reverse_complement

Takes a string of nucleotide bases and returns the reverse complement of that string. Accepts inputs of ```String``` and ```SubString{String}``` (input from a FastqRecord)

Supported keyword arguments include:

* ```sequence::Union{String, SubString{String}}```: A string sequence of nucleotide bases or sequence entry from a FastqRecord

# Example:
```julia
reverse_complement("ATCGT")
"ACGAT"
```
"""
function reverse_complement(sequence::Union{String, SubString{String}})
    rc = reverse(replace(sequence,
                "C" => "G",
                "G" => "C",
                "A" => "T",
                "T" => "A"))
end

"""
    Nucleotide_Essentials.demultiplex_se

Compares a list of provided barcodes with the provided multiplexed reads and separates the reads into individual .fastq files. If a barcode is found within the read, the barcode is removed from the sequence. The quality data of the reads is preserved and written to the outputted .fastq file. If a barcode is not found, the sequnce and quality is written to the unassigned .fastq file unchanged.

The mapping file must be either a .csv or .txt file with two columns. The first column heading must be `SampleID` and the second column heading must be `BarcodeSequence`.

EXAMPLE MAPPING FILE:

| SampleID | BarcodeSequence |
|----------|-----------------|
| Sample1  | Barcode1        |
| Sample2  | Barcode2        |
| Sample3  | Barcode3        |
| Sample4  | Barcode4        |
| Sample5  | Barcode5        |
| Sample6  | Barcode6        |
| Sample7  | Barcode7        |
| Sample8  | Barcode8        |

Supported keyword arguments include:

* ```R1::String```: Path to multiplexed reads
* ```Map::String```: Path to the mapping file
* ```mismatch::Int64=0``` (optional): Number of allowed mismatches in barcode. Potential options include 0 or 1. If 1 mismatch, computation time will significantly increase. Default is to allow for 0 mismatches (exact matches only).
* ```debug::Bool=false``` (optional): If true, a log file will be created and debugging data will be printed while the function is running (default is false).

# Example:
```julia
demultiplex_se("multiplexreads.fastq", "mapping_file.fastq")
```
"""
function demultiplex_se(
    R1::String,
    Map::String,
    mismatch::Int64=0,
    debug::Bool=false
    )

    @debug("Making output directory")
    Dir = string(rsplit(abspath(R1), "/", limit = 2)[1])
    Out = "/demultiplex_output"
    mkdir(string(Dir, Out))

    if debug === true
        # Set up location for the log
        ENV["JULIA_DEBUG"] = Main
        if isdir(string(Dir, "/logs")) === false
            mkdir(string(Dir, "/logs"))
        end
        io = open(string(Dir, "/logs/log_demultiplex.txt"), "w+")
        logger = ConsoleLogger(io)
        global_logger(logger)
    else
        Logging.disable_logging(Logging.Info)
    end

    pre_dmx = readFastq(R1) # these are our forward multiplexed reads
    seqs = pre_dmx.sequence
    qual = pre_dmx.quality
    og_name = pre_dmx.ID

    # set up our mapping file - this is our mapping file in txt format
    if endswith(Map, ".csv") == true
        @debug("Reading in our mapping file - it is a .csv file")
        mapping = CSV.read(Map, DataFrame)
        barcodes = Array(mapping[!, :BarcodeSequence])
        barcode_lengths = []
        for b in barcodes
            hold_length = length(b)
            push!(barcode_lengths, hold_length)
        end
        barcode_length = median(barcode_lengths)
        IDs = Array(mapping[!, :SampleID])

        # Make sure the number of barcodes matches the number of sample IDs
        if length(unique(barcodes)) != length(unique(IDs))
            error("The number of barcodes and unique samples should match, check mapping file for errors")
        end
    end

    # set up our mapping file - this is our mapping file in txt format
    if endswith(Map, ".txt") == true
        @debug("Reading in our mapping file - it is a .txt file")
        mapping = CSV.read(Map, DataFrame)
        barcodes = Array(mapping[!, :BarcodeSequence])
        barcode_lengths = []
        for b in barcodes
            hold_length = length(b)
            push!(barcode_lengths, hold_length)
        end
        barcode_length = median(barcode_lengths)
        IDs = Array(mapping[!, :SampleID])

        # Make sure the number of barcodes matches the number of sample IDs
        if length(unique(barcodes)) != length(unique(IDs))
            error("The number of barcodes and unique samples should match, check mapping file for errors")
        end
    end

    # dmuxing
    @debug("Beginning demultiplexing algorithm")
    for s in 1:length(seqs)
        # demultiplexing not allowing for any mismatches
        if mismatch == 0
            @debug(string("Beginning demultiplexing for sequence #", s, " of ", length(seqs)))
            lookat = Integer(barcode_length + 5)
            this_seq = first(seqs[s], lookat)
            for b in 1:length(barcodes)
                @debug(string("Querying against barcode #", b))
                if occursin(barcodes[b], this_seq)
                    sample = string(IDs[b])
                    file_end = string("_R1.fastq")
                    if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                        @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                        out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                        write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                        break
                    else
                        @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                        out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                        iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                        write(iostream, out_string2)
                        close(iostream)
                        break
                    end
                elseif b == length(barcodes)
                    @debug("Beginning last barcode")
                    if occursin(barcodes[b], this_seq)
                        sample = string(IDs[b])
                        file_end = string("_R1.fastq")
                        if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                            break
                        else
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                            write(iostream, out_string2)
                            close(iostream)
                            break
                        end
                    else
                        if isfile(string(Dir, "/", Out, "/unassigned_R1.fastq")) == false
                            @debug("No unassigned Fastq file detected - Creating file")
                            @debug("Reached last barcode - writing to unassigned Fastq file.")
                            out_string3 = string("@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                            write(string(Dir, "/", Out, "/unassigned_R1.fastq"), out_string3)
                        else
                            @debug("Reached last barcode - writing to unassigned Fastq file.")
                            out_string4 = string("\n@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                            iostream = open(string(Dir, "/", Out, "/unassigned_R1.fastq"), "a")
                            write(iostream, out_string4)
                            close(iostream)
                        end
                    end
                else
                    continue
                end
            end

            # demultiplexing allowing for mismatches
        else
            @debug(string("Beginning demultiplexing for sequence #", s, " of ", length(seqs)))
            lookat = Integer(barcode_length + 5)
            this_seq = first(seqs[s], lookat)
            for b in 1:length(barcodes)
                @debug(string("Building barcode pool for barcode #", b, ": ", barcodes[b]))
                temp_barcode_pool = potential_mismatches(barcodes[b], mismatch)
                @debug("Finished building potential barcode pool")
                for drop in 1:length(temp_barcode_pool)
                    @debug(string("Querying for variant #", drop, " from barcode #", b))
                    if occursin(temp_barcode_pool[drop], this_seq)
                        sample = string(IDs[b])
                        file_end = string("_R1.fastq")
                        if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                            break
                        else
                            @debug("Else after if 1.1")
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                            write(iostream, out_string2)
                            close(iostream)
                            break
                        end
                    elseif drop == length(temp_barcode_pool)
                        @debug("Beginning last variant")
                        if occursin(temp_barcode_pool[drop], this_seq)
                            sample = string(IDs[b])
                            file_end = string("_R1.fastq")
                            if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                                @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                                out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                                write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                                break
                            else
                                @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                                out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                                iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                                write(iostream, out_string2)
                                close(iostream)
                                break
                            end
                        else
                            if b == length(barcodes)
                                @debug("Reached last barcode - writing to unassigned Fastq file.")
                                if isfile(string(Dir, "/", Out, "/unassigned_R1.fastq")) == false
                                    @debug("No unassigned Fastq file detected - Creating file")
                                    out_string3 = string("@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                                    write(string(Dir, "/", Out, "/unassigned_R1.fastq"), out_string3)
                                else
                                    @debug("Reached last barcode - writing to unassigned Fastq file.")
                                    out_string4 = string("\n@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                                    iostream = open(string(Dir, "/", Out, "/unassigned_R1.fastq"), "a")
                                    write(iostream, out_string4)
                                    close(iostream)
                                end
                            else
                                @debug("Reached last variant in barcode pool - moving on to the next pool")
                                continue
                            end
                        end
                    else
                        @debug("I'm not sure?")
                        continue
                    end
                end
            end
        end
    end

    # output statistics on demultiplexed samples
    sampleholds = []
    reads = []
    files = glob("*.fastq", string(Dir, "/", Out))
    for file in files
        push!(sampleholds, rsplit(file, "/", limit=2)[2])
        push!(reads, (countlines(file) / 4))
    end
    stats = DataFrame(Filename=sampleholds, AssignedReads=reads)
    CSV.write(string(Dir, "/", Out, "/demultiplex_stats.csv"), stats)

    if debug === true
        # stop logger
        close(io)
    end
end # end function demultiplex_se

"""
    Nucleotide_Essentials.demultiplex_pe

Compares a list of provided barcodes with the provided paired-end multiplexed reads and separates the reads into individual .fastq files. If a barcode is found within R1 reads, the barcode is removed from the sequence. The quality data of the reads is preserved and written to the outputted .fastq file. If a barcode is not found, the sequnce and quality is written to the R1 unassigned .fastq file unchanged. If a barcode is found within R2 reads, the barcode is removed from the sequence. The quality data of the reads is preserved and written to the outputted .fastq file. If a barcode is not found, the sequnce and quality is written to the R2 unassigned .fastq file unchanged.

*Dual-indexed reads are not _yet_ supported*

The mapping file must be either a .csv or .txt file with two columns. The first column heading must be `SampleID` and the second column heading must be `BarcodeSequence`.

EXAMPLE MAPPING FILE:

| SampleID | BarcodeSequence |
|----------|-----------------|
| Sample1  | Barcode1        |
| Sample2  | Barcode2        |
| Sample3  | Barcode3        |
| Sample4  | Barcode4        |
| Sample5  | Barcode5        |
| Sample6  | Barcode6        |
| Sample7  | Barcode7        |
| Sample8  | Barcode8        |

Supported keyword arguments include:

* ```R1::String```: Path to forward multiplexed reads
* ```R2::String```: Path to reverse multiplexed reads
* ```Map::String```: Path to the mapping file
* ```mismatch::Int64=0``` (optional): Number of allowed mismatches in barcode. Potential options include 0 or 1. If 1 mismatch, computation time will significantly increase. Default is to allow for 0 mismatches (exact matches only).
* ```debug::Bool=false``` (optional): If true, a log file will be created and debugging data will be printed while the function is running (default is false).

# Example:
```julia
demultiplex_pe("forward_multiplexreads.fastq", "reverse_multiplexreads.fastq", "mapping_file.fastq")
```
"""
function demultiplex_pe(
    R1::String,
    R2::String,
    Map::String,
    mismatch::Int64=0,
    debug::Bool=false
    )

    @debug("Making output directory")
    Dir = string(rsplit(abspath(R1), "/", limit = 2)[1])
    Out = "/demultiplex_output"
    mkdir(string(Dir, Out))

    if debug === true
        # Set up location for the log
        ENV["JULIA_DEBUG"] = Main
        if isdir(string(Dir, "/logs")) === false
            mkdir(string(Dir, "/logs"))
        end
        io = open(string(Dir, "/logs/log_demultiplex.txt"), "w+")
        logger = ConsoleLogger(io)
        global_logger(logger)
    else
        Logging.disable_logging(Logging.Info)
    end

    pre_dmx = readFastq(R1) # these are our forward multiplexed reads
    pre_dmx2 = readFastq(R2) # these are our forward multiplexed reads
    seqs = pre_dmx.sequence
    seqs2 = pre_dmx2.sequence
    qual = pre_dmx.quality
    qual2 = pre_dmx2.quality
    og_name = pre_dmx.ID
    og_name2 = pre_dmx2.ID

    # set up our mapping file - this is our mapping file in txt format
    if endswith(Map, ".csv") == true
        @debug("Reading in our mapping file - it is a .csv file")
        mapping = CSV.read(Map, DataFrame)
        barcodes = Array(mapping[!, :BarcodeSequence])
        barcode_lengths = []
        for b in barcodes
            hold_length = length(b)
            push!(barcode_lengths, hold_length)
        end
        barcode_length = median(barcode_lengths)
        IDs = Array(mapping[!, :SampleID])

        # Make sure the number of barcodes matches the number of sample IDs
        if length(unique(barcodes)) != length(unique(IDs))
            error("The number of barcodes and unique samples should match, check mapping file for errors")
        end
    end

    # set up our mapping file - this is our mapping file in txt format
    if endswith(Map, ".txt") == true
        @debug("Reading in our mapping file - it is a .txt file")
        mapping = CSV.read(Map, DataFrame)
        barcodes = Array(mapping[!, :BarcodeSequence])
        barcode_lengths = []
        for b in barcodes
            hold_length = length(b)
            push!(barcode_lengths, hold_length)
        end
        barcode_length = median(barcode_lengths)
        IDs = Array(mapping[!, :SampleID])

        # Make sure the number of barcodes matches the number of sample IDs
        if length(unique(barcodes)) != length(unique(IDs))
            error("The number of barcodes and unique samples should match, check mapping file for errors")
        end
    end

    # dmuxing forward reads
    @debug("Beginning demultiplexing algorithm for forward reads")
    for s in 1:length(seqs)
        # demultiplexing not allowing for any mismatches
        if mismatch == 0
            @debug(string("Beginning demultiplexing for sequence #", s, " of ", length(seqs)))
            lookat = Integer(barcode_length + 5)
            this_seq = first(seqs[s], lookat)
            for b in 1:length(barcodes)
                @debug(string("Querying against barcode #", b))
                if occursin(barcodes[b], this_seq)
                    sample = string(IDs[b])
                    file_end = string("_R1.fastq")
                    if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                        @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                        out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                        write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                        break
                    else
                        @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                        out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                        iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                        write(iostream, out_string2)
                        close(iostream)
                        break
                    end
                elseif b == length(barcodes)
                    @debug("Beginning last barcode")
                    if occursin(barcodes[b], this_seq)
                        sample = string(IDs[b])
                        file_end = string("_R1.fastq")
                        if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                            break
                        else
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                            write(iostream, out_string2)
                            close(iostream)
                            break
                        end
                    else
                        if isfile(string(Dir, "/", Out, "/unassigned_R1.fastq")) == false
                            @debug("No unassigned Fastq file detected - Creating file")
                            @debug("Reached last barcode - writing to unassigned Fastq file.")
                            out_string3 = string("@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                            write(string(Dir, "/", Out, "/unassigned_R1.fastq"), out_string3)
                        else
                            @debug("Reached last barcode - writing to unassigned Fastq file.")
                            out_string4 = string("\n@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                            iostream = open(string(Dir, "/", Out, "/unassigned_R1.fastq"), "a")
                            write(iostream, out_string4)
                            close(iostream)
                        end
                    end
                else
                    continue
                end
            end

            # demultiplexing allowing for mismatches
        else
            @debug(string("Beginning demultiplexing for sequence #", s, " of ", length(seqs)))
            lookat = Integer(barcode_length + 5)
            this_seq = first(seqs[s], lookat)
            for b in 1:length(barcodes)
                @debug(string("Building barcode pool for barcode #", b, ": ", barcodes[b]))
                temp_barcode_pool = potential_mismatches(barcodes[b], mismatch)
                @debug("Finished building potential barcode pool")
                for drop in 1:length(temp_barcode_pool)
                    @debug(string("Querying for variant #", drop, " from barcode #", b))
                    if occursin(temp_barcode_pool[drop], this_seq)
                        sample = string(IDs[b])
                        file_end = string("_R1.fastq")
                        if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                            break
                        else
                            @debug("Else after if 1.1")
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                            iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                            write(iostream, out_string2)
                            close(iostream)
                            break
                        end
                    elseif drop == length(temp_barcode_pool)
                        @debug("Beginning last variant")
                        if occursin(temp_barcode_pool[drop], this_seq)
                            sample = string(IDs[b])
                            file_end = string("_R1.fastq")
                            if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                                @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                                out_string = string("@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                                write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                                break
                            else
                                @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                                out_string2 = string("\n@", IDs[b], "\n", chop(seqs[s], head = Integer(barcode_length)), "\n+\n", chop(qual[s], head = Integer(barcode_length)))
                                iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                                write(iostream, out_string2)
                                close(iostream)
                                break
                            end
                        else
                            if b == length(barcodes)
                                @debug("Reached last barcode - writing to unassigned Fastq file.")
                                if isfile(string(Dir, "/", Out, "/unassigned_R1.fastq")) == false
                                    @debug("No unassigned Fastq file detected - Creating file")
                                    out_string3 = string("@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                                    write(string(Dir, "/", Out, "/unassigned_R1.fastq"), out_string3)
                                else
                                    @debug("Reached last barcode - writing to unassigned Fastq file.")
                                    out_string4 = string("\n@", og_name[s], "\n", seqs[s], "\n+\n", qual[s])
                                    iostream = open(string(Dir, "/", Out, "/unassigned_R1.fastq"), "a")
                                    write(iostream, out_string4)
                                    close(iostream)
                                end
                            else
                                @debug("Reached last variant in barcode pool - moving on to the next pool")
                                continue
                            end
                        end
                    else
                        @debug("I'm not sure?")
                        continue
                    end
                end
            end
        end
    end

    # dmuxing reverse reads
    @debug("Beginning demultiplexing algorithm for reverse reads")
    for s in 1:length(seqs2)
        # demultiplexing not allowing for any mismatches
        if mismatch == 0
            @debug(string("Beginning demultiplexing for sequence #", s, " of ", length(seqs2)))
            lookat = Integer(barcode_length + 5)
            this_line = reverse_complement(seqs2[s])
            this_seq = (first(this_line, lookat))
            for b in 1:length(barcodes)
                @debug(string("Querying against barcode #", b))
                if occursin(barcodes[b], this_seq)
                    sample = string(IDs[b])
                    file_end = string("_R2.fastq")
                    if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                        @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                        out_string = string("@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                        write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                        break
                    else
                        @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                        out_string2 = string("\n@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                        iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                        write(iostream, out_string2)
                        close(iostream)
                        break
                    end
                elseif b == length(barcodes)
                    @debug("Beginning last barcode")
                    if occursin(barcodes[b], this_seq)
                        sample = string(IDs[b])
                        file_end = string("_R2.fastq")
                        if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string = string("@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                            write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                            break
                        else
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string2 = string("\n@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                            iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                            write(iostream, out_string2)
                            close(iostream)
                            break
                        end
                    else
                        if isfile(string(Dir, "/", Out, "/unassigned_R2.fastq")) == false
                            @debug("No unassigned Fastq file detected - Creating file")
                            @debug("Reached last barcode - writing to unassigned Fastq file.")
                            out_string3 = string("@", og_name2[s], "\n", seqs2[s], "\n+\n", qual2[s])
                            write(string(Dir, "/", Out, "/unassigned_R2.fastq"), out_string3)
                        else
                            @debug("Reached last barcode - writing to unassigned Fastq file.")
                            out_string4 = string("\n@", og_name2[s], "\n", seqs2[s], "\n+\n", qual2[s])
                            iostream = open(string(Dir, "/", Out, "/unassigned_R2.fastq"), "a")
                            write(iostream, out_string4)
                            close(iostream)
                        end
                    end
                else
                    continue
                end
            end

            # demultiplexing allowing for mismatches
        else
            @debug(string("Beginning demultiplexing for sequence #", s, " of ", length(seqs2)))
            lookat = Integer(barcode_length + 5)
            this_line = reverse_complement(seqs2[s])
            this_seq = (first(this_line, lookat))
            for b in 1:length(barcodes)
                @debug(string("Building barcode pool for barcode #", b, ": ", barcodes[b]))
                temp_barcode_pool = potential_mismatches(barcodes[b], mismatch)
                @debug("Finished building potential barcode pool")
                for drop in 1:length(temp_barcode_pool)
                    @debug(string("Querying for variant #", drop, " from barcode #", b))
                    if occursin(temp_barcode_pool[drop], this_seq)
                        sample = string(IDs[b])
                        file_end = string("_R2.fastq")
                        if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string = string("@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                            write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                            break
                        else
                            @debug("Else after if 1.1")
                            @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                            out_string2 = string("\n@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                            iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                            write(iostream, out_string2)
                            close(iostream)
                            break
                        end
                    elseif drop == length(temp_barcode_pool)
                        @debug("Beginning last variant")
                        if occursin(temp_barcode_pool[drop], this_seq)
                            sample = string(IDs[b])
                            file_end = string("_R2.fastq")
                            if isfile(string(Dir, "/", Out, "/", sample, file_end)) == false
                                @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                                out_string = string("@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                                write(string(Dir, "/", Out, "/", sample, file_end), out_string)
                                break
                            else
                                @debug(string("Barcode found - writing to fastq file for ", IDs[b]))
                                out_string2 = string("\n@", IDs[b], "\n", (chop(seqs2[s], tail = Integer(barcode_length))), "\n+\n", (chop(qual2[s], tail = Integer(barcode_length))))
                                iostream = open(string(Dir, "/", Out, "/", sample, file_end), "a")
                                write(iostream, out_string2)
                                close(iostream)
                                break
                            end
                        else
                            if b == length(barcodes)
                                @debug("Reached last barcode - writing to unassigned Fastq file.")
                                if isfile(string(Dir, "/", Out, "/unassigned_R2.fastq")) == false
                                    @debug("No unassigned Fastq file detected - Creating file")
                                    out_string3 = string("@", og_name2[s], "\n", seqs2[s], "\n+\n", qual2[s])
                                    write(string(Dir, "/", Out, "/unassigned_R2.fastq"), out_string3)
                                else
                                    @debug("Reached last barcode - writing to unassigned Fastq file.")
                                    out_string4 = string("\n@", og_name2[s], "\n", seqs2[s], "\n+\n", qual2[s])
                                    iostream = open(string(Dir, "/", Out, "/unassigned_R2.fastq"), "a")
                                    write(iostream, out_string4)
                                    close(iostream)
                                end
                            else
                                @debug("Reached last variant in barcode pool - moving on to the next pool")
                                continue
                            end
                        end
                    else
                        @debug("I'm not sure?")
                        continue
                    end
                end
            end
        end
    end

    # output statistics on demultiplexed samples
    sampleholds = []
    reads = []
    files = glob("*.fastq", string(Dir, "/", Out))
    for file in files
        push!(sampleholds, rsplit(file, "/", limit=2)[2])
        push!(reads, (countlines(file) / 4))
    end
    stats = DataFrame(Filename=sampleholds, AssignedReads=reads)
    CSV.write(string(Dir, "/", Out, "/demultiplex_stats.csv"), stats)

    if debug === true
        # stop logger
        close(io)
    end
end # end function demultiplex_pe


text:
# Custom Layers needed for ViT
"""
    MultiHeadAttention(in_planes::Int, number_heads::Int; qkv_bias::Bool=false, attention_dropout_rate::T=0.0f0, projection_dropout_rate::T=0.0f0) where {T}

Multi-head self-attention layer
"""
struct MultiHeadAttention{Q, A, P} <:
       Lux.AbstractExplicitContainerLayer{(:qkv_layer, :attention_dropout, :projection)}
    number_heads::Int
    qkv_layer::Q
    attention_dropout::A
    projection::P
end

function MultiHeadAttention(in_planes::Int,
                            number_heads::Int;
                            qkv_bias::Bool=false,
                            attention_dropout_rate::T=0.0f0,
                            projection_dropout_rate::T=0.0f0) where {T}
    @assert in_planes % number_heads==0 "`in_planes` should be divisible by `number_heads`"
    qkv_layer = Dense(in_planes, in_planes * 3; bias=qkv_bias)
    attention_dropout = Dropout(attention_dropout_rate)
    projection = Chain(Dense(in_planes, in_planes), Dropout(projection_dropout_rate))

    return MultiHeadAttention(number_heads, qkv_layer, attention_dropout, projection)
end

function (m::MultiHeadAttention)(x::AbstractArray{T, 3}, ps, st) where {T}
    nfeatures, seq_len, batch_size = size(x)

    x_reshaped = reshape(x, nfeatures, seq_len * batch_size)
    qkv, st_qkv = m.qkv_layer(x_reshaped, ps.qkv_layer, st.qkv_layer)
    qkv_reshaped = reshape(qkv, nfeatures ÷ m.number_heads, m.number_heads, seq_len,
                           3 * batch_size)
    query, key, value = fast_chunk(qkv_reshaped, Val(3), Val(4))

    scale = convert(T, sqrt(size(query, 1) / m.number_heads))
    key_reshaped = reshape(permutedims(key, (2, 1, 3, 4)), m.number_heads,
                           nfeatures ÷ m.number_heads, seq_len * batch_size)
    query_reshaped = reshape(query, nfeatures ÷ m.number_heads, m.number_heads,
                             seq_len * batch_size)

    attention = softmax(batched_mul(query_reshaped, key_reshaped) .* scale)
    attention, st_attention = m.attention_dropout(attention, ps.attention_dropout,
                                                  st.attention_dropout)

    value_reshaped = reshape(value, nfeatures ÷ m.number_heads, m.number_heads,
                             seq_len * batch_size)
    pre_projection = reshape(batched_mul(attention, value_reshaped),
                             (nfeatures, seq_len, batch_size))
    y, st_projection = m.projection(reshape(pre_projection, size(pre_projection, 1), :),
                                    ps.projection, st.projection)

    st_ = (qkv_layer=st_qkv, attention=st_attention, projection=st_projection)
    return reshape(y, :, seq_len, batch_size), st_
end

"""
    ClassTokens(dim; init=Lux.zeros32)

Appends class tokens to an input with embedding dimension `dim` for use in many vision transformer namels.
"""
struct ClassTokens{I} <: Lux.AbstractExplicitLayer
    dim::Int
    init::I
end

ClassTokens(dim::Int; init=Lux.zeros32) = ClassTokens(dim, init)

Lux.initialparameters(rng::AbstractRNG, c::ClassTokens) = (token=c.init(rng, c.dim, 1, 1),)

_fill_like(y::AbstractArray{T, 3}) where {T} = fill!(similar(y, 1, 1, size(y, 3)), one(T))
ChainRulesCore.@non_differentiable _fill_like(y)

function (m::ClassTokens)(x::AbstractArray{T, 3}, ps, st) where {T}
    # Generic Alternative: Repeat is extremely inefficient on GPUs and even in general
    tokens = ps.token .* _fill_like(x)
    return hcat(tokens, x), st
end

"""
    ViPosEmbedding(embedsize, npatches; init = (rng, dims...) -> randn(rng, Float32, dims...))

Positional embedding layer used by many vision transformer-like namels.
"""
struct ViPosEmbedding{I} <: Lux.AbstractExplicitLayer
    embedding_size::Int
    number_patches::Int
    init::I
end

function ViPosEmbedding(embedding_size::Int, number_patches::Int;
                        init=(rng, dims...) -> randn(rng, Float32, dims...))
    return ViPosEmbedding(embedding_size, number_patches, init)
end

function Lux.initialparameters(rng::AbstractRNG, v::ViPosEmbedding)
    return (vectors=v.init(rng, v.embedding_size, v.number_patches),)
end

(v::ViPosEmbedding)(x, ps, st) = x .+ ps.vectors, st

# Helper Functions
"""
   transformer\\_encoder(in\\_planes, depth, number\\_heads; mlp\\_ratio = 4.0f0, dropout = 0.0f0)

Transformer as used in the base ViT architecture. ([reference](https://arxiv.org/abs/2010.11929)).

## Arguments
- `in_planes`: number of input channels
- `depth`: number of attention blocks
- `number_heads`: number of attention heads
- `mlp_ratio`: ratio of MLP layers to the number of input channels
- `dropout_rate`: dropout rate
"""
function transformer_encoder(in_planes, depth, number_heads; mlp_ratio=4.0f0,
                             dropout_rate=0.0f0)
    # Currently Lux doesn't have a LayerNorm implementation. Relying on layer_norm
    hidden_planes = floor(Int, mlp_ratio * in_planes)
    layers = [Chain(SkipConnection(Chain(layer_norm(in_planes),
                                         MultiHeadAttention(in_planes,
                                                            number_heads;
                                                            attention_dropout_rate=dropout_rate,
                                                            projection_dropout_rate=dropout_rate)),
                                   +),
                    SkipConnection(Chain(layer_norm(in_planes),
                                         Chain(Dense(in_planes => hidden_planes, gelu),
                                               Dropout(dropout_rate),
                                               Dense(hidden_planes => in_planes),
                                               Dropout(dropout_rate))),
                                   +)) for _ in 1:depth]
    return Chain(layers...)
end

function patch_embedding(imsize::Tuple{<:Int, <:Int}=(224, 224);
                         in_channels=3,
                         patch_size::Tuple{<:Int, <:Int}=(16, 16),
                         embed_planes=768,
                         norm_layer=in_planes -> NoOpLayer(),
                         flatten=true)
    im_width, im_height = imsize
    patch_width, patch_height = patch_size

    @assert (im_width % patch_width == 0) && (im_height % patch_height == 0)
    "Image dimensions must be divisible by the patch size."

    return Chain(Conv(patch_size, in_channels => embed_planes; stride=patch_size),
                 flatten ? flatten_spatial : identity,
                 norm_layer(embed_planes))
end

# ViT Implementation
function vision_transformer(;
                            imsize::Tuple{<:Int, <:Int}=(256, 256),
                            in_channels::Int=3,
                            patch_size::Tuple{<:Int, <:Int}=(16, 16),
                            embed_planes::Int=768,
                            depth::Int=6,
                            number_heads=16,
                            mlp_ratio=4.0f0,
                            dropout_rate=0.1f0,
                            embedding_dropout_rate=0.1f0,
                            pool::Symbol=:class,
                            num_classes::Int=1000,
                            kwargs...)
    @assert pool in (:class, :mean) "Pool type must be either :class (class token) or :mean (mean pooling)"
    number_patches = prod(imsize .÷ patch_size)

    return Chain(Chain(patch_embedding(imsize; in_channels, patch_size, embed_planes),
                       ClassTokens(embed_planes),
                       ViPosEmbedding(embed_planes, number_patches + 1),
                       Dropout(embedding_dropout_rate),
                       transformer_encoder(embed_planes, depth, number_heads; mlp_ratio,
                                           dropout_rate),
                       (pool == :class) ? SelectDim(2, 1) : seconddimmean),
                 Chain(layer_norm(embed_planes), Dense(embed_planes, num_classes, tanh)))
end

const VIT_CONFIGS = Dict(:tiny => (depth=12, embed_planes=192, number_heads=3),
                         :small => (depth=12, embed_planes=384, number_heads=6),
                         :base => (depth=12, embed_planes=768, number_heads=12),
                         :large => (depth=24, embed_planes=1024, number_heads=16),
                         :huge => (depth=32, embed_planes=1280, number_heads=16),
                         :giant => (depth=40, embed_planes=1408, number_heads=16,
                                    mlp_ratio=48 / 11),
                         :gigantic => (depth=48, embed_planes=1664, number_heads=16,
                                       mlp_ratio=64 / 13))

function vision_transformer(name::Symbol; kwargs...)
    assert_name_present_in(name, keys(VIT_CONFIGS))
    model = vision_transformer(; VIT_CONFIGS[name]..., kwargs...)
    return initialize_model(name, model; kwargs...)
end

ext: jl