from __future__ import print_function, division
import sys
import os
import signal

if sys.version_info < (3, 0):
    from StringIO import StringIO
else:
    from io import StringIO

try:
    from Bio import AlignIO, SeqIO
except ImportError:
    sys.stderr.write("Error! BioPython is not detected!\n")
    sys.exit(1)

import random
import subprocess
import re
from collections import namedtuple, defaultdict
import argparse
import json

'''
BLCA Core annotation tool
'''


class SamEntry(object):
    def __init__(self, raw_row):
        self.qname = raw_row[0]
        self.flag = raw_row[1]
        self.rname = raw_row[2]
        self.pos = raw_row[3]
        self.mapq = raw_row[4]
        self.cigar = raw_row[5]
        self.rnext = raw_row[6]
        self.pnext = raw_row[7]
        self.tlen = raw_row[8]
        self.seq = raw_row[9]
        self.qual = raw_row[10]
        self.alignment_scores = [int(score.split(':')[-1]) for score in raw_row[11:13]]
        # sometimes it is a column early  because there is no second best match
        if raw_row[17][:2] == 'MD':
            self.md_z_flags = raw_row[17].split(':')[-1]
        else:
            self.md_z_flags = raw_row[18].split(':')[-1]

        total_match_count, total_count = self.calulate_match_count(self.md_z_flags)
        self.identity_ratio_old = total_match_count / total_count
        self.identity_ratio = total_match_count / float(len(self.seq))

        self.soft_clipped_ratio = self.cigar_total_soft_clipping() / total_count

        self.total_match = total_match_count

    def calulate_match_count(self, md_z_flags):
        tokenizer = re.compile(r'(\d+)|(\^[A-Z])|([A-Z])')
        total_match_count = 0.0
        total_mismatch_count = 0.0
        for item in tokenizer.finditer(md_z_flags):
            match_count, deletion, mismatch = item.groups()
            if match_count:
                total_match_count += int(match_count)
            if deletion:
                total_mismatch_count += 1
            if mismatch:
                total_mismatch_count += 1

        return total_match_count, (total_match_count + total_mismatch_count)

    # I am not completely sure if this is the logic that you want to use for checking unmapped at the ends
    # This function returns the max S at either end of the CIGAR score, 0  if there is no S at both ends
    def cigar_max_s(self):
        beginning_s, ending_s = self.get_soft_clipping()
        return max(beginning_s, ending_s)

    def cigar_total_soft_clipping(self):
        beginning_s, ending_s = self.get_soft_clipping()
        return beginning_s + ending_s

    def get_soft_clipping(self):
        tokenizer = re.compile(r'(?:\d+)|(?:[A-Z=])')
        cigar_elements = tokenizer.findall(self.cigar)

        if cigar_elements[1] == 'S':
            beginning_s = int(cigar_elements[0])
        else:
            beginning_s = 0

        if cigar_elements[-1] == 'S':
            ending_s = int(cigar_elements[-2])
        else:
            ending_s = 0

        return (beginning_s, ending_s)


class NotAvailableHandler(object):
    def __init__(self):
        self.count = 0
        self.na_forms = {'na', 'NA', 'Not Available', 'not available', 'Not available', 'nan'}
        self.encoded_na_form = 'NA;;'

    def encode_if_na(self, taxon):
        if taxon not in self.na_forms:
            return taxon
        self.count += 1
        return 'NA;;{}'.format(self.count)

    def decode_if_na(self, taxon):
        if not taxon.startswith(self.encoded_na_form):
            return taxon

        return 'NA'


parser = argparse.ArgumentParser(description='Bayesian-based LCA taxonomic classification method')

##### Required arguments #####
required = parser.add_argument_group('required arguments')
required.add_argument("-i", "--sam", help="Input SAM file", type=str, required=True)
required.add_argument('-q', '--reference', help="Reference fasta file", type=str, required=True)
required.add_argument("-r", "--tax", help="reference taxonomy file for the Database", type=str, required=True)

##### Taxonomy filtering arguments #####
taxoptions = parser.add_argument_group('taxonomy profiling options [filtering of hits]')
taxoptions.add_argument("-n", "--nper", help="number of times to bootstrap. Default: 100", type=int, default=100)
taxoptions.add_argument("-b", "--iset", help="minimum identity score to include", type=float, default=0.8)
taxoptions.add_argument('-l', '--length', help="minimum length of hit to include relative to query", type=float,
                        default=0.5)
taxoptions.add_argument('-s', '--softclipping', help='maximum soft clipped ratio to include', type=float, default=0.2)
##### Alignment control arguments #####
alignoptions = parser.add_argument_group('alignment control arguments')
alignoptions.add_argument("-m", "--match", default=1.0, help="alignment match score. Default: 1", type=float)
alignoptions.add_argument("-f", "--mismatch", default=-2.5, help="alignment mismatch penalty. Default: -2.5",
                          type=float)
