Skip to content

Commit

Permalink
Fix: Replace biomformat with external biom in remove_bimera
Browse files Browse the repository at this point in the history
The `biomformat` library is not well maintained and is only capable of
creating biom files of v1.0.0. This does not work properly when there is
only one column in the OTU table. Hence, we decided to use the external
`biom` command to create the table
  • Loading branch information
dileep-kishore committed Nov 23, 2021
1 parent c3b5932 commit 6e3c609
Showing 1 changed file with 10 additions and 11 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -3,27 +3,26 @@
suppressWarnings(library(dada2))

chimera.method <- "${chimera_method}"
multithread <- ${ncpus}
multithread <- as.integer("${ncpus}")

# load table and convert to matrix
seq.table.file <- "${seqtable_file}"
seq.table <- read.csv(seq.table.file, sep="\\t", check.names=FALSE)
rownames(seq.table) <- seq.table[,1]
seq.table <- seq.table[ -c(1) ]
seq.table <- read.csv(seq.table.file, sep = "\\t", check.names = FALSE)
rownames(seq.table) <- seq.table[, 1]
seq.table <- seq.table[-c(1)]
seq.mat <- data.matrix(seq.table)
mode(seq.mat) <- "integer"

# chimera removal
seq.table.nochim <- removeBimeraDenovo(t(seq.mat), method=chimera.method, multithread=multithread, verbose=TRUE)
seq.table.nochim <- removeBimeraDenovo(t(seq.mat), method = chimera.method, multithread = multithread, verbose = TRUE)

# Export sequences
uniquesToFasta(
getUniques(seq.table.nochim),
fout="unhashed_rep_seqs.fasta",
ids=paste0("Seq", seq(length(getUniques(seq.table.nochim))))
getUniques(seq.table.nochim),
fout = "unhashed_rep_seqs.fasta",
ids = paste0("Seq", seq(length(getUniques(seq.table.nochim))))
)
write.table(t(seq.table.nochim), file = "seq_table_nochim.tsv", col.names = NA, sep = "\\t", quote = FALSE)

# Export to biom
library(biomformat)
st.biom <- make_biom(t(seq.table.nochim))
write_biom(st.biom, "unhashed_otu_table.biom")
system("biom convert -i seq_table_nochim.tsv -o unhashed_otu_table.biom --to-hdf5")

0 comments on commit 6e3c609

Please sign in to comment.