require(dada2, quietly=TRUE)

#Parameters
if (is.null(snakemake@params$maxN)) {snakemake@params$maxN <- 0}
if (is.null(snakemake@params$truncQ)) {snakemake@params$truncQ <- 2}
if (is.null(snakemake@params$truncLen)) {snakemake@params$truncLen <- 0}
if (is.null(snakemake@params$maxEE_R1)) {snakemake@params$maxEE_R1 <- Inf}
if (is.null(snakemake@params$maxEE_R2)) {snakemake@params$maxEE_R2 <- Inf}
if (is.null(snakemake@params$rm_phix)) {snakemake@params$rm_phix <- TRUE}
if (is.null(snakemake@params$minLen)) {snakemake@params$minLen <- 0}

maxEE <- c(snakemake@params$maxEE_R1, snakemake@params$maxEE_R2)

# Performance
threads <- snakemake@threads

out <- filterAndTrim(fwd = snakemake@input$R1, filt = snakemake@output$R1,
              rev = snakemake@input$R2, filt.rev = snakemake@output$R2,
              maxN = snakemake@params$maxN, truncQ = snakemake@params$truncQ,
              maxEE = maxEE, rm.phix = snakemake@params$rm_phix,
              minLen = snakemake@params$minLen, multithread = threads,
              verbose=TRUE)

saveRDS(out, snakemake@output$rds)

args = commandArgs(trailingOnly=TRUE)

itsx_in <- snakemake@input$its_in
itsx_out <- snakemake@input$its_out
sq.nc <- snakemake@input$seqtab
out <- snakemake@output$seqtab
out_mock <- snakemake@output$mock
threads <- snakemake@threads
mock <- snakemake@params$mock

mock_samples <- strsplit(mock, ",")[[1]]

library(dada2); packageVersion("dada2")
library(dplyr)
library(Biostrings)

dna <- readDNAStringSet(itsx_in)
dna_itsx <- as.character(readDNAStringSet(itsx_out))
seqtab.nc <- readRDS(sq.nc)
colnames(seqtab.nc) <- names(dna)
seqtab.nc <- seqtab.nc[,colnames(seqtab.nc)%in%names(dna_itsx)]

colnames(seqtab.nc) <- dna_itsx
seqtab.nc <- t(seqtab.nc)

#Remove sequences below 50bp
if (any(nchar(getSequences(t(seqtab.nc)))<50)) {
  dna_itsx <- dna_itsx[!(nchar(dna_itsx)<50)]
  seqtab.nc <- seqtab.nc[nchar(rownames(seqtab.nc))>50,]
}


#At this point we can summarise all identical sequences 
seqtab.nc2 <- cbind.data.frame(sequence=rownames(seqtab.nc), seqtab.nc)
seqtab.nc3 <- group_by(seqtab.nc2, sequence) %>% 
  summarise_each(funs(sum))

if (length(intersect(colnames(seqtab.nc3), mock_samples)) > 0) {
  mock_cols <- append(c("sequence"), intersect(colnames(seqtab.nc3), mock_samples))
  mock <- seqtab.nc3 %>% select(all_of(mock_cols))
  seqtab.nc4 <- seqtab.nc3 %>% select(!any_of(mock_cols))
  saveRDS(mock, out_mock)
} else {
  seqtab.nc4 <- seqtab.nc3[,colnames(seqtab.nc3)!="sequence"]
}

rownames(seqtab.nc4) <- seqtab.nc3$sequence
seqtab.nc4 <- data.matrix(t(seqtab.nc4))

saveRDS(seqtab.nc4, out)

from argparse import ArgumentParser
from Bio.Seq import reverse_complement


def revcomp(seq):
    return reverse_complement(seq)


def main(args):
    r = revcomp(args.seq)
    print(r, end="")


if __name__ == "__main__":
    parser = ArgumentParser()
    parser.add_argument("seq", type=str,
                        help="Input fasta string")
    args = parser.parse_args()
    main(args)

shell:
    """
    A=$(cat {input.fwd_rc})
    a=$(cat {input.rev_rc})

    cutadapt -g {params.FWD} -a $a -G {params.REV} -A $A -j {threads} \
     -n {params.n} -o {output.R1} -p {output.R2} --minimum-length {params.min_len} \
     {input.R1} {input.R2} > {log} 2>&1
    """

shell:
     """
     ITSx -i {input} -o {params.prefix}/itsx --cpu {threads} --multi_thread T \
        --preserve T --partial 50 --minlen 50 > {log} 2>&1
     """

shell:
    """
    swarm -t {threads} -d 3 -z --output-file {output.txt} \
        --seeds {output.seeds} {input} > {log} 2>&1
    """