alignoptions.add_argument("-g", "--ngap", default=-2.0, help="alignment gap penalty. Default: -2", type=float)
##### Other arguments #####
optional = parser.add_argument_group('other arguments')
optional.add_argument('-p', '--muscle', help='Path to call muscle default: muscle', default='muscle')

optional.add_argument("-o","--outfile",help="output file name. Default: <fasta>.blca.out",type=str)
optional.add_argument("-v","--votesfile",help="votes file name. Default: <fasta>.blca.votes.jsonlines",type=str)
optional.add_argument('--continue_mode', help="continue from a previous run by appending to the same output file", action='store_true')
optional.add_argument("--muscle_use_diags", help="pass the diag argument to muscle", action='store_true')
optional.add_argument("--muscle_max_iterations", help="set the max number of iterations for muscle", type=int, default=16)

##### parse arguments #####
args = parser.parse_args()

### bootstrap times ###
nper = args.nper  # number of bootstrap to permute
### Filter hits per query ###
iset = args.iset  # identify threshold
### Alignment options ###
ngap = args.ngap  # gap penalty
match = args.match  # match score
mismatch = args.mismatch  # mismatch penalty
min_length = args.length
sam_file_name = args.sam
outfile_name = args.outfile or (sam_file_name + '.blca.out')
votesfile_name = args.votesfile or (sam_file_name + '.blca.votes.jsonlines')
reference_fasta = args.reference
tax = args.tax
muscle_path = args.muscle
max_soft_clipping_allowed = args.softclipping
continue_mode = args.continue_mode
muscle_use_diags = args.muscle_use_diags
muscle_max_iterations = args.muscle_max_iterations


levels = ["superkingdom", "phylum", "class", "order", "family", "genus", "species"]


def check_program(prgname):
    '''Check whether a program has been installed and put in the PATH'''
    path = os.popen("which " + prgname).read().rstrip()
    if len(path) > 0 and os.path.exists(path):
        print(prgname + " is located in your PATH!")
    else:
        print("ERROR: " + prgname + " is NOT in your PATH, please set up " + prgname + "!")
        sys.exit(1)


def get_dic_from_aln(aln):
    '''Read in alignment and convert it into a dictionary'''
    alignment = AlignIO.read(aln, "clustal")
    alndic = {}
    for r in alignment:
        alndic[r.id] = list(r.seq)
    return alndic


def pairwise_score(alndic, query, match, mismatch, ngap):
    '''Calculate pairwise alignment score given a query'''
    nt = ["A", "C", "T", "G", "g", "a", "c", "t"]
    hitscore = {}
    for k, v in alndic.items():
        if k != query:
            hitscore[k] = 0
            for i in range(len(v)):
                if (alndic[query][i] in nt) and (v[i] in nt) and (alndic[query][i] == v[i]):
                    hitscore[k] += float(match)
                elif (alndic[query][i] not in nt) and (v[i] not in nt) and (alndic[query][i] == v[i]):
                    hitscore[k] += float(0)
                elif ((alndic[query][i] not in nt) or (v[i] not in nt)) and (alndic[query][i] != v[i]):
                    hitscore[k] += float(mismatch)
                elif (alndic[query][i] in nt) and (v[i] in nt) and (alndic[query][i] != v[i]):
                    hitscore[k] += float(ngap)
    total = float(sum(hitscore.values()))
    if total <= 0:
        total = 1
    for k, v in hitscore.items():
        hitscore[k] = v / total
    return hitscore


def random_aln_score(alndic, query, match, mismatch, ngap):
    '''Randomize the alignment, and calculate the score'''
    nt = ["A", "C", "T", "G", "g", "a", "c", "t"]
    idx = []
    for i in range(len(list(alndic.values())[0])):
        idx.append(random.choice(range(len(list(alndic.values())[0]))))

    hitscore = {}
    for k, v in alndic.items():
        if k != query:
            hitscore[k] = 0
            for i in idx:
                if (alndic[query][i] in nt) and (v[i] in nt) and (alndic[query][i] == v[i]):
                    hitscore[k] += float(match)
                elif (alndic[query][i] not in nt) and (v[i] not in nt) and (alndic[query][i] == v[i]):
                    hitscore[k] += float(0)
                elif ((alndic[query][i] not in nt) or (v[i] not in nt)) and (alndic[query][i] != v[i]):
                    hitscore[k] += float(mismatch)
                elif (alndic[query][i] in nt) and (v[i] in nt) and (alndic[query][i] != v[i]):
                    hitscore[k] += float(ngap)
    return hitscore


def get_gap_pos(query, alndic):
    '''Get the gap position in the alignment'''
    for i in range(len(alndic[query])):
        if alndic[query][i] != "-":
            start = i
            break
    for i in range(len(alndic[query]) - 1, 0, -1):
        if alndic[query][i] != "-":
            end = i
            break
    return start, end


