biom export

dada2.smk

@@ -13,7 +13,8 @@ rule makeSequenceTable:
 		expand("work/dada/{sample}.rds", sample=SAMPLES)
 	output:
 		table = "work/dada/seqtab.rds",
-		fasta = "work/dada/seqtab.fasta"
+		fasta = "work/dada/seqtab.fasta",
+		biom = "work/dada/dada2.biom"
 	threads:
 		config["THREADS"]
 	script:

makeSequenceTable.R

 library(dada2)
+library(biomformat)
 
-seqtab <- makeSequenceTable(lapply(snakemake@input, readRDS), orderBy = "abundance")
+dada_list <- lapply(snakemake@input, readRDS)
+names(dada_list) <- lapply(snakemake@input, function(x){basename(tools::file_path_sans_ext(x))})
+
+seqtab <- makeSequenceTable(dada_list, orderBy = "abundance")
 seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread = snakemake@threads, verbose = TRUE)
+
 saveRDS(seqtab.nochim, snakemake@output$table)
+
 uniquesToFasta(seqtab.nochim, snakemake@output$fasta)
+
+biom <- make_biom(seqtab.nochim)
+write_biom(biom, snakemake@output$biom)