def cut_gap(alndic, start, end):
    '''Given a start and end gap position, truncate the alignmnet'''
    trunc_alndic = {}
    for k_truc, v_truc in alndic.items():
        trunc_alndic[k_truc] = v_truc[start:end]
    return trunc_alndic


def read_tax_acc(taxfile, not_available_handler):
    tx = open(taxfile)
    acctax = {}
    for l in tx:
        lne = l.rstrip().strip(";").split("\t")
        if len(lne) != 2:
            continue
        if (levels[0] + ':') not in l:
            taxons = [not_available_handler.encode_if_na(taxon) for taxon in lne[1].split(';')]
            acctax[lne[0].split('.')[0]] = dict(zip(levels, taxons))
        else:
            pairs = [x.split(":", 1) for x in lne[1].split(";")]
            encoded = [(level, not_available_handler.encode_if_na(taxon)) for level, taxon in pairs]
            acctax[lne[0].split(".")[0]] = dict(encoded)
    tx.close()
    return acctax


################################################################
##
## 	Running Script Start
##
################################################################

## check whether muscle is located in the path
# check_program("muscle")

### read in pre-formatted lineage information ###
na_handler = NotAvailableHandler()
acc2tax = read_tax_acc(tax, na_handler)
print("> 1 > Read in taxonomy information!")

reference_sequences = {}
with open(reference_fasta) as f:
    for r in SeqIO.parse(f, "fasta"):
        reference_sequences[r.id] = str(r.seq)

print("> 2 > Read in reference db")

SequenceInfo = namedtuple('SequenceInfo', ['seq', 'hits'])
### read in input fasta file ###
input_sequences = {}
possible_rejects = set()
with open(sam_file_name) as sam_file:
    for line in sam_file:
        pieces = line.strip().split('\t')
        entry = SamEntry(pieces)

        if entry.identity_ratio_old < iset:
            possible_rejects.add(entry.qname)
        if entry.soft_clipped_ratio > max_soft_clipping_allowed:
            possible_rejects.add(entry.qname)
        elif entry.rname not in reference_sequences:
            possible_rejects.add(entry.qname)
        elif len(reference_sequences[entry.rname]) / float(len(entry.seq)) < min_length:
            possible_rejects.add(entry.qname)
        elif entry.qname not in input_sequences:
            input_sequences[entry.qname] = SequenceInfo(seq=entry.seq, hits=[entry.rname])
        else:
            input_sequences[entry.qname].hits.append(entry.rname)

rejects = possible_rejects.difference(set(input_sequences))

print("> 3 > Read in bowtie2 output!")

already_assigned = set()
if continue_mode:
    for line in open(outfile_name):
        already_assigned.add(line.split('\t')[0])

    outfile = open(outfile_name, 'a')
else:
    outfile = open(outfile_name, 'w')

for seqn, info in input_sequences.items():
    if seqn in already_assigned:
        continue

print("> 3 > Read in bowtie2 output!")

count = 0
outfile = open(outfile_name, 'w')
votesfile = open(votesfile_name, 'w')

for seqn, info in input_sequences.items():
    count += 1

    if seqn in acc2tax:
        print("[WARNING] Your sequence " + seqn + " has the same ID as the reference database! Please correct it!")
        print("...Skipping sequence " + seqn + " ......")
        outfile.write(seqn + "\tSkipped\n")
        continue

    ### Get all the hits list belong to the same query ###
    ### Add query fasta sequence to extracted hit fasta ###
    fifsa = []
    for hit in info.hits:
        if hit not in reference_sequences:
            print("Missing reference sequence for " + hit)
            continue
        fifsa.append(">{}\n{}\n".format(hit, reference_sequences[hit]))
    fifsa.append(">" + seqn + "\n" + info.seq)
    fifsa = "\n".join(fifsa)
    # os.system("rm " + seqn + ".dblist")

    ### Run muscle ###
    muscle_call = [muscle_path, '-quiet', '-clw', '-maxiters',  str(muscle_max_iterations)]
    if muscle_use_diags:
        muscle_call.append('-diags')
    proc = subprocess.Popen(muscle_call, stdout=subprocess.PIPE, stdin=subprocess.PIPE, stderr=subprocess.PIPE)
    outs, errs = proc.communicate(fifsa.encode('utf-8'))

    if proc.returncode == -signal.SIGSEGV:
        print("Error: segmentation fault")

    alndic = get_dic_from_aln(StringIO(outs.decode('utf-8')))
    # os.system("rm " + seqn + ".hits.fsa")
    # os.system("rm " + seqn + ".muscle")
    #    	print "Processing:",k1

    ### get gap position and truncate the alignment###
    ## SAARA: Was getting here a keyerror, but wanted to see if otherwise working so added the exception
    # try:
    #     start, end = get_gap_pos(seqn, alndic)
    # except KeyError:
    #     continue
    start, end = get_gap_pos(seqn, alndic)
    trunc_alndic = cut_gap(alndic, start, end)
    orgscore = pairwise_score(trunc_alndic, seqn, match, mismatch, ngap)
    ### start bootstrap ###
    perdict = {}  # record alignmet score for each iteration
    pervote = {}  # record vote after nper bootstrap

    for j in range(nper):
        random_scores = random_aln_score(trunc_alndic, seqn, match, mismatch, ngap)
        perdict[j] = random_scores
        max_score = max(random_scores.values())
        hits_with_max_score = [k3 for k3, v3 in random_scores.items() if v3 == max_score]
        vote_share = 1.0 / len(hits_with_max_score)
        for hit in hits_with_max_score:
            if hit in pervote:
                pervote[hit] += vote_share
            else:
                pervote[hit] = vote_share

    ### normalize vote by total votes ###
    ttlvote = sum(pervote.values())
    for k4, v4 in pervote.items():
        pervote[k4] = v4 / ttlvote * 100
    ###

    votes_by_level = {}
    for level in levels:
        votes_by_level[level] = defaultdict(int)

    for hit in orgscore.keys():
        short_hit_name = hit.split(".")[0]
        if short_hit_name not in acc2tax:
            print("Missing taxonomy info for ", short_hit_name)
            continue
        hit_taxonomy = acc2tax[short_hit_name]
        for level in levels:
            # deal with missing values in the taxonomy
            if level not in hit_taxonomy:
                hit_taxonomy[level] = na_handler.encode_if_na("NA")

            if hit in pervote:
                votes_by_level[level][hit_taxonomy[level]] += pervote[hit]
            else:
                votes_by_level[level][hit_taxonomy[level]] += 0

    votesfile.write(json.dumps({"seqn": seqn, "votes": votes_by_level}) + "\n")

    try:
        outfile.write(seqn + "\t")
        for level in levels:
            levels_votes = votes_by_level[level]
            outfile.write(level + ":" + na_handler.decode_if_na(max(levels_votes, key=levels_votes.get)) + ";")
        outfile.write("\t")
        for level in levels:
            levels_votes = votes_by_level[level]
            outfile.write(level + ":" + str(max(levels_votes.values())) + ";")
        outfile.write("\t" + ";".join(info.hits))
        outfile.write("\n")
    except ValueError as e:
        print(f"ValueError for {seqn}, possibly due to missing votes")
        outfile.write(seqn + "\tUnclassified\n")

for seqn in rejects:
    outfile.write(seqn + "\tUnclassified\n")

outfile.close()
votesfile.close()

library(plyr)
library(dplyr)
library(yaml)
library(dada2)
library(Biostrings)
library(ggplot2)
library(tidyr)
library(tibble)

args <- commandArgs(trailingOnly = T)
config <- read_yaml(args[2])

# Set the different paths for all the supplied libraries
# NOTICE: only forward files given as input
filtFs <- args[3:length(args)]
filtFs_single <- gsub("_1P", "_1U", filtFs)

filtRs <- gsub("_1P", "_2P", filtFs)
filtRs_single <- gsub("_2P", "_2U", filtRs)

outpath <- args[1]

# The parameters nti/ntj, need the chosen nucleotides in a vector format:
nti <- strsplit(config$DADA2$plotERRORS$nti, "")[[1]]
if (length(nti) == 0) {
  message("Default value given to nti, all nucleotide transitions will be shown in the error plot")
  nti <- c("A", "C", "G", "T")
}

ntj <- strsplit(config$DADA2$plotERRORS$ntj, "")[[1]]
if (length(nti) == 0) {
  message("Default value given to ntj, all nucleotide transitions will be shown in the error plot")
  ntj <- c("A", "C", "G", "T")
}

#For the analysis of Novaseq data, use a different error model, because of binned quality values
#https://github.com/ErnakovichLab/dada2_ernakovichlab/tree/split_for_premise
#https://github.com/benjjneb/dada2/issues/1307 and issue 791

loessErrfun_mod4 <- function(trans) {
  qq <- as.numeric(colnames(trans))
  est <- matrix(0, nrow=0, ncol=length(qq))
  for(nti in c("A","C","G","T")) {
    for(ntj in c("A","C","G","T")) {
      if(nti != ntj) {
        errs <- trans[paste0(nti,"2",ntj),]
        tot <- colSums(trans[paste0(nti,"2",c("A","C","G","T")),])
        rlogp <- log10((errs+1)/tot)  # 1 psuedocount for each err, but if tot=0 will give NA
        rlogp[is.infinite(rlogp)] <- NA
        df <- data.frame(q=qq, errs=errs, tot=tot, rlogp=rlogp)

        # original
        # ###! mod.lo <- loess(rlogp ~ q, df, weights=errs) ###!
        # mod.lo <- loess(rlogp ~ q, df, weights=tot) ###!
        # #        mod.lo <- loess(rlogp ~ q, df)

        # jonalim's solution
        # https://github.com/benjjneb/dada2/issues/938
        mod.lo <- loess(rlogp ~ q, df, weights = log10(tot),degree = 1, span = 0.95)

        pred <- predict(mod.lo, qq)
        maxrli <- max(which(!is.na(pred)))
        minrli <- min(which(!is.na(pred)))
        pred[seq_along(pred)>maxrli] <- pred[[maxrli]]
        pred[seq_along(pred)<minrli] <- pred[[minrli]]
        est <- rbind(est, 10^pred)
      } # if(nti != ntj)
    } # for(ntj in c("A","C","G","T"))
  } # for(nti in c("A","C","G","T"))

  # HACKY
  MAX_ERROR_RATE <- 0.25
  MIN_ERROR_RATE <- 1e-7
  est[est>MAX_ERROR_RATE] <- MAX_ERROR_RATE
  est[est<MIN_ERROR_RATE] <- MIN_ERROR_RATE

  # enforce monotonicity
  # https://github.com/benjjneb/dada2/issues/791
  estorig <- est
  est <- est %>%
    data.frame() %>%
    mutate_all(funs(case_when(. < X40 ~ X40,
                              . >= X40 ~ .))) %>% as.matrix()
  rownames(est) <- rownames(estorig)
  colnames(est) <- colnames(estorig)

  # Expand the err matrix with the self-transition probs
  err <- rbind(1-colSums(est[1:3,]), est[1:3,],
               est[4,], 1-colSums(est[4:6,]), est[5:6,],
               est[7:8,], 1-colSums(est[7:9,]), est[9,],
               est[10:12,], 1-colSums(est[10:12,]))
  rownames(err) <- paste0(rep(c("A","C","G","T"), each=4), "2", c("A","C","G","T"))
  colnames(err) <- colnames(trans)
  # Return
  return(err)
}

# Get sample names
sample.names <- gsub("_1P.fastq.gz", "", basename(filtFs))
message(paste0("Sample ", sample.names, " will be analyzed", collapse = "\n"))

# Assign names to files
names(filtFs) <- sample.names
names(filtRs) <- sample.names
names(filtFs_single) <- sample.names
names(filtRs_single) <- sample.names

allfiles <- list(filtFs, filtRs, filtFs_single, filtRs_single)
names(allfiles) <- c("filtFs", "filtRs", "filtFs_single", "filtRs_single")
errs <- list()
dereps <- list()
dadas <- list()
seqtab <- list()
files_exist <- list()


for (i in 1:4) {

  message("Check that files are not empty")

  files_loop <- c()

  for (j in 1:length(allfiles[[i]])) {
    info <- file.info(allfiles[[i]][j])
    if (!is.na(info$size) & info$size > 20) { # 20 is the size of an empty .gz file
      files_loop <- c(files_loop, allfiles[[i]][j])
    }
  }

  if (is.null(files_loop)) {
    files_exist[i] <- list(NULL)
  } else {
    files_exist[[i]] <- files_loop
  }
}

# Update sample names with existing paired files
sample.names<-names(files_exist[[1]])

# Loop through all file types (forward, reverse, unpaired forward, unpaired reverse) for learning errors and dereplicating
for (i in 1:4) {

  if (length(files_exist[[i]])!=0) {        #any(file.exists(files_exist[[i]]))

    message(paste("learning error rates of files:", i, ": (1) forward paired (2) reverse paired (3) forward single and (4) reverse single " , sep=" "))

    #Add here a more generic character matching
    if (config$meta$sequencing$seq_meth=="NovaSeq6000"){ 
    message("learning error rates using a modified error model for NovaSeq data")

    errs[[i]] <- learnErrors(files_exist[[i]], #allfiles[[i]][file.exists(allfiles[[i]])]
                            multithread = config$DADA2$learnERRORS$multithread,
                            nbases = as.numeric(config$DADA2$learnERRORS$nbases),
                            randomize = config$DADA2$learnERRORS$randomize,
                            MAX_CONSIST = as.numeric(config$DADA2$learnERRORS$MAX_CONSIST),
                            OMEGA_C = as.numeric(config$DADA2$learnERRORS$OMEGA_C),
                            verbose = config$DADA2$learnERRORS$verbose,
                            errorEstimationFunction = loessErrfun_mod4,
                            )


    } else {


    errs[[i]] <- learnErrors(files_exist[[i]], #allfiles[[i]][file.exists(allfiles[[i]])]
                            multithread = config$DADA2$learnERRORS$multithread,
                            nbases = as.numeric(config$DADA2$learnERRORS$nbases),
                            randomize = config$DADA2$learnERRORS$randomize,
                            MAX_CONSIST = as.numeric(config$DADA2$learnERRORS$MAX_CONSIST),
                            OMEGA_C = as.numeric(config$DADA2$learnERRORS$OMEGA_C),
                            verbose = config$DADA2$learnERRORS$verbose)

    }

    message("Making error estimation plots of reads")
    png(filename = paste0(outpath, "06-report/dada2/error_profile_", names(allfiles)[i], ".png"))
      p_ERR <- plotErrors(
        errs[[i]],
        nti = nti,
        ntj = ntj,
        obs = config$DADA2$plotERRORS$obs,
        err_out = config$DADA2$plotERRORS$err_out,
        err_in = config$DADA2$plotERRORS$err_in,
        nominalQ = config$DADA2$plotERRORS$nominalQ)
      print(p_ERR)
    dev.off()

    message("Running dereplication of reads")
    #print(paste("files that exist:", files_exist[[i]]))
    dereps[[i]] <- derepFastq(files_exist[[i]], #allfiles[[i]][file.exists(allfiles[[i]])]
                              n = as.numeric(config$DADA2$derepFastq$num),
                              config$DADA2$learnERRORS$verbose)

    message("Running dada on reads")
    dadas[[i]] <- dada(dereps[[i]],
                      errs[[i]],
                      selfConsist = config$DADA2$dada$selfConsist,
                      pool = config$DADA2$dada$pool,
                      priors = config$DADA2$dada$priors,
                      multithread = config$DADA2$learnERRORS$multithread,
                      verbose = config$DADA2$learnERRORS$verbose)

    # Convert to list in case there's only one file
    if (is(dadas[[i]], "dada")) {
      message("Convert to list in case there's only one file")
      dadas[[i]] <- list(dadas[[i]])
      names(dadas[[i]]) <- names(files_exist[[i]])#allfiles[[i]][file.exists(allfiles[[i]])])
    }
    if (is(dereps[[i]], "derep")) {
      dereps[[i]] <- list(dereps[[i]])
      names(dereps[[i]]) <- names(files_exist[[i]])#allfiles[[i]][file.exists(allfiles[[i]])])
    }

  # If no files found for the paired reads (should not be the case!)
  } else if (!grepl("single", names(allfiles)[i])) {

    message("Error: no paired reads to process")
    stop()

  # If no files found for the unpaired reads, continue with the workflow
  # It has to be made sure in the snakefile that this step is run despite not requiring output files
  } else {

    message("No further unpaired reads to process")

  }
}

# Merge forward and reverse paired reads
message("Attempting merge")

# Define function for merging and formatting seqtabs: next steps require an integer matrix
merge_format_seqtab <- function(seqtab1, seqtab2) {
  df1 <- reshape2::melt(seqtab1, varnames = c("sample", "sequence"))
  df2 <- reshape2::melt(seqtab2, varnames = c("sample", "sequence"))
  df <- bind_rows(df1, df2) %>%
    group_by(sample, sequence) %>%
    summarize(value = sum(value)) %>%
    ungroup()
  m <- reshape2::acast(df, sample ~ sequence, value.var = "value")
  mode(m) <- "integer"
  return(m)
}

if (config$DADA2$mergePairs$include) {
  message("merging pairs")
  mergers <- mergePairs(dadas[[1]],
                        dereps[[1]],
                        dadas[[2]],
                        dereps[[2]],
                        minOverlap = as.numeric(config$DADA2$mergePairs$minOverlap),
                        maxMismatch = as.numeric(config$DADA2$mergePairs$maxMismatch),
                        returnRejects = config$DADA2$mergePairs$returnRejects,
                        propagateCol = config$DADA2$mergePairs$propagateCol,
                        justConcatenate = config$DADA2$mergePairs$justConcatenate,
                        trimOverhang = config$DADA2$mergePairs$trimOverhang,
                        verbose = config$DADA2$learnERRORS$verbose)
  if (is.data.frame(mergers)) {
    mergers <- list(mergers) %>% setNames(sample.names)
  }

  # Create read/ASV mapping
  mergers_all <- mergePairs(dadas[[1]], dereps[[1]], dadas[[2]], dereps[[2]],
    minOverlap = as.numeric(config$DADA2$mergePairs$minOverlap), maxMismatch = as.numeric(config$DADA2$mergePairs$maxMismatch),
    returnRejects = TRUE, propagateCol = config$DADA2$mergePairs$propagateCol,
    justConcatenate = config$DADA2$mergePairs$justConcatenate, trimOverhang = config$DADA2$mergePairs$trimOverhang,
    verbose = config$DADA2$learnERRORS$verbose)
  if (is.data.frame(mergers_all)) {
    mergers_all <- list(mergers_all) %>% setNames(sample.names)
  }

  mapping <- lapply(names(mergers_all), function(name) {
    merger_to_dada <- lapply(mergers_all[[name]]$forward, function(x) {
      which(dadas[[1]][[name]]$map == x)
    })
    merger_to_read <- lapply(merger_to_dada, function(x) {
      which(dereps[[1]][[name]]$map %in% x)
    })
    names(merger_to_read) <- mergers_all[[name]]$sequence
    merger_to_read
  })
  names(mapping) <- names(mergers_all)
  for (name in names(mapping)) {
    fq <- microseq::readFastq(filtFs[[name]])
    headers <- sub(" .*", "", fq$Header)
    mapping[[name]] <- mapping[[name]] %>%
      enframe(name = "sequence", value = "read") %>%
      unnest() %>%
      filter(sequence != "") %>%
      mutate(read = headers[read]) %>%
      arrange(sequence, read)
  }

  seqtab <- makeSequenceTable(mergers)

  # When merging is done with returnRejects=TRUE, the abundance of the rejected merges is returned, but not the sequence
  # We want to collect also these single sequences and add them to the seqtab (to avoid loosing ANY data)
  if (config$DADA2$mergePairs$returnRejects == TRUE) {

    unmerged_f <- list()
    unmerged_r <- list()
    concatenated <- list()

    for (i in 1:length(sample.names)) {
      unmerged_f[[i]] <- dadas[[1]][[sample.names[i]]]$sequence[mergers[[sample.names[i]]]$forward[!mergers[[sample.names[i]]]$accept]]
      unmerged_r[[i]] <- dadas[[2]][[sample.names[i]]]$sequence[mergers[[sample.names[i]]]$reverse[!mergers[[sample.names[i]]]$accept]]
      # Here for the rejected reads (!merger$sample$accept) the indices are collected (merger$sample$forward, merger$sample$reverse)
      # The sequences are sourced from the original dada-file (dadaF$sample$denoised, dadaR$sample$denoised)
      # It seems that concatenating these reads and keeping them for further analyses can result in better taxonomic coverage (Dacey et al. 2021 https://doi.org/10.1186/s12859-021-04410-2)
      # Abundances for these reads is taken from the merged abundances.
      # reverse complement reverse reads so that the following taxonomic assignment will work optimally.
      unmerged_r[[i]] <- sapply(sapply(sapply(unmerged_r[[i]], DNAString), Biostrings::reverseComplement), toString)
      sequence <- paste0(unmerged_f[[i]], unmerged_r[[i]])
      abundance <- mergers[[sample.names[i]]]$abundance[!mergers[[sample.names[i]]]$accept]
      concatenated[[i]] <- tibble(sequence, abundance)
    }
    names(concatenated) <- sample.names
    #names(unmerged_r)=sample.names

    # Make sequence table
    seqtab_unmerged <- makeSequenceTable(concatenated)

    # The merged returnrejects = T seqtab also contains a column with an empty header,
    # This is all rejected (non-merged) abundances combined.
    # This column messes with future steps, so we want to remove it.
    # We have instead collected the abundances of the unmerged reads to add to the table with sequences.
    #seqtab <- seqtab[,-which(colnames(seqtab) == "")]

    # Merge with paired reads and format for the next steps (integer matrix)
    seqtab <- merge_format_seqtab(seqtab, seqtab_unmerged)

  }

} else {

  message("no merging of paired reads")
  # Combine forward and reverse sequences to one table
  seqtab1 <- makeSequenceTable(dadas[[1]])
  seqtab2 <- makeSequenceTable(dadas[[2]])
  # Reverse complement reverse reads so that the following taxonomic assignment will work optimally.
  colnames(seqtab2) <- sapply(sapply(sapply(colnames(seqtab2), DNAString), Biostrings::reverseComplement), toString)
  seqtab <- cbind(seqtab1, seqtab2)

}

# Add ASVs from single reads to full table, and format table to the right format to continue with the pipeline
# It has to be an integer matrix with samples as rownames and sequences as column names
if (length(files_exist[[3]])!=0) {
  message("Adding ASVs from unpaired forward reads to ASV-table")
  seqtab3 <- makeSequenceTable(dadas[[3]])
  seqtab <- merge_format_seqtab(seqtab, seqtab3)
}

if (length(files_exist[[4]])!=0) {
  message("Adding ASVs from unpaired reverse reads to ASV-table")
  seqtab4 <- makeSequenceTable(dadas[[4]])
  # Here also the sequences from the reverse reads are reverse complemented before they are added to the sequence table
  colnames(seqtab4) <- sapply(sapply(sapply(colnames(seqtab4), DNAString), reverseComplement), toString)
  seqtab <- merge_format_seqtab(seqtab, seqtab4)
}

# Chimeras removed from the full combined table as per: https://github.com/benjjneb/dada2/issues/1235:
# We think the best way (in most cases, using current common techs -- such is the challenge of recommendations) is to combine the tables from multiple runs and then remove chimeras on that combined table.
# Look into details of how chimera removal is done to understand if this is smart
# In their case they are looking at equal length reads, possibly the single reads and paired reads should not be combined for chimera removal?

message("removing chimeras")
seqtab.nochim <- removeBimeraDenovo(seqtab,
                                    method = config$DADA2$removeBimeraDenovo$method,
                                    multithread = config$DADA2$learnERRORS$multithread,
                                    verbose = config$DADA2$learnERRORS$verbose)

print(dim(seqtab.nochim))
# ASVs denominated by the actual sequence, we want to simplify the names.
new.names <- c(paste("asv.", 1:length(colnames(seqtab.nochim)), sep = ""))
message(head(new.names))

# Save fasta, before changing names in the seqtab table
message("Making fasta table")
uniquesToFasta(seqtab.nochim, fout = paste0(outpath, "03-dada2/rep-seqs.fna"), ids = new.names)

# Export read/asv mapping
if (exists("mapping")) {
  message("Exporting read/asv mapping")
  asv_sequences <- data.frame(asv = new.names, sequence = colnames(seqtab.nochim))
  for (name in names(mapping)) {
    mapping[[name]] <- mapping[[name]] %>%
      left_join(asv_sequences, by = "sequence") %>%
      select(read, asv) %>%
      filter(!is.na(asv))
    dir.create(file.path(paste0(outpath, "03-dada2/mapping/")))
    dir.create(file.path(paste0(outpath, "03-dada2/mapping/", name)))
    write.table(mapping[[name]], paste0(outpath, "03-dada2/mapping/", name, "/", name, "_mapping.txt"), sep = "\t", col.names = TRUE, row.names = FALSE, quote = FALSE)
  }
}

# Finally, change names in otu table
message("Changing sequence names")
colnames(seqtab.nochim) <- new.names
#Replace any NAs that are possibly in the seqtab with 0s
seqtab.nochim[is.na(seqtab.nochim)] <- 0


# This show sequence length distributions (see if you should include this)
#seq_hist <- table(nchar(getSequences(seqtab)))
#fname_seqh <- paste(args[6],"seq_hist.txt",sep="")
#write.table(seq_hist, file = fname_seqh  , sep = "\t", quote=FALSE, col.names = FALSE)

# Collect results of how many reads are available at each step in a table:
getN <- function(x) sum(getUniques(x))

# Make a table with all information on the reads retained from the run, if paired reads were merged:
message("Making summary table")
if (config$DADA2$mergePairs$include) {
  track <- cbind(sapply(dadas[[1]], getN), sapply(dadas[[2]], getN), sapply(mergers, getN), rowSums(seqtab.nochim), rowSums(seqtab.nochim != 0))
  colnames(track) <- c("denoisedF", "denoisedR", "merged", "nonchim", "ASVs")
# If paired reads were not merged:
} else {
  track <- cbind(sapply(dadas[[1]], getN), sapply(dadas[[2]], getN), rowSums(seqtab.nochim), rowSums(seqtab.nochim != 0))
  colnames(track) <- c("denoisedF", "denoisedR", "nonchim", "ASVs")
}

# Add also information on the reads that came from possibly evaluated single reads
if (length(files_exist[[3]])!=0) {
  message("Adding info from unpaired forward reads to the summary table")
  denoisedF_single <- sapply(dadas[[3]], getN)
  track <- merge(track, data.frame(denoisedF_single), by = 0, all = TRUE)
  rownames(track) <- track$Row.names
  track <- subset(track, select = -Row.names)
  # Reorder columns:
  track <- track %>% relocate(denoisedF_single, .before = nonchim)
}

if (length(files_exist[[4]])!=0) {
  message("Adding info from unpaired reverse reads to the summary table")
  denoisedR_single <- sapply(dadas[[4]], getN)
  track <- merge(track, data.frame(denoisedR_single), by = 0, all = TRUE)
  rownames(track) <- track$Row.names
  track <- subset(track, select = -Row.names)
  track <- track %>% relocate(denoisedR_single, .before = nonchim)
}

rownames(track)
#rownames(track) <- sample.names
#message(track)

# Read results of filtering step and append the results of ASV step:
out <- read.table(paste0(outpath, "06-report/dada2/dada2_filtering_stats.txt"), header = TRUE)
#track <- cbind(out, track)
track <- merge(out, track, by = 0, all = TRUE)
rownames(track) <- track$Row.names
track <- subset(track, select = -Row.names)


# Write output tables of reads and the otu_table:
write.table(track, paste0(outpath, "06-report/dada2/dada2_stats.txt"), row.names = TRUE, col.names = TRUE, quote = FALSE)
write.table(t(seqtab.nochim), paste0(outpath,"03-dada2/seqtab-nochim.txt"), sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE)

# Fix for https://github.com/tidyverse/ggplot2/issues/2787
if (file.exists("Rplots.pdf")) {
  file.remove("Rplots.pdf")